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_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, iterlog, timelog, myrank
    logical :: writelog
 
    contact_elim = hecmw_mat_get_contact_elim(hecmat)
    solver_type = hecmw_mat_get_solver_type(hecmat)
    iterlog = hecmw_mat_get_iterlog(hecmat)
    timelog = hecmw_mat_get_timelog(hecmat)
    myrank = hecmw_comm_get_rank()
    writelog = (myrank==0 .and. (iterlog>0 .or. timelog>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
 
    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)
      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)
 
    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)
      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
 
 
  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