hecmw_mat_resid_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_hecmw_mat_resid_contact
  use hecmw_util
  use hecmw_solver_las
  use hecmw_matrix_misc
  use hecmw_solver_misc
  use m_hecmw_comm_f
  implicit none
 
  private
  public :: hecmw_get_residual_contact
  public :: hecmw_get_rel_resid_l2_contact
  public :: hecmw_get_resid_max_contact
 
contains
 
  subroutine hecmw_get_residual_contact(hecMESH, hecMAT, hecLagMAT, r, rlag)
    implicit none
    type(hecmwst_local_mesh) :: hecmesh
    type(hecmwst_matrix)     :: hecmat
    type(hecmwst_matrix_lagrange), intent(in) :: heclagmat
    real(kind=kreal), intent(out) :: r(:)
    real(kind=kreal), intent(out) :: rlag(:)
    real(kind=kreal)   :: tcomm, sum
    integer(kind=kint) :: i, idof, ii, i0, ls, le, l, loc0, ll, j, loc, l0, k
    integer(kind=kint) :: ndof, npndof
    !! r    = b     - K    x  - Ulag lag
    !! rlag = blag  - Llag x
    ndof=hecmat%NDOF
    npndof=hecmat%NP*ndof
    do i=1,npndof
      r(i)=0.d0
    enddo
    ! r = b - K x
    tcomm = 0.d0
    call hecmw_matresid(hecmesh, hecmat, hecmat%X, hecmat%B, r, tcomm)
    if (heclagmat%num_lagrange > 0) then
      ! r = r - Ulag lag
      do i=1,hecmat%N
        i0=(i-1)*ndof
        do idof=1,ndof
          ii=i0+idof
          ls=heclagmat%indexU_lagrange(i-1)+1
          le=heclagmat%indexU_lagrange(i)
          sum=0.d0
          do l=ls,le
            loc=(l-1)*ndof+idof
            ll=heclagmat%itemU_lagrange(l)
            !sum=sum+hecLagMAT%AU_lagrange(loc)*hecLagMAT%Lagrange(ll)
            sum=sum+heclagmat%AU_lagrange(loc)*hecmat%X(npndof+ll)
          enddo
          r(ii)=r(ii)-sum
        enddo
      enddo
      ! rlag = blag - Llag x
      do i=1,heclagmat%num_lagrange
        ls=heclagmat%indexL_lagrange(i-1)+1
        le=heclagmat%indexL_lagrange(i)
        sum=0.d0
        do l=ls,le
          loc0=(l-1)*ndof
          j=heclagmat%itemL_lagrange(l)
          l0=(j-1)*ndof
          do k=1,ndof
            loc=loc0+k
            ll=l0+k
            sum=sum+heclagmat%AL_lagrange(loc)*hecmat%X(ll)
          enddo
        enddo
        rlag(i)=hecmat%B(npndof+i)-sum
      enddo
    end if
  end subroutine hecmw_get_residual_contact
 
  function hecmw_get_rel_resid_l2_contact(hecMESH, hecMAT, hecLagMAT)
    implicit none
    real(kind=kreal)         :: hecmw_get_rel_resid_l2_contact
    type(hecmwst_local_mesh) :: hecmesh
    type(hecmwst_matrix)     :: hecmat
    type(hecmwst_matrix_lagrange), intent(in) :: heclagmat
    real(kind=kreal), allocatable :: r(:)
    real(kind=kreal), allocatable :: rlag(:)
    real(kind=kreal)   :: bnorm2, rnorm2
    real(kind=kreal)   :: rlagnorm2
    real(kind=kreal)   :: tcomm
    integer(kind=kint) :: i
    allocate(r(hecmat%NDOF*hecmat%NP),rlag(heclagmat%num_lagrange))
    ! |b|^2
    call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, hecmat%B, hecmat%B, bnorm2, tcomm)
    call hecmw_get_residual_contact(hecmesh, hecmat, heclagmat, r, rlag)
    ! |r|^2
    call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, r, r, rnorm2, tcomm)
    ! |rlag|^2
    rlagnorm2=0.d0
    do i=1,heclagmat%num_lagrange
      rlagnorm2=rlagnorm2+rlag(i)*rlag(i)
    enddo
    deallocate(r,rlag)
    call hecmw_allreduce_r1(hecmesh, rlagnorm2, hecmw_sum)
    ! |r_total|^2 = |r|^2 + |rlag|^2
    if (heclagmat%num_lagrange > 0) rnorm2=rnorm2+rlagnorm2
    ! |r_total| / |b|
    hecmw_get_rel_resid_l2_contact = sqrt(rnorm2 / bnorm2)
  end function hecmw_get_rel_resid_l2_contact
 
  function hecmw_get_resid_max_contact(hecMESH, hecMAT, hecLagMAT)
    implicit none
    real(kind=kreal)         :: hecmw_get_resid_max_contact
    type(hecmwst_local_mesh) :: hecmesh
    type(hecmwst_matrix)     :: hecmat
    type(hecmwst_matrix_lagrange), intent(in) :: heclagmat
    real(kind=kreal), allocatable :: r(:)
    real(kind=kreal), allocatable :: rlag(:)
    real(kind=kreal) :: rmax, rlagmax
    real(kind=kreal) :: tcomm
    allocate(r(hecmat%NDOF*hecmat%NP),rlag(heclagmat%num_lagrange))
    call hecmw_get_residual_contact(hecmesh, hecmat, heclagmat, r, rlag)
    rmax = maxval(dabs(r))
    if (heclagmat%num_lagrange > 0) then
      rlagmax = maxval(dabs(rlag))
      if (rlagmax > rmax) rmax = rlagmax
    endif
    deallocate(r,rlag)
    call hecmw_allreduce_r1(hecmesh, rmax, hecmw_max)
    hecmw_get_resid_max_contact = rmax
  end function hecmw_get_resid_max_contact
 
end module m_hecmw_mat_resid_contact