solve_LINEQ_mf_contact.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2026 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module m_solve_lineq_mf_contact
  use hecmw_util
  use m_hecmw_comm_f
  use hecmw_matrix_misc
  use hecmw_matrix_ass
  use hecmw_matrix_dump
  use hecmw_solver_direct_mf
  use hecmw_ebc_defer
 
  private
  public :: solve_lineq_mf_contact_init
  public :: solve_lineq_mf_contact
 
  logical, save :: INITIALIZED = .false.
  logical, save :: NEED_ANALYSIS = .true.
  logical, save :: IS_SYMMETRIC = .true.
  type (hecmwST_matrix), save :: mfMAT
 
  integer(kind=kint), parameter :: MFC_IR_STEPS = 5
  real(kind=kreal), parameter :: mfc_ir_tol = 1.0d-8
 
contains
 
  subroutine solve_lineq_mf_contact_init(hecMESH,hecMAT,hecLagMAT,is_sym)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(inout) :: hecmat
    type (hecmwst_matrix_lagrange), intent(in) :: heclagmat
    logical, intent(in) :: is_sym
 
    if (initialized) then
      call hecmw_mat_finalize(mfmat)
      initialized = .false.
    endif
    call hecmw_mat_init(mfmat)
    is_symmetric = is_sym
    need_analysis = .true.
    initialized = .true.
  end subroutine solve_lineq_mf_contact_init
 
  subroutine solve_lineq_mf_contact(hecMESH,hecMAT,hecLagMAT,hecEBC,istat)
    implicit none
    type (hecmwst_local_mesh), intent(inout) :: hecmesh
    type (hecmwst_matrix), intent(inout) :: hecmat
    type (hecmwst_matrix_lagrange), intent(inout) :: heclagmat
    type (hecmwst_ebc), intent(inout) :: hecebc
    integer(kind=kint), intent(out) :: istat
    integer(kind=kint) :: ntdf, mpc_method, iterlog, timelog, irmax, it
    real(kind=kreal) :: bnrm, rnrm, rprev, irtol, t1, t2
    real(kind=kreal), allocatable :: bb(:), xx(:), rr(:), xp(:)
    ! the unit hecmw_solve passes to the direct solvers for error messages
    integer(kind=kint), parameter :: imsg = 51
 
    mpc_method = hecmw_mat_get_mpc_method(hecmat)
    if (mpc_method < 1 .or. 3 < mpc_method) then
      mpc_method = 1
      call hecmw_mat_set_mpc_method(hecmat,mpc_method)
    endif
    if (mpc_method /= 1) then
      write(*,*) 'ERROR: MPCMETHOD other than penalty is not available for DIRECTmf solver', &
          ' in contact analysis without elimination'
      stop
    endif
    call hecmw_mat_ass_equation(hecmesh, hecmat)
    call hecmw_mat_ass_equation_rhs(hecmesh, hecmat)
    call hecmw_ebc_apply(hecmesh, hecmat, hecebc)
 
    call hecmw_mat_dump(hecmat, hecmesh)
 
    if (need_analysis) then
      call mf_contact_set_profile(hecmat, heclagmat)
      need_analysis = .false.
    endif
 
    call mf_contact_set_values(hecmat, heclagmat)
    ntdf = hecmat%NP*hecmat%NDOF + heclagmat%num_lagrange
    mfmat%B(1:ntdf) = hecmat%B(1:ntdf)
    mfmat%X(1:ntdf) = 0.0d0
    mfmat%Iarray(97) = 1
    mfmat%symmetric = is_symmetric
 
    call hecmw_solve_direct_mf(hecmesh, mfmat, imsg)
    istat = 0
 
    ! refine on the saddle-point system itself: the refinement of the driver measures its
    ! norms through the mesh, which does not cover the Lagrange rows
    irmax = hecmat%Iarray(44)
    if (irmax == 0) irmax = mfc_ir_steps
    if (irmax > 0) then
      iterlog = hecmw_mat_get_iterlog(hecmat)
      timelog = hecmw_mat_get_timelog(hecmat)
      irtol = hecmat%Rarray(42)
      if (.not. (irtol > 0.0d0)) irtol = mfc_ir_tol
      allocate(bb(ntdf), xx(ntdf), rr(ntdf), xp(ntdf))
      bb(1:ntdf) = mfmat%B(1:ntdf)
      xx(1:ntdf) = mfmat%X(1:ntdf)
      bnrm = sqrt(dot_product(bb(1:ntdf), bb(1:ntdf)))
      t1 = hecmw_wtime()
      rprev = huge(rprev)
      do it = 1, irmax
        call mf_contact_resid(ntdf, xx, bb, rr)
        rnrm = sqrt(dot_product(rr(1:ntdf), rr(1:ntdf))) / max(bnrm, tiny(bnrm))
        if (iterlog > 0) write(*,'(a,i0,a,1pe11.4)') '[DIRECTmf]: refinement ', it - 1, ': residual = ', rnrm
        if (rnrm <= irtol) exit
        ! a step that grew the residual would hand back a worse solution than it got,
        ! so the step is undone and the sweep stops with a warning
        if (rnrm > rprev) then
          xx(1:ntdf) = xp(1:ntdf)
          write(*,'(a,1pe11.4)') '[DIRECTmf]: WARNING: refinement diverged, reverted to residual = ', rprev
          exit
        endif
        rprev = rnrm
        xp(1:ntdf) = xx(1:ntdf)
        mfmat%B(1:ntdf) = rr(1:ntdf)
        call hecmw_solve_direct_mf(hecmesh, mfmat, imsg)
        xx(1:ntdf) = xx(1:ntdf) + mfmat%X(1:ntdf)
      enddo
      t2 = hecmw_wtime()
      if (timelog > 0) write(*,'(a,f10.3,a)') '[DIRECTmf]: refinement done (', t2 - t1, ' sec)'
      mfmat%X(1:ntdf) = xx(1:ntdf)
      deallocate(bb, xx, rr, xp)
    endif
 
    hecmat%X(1:ntdf) = mfmat%X(1:ntdf)
 
    call hecmw_mat_dump_solution(hecmat)
  end subroutine solve_lineq_mf_contact
 
  subroutine mf_contact_set_profile(hecMAT, hecLagMAT)
    implicit none
    type (hecmwst_matrix), intent(in) :: hecmat
    type (hecmwst_matrix_lagrange), intent(in) :: heclagmat
    integer(kind=kint) :: np, ndof, nlag, ntdf, npl, npu
    integer(kind=kint) :: i, j, k, l, cl, cu
 
    np = hecmat%NP
    ndof = hecmat%NDOF
    nlag = heclagmat%num_lagrange
    ntdf = np*ndof + nlag
    npl = hecmat%NPL*ndof*ndof + np*ndof*(ndof-1)/2
    npu = hecmat%NPU*ndof*ndof + np*ndof*(ndof-1)/2
    if (nlag > 0) then
      npl = npl + heclagmat%numL_lagrange*ndof
      npu = npu + heclagmat%numU_lagrange*ndof
    endif
 
    call hecmw_mat_finalize(mfmat)
    mfmat%N = ntdf
    mfmat%NP = ntdf
    mfmat%NDOF = 1
    mfmat%NPL = npl
    mfmat%NPU = npu
    allocate(mfmat%indexL(0:ntdf), mfmat%indexU(0:ntdf))
    allocate(mfmat%itemL(npl), mfmat%itemU(npu))
    allocate(mfmat%D(ntdf), mfmat%AL(npl), mfmat%AU(npu))
    allocate(mfmat%B(ntdf), mfmat%X(ntdf))
 
    cl = 0
    cu = 0
    mfmat%indexL(0) = 0
    mfmat%indexU(0) = 0
    do i = 1, np
      do j = 1, ndof
        do l = hecmat%indexL(i-1)+1, hecmat%indexL(i)
          do k = 1, ndof
            cl = cl + 1
            mfmat%itemL(cl) = (hecmat%itemL(l)-1)*ndof + k
          enddo
        enddo
        do k = 1, j-1
          cl = cl + 1
          mfmat%itemL(cl) = (i-1)*ndof + k
        enddo
        do k = j+1, ndof
          cu = cu + 1
          mfmat%itemU(cu) = (i-1)*ndof + k
        enddo
        do l = hecmat%indexU(i-1)+1, hecmat%indexU(i)
          do k = 1, ndof
            cu = cu + 1
            mfmat%itemU(cu) = (hecmat%itemU(l)-1)*ndof + k
          enddo
        enddo
        if (nlag > 0) then
          do l = heclagmat%indexU_lagrange(i-1)+1, heclagmat%indexU_lagrange(i)
            cu = cu + 1
            mfmat%itemU(cu) = np*ndof + heclagmat%itemU_lagrange(l)
          enddo
        endif
        mfmat%indexL((i-1)*ndof+j) = cl
        mfmat%indexU((i-1)*ndof+j) = cu
      enddo
    enddo
    do i = 1, nlag
      do l = heclagmat%indexL_lagrange(i-1)+1, heclagmat%indexL_lagrange(i)
        do k = 1, ndof
          cl = cl + 1
          mfmat%itemL(cl) = (heclagmat%itemL_lagrange(l)-1)*ndof + k
        enddo
      enddo
      mfmat%indexL(np*ndof+i) = cl
      mfmat%indexU(np*ndof+i) = cu
    enddo
 
    mfmat%Iarray = hecmat%Iarray
    mfmat%Rarray = hecmat%Rarray
    if (mfmat%Iarray(43) /= 0) then
      ! the user asked for BLR, so the fallback is reported regardless of LOGLEVEL
      if (hecmw_comm_get_rank() == 0) write(*,*) &
        '[DIRECTmf]: BLR is not available in contact analysis without elimination; disabled'
    endif
    ! the refinement stays with the caller (mf_contact_resid); the driver would measure it
    ! through the mesh
    mfmat%Iarray(43) = 0
    mfmat%Iarray(44) = -1
    mfmat%Iarray(48) = 0
    mfmat%Iarray(98) = 1
  end subroutine mf_contact_set_profile
 
  subroutine mf_contact_set_values(hecMAT, hecLagMAT)
    implicit none
    type (hecmwst_matrix), intent(in) :: hecmat
    type (hecmwst_matrix_lagrange), intent(in) :: heclagmat
    integer(kind=kint) :: np, ndof, nlag
    integer(kind=kint) :: i, j, k, l, cl, cu
 
    np = hecmat%NP
    ndof = hecmat%NDOF
    nlag = heclagmat%num_lagrange
 
    cl = 0
    cu = 0
    do i = 1, np
      do j = 1, ndof
        do l = hecmat%indexL(i-1)+1, hecmat%indexL(i)
          do k = 1, ndof
            cl = cl + 1
            mfmat%AL(cl) = hecmat%AL(((l-1)*ndof+j-1)*ndof+k)
          enddo
        enddo
        do k = 1, j-1
          cl = cl + 1
          mfmat%AL(cl) = hecmat%D(((i-1)*ndof+j-1)*ndof+k)
        enddo
        mfmat%D((i-1)*ndof+j) = hecmat%D(((i-1)*ndof+j-1)*ndof+j)
        do k = j+1, ndof
          cu = cu + 1
          mfmat%AU(cu) = hecmat%D(((i-1)*ndof+j-1)*ndof+k)
        enddo
        do l = hecmat%indexU(i-1)+1, hecmat%indexU(i)
          do k = 1, ndof
            cu = cu + 1
            mfmat%AU(cu) = hecmat%AU(((l-1)*ndof+j-1)*ndof+k)
          enddo
        enddo
        if (nlag > 0) then
          do l = heclagmat%indexU_lagrange(i-1)+1, heclagmat%indexU_lagrange(i)
            cu = cu + 1
            mfmat%AU(cu) = heclagmat%AU_lagrange((l-1)*ndof+j)
          enddo
        endif
      enddo
    enddo
    do i = 1, nlag
      mfmat%D(np*ndof+i) = 0.0d0
      do l = heclagmat%indexL_lagrange(i-1)+1, heclagmat%indexL_lagrange(i)
        do k = 1, ndof
          cl = cl + 1
          mfmat%AL(cl) = heclagmat%AL_lagrange((l-1)*ndof+k)
        enddo
      enddo
    enddo
  end subroutine mf_contact_set_values
 
  subroutine mf_contact_resid(ntdf, x, b, r)
    implicit none
    integer(kind=kint), intent(in) :: ntdf
    real(kind=kreal), intent(in) :: x(:), b(:)
    real(kind=kreal), intent(out) :: r(:)
    integer(kind=kint) :: i, l
    real(kind=kreal) :: s
 
    do i = 1, ntdf
      s = mfmat%D(i)*x(i)
      do l = mfmat%indexL(i-1)+1, mfmat%indexL(i)
        s = s + mfmat%AL(l)*x(mfmat%itemL(l))
      enddo
      do l = mfmat%indexU(i-1)+1, mfmat%indexU(i)
        s = s + mfmat%AU(l)*x(mfmat%itemU(l))
      enddo
      r(i) = b(i) - s
    enddo
  end subroutine mf_contact_resid
 
end module m_solve_lineq_mf_contact