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
 
  ! 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
 
  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             
 
    integer :: num_substep
    integer :: max_iter
    integer :: max_contiter
    integer :: amp_id
    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
 
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%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
    stepinfo%converg_lag = 1.d-4
    stepinfo%converg_ddisp = 1.d-8
    stepinfo%maxres = 1.d+10
    stepinfo%timepoint_id = 0
    stepinfo%AincParam_id = 0
  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 fstr_print_steps( nfile, steps )
    integer, intent(in)         :: nfile
    type(step_info), intent(in) :: steps(:)
    integer :: i, j, nstep, nbc
    nstep = size(steps)
 
    write( nfile, * ) "-----Information of steps:",nstep
 
    do i=1,nstep
      write( nfile, * ) "  -----Step:",i
      write(nfile,*) steps(i)%solution, steps(i)%elapsetime, steps(i)%converg, &
        steps(i)%num_substep, steps(i)%max_iter
      if( associated( steps(i)%Boundary ) ) then
        nbc = size( steps(i)%Boundary )
        write(nfile,*) "  Boundary conditions"
        write(nfile,*) ( steps(i)%Boundary(j),j=1,nbc )
      endif
      if( associated( steps(i)%Load ) ) then
        nbc = size( steps(i)%Load )
        write(nfile,*) "  External load conditions"
        write(nfile,*) ( steps(i)%Load(j),j=1,nbc )
      endif
      if( associated( steps(i)%Contact ) ) then
        nbc = size( steps(i)%Contact )
        write(nfile,*) "  Contact conditions"
        write(nfile,*) ( steps(i)%Contact(j),j=1,nbc )
      endif
      if( associated( steps(i)%ElemActivation ) ) then
        nbc = size( steps(i)%ElemActivation )
        write(nfile,*) "  ElemActivation conditions"
        write(nfile,*) ( steps(i)%ElemActivation(j),j=1,nbc )
      endif
    enddo
  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
 
end module