m_step.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_step
  use hecmw
  use fstr_ctrl_util_f
 
  implicit none
 
  integer, parameter :: stepstatic = 1
  integer, parameter :: stepvisco  = 2
  integer, parameter :: stepfixedinc = 1
  integer, parameter :: stepautoinc  = 2
  integer, parameter :: stepampramp = 1  ! prescribed disp. ramped over the step (default for static)
  integer, parameter :: stepampstep = 2  ! prescribed disp. applied in full at step start (default for visco)
 
  ! statistics of newton iteration
  integer(kind=kint),parameter :: knstmaxit  = 1 ! maximum number of newton iteration
  integer(kind=kint),parameter :: knstsumit  = 2 ! total number of newton iteration
  integer(kind=kint),parameter :: knstciter  = 3 ! number of contact iteration
  integer(kind=kint),parameter :: knstdresn  = 4 ! reason of not to converged
 
  ! indices of the convergence criteria: (quantity, DOF group, norm)
  integer(kind=kint),parameter :: kcnvresidual    = 1
  integer(kind=kint),parameter :: kcnvcorrection  = 2
  integer(kind=kint),parameter :: kcnvtranslation = 1
  integer(kind=kint),parameter :: kcnvrotation    = 2
  integer(kind=kint),parameter :: kcnvlagrange    = 3
  integer(kind=kint),parameter :: kcnvl2          = 1
  integer(kind=kint),parameter :: kcnvmax         = 2
 
  type step_info
    integer             :: solution
    integer             :: inc_type
    character( len=80 ) :: control
    character( len=80 ) :: convcontrol
    real(kind=kreal)    :: converg            
    real(kind=kreal)    :: converg_lag        
    real(kind=kreal)    :: converg_ddisp      
    real(kind=kreal)    :: maxres             
    logical             :: cnv_check(2,3,2)
    real(kind=kreal)    :: cnv_tol(2,3,2)     
 
    integer :: num_substep
    integer :: max_iter
    integer :: max_contiter
    integer :: amp_id
    integer :: amp_default_type
    real(kind=kreal) :: initdt                
    real(kind=kreal) :: elapsetime            
    real(kind=kreal) :: mindt                 
    real(kind=kreal) :: maxdt                 
    real(kind=kreal) :: starttime             
    integer, pointer :: boundary(:)=>null()   
    integer, pointer :: load(:)=>null()       
    integer, pointer :: contact(:)=>null()    
    integer, pointer :: elemactivation(:)=>null()      
    integer :: timepoint_id
    integer :: aincparam_id
  end type
 
  type tparamautoinc
    character(HECMW_NAME_LEN)  :: name
    real(kind=kreal)     :: ainc_rs            
    real(kind=kreal)     :: ainc_rl            
    integer( kind=kint ) :: nrbound_s(10)
    integer( kind=kint ) :: nrbound_l(10)
    integer( kind=kint ) :: nrtimes_s
    integer( kind=kint ) :: nrtimes_l
    real(kind=kreal)     :: ainc_rc            
    integer( kind=kint ) :: cbbound
  end type
 
  type tparamconverg
    character(HECMW_NAME_LEN) :: name
    logical              :: given(2,3,2)
    logical              :: check(2,3,2)
    real(kind=kreal)     :: tol(2,3,2)        
  end type
 
contains
 
  subroutine init_stepinfo( stepinfo )
    type( step_info ), intent(out) :: stepinfo
    stepinfo%solution = stepstatic
    stepinfo%inc_type = stepfixedinc
    stepinfo%num_substep = 1
    stepinfo%max_iter = 50
    stepinfo%max_contiter = 10
    stepinfo%amp_id = -1
    stepinfo%amp_default_type = stepampramp
    stepinfo%initdt = 1.d0
    stepinfo%mindt = 1.d-4
    stepinfo%maxdt = 1.d0
    stepinfo%elapsetime = 1.d0
    stepinfo%starttime = 0.d0
    stepinfo%converg = 1.d-3
    ! The multiplier criterion guards against a gap constraint left unsatisfied, which the force residual cannot
    ! see; the accuracy of the solution is carried by the residual criterion, so this one is left loose enough
    ! that it does not decide the iteration count on its own.
    stepinfo%converg_lag = 1.d-2
    stepinfo%converg_ddisp = 1.d-8
    stepinfo%maxres = 1.d+10
    stepinfo%timepoint_id = 0
    stepinfo%AincParam_id = 0
    call setup_stepinfo_converg( stepinfo )
  end subroutine
 
  subroutine setup_stepinfo_converg( stepinfo, cnvparam )
    type( step_info ), intent(inout)            :: stepinfo
    type( tparamconverg ), intent(in), optional :: cnvparam
 
    integer :: iq, ig, inorm
 
    stepinfo%cnv_check(:,:,:) = .false.
    stepinfo%cnv_tol(:,:,:) = 0.d0
    stepinfo%cnv_check(kcnvresidual, kcnvtranslation:kcnvrotation, kcnvl2) = .true.
    stepinfo%cnv_tol(kcnvresidual, kcnvtranslation:kcnvrotation, kcnvl2) = stepinfo%converg
    stepinfo%cnv_check(kcnvcorrection, kcnvtranslation:kcnvrotation, kcnvl2) = .true.
    stepinfo%cnv_tol(kcnvcorrection, kcnvtranslation:kcnvrotation, kcnvl2) = stepinfo%converg_ddisp
    stepinfo%cnv_check(kcnvcorrection, kcnvlagrange, kcnvl2) = .true.
    stepinfo%cnv_tol(kcnvcorrection, kcnvlagrange, kcnvl2) = stepinfo%converg_lag
    if( .not. present(cnvparam) ) return
    do inorm = kcnvl2, kcnvmax
      do ig = kcnvtranslation, kcnvlagrange
        do iq = kcnvresidual, kcnvcorrection
          if( .not. cnvparam%given(iq,ig,inorm) ) cycle
          stepinfo%cnv_check(iq,ig,inorm) = cnvparam%check(iq,ig,inorm)
          stepinfo%cnv_tol(iq,ig,inorm) = cnvparam%tol(iq,ig,inorm)
        end do
      end do
    end do
  end subroutine
 
  subroutine setup_stepinfo_starttime( stepinfos )
    type( step_info ), pointer, intent(inout) :: stepinfos(:)
 
    integer :: i
 
    stepinfos(1)%starttime = 0.d0
    do i=1,size(stepinfos)-1
      stepinfos(i+1)%starttime = stepinfos(i)%starttime + stepinfos(i)%elapsetime
    end do
  end subroutine
 
  logical function isboundaryactive( bnd, stepinfo )
    integer, intent(in)           :: bnd
    type( step_info ), intent(in) :: stepinfo
    isboundaryactive = .false.
    if( .not. associated( stepinfo%Boundary ) ) return
    if( any( stepinfo%Boundary== bnd ) ) isboundaryactive = .true.
  end function
 
  logical function isloadactive( bnd, stepinfo )
    integer, intent(in)           :: bnd
    type( step_info ), intent(in) :: stepinfo
    isloadactive = .false.
    if( .not. associated( stepinfo%Load ) ) return
    if( any( stepinfo%Load == bnd ) ) isloadactive = .true.
  end function
 
  logical function iscontactactive( bnd, stepinfo )
    integer, intent(in)           :: bnd
    type( step_info ), intent(in) :: stepinfo
    iscontactactive = .false.
    if( .not. associated( stepinfo%Contact ) ) return
    if( any( stepinfo%Contact== bnd ) ) iscontactactive = .true.
  end function
 
  logical function iselemactivationactive( bnd, stepinfo )
    integer, intent(in)           :: bnd
    type( step_info ), intent(in) :: stepinfo
    iselemactivationactive = .false.
    if( .not. associated( stepinfo%ElemActivation ) ) return
    if( any( stepinfo%ElemActivation== bnd ) ) iselemactivationactive = .true.
  end function
 
  subroutine free_stepinfo( step )
    type(step_info), intent(inout) :: step
    if( associated( step%Boundary ) ) deallocate( step%Boundary )
    if( associated( step%Load ) )     deallocate( step%Load )
    if( associated( step%Contact ) )  deallocate( step%Contact )
    if( associated( step%ElemActivation ) )    deallocate( step%ElemActivation )
  end subroutine
 
  subroutine init_aincparam( aincparam )
    type( tparamautoinc ), intent(out) :: aincparam
 
    aincparam%name      = ''
    aincparam%ainc_Rs   = 0.5d0
    aincparam%ainc_Rl   = 1.75d0
    aincparam%NRbound_s = 0
    aincparam%NRbound_s(knstmaxit) = 20
    aincparam%NRbound_s(knstsumit) = 100
    aincparam%NRbound_s(knstciter) = 10
    aincparam%NRbound_l = 0
    aincparam%NRbound_l(knstmaxit) = 10
    aincparam%NRbound_l(knstsumit) = 50
    aincparam%NRbound_l(knstciter) = 5
    aincparam%NRtimes_s = 2
    aincparam%NRtimes_l = 2
    aincparam%ainc_Rc   = 0.25d0
    aincparam%CBbound   = 5
  end subroutine
 
  subroutine init_convergparam( cnvparam )
    type( tparamconverg ), intent(out) :: cnvparam
 
    cnvparam%name  = ''
    cnvparam%given = .false.
    cnvparam%check = .false.
    cnvparam%tol   = 0.d0
  end subroutine
 
end module