hecmw_solver_contact.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module m_solve_lineq_contact
 
  use hecmw_util
  use m_solve_lineq_mkl_contact
  use m_solve_lineq_direct_serial_lag
  use m_solve_lineq_mumps_contact
  use m_solve_lineq_mf_contact
  use m_solve_lineq_contact_elim
  use hecmw_matrix_misc
  use m_hecmw_comm_f
  use hecmw_ebc_defer
 
  implicit none
 
  private
  public :: solve_lineq_contact_init
  public :: solve_lineq_contact
 
contains
 
  subroutine solve_lineq_contact_init(hecMESH,hecMAT,hecLagMAT,is_sym)
    type (hecmwst_local_mesh)                :: hecmesh
    type (hecmwst_matrix)                    :: hecmat
    type (hecmwst_matrix_lagrange)           :: heclagmat
    logical :: is_sym
    integer(kind=kint) :: contact_elim, solver_type, loglevel, myrank
    logical :: writelog
 
    contact_elim = hecmw_mat_get_contact_elim(hecmat)
    solver_type = hecmw_mat_get_solver_type(hecmat)
    loglevel = max(0, hecmw_mat_get_loglevel(hecmat)) ! LOGLEVEL is -1 when not given
    myrank = hecmw_comm_get_rank()
    writelog = (myrank==0 .and. loglevel>0)
 
    if( contact_elim==0 )then ! auto
      if( solver_type==1 )then ! iterative
        contact_elim=1
      else ! direct
        contact_elim=-1
      endif
      call hecmw_mat_set_contact_elim(hecmat,contact_elim)
    endif
 
    ! without a Lagrange multiplier anywhere (ALAGRANGE, or SLAGRANGE before any activation)
    ! there is no saddle point part to keep, so the single-process limitation of the
    ! no-elimination DIRECTmf path does not apply: the elimination path solves the plain
    ! distributed matrix as is. The stored CONTACT_ELIM is not changed, so the limitation
    ! still guards the solves after multipliers appear.
    if( contact_elim/=1 .and. solver_type==7 .and. hecmw_comm_get_size() > 1 )then
      if( .not. any_lagrange(hecmesh, heclagmat) ) contact_elim = 1
    endif
 
    if( contact_elim==1 )then
      if( writelog ) write(*,*) 'solve contact with elimination'
      call solve_lineq_contact_elim_init(hecmesh,hecmat,heclagmat,is_sym)
    else
      if( writelog ) write(*,*) 'solve contact without elimination'
      if( solver_type==1 )then
        write(*,*) 'ERROR: iterative solver without elimination not available in contact analysis'
        call hecmw_abort(hecmw_comm_get_comm())
      elseif( solver_type==2 )then
        call solve_lineq_serial_lag_hecmw_init(hecmat,heclagmat,is_sym)
      else if( solver_type==3 )then
        call solve_lineq_mkl_contact_init(hecmesh,is_sym)
      elseif( solver_type==5 ) then
        call solve_lineq_mumps_contact_init(hecmesh,hecmat,heclagmat,is_sym)
      elseif( solver_type==7 ) then
        if( hecmw_comm_get_size() > 1) then
          write(*,*) 'ERROR: !SOLVER,METHOD=DIRECTmf not available in parallel contact analysis',&
              ' without elimination; specify CONTACT_ELIM=1 or use MUMPS or DIRECTmkl instead'
          call hecmw_abort(hecmw_comm_get_comm())
        endif
        call solve_lineq_mf_contact_init(hecmesh,hecmat,heclagmat,is_sym)
      else
        write(*,*) 'ERROR: specified solver not available in contact analysis without elimination;',&
            ' please use MUMPS or DIRECTmkl instead'
        call hecmw_abort(hecmw_comm_get_comm())
      endif
    endif
  end subroutine solve_lineq_contact_init
 
 
  subroutine solve_lineq_contact(hecMESH,hecMAT,hecLagMAT,conMAT,hecEBC,istat,rf,is_contact_active)
 
    type (hecmwst_local_mesh)                :: hecmesh
    type (hecmwst_matrix)                    :: hecmat
    type (hecmwst_matrix_lagrange)           :: heclagmat
    type (hecmwst_matrix)                    :: conmat
    type (hecmwst_ebc)                       :: hecebc
    integer(kind=kint), intent(out)          :: istat
    real(kind=kreal), optional               :: rf
    logical                                  :: is_contact_active
 
    real(kind=kreal)                         :: factor
    real(kind=kreal) :: t1, t2
    integer(kind=kint) :: ndof
    integer(kind=kint) :: contact_elim, solver_type
 
    contact_elim = hecmw_mat_get_contact_elim(hecmat)
    solver_type = hecmw_mat_get_solver_type(hecmat)
 
    ! the same bypass as the init: no multiplier anywhere means no saddle point part
    if( contact_elim/=1 .and. solver_type==7 .and. hecmw_comm_get_size() > 1 )then
      if( .not. any_lagrange(hecmesh, heclagmat) ) contact_elim = 1
    endif
 
    factor = 1.0d0
    if( present(rf) )factor = rf
 
    t1 = hecmw_wtime()
 
    istat = 0
    ! the direct-solver branches below report failure through istat only and never
    ! touch these flags; clear them so a stale verdict from an earlier solve is not
    ! mistaken for this one's. The elimination branch overwrites them with the real result.
    call hecmw_mat_set_flag_converged(hecmat, 1)
    call hecmw_mat_set_flag_diverged(hecmat, 0)
    if( contact_elim==1 )then
      call solve_lineq_contact_elim(hecmesh,hecmat,heclagmat,hecebc,istat,conmat,is_contact_active)
    else
      if( solver_type==1 )then
        write(*,*) 'ERROR: iterative solver without elimination not available in contact analysis'
        call hecmw_abort(hecmw_comm_get_comm())
      elseif( solver_type==2 )then
        if( hecmw_comm_get_size() > 1) then
          write(*,*) 'ERROR: !SOLVER,METHOD=DIRECT not available in parallel contact analysis;',&
              ' please use MUMPS or DIRECTmkl instead'
          call hecmw_abort(hecmw_comm_get_comm())
        else
          call add_conmat_to_hecmat(hecmat,conmat,heclagmat)
          call solve_lineq_serial_lag_hecmw(hecmesh,hecmat,heclagmat,hecebc)
        endif
      elseif( solver_type==3 )then
        if( hecmw_comm_get_size() > 1) then
          call solve_lineq_mkl_contact(hecmesh,hecmat,heclagmat,hecebc,istat,conmat)
        else
          call add_conmat_to_hecmat(hecmat,conmat,heclagmat)
          call solve_lineq_mkl_contact(hecmesh,hecmat,heclagmat,hecebc,istat)
        endif
      elseif( solver_type==5 ) then
        call solve_lineq_mumps_contact(hecmesh,hecmat,heclagmat,hecebc,istat,conmat)
      elseif( solver_type==7 ) then
        if( hecmw_comm_get_size() > 1) then
          write(*,*) 'ERROR: !SOLVER,METHOD=DIRECTmf not available in parallel contact analysis',&
              ' without elimination; specify CONTACT_ELIM=1 or use MUMPS or DIRECTmkl instead'
          call hecmw_abort(hecmw_comm_get_comm())
        else
          call add_conmat_to_hecmat(hecmat,conmat,heclagmat)
          call solve_lineq_mf_contact(hecmesh,hecmat,heclagmat,hecebc,istat)
        endif
      else
        write(*,*) 'ERROR: specified solver not available in contact analysis without elimination;',&
            ' please use MUMPS or DIRECTmkl instead'
        call hecmw_abort(hecmw_comm_get_comm())
      endif
    endif
 
    ndof = hecmat%NDOF
    call hecmw_update_r(hecmesh,hecmat%X,hecmesh%n_node, ndof)
 
    t2 = hecmw_wtime()
    if (hecmw_mat_get_timelog(hecmat) .ge. 1) then
      if ( hecmw_comm_get_rank() ==0) write(*,*) ' solve time :', t2 - t1
    endif
 
    hecmat%X=factor*hecmat%X
 
  end subroutine solve_lineq_contact
 
 
  logical function any_lagrange(hecMESH, hecLagMAT)
    type (hecmwst_local_mesh)      :: hecmesh
    type (hecmwst_matrix_lagrange) :: heclagmat
    integer(kind=kint) :: n
 
    n = heclagmat%num_lagrange
    call hecmw_allreduce_i1(hecmesh, n, hecmw_max)
    any_lagrange = n > 0
  end function any_lagrange
 
 
  subroutine add_conmat_to_hecmat(hecMAT,conMAT,hecLagMat)
    type (hecmwst_matrix),          intent(inout) :: hecmat
    type (hecmwst_matrix),          intent(in)    :: conmat
    type (hecmwst_matrix_lagrange), intent(in)    :: heclagmat
 
 
    integer(kind=kint) :: ndof,ndof2,i
 
    ndof = hecmat%NDOF
    ndof2 = ndof*ndof
 
    do i=1,hecmat%NP*ndof + heclagmat%num_lagrange
      hecmat%B(i) = hecmat%B(i) + conmat%B(i)
    enddo
 
    do i=1,hecmat%NP*ndof2
      hecmat%D(i) = hecmat%D(i) + conmat%D(i)
    enddo
 
    do i=1,hecmat%NPL*ndof2
      hecmat%AL(i) = hecmat%AL(i) + conmat%AL(i)
    enddo
 
    do i=1,hecmat%NPU*ndof2
      hecmat%AU(i) = hecmat%AU(i) + conmat%AU(i)
    enddo
  end subroutine add_conmat_to_hecmat
 
 
end module m_solve_lineq_contact