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_convergence_measure
  public :: fstr_convergence_state
  public :: fstr_check_convergence
  public :: fstr_check_convergence_main
  public :: fstr_check_linear_solver
 
  type fstr_convergence_measure
    real(kind=kreal)   :: value = 0.0d0      
    real(kind=kreal)   :: tol   = 0.0d0      
    logical            :: check = .false.    
    logical            :: ok    = .false.    
    logical            :: absolute = .false. 
    integer(kind=kint) :: node  = 0
    integer(kind=kint) :: dof   = 0
  end type fstr_convergence_measure
 
  type fstr_convergence_state
    type(fstr_convergence_measure) :: m(2,3,2)
  end type fstr_convergence_state
 
  real(kind=kreal), parameter :: fref_floor = 1.0d-8
 
contains
 
  subroutine fstr_check_convergence( hecMESH, hecMAT, fstrSOLID, fstrPR, ndof, iter, sub_step, cstep, &
      residual_vec, cnvstat, iterStatus, hecLagMAT )
    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(:)
    type(fstr_convergence_state), intent(out) :: cnvstat
    integer(kind=kint), intent(out)           :: iterstatus
    type(hecmwst_matrix_lagrange), intent(in), optional :: heclagmat
 
    real(kind=kreal)    :: res_for_check
    logical             :: do_failure_check
 
    ! --- core convergence check (customizable) ---
    call fstr_check_convergence_main( hecmesh, hecmat, fstrsolid, fstrpr, ndof, iter, cstep, &
        residual_vec, cnvstat, iterstatus, do_failure_check, res_for_check, heclagmat )
 
    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
      write(ilog,'(a,i5,a,i5)') '### Fail to Converge  : at total_step=', cstep, '  sub_step=', sub_step
      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, cnvstat, iterStatus, do_failure_check, res_for_check, hecLagMAT )
    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(:)
    type(fstr_convergence_state), intent(out) :: cnvstat
    integer(kind=kint), intent(out)           :: iterstatus
    logical, intent(out)                      :: do_failure_check
    real(kind=kreal), intent(out)             :: res_for_check
    type(hecmwst_matrix_lagrange), intent(in), optional :: heclagmat
 
    logical :: has_dx
 
    iterstatus = kitrcontinue
    do_failure_check = .false.
    res_for_check = 0.0d0
 
    has_dx = .not. ( fstrpr%solution_type == kstdynamic .and. iter == 1 )
 
    call fstr_evaluate_convergence( hecmesh, hecmat, fstrsolid, ndof, cstep, has_dx, residual_vec, cnvstat, heclagmat )
 
    if( hecmesh%my_rank == 0 ) call fstr_print_convergence_state( iter, cnvstat )
 
    if( .not. has_dx ) return
 
    if( fstr_decide_convergence( cnvstat ) ) then
      iterstatus = kitrconverged
      return
    endif
 
    do_failure_check = .true.
    res_for_check = max( cnvstat%m(kcnvresidual, kcnvtranslation, kcnvl2)%value, &
                         cnvstat%m(kcnvresidual, kcnvrotation, kcnvl2)%value )
 
  end subroutine fstr_check_convergence_main
 
  logical function fstr_decide_convergence( cnvstat )
    implicit none
    type(fstr_convergence_state), intent(in) :: cnvstat
 
    type(fstr_convergence_measure) :: lag(2), dx(4), res(4)
 
    lag = [ cnvstat%m(kcnvcorrection, kcnvlagrange, :) ]
    dx  = [ cnvstat%m(kcnvcorrection, kcnvtranslation:kcnvrotation, :) ]
    res = [ cnvstat%m(kcnvresidual,   kcnvtranslation:kcnvrotation, :) ]
 
    fstr_decide_convergence = cnv_all_ok(lag) .and. &
      ( ( cnv_any_check(dx)  .and. cnv_all_ok(dx)  ) .or. &
        ( cnv_any_check(res) .and. cnv_all_ok(res) ) )
 
  end function fstr_decide_convergence
 
  logical function cnv_all_ok( ms )
    implicit none
    type(fstr_convergence_measure), intent(in) :: ms(:)
 
    cnv_all_ok = all( ms%ok .or. .not. ms%check )
 
  end function cnv_all_ok
 
  logical function cnv_any_check( ms )
    implicit none
    type(fstr_convergence_measure), intent(in) :: ms(:)
 
    cnv_any_check = any( ms%check )
 
  end function cnv_any_check
 
  subroutine fstr_evaluate_convergence( hecMESH, hecMAT, fstrSOLID, ndof, cstep, has_dx, residual_vec, cnvstat, &
      hecLagMAT )
    implicit none
 
    type(hecmwst_local_mesh), intent(in)      :: hecmesh
    type(hecmwst_matrix), intent(in)          :: hecmat
    type(fstr_solid), intent(in)              :: fstrsolid
    integer(kind=kint), intent(in)            :: ndof
    integer(kind=kint), intent(in)            :: cstep
    logical, intent(in)                       :: has_dx
    real(kind=kreal), intent(in)              :: residual_vec(:)
    type(fstr_convergence_state), intent(out) :: cnvstat
    type(hecmwst_matrix_lagrange), intent(in), optional :: heclagmat
 
    integer(kind=kint), parameter :: nsum = 13
    real(kind=kreal)   :: sq(nsum)
    integer(kind=kint) :: i, npndof, num_lagrange
    real(kind=kreal)   :: dx_t, dx_r, fref_t
    logical            :: has_rot, has_lag
 
    num_lagrange = 0
    if( present(heclagmat) ) num_lagrange = heclagmat%num_lagrange
    npndof = hecmat%NP*ndof
 
    has_rot = ( ndof == 6 )
 
    ! --- squares of the norms, reduced in a single collective ---
    sq(:) = 0.0d0
    call fstr_get_sqnorm_dofgroup( hecmesh, ndof, residual_vec,      sq(1), sq(2) )
    call fstr_get_sqnorm_dofgroup( hecmesh, ndof, fstrsolid%QFORCE,  sq(3), sq(4) )
    call fstr_get_sqnorm_dofgroup( hecmesh, ndof, fstrsolid%DFORCE,  sq(5), sq(6) )
    if( has_dx ) then
      call fstr_get_sqnorm_dofgroup( hecmesh, ndof, hecmat%X,         sq(7), sq(8) )
      call fstr_get_sqnorm_dofgroup( hecmesh, ndof, fstrsolid%dunode, sq(9), sq(10) )
      do i = 1, num_lagrange
        sq(11) = sq(11) + hecmat%X(npndof+i)*hecmat%X(npndof+i)
        sq(12) = sq(12) + heclagmat%Lagrange(i)*heclagmat%Lagrange(i)
      enddo
    endif
    sq(13) = dble(num_lagrange)
    call hecmw_allreduce_r( hecmesh, sq, nsum, hecmw_sum )
 
    ! Lagrange rows may be absent from a subdomain while present in the model, so the criteria must be selected from
    ! the reduced count.
    has_lag = ( sq(13) > 0.5d0 )
 
    cnvstat%m(:,:,:)%tol = fstrsolid%step_ctrl(cstep)%cnv_tol(:,:,:)
    cnvstat%m(:,:,:)%check = fstrsolid%step_ctrl(cstep)%cnv_check(:,:,:)
    if( .not. has_rot ) cnvstat%m(:,kcnvrotation,:)%check = .false.
    if( .not. has_lag ) cnvstat%m(:,kcnvlagrange,:)%check = .false.
    if( .not. has_dx ) cnvstat%m(kcnvcorrection,:,:)%check = .false.
 
    ! --- residual relative to the reference force ---
    fref_t = max( sqrt(sq(3)), sqrt(sq(5)) )
    call fstr_set_relative_measure( cnvstat%m(kcnvresidual,kcnvtranslation,kcnvl2), &
        sqrt(sq(1)), fref_t )
    if( has_rot ) then
      call fstr_set_relative_measure( cnvstat%m(kcnvresidual,kcnvrotation,kcnvl2), &
          sqrt(sq(2)), max( sqrt(sq(4)), sqrt(sq(6)) ) )
    endif
 
    ! --- correction relative to the accumulated increment ---
    ! A vanishing denominator leaves the ratio at 1, i.e. not converged.
    if( has_dx ) then
      dx_t = 1.0d0
      if( sq(9) > 0.0d0 ) dx_t = sqrt( sq(7)/sq(9) )
      dx_r = 1.0d0
      if( sq(10) > 0.0d0 ) dx_r = sqrt( sq(8)/sq(10) )
      cnvstat%m(kcnvcorrection,kcnvtranslation,kcnvl2)%value = dx_t
      if( has_rot ) cnvstat%m(kcnvcorrection,kcnvrotation,kcnvl2)%value = dx_r
      ! The multiplier is a contact force, so the reference force bounds its norm from below.
      if( has_lag ) then
        call fstr_set_relative_measure( cnvstat%m(kcnvcorrection,kcnvlagrange,kcnvl2), &
            sqrt(sq(11)), max( sqrt(sq(12)), fref_t ) )
      endif
    endif
 
    ! The flags are the same on every subdomain, so the collectives of the max norms are skipped consistently.
    if( any( cnvstat%m(:,:,kcnvmax)%check ) ) then
      call fstr_evaluate_convergence_max( hecmesh, hecmat, fstrsolid, ndof, has_dx, residual_vec, cnvstat, &
          heclagmat )
    endif
 
    ! A NaN component fails every comparison above: it never wins a maximum and never passes the test on a denominator,
    ! so a ratio would keep its default. The sums are the only place it survives, and every criterion takes it over.
    do i = 1, nsum-1
      if( sq(i) /= sq(i) ) cnvstat%m(:,:,:)%value = sq(i)
    enddo
 
    cnvstat%m(:,:,:)%ok = cnvstat%m(:,:,:)%check .and. ( cnvstat%m(:,:,:)%value < cnvstat%m(:,:,:)%tol )
 
  end subroutine fstr_evaluate_convergence
 
  subroutine fstr_evaluate_convergence_max( hecMESH, hecMAT, fstrSOLID, ndof, has_dx, residual_vec, cnvstat, &
      hecLagMAT )
    implicit none
 
    type(hecmwst_local_mesh), intent(in)        :: hecmesh
    type(hecmwst_matrix), intent(in)            :: hecmat
    type(fstr_solid), intent(in)                :: fstrsolid
    integer(kind=kint), intent(in)              :: ndof
    logical, intent(in)                         :: has_dx
    real(kind=kreal), intent(in)                :: residual_vec(:)
    type(fstr_convergence_state), intent(inout) :: cnvstat
    type(hecmwst_matrix_lagrange), intent(in), optional :: heclagmat
 
    integer(kind=kint), parameter :: nmax = 12
    ! entries of vmax located, and the criteria they belong to
    integer(kind=kint), parameter :: nloc = 4
    integer(kind=kint), parameter :: imax_loc(nloc) = [ 1, 2, 7, 8 ]
    integer(kind=kint), parameter :: iq_loc(nloc) = [ kcnvresidual, kcnvresidual, kcnvcorrection, kcnvcorrection ]
    integer(kind=kint), parameter :: ig_loc(nloc) = [ kcnvtranslation, kcnvrotation, kcnvtranslation, kcnvrotation ]
    real(kind=kreal)   :: vmax(nmax), gmax(nmax)
    integer(kind=kint) :: loc(2,nloc), gnode(nloc), gdof(nloc)
    integer(kind=kint) :: i, k, npndof, num_lagrange
    real(kind=kreal)   :: dx_t, dx_r, fref_t
 
    num_lagrange = 0
    if( present(heclagmat) ) num_lagrange = heclagmat%num_lagrange
    npndof = hecmat%NP*ndof
 
    vmax(:) = 0.0d0
    loc(:,:) = 0
    call fstr_get_maxabs_dofgroup( hecmesh, ndof, residual_vec,      vmax(1), vmax(2), loc(:,1), loc(:,2) )
    call fstr_get_maxabs_dofgroup( hecmesh, ndof, fstrsolid%QFORCE,  vmax(3), vmax(4) )
    call fstr_get_maxabs_dofgroup( hecmesh, ndof, fstrsolid%DFORCE,  vmax(5), vmax(6) )
    if( has_dx ) then
      call fstr_get_maxabs_dofgroup( hecmesh, ndof, hecmat%X,         vmax(7), vmax(8), loc(:,3), loc(:,4) )
      call fstr_get_maxabs_dofgroup( hecmesh, ndof, fstrsolid%dunode, vmax(9), vmax(10) )
      do i = 1, num_lagrange
        vmax(11) = max( vmax(11), abs(hecmat%X(npndof+i)) )
        vmax(12) = max( vmax(12), abs(heclagmat%Lagrange(i)) )
      enddo
    endif
    gmax(:) = vmax(:)
    call hecmw_allreduce_r( hecmesh, gmax, nmax, hecmw_max )
 
    do k = 1, nloc
      if( gmax(imax_loc(k)) > vmax(imax_loc(k)) ) loc(1,k) = 0
    enddo
    gnode(:) = loc(1,:)
    call hecmw_allreduce_i( hecmesh, gnode, nloc, hecmw_max )
    do k = 1, nloc
      gdof(k) = 0
      if( loc(1,k) == gnode(k) ) gdof(k) = loc(2,k)
    enddo
    call hecmw_allreduce_i( hecmesh, gdof, nloc, hecmw_max )
 
    ! --- residual relative to the reference force ---
    fref_t = max( gmax(3), gmax(5) )
    call fstr_set_relative_measure( cnvstat%m(kcnvresidual,kcnvtranslation,kcnvmax), &
        gmax(1), fref_t )
    if( ndof == 6 ) then
      call fstr_set_relative_measure( cnvstat%m(kcnvresidual,kcnvrotation,kcnvmax), &
          gmax(2), max( gmax(4), gmax(6) ) )
    endif
 
    ! --- correction relative to the accumulated increment, with the same rules as the L2 norms ---
    if( has_dx ) then
      dx_t = 1.0d0
      if( gmax(9) > 0.0d0 ) dx_t = gmax(7) / gmax(9)
      dx_r = 1.0d0
      if( gmax(10) > 0.0d0 ) dx_r = gmax(8) / gmax(10)
      cnvstat%m(kcnvcorrection,kcnvtranslation,kcnvmax)%value = dx_t
      if( ndof == 6 ) cnvstat%m(kcnvcorrection,kcnvrotation,kcnvmax)%value = dx_r
      call fstr_set_relative_measure( cnvstat%m(kcnvcorrection,kcnvlagrange,kcnvmax), &
          gmax(11), max( gmax(12), fref_t ) )
    endif
 
    do k = 1, nloc
      cnvstat%m(iq_loc(k),ig_loc(k),kcnvmax)%node = gnode(k)
      cnvstat%m(iq_loc(k),ig_loc(k),kcnvmax)%dof = gdof(k)
    enddo
 
  end subroutine fstr_evaluate_convergence_max
 
  subroutine fstr_set_relative_measure( ms, val, fref )
    implicit none
    type(fstr_convergence_measure), intent(inout) :: ms
    real(kind=kreal), intent(in)                  :: val   
    real(kind=kreal), intent(in)                  :: fref  
 
    ms%absolute = ( fref < fref_floor )
    if( ms%absolute ) then
      ms%value = val
    else
      ms%value = val / fref
    endif
 
  end subroutine fstr_set_relative_measure
 
  subroutine fstr_get_sqnorm_dofgroup( hecMESH, ndof, vec, sq_t, sq_r )
    implicit none
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in)       :: ndof
    real(kind=kreal), intent(in)         :: vec(:)
    real(kind=kreal), intent(out)        :: sq_t
    real(kind=kreal), intent(out)        :: sq_r
 
    integer(kind=kint) :: i, idof, idx
 
    sq_t = 0.0d0
    sq_r = 0.0d0
    do i = 1, hecmesh%nn_internal
      idx = ndof*(i-1)
      do idof = 1, min(ndof,3)
        sq_t = sq_t + vec(idx+idof)*vec(idx+idof)
      enddo
      if( ndof /= 6 ) cycle
      do idof = 4, 6
        sq_r = sq_r + vec(idx+idof)*vec(idx+idof)
      enddo
    enddo
 
  end subroutine fstr_get_sqnorm_dofgroup
 
  subroutine fstr_get_maxabs_dofgroup( hecMESH, ndof, vec, vmax_t, vmax_r, loc_t, loc_r )
    implicit none
    type(hecmwst_local_mesh), intent(in)      :: hecmesh
    integer(kind=kint), intent(in)            :: ndof
    real(kind=kreal), intent(in)              :: vec(:)
    real(kind=kreal), intent(out)             :: vmax_t
    real(kind=kreal), intent(out)             :: vmax_r
    integer(kind=kint), intent(out), optional :: loc_t(2)
    integer(kind=kint), intent(out), optional :: loc_r(2)
 
    integer(kind=kint) :: i, idof, idx, gid
    integer(kind=kint) :: lt(2), lr(2)
    real(kind=kreal)   :: a
 
    vmax_t = 0.0d0
    vmax_r = 0.0d0
    lt(:) = 0
    lr(:) = 0
    ! Ties between nodes go to the largest node ID, as in the reduction across subdomains, so that the location does
    ! not depend on the partitioning.
    do i = 1, hecmesh%nn_internal
      idx = ndof*(i-1)
      gid = hecmesh%global_node_ID(i)
      do idof = 1, min(ndof,3)
        a = abs(vec(idx+idof))
        if( a > vmax_t .or. ( a >= vmax_t .and. a > 0.0d0 .and. gid > lt(1) ) ) then
          vmax_t = a
          lt(1) = gid
          lt(2) = idof
        endif
      enddo
      if( ndof /= 6 ) cycle
      do idof = 4, 6
        a = abs(vec(idx+idof))
        if( a > vmax_r .or. ( a >= vmax_r .and. a > 0.0d0 .and. gid > lr(1) ) ) then
          vmax_r = a
          lr(1) = gid
          lr(2) = idof
        endif
      enddo
    enddo
    if( present(loc_t) ) loc_t(:) = lt(:)
    if( present(loc_r) ) loc_r(:) = lr(:)
 
  end subroutine fstr_get_maxabs_dofgroup
 
  subroutine fstr_print_convergence_state( iter, cnvstat )
    implicit none
    integer(kind=kint), intent(in)           :: iter
    type(fstr_convergence_state), intent(in) :: cnvstat
 
    character(len=512) :: line
    integer(kind=kint) :: iq, ig
 
    write(line,'(a,i8)') " iter:", iter
    do iq = kcnvresidual, kcnvcorrection
      do ig = kcnvtranslation, kcnvlagrange
        if( .not. cnvstat%m(iq,ig,kcnvl2)%check ) cycle
        write(line(len_trim(line)+1:),'(a,a,a,1pe11.4)') ", ", &
            trim(fstr_convergence_label( iq, ig, kcnvl2, cnvstat%m(iq,ig,kcnvl2)%absolute )), ":", &
            cnvstat%m(iq,ig,kcnvl2)%value
      enddo
    enddo
    write(*,'(a)') trim(line)
 
    line = ""
    do iq = kcnvresidual, kcnvcorrection
      do ig = kcnvtranslation, kcnvlagrange
        if( .not. cnvstat%m(iq,ig,kcnvmax)%check ) cycle
        if( len_trim(line) > 0 ) line = trim(line) // ","
        write(line(len_trim(line)+1:),'(a,a,a,1pe11.4)') " ", &
            trim(fstr_convergence_label( iq, ig, kcnvmax, cnvstat%m(iq,ig,kcnvmax)%absolute )), ":", &
            cnvstat%m(iq,ig,kcnvmax)%value
        if( cnvstat%m(iq,ig,kcnvmax)%node > 0 ) then
          write(line(len_trim(line)+1:),'(a,i0,a,i0,a)') " (node ", cnvstat%m(iq,ig,kcnvmax)%node, &
              ", dof ", cnvstat%m(iq,ig,kcnvmax)%dof, ")"
        endif
      enddo
    enddo
    if( len_trim(line) > 0 ) write(*,'(a,a)') repeat(" ", 15), trim(line)
 
  end subroutine fstr_print_convergence_state
 
  function fstr_convergence_label( iq, ig, inorm, absolute ) result( label )
    implicit none
    integer(kind=kint), intent(in) :: iq
    integer(kind=kint), intent(in) :: ig
    integer(kind=kint), intent(in) :: inorm
    logical, intent(in)            :: absolute
    character(len=32)              :: label
 
    character(len=5), parameter  :: res_name(3)  = [ 'force', 'mom  ', 'lag  ' ]
    character(len=10), parameter :: corr_name(3) = [ 'disp.corr.', 'rot.corr. ', 'lag.corr. ' ]
 
    if( iq == kcnvresidual ) then
      label = 'res(' // trim(res_name(ig))
      if( inorm == kcnvmax ) label = trim(label) // ',max'
      if( absolute ) label = trim(label) // ',abs'
      label = trim(label) // ')'
    else
      label = corr_name(ig)
      if( inorm == kcnvmax ) label = trim(label) // '(max)'
    endif
 
  end function fstr_convergence_label
 
  subroutine fstr_check_linear_solver( hecMESH, hecMAT, fstrSOLID, cstep, sub_step, iterStatus, istat )
    implicit none
 
    type(hecmwst_local_mesh), intent(in)     :: hecmesh
    type(hecmwst_matrix), intent(in)         :: hecmat
    type(fstr_solid), intent(inout)          :: fstrsolid
    integer(kind=kint), intent(in)           :: cstep
    integer(kind=kint), intent(in)           :: sub_step
    integer(kind=kint), intent(out)          :: iterstatus
    integer(kind=kint), intent(in), optional :: istat
 
    iterstatus = kitrcontinue
 
    ! direct solvers report a genuine failure only through istat
    if( present(istat) ) then
      if( istat /= 0 ) then
        iterstatus = kitrdiverged
        fstrsolid%NRstat_i(knstdresn) = 4
      endif
    endif
 
    if( iterstatus == kitrcontinue ) then
      if( hecmw_mat_get_flag_diverged(hecmat) /= kno ) then
        ! broke down: indefinite preconditioner/matrix, or NaN
        iterstatus = kitrdiverged
        fstrsolid%NRstat_i(knstdresn) = 6
      else if( hecmw_mat_get_solver_type(hecmat) == 1 .and. &
               hecmw_mat_get_flag_converged(hecmat) == kno ) then
        ! ran out of iterations. Restricted to iterative solvers because a direct solver leaves the same flag unset
        ! whenever its residual exceeds the tolerance, and hecmw_solve only warns and keeps going in that case
        iterstatus = kitrdiverged
        fstrsolid%NRstat_i(knstdresn) = 5
      endif
    endif
 
    if( iterstatus == kitrcontinue ) return
 
    if( hecmesh%my_rank == 0 ) then
      write(*,'(a,i5,a,i5)') '     ### Fail to Converge  : at total_step=', cstep, '  sub_step=', sub_step
    endif
    fstrsolid%CutBack_stat = fstrsolid%CutBack_stat + 1
 
  end subroutine fstr_check_linear_solver
 
end module m_fstr_iterationcontrol