fstr_CheckConvergence.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_fstr_iterationcontrol
  use m_fstr
  implicit none
 
  private
  public :: fstr_check_convergence
  public :: fstr_check_convergence_main
 
contains
 
  subroutine fstr_check_convergence( &
      hecMESH, hecMAT, fstrSOLID, fstrPR, &
      ndof, iter, sub_step, cstep, &
      residual_vec, nresid, &
      resb, res_prev, &
      n_node_global, &
      iterStatus, &
      maxDLag, converg_dlag )
    implicit none
 
    type(hecmwst_local_mesh), intent(in)    :: hecmesh
    type(hecmwst_matrix), intent(in)        :: hecmat
    type(fstr_solid), intent(inout)         :: fstrsolid
    type(fstr_param), intent(in)            :: fstrpr
    integer(kind=kint), intent(in)          :: ndof
    integer(kind=kint), intent(in)          :: iter
    integer(kind=kint), intent(in)          :: sub_step
    integer(kind=kint), intent(in)          :: cstep
    real(kind=kreal), intent(in)            :: residual_vec(:)
    integer(kind=kint), intent(in)          :: nresid
    real(kind=kreal), intent(inout)         :: resb
    real(kind=kreal), intent(inout)         :: res_prev
    integer(kind=kint), intent(in)          :: n_node_global
    integer(kind=kint), intent(out)         :: iterstatus
    real(kind=kreal), intent(in), optional  :: maxdlag
    real(kind=kreal), intent(in), optional  :: converg_dlag
 
    real(kind=kreal)    :: res_for_check
    logical             :: do_failure_check, is_dynamic, is_contact
 
    is_dynamic = (fstrpr%solution_type == kstdynamic)
    is_contact = (n_node_global > 0)
 
    ! --- core convergence check (customizable) ---
    call fstr_check_convergence_main( &
        hecmesh, hecmat, fstrsolid, fstrpr, &
        ndof, iter, cstep, &
        residual_vec, nresid, &
        resb, res_prev, &
        n_node_global, &
        iterstatus, &
        do_failure_check, res_for_check, &
        maxdlag, converg_dlag )
 
    if( iterstatus == kitrconverged ) return
    if( .not. do_failure_check ) return
 
    ! --- common divergence / NaN classification ---
    if( res_for_check /= res_for_check ) then
      iterstatus = kitrfloatingerror
    else if( iter == fstrsolid%step_ctrl(cstep)%max_iter .or. &
             res_for_check > fstrsolid%step_ctrl(cstep)%maxres ) then
      iterstatus = kitrdiverged
    endif
 
    if( iterstatus == kitrcontinue ) return
 
    ! --- common failure handling: log + stats update ---
    if( hecmesh%my_rank == 0 ) then
      ! Static non-contact path (legacy fstr_check_iteration_converged) also
      ! wrote to ILOG; preserve that behaviour for backward compatibility.
      if( .not. is_dynamic .and. .not. is_contact ) then
        write(ilog,'(a,i5,a,i5)') '### Fail to Converge  : at total_step=', cstep, '  sub_step=', sub_step
      endif
      write(   *,'(a,i5,a,i5)') '     ### Fail to Converge  : at total_step=', cstep, '  sub_step=', sub_step
    endif
    fstrsolid%NRstat_i(knstmaxit) = max(fstrsolid%NRstat_i(knstmaxit), iter)
    fstrsolid%NRstat_i(knstsumit) = fstrsolid%NRstat_i(knstsumit) + iter
    fstrsolid%CutBack_stat = fstrsolid%CutBack_stat + 1
    if( iterstatus == kitrdiverged .and. &
        iter == fstrsolid%step_ctrl(cstep)%max_iter ) then
      fstrsolid%NRstat_i(knstdresn) = 1
    else
      ! kitrDiverged due to maxres, or kitrFloatingError due to NaN
      fstrsolid%NRstat_i(knstdresn) = 2
    endif
 
  end subroutine fstr_check_convergence
 
  subroutine fstr_check_convergence_main( &
      hecMESH, hecMAT, fstrSOLID, fstrPR, &
      ndof, iter, cstep, &
      residual_vec, nresid, &
      resb, res_prev, &
      n_node_global, &
      iterStatus, &
      do_failure_check, res_for_check, &
      maxDLag, converg_dlag )
    implicit none
 
    type(hecmwst_local_mesh), intent(in)    :: hecmesh
    type(hecmwst_matrix), intent(in)        :: hecmat
    type(fstr_solid), intent(inout)         :: fstrsolid
    type(fstr_param), intent(in)            :: fstrpr
    integer(kind=kint), intent(in)          :: ndof
    integer(kind=kint), intent(in)          :: iter
    integer(kind=kint), intent(in)          :: cstep
    real(kind=kreal), intent(in)            :: residual_vec(:)
    integer(kind=kint), intent(in)          :: nresid
    real(kind=kreal), intent(inout)         :: resb
    real(kind=kreal), intent(inout)         :: res_prev
    integer(kind=kint), intent(in)          :: n_node_global
    integer(kind=kint), intent(out)         :: iterstatus
    logical, intent(out)                    :: do_failure_check
    real(kind=kreal), intent(out)           :: res_for_check
    real(kind=kreal), intent(in), optional  :: maxdlag
    real(kind=kreal), intent(in), optional  :: converg_dlag
 
    real(kind=kreal) :: res_sq, res_nrm, qnrm, rres, xnrm, dunrm, rxnrm
    real(kind=kreal) :: relres, res_normalized
    logical :: is_dynamic, is_contact
 
    iterstatus = kitrcontinue
    do_failure_check = .false.
    res_for_check = 0.0d0
 
    is_dynamic = (fstrpr%solution_type == kstdynamic)
    is_contact = (n_node_global > 0)
 
    ! =========================================================================
    ! PATH 1: Static, non-contact  (QFORCE normalization + disp correction)
    ! =========================================================================
    if( .not. is_dynamic .and. .not. is_contact ) then
      call hecmw_innerproduct_r(hecmesh, ndof, hecmat%B, hecmat%B, res_sq)
      res_nrm = sqrt(res_sq)
 
      call hecmw_innerproduct_r(hecmesh, ndof, hecmat%X, hecmat%X, xnrm)
      xnrm = sqrt(xnrm)
 
      call hecmw_innerproduct_r(hecmesh, ndof, fstrsolid%QFORCE, fstrsolid%QFORCE, qnrm)
      qnrm = sqrt(qnrm)
      if( qnrm < 1.0d-8 ) qnrm = 1.0d0
 
      if( iter == 1 ) then
        dunrm = xnrm
      else
        call hecmw_innerproduct_r(hecmesh, ndof, fstrsolid%dunode, fstrsolid%dunode, dunrm)
        dunrm = sqrt(dunrm)
      endif
 
      rres = res_nrm / qnrm
      rxnrm = xnrm / dunrm
 
      if( hecmesh%my_rank == 0 ) then
        if( qnrm == 1.0d0 ) then
          write(*,"(a,i8,a,1pe11.4,a,1pe11.4)") " iter:", iter, ", residual(abs):", rres, ", disp.corr.:", rxnrm
        else
          write(*,"(a,i8,a,1pe11.4,a,1pe11.4)") " iter:", iter, ", residual:", rres, ", disp.corr.:", rxnrm
        endif
      endif
 
      ! Convergence (skip if linear solver flagged divergence)
      if( hecmw_mat_get_flag_diverged(hecmat) == kno ) then
        if( rres < fstrsolid%step_ctrl(cstep)%converg .or. &
            rxnrm < fstrsolid%step_ctrl(cstep)%converg_ddisp ) then
          iterstatus = kitrconverged
          return
        endif
      endif
 
      do_failure_check = .true.
      res_for_check = rres
 
    ! =========================================================================
    ! PATH 2: Static, contact  (self-normalization by n_node_global)
    ! =========================================================================
    else if( .not. is_dynamic .and. is_contact ) then
      res_sq = dot_product(residual_vec(1:nresid), residual_vec(1:nresid))
      call hecmw_allreduce_r1(hecmesh, res_sq, hecmw_sum)
      res_normalized = sqrt(res_sq) / dble(n_node_global)
 
      if( iter == 1 ) resb = res_normalized
      if( resb == 0.0d0 ) then
        resb = 1.0d0
      else
        relres = dabs(res_prev - res_normalized) / resb
      endif
 
      if( hecmesh%my_rank == 0 ) then
        if( iter == 1 ) then
          write(*, '(a,i3,a,2e15.7)') ' - Residual(', iter, ') =', res_normalized, 0.0d0
        else
          write(*, '(a,i3,a,2e15.7)') ' - Residual(', iter, ') =', res_normalized, relres
        endif
      endif
 
      ! Convergence (iter==1: absolute only; iter>1: also relative)
      if( iter > 1 ) then
        if( res_normalized < fstrsolid%step_ctrl(cstep)%converg .or. &
            relres < fstrsolid%step_ctrl(cstep)%converg_ddisp ) then
          iterstatus = kitrconverged
          res_prev = res_normalized
          return
        endif
      else
        if( res_normalized < fstrsolid%step_ctrl(cstep)%converg ) then
          iterstatus = kitrconverged
          res_prev = res_normalized
          return
        endif
      endif
 
      res_prev = res_normalized
 
      do_failure_check = .true.
      res_for_check = res_normalized
 
    ! =========================================================================
    ! PATH 3: Dynamic, contact  (self-normalization + maxDLag condition)
    ! =========================================================================
    else if( is_dynamic .and. is_contact ) then
      res_sq = dot_product(residual_vec(1:nresid), residual_vec(1:nresid))
      call hecmw_allreduce_r1(hecmesh, res_sq, hecmw_sum)
 
      if( iter == 1 ) resb = res_sq
      res_normalized = dsqrt(res_sq / resb)
 
      ! Only check when nlgeom and ndof /= 4
      if( fstrpr%nlgeom .and. ndof /= 4 ) then
        if( hecmesh%my_rank == 0 ) then
          write(*,'(a,i5,a,1pe12.4)') "iter: ", iter, ", res: ", res_normalized
          write(ista,'(''iter='',I5,''- Residual'',E15.7)') iter, res_normalized
          if( present(maxdlag) ) then
            write(*,'(a,1e15.7)') ' - MaxDLag =', maxdlag
            write(ista,'(a,1e15.7)') ' - MaxDLag =', maxdlag
          endif
        endif
        if( present(maxdlag) .and. present(converg_dlag) ) then
          if( res_normalized < fstrsolid%step_ctrl(cstep)%converg .and. &
              maxdlag < converg_dlag ) then
            iterstatus = kitrconverged
            return
          endif
        else
          if( res_normalized < fstrsolid%step_ctrl(cstep)%converg ) then
            iterstatus = kitrconverged
            return
          endif
        endif
 
        do_failure_check = .true.
        res_for_check = res_normalized
      endif
 
    endif
 
  end subroutine fstr_check_convergence_main
 
end module m_fstr_iterationcontrol