fstr_Ctrl_TimeInc.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_timeinc
  use m_fstr
  use m_step
  use m_timepoint
 
  implicit none
 
  real(kind=kreal), private :: current_time = 0.d0 
  real(kind=kreal), private :: time_inc = 0.d0 
  real(kind=kreal), private :: time_inc_base = 0.d0 
 
contains
 
  !! basic functions to get time/timeinc
  real(kind=kreal) function fstr_get_time()
    fstr_get_time = current_time
  end function
 
  real(kind=kreal) function fstr_get_timeinc()
    fstr_get_timeinc = time_inc
  end function
 
  real(kind=kreal) function fstr_get_timeinc_base()
    fstr_get_timeinc_base = time_inc_base
  end function
 
  !! basic functions to set/update time/timeinc
  subroutine fstr_set_time( time )
    real(kind=kreal), intent(in) :: time  
    current_time = time
  end subroutine
 
  subroutine fstr_set_timeinc( dtime )
    real(kind=kreal), intent(in) :: dtime  
    time_inc = dtime
  end subroutine
 
  subroutine fstr_set_timeinc_base( dtime_base )
    real(kind=kreal), intent(in) :: dtime_base  
    time_inc_base = dtime_base
  end subroutine
 
  subroutine fstr_proceed_time()
    current_time = current_time + time_inc
  end subroutine
 
  !! functions to print status
  subroutine fstr_timeinc_printstatus_init
    write(ista,'(A10,"-+-",A60,"-+-",A40)') repeat("-",10),repeat("-",60),repeat("-",40)
    write(ista,'(2A5," | ",2A5,2A7,3A12," | ",A40)') '     ', '     ', ' ', ' # of',  'MAX #', 'TOT #',  '','', '', ''
    write(ista,'(2A5," | ",2A5,2A7,3A12," | ",A7,A33)') &
      & 'STEP', 'SUB', 'STAT', ' CONT', 'NEWTON', 'NEWTON', 'START', 'TIME','END', 'MESSAGE',''
    write(ista,'(2A5," | ",2A5,2A7,3A12," | ",A40)') &
      & '', 'STEP', '  ', 'ITER',  'ITER', 'ITER',  'TIME', 'INC', 'TIME',''
    write(ista,'(A10,"-+-",A60,"-+-",A40)') repeat("-",10),repeat("-",60),repeat("-",40)
  end subroutine
 
  subroutine fstr_timeinc_printstatus( stepinfo, fstrPARAM, totstep, substep, NRstatI, NRstatR, AutoINC_stat, Cutback_stat )
    type(step_info), intent(in)       ::  stepinfo
    type(fstr_param), intent(in)      ::  fstrPARAM
    integer(kind=kint), intent(in)    ::  totstep
    integer(kind=kint), intent(in)    ::  substep
    integer(kind=kint), intent(in)    ::  NRstatI(:)  !previous Newton Raphson Iteration Log
    real(kind=kreal), intent(in)      ::  nrstatr(:)  !previous Newton Raphson Iteration Log
    integer(kind=kint), intent(inout) ::  AutoINC_stat
    integer(kind=kint), intent(in)    ::  Cutback_stat
 
    character(len=5) :: cstep, csstep, cstate, ccont
    character(len=7) :: cmaxn, ctotn
    character(len=12) :: ctime, cdtime, etime
    character(len=40) :: message
    type(tparamautoinc) :: pAinc
 
    painc = fstrparam%ainc(stepinfo%AincParam_id)
 
    write(cstep,'(I5)') totstep
    write(csstep,'(I5)') substep
    write(ccont,'(I5)') nrstati(knstciter)
    write(cmaxn,'(I7)') nrstati(knstmaxit)
    write(ctotn,'(I7)') nrstati(knstsumit)
    write(ctime,'(1pE12.4)') current_time
    write(cdtime,'(1pE12.4)') time_inc
 
    if( cutback_stat > 0 ) then
      write(etime,'(1pE12.4)') current_time
      write(cstate,'(I4,A)') cutback_stat,'F'
      if( nrstati(knstdresn) == 1 ) write(message,'(A)') 'Failed to converge due to MAXITER.'
      if( nrstati(knstdresn) == 2 ) write(message,'(A)') 'Failed to converge due to MAXRES.'
      if( nrstati(knstdresn) == 3 ) write(message,'(A)') 'Failed to converge due to MAXCONTITER.'
      if( nrstati(knstdresn) == 4 ) write(message,'(A)') 'Failed to converge due to MatSolveError.'
      if( cutback_stat == painc%CBbound ) write(message,'(A)') '# of successive cutback reached max.'
    else
      write(etime,'(1pE12.4)') current_time+time_inc
      write(cstate,'(A5)') 'S'
      write(message,'(A)') ''
    endif
 
    write(ista,'(2A5," | ",2A5,2A7,3A12," | ",A40)') cstep, csstep, cstate, ccont, &
      &  cmaxn, ctotn, ctime, cdtime, etime, message
 
  end subroutine
 
  subroutine fstr_timeinc_printstatus_final(success_flag)
    logical, intent(in) :: success_flag
    write(ista,'(A10,"-+-",A60,"-+-",A40)') repeat("-",10),repeat("-",60),repeat("-",40)
    if(success_flag) then
      write(ista,'(A)') 'FSTR_SOLVE_NLGEOM HAS COMPLETED SUCCESSFULLY'
    else
      write(ista,'(A)') 'FSTR_SOLVE_NLGEOM HAS NOT COMPLETED SUCCESSFULLY'
    end if
  end subroutine
 
  !! true if step finished
  logical function fstr_timeinc_isstepfinished( stepinfo )
    type(step_info), intent(in)  :: stepinfo
 
    real(kind=kreal) :: endtime, tolerance
 
    endtime = stepinfo%starttime + stepinfo%elapsetime
    ! roundoff error of current_time scales with the absolute time, which can
    ! exceed a tolerance based only on elapsetime when elapsetime << |endtime|
    tolerance = 1.d-12 * max(stepinfo%elapsetime, dabs(endtime))
    fstr_timeinc_isstepfinished = (endtime-current_time < tolerance)
  end function
 
  !! true if current time is time point
  logical function fstr_timeinc_istimepoint( stepinfo, fstrPARAM )
    type(step_info), intent(in)  :: stepinfo
    type(fstr_param), intent(in) :: fstrparam
 
    fstr_timeinc_istimepoint = .false.
    if( stepinfo%inc_type == stepfixedinc ) return
    if( stepinfo%timepoint_id == 0 ) return
 
    fstr_timeinc_istimepoint = is_at_timepoints(current_time,stepinfo%starttime, &
      &  fstrparam%timepoints(stepinfo%timepoint_id))
  end function
 
  !! set time_inc from stepinfo and Newton Raphson iteration status
  subroutine fstr_timeinc_settimeincrement( stepinfo, fstrPARAM, substep, NRstatI, NRstatR, AutoINC_stat, Cutback_stat )
    type(step_info), intent(in)       ::  stepinfo
    type(fstr_param), intent(in)      ::  fstrPARAM
    integer(kind=kint), intent(in)    ::  substep
    integer(kind=kint), intent(in)    ::  NRstatI(:)  !previous Newton Raphson Iteration Log
    real(kind=kreal), intent(in)      ::  nrstatr(:)  !previous Newton Raphson Iteration Log
    integer(kind=kint), intent(inout) ::  AutoINC_stat
    integer(kind=kint), intent(in)    ::  Cutback_stat
 
    real(kind=kreal) :: timeinc0
    real(kind=kreal) :: endtime, remain
    logical          :: to_be_decreased, to_be_increased
    type(tparamautoinc) :: pAinc
 
    painc = fstrparam%ainc(stepinfo%AincParam_id)
 
    if( stepinfo%inc_type == stepfixedinc ) then
      timeinc0 = stepinfo%initdt
    else ! INCTYPE==AUTO
      if( cutback_stat > 0 ) then
        timeinc0 = painc%ainc_Rc*time_inc
        if(myrank == 0) write(*,'(2(A,E10.3))') 'time increment is decreased from ', time_inc_base, ' to ', timeinc0
        autoinc_stat = -1
      else
        if( substep == 1 ) then
          timeinc0 = stepinfo%initdt
          autoinc_stat = 0
        else
          timeinc0 = time_inc_base
 
          !decrease condition
          to_be_decreased = .false.
          if( nrstati(knstmaxit) > painc%NRbound_s(knstmaxit) ) to_be_decreased = .true.
          if( nrstati(knstsumit) > painc%NRbound_s(knstsumit) ) to_be_decreased = .true.
          if( nrstati(knstciter) > painc%NRbound_s(knstciter) ) to_be_decreased = .true.
 
          !increase condition
          to_be_increased = .true.
          if( nrstati(knstmaxit) > painc%NRbound_l(knstmaxit) ) to_be_increased = .false.
          if( nrstati(knstsumit) > painc%NRbound_l(knstsumit) ) to_be_increased = .false.
          if( nrstati(knstciter) > painc%NRbound_l(knstciter) ) to_be_increased = .false.
 
          ! count # of times that increase/decrease condition has been satisfied
          ! AutoINC_stat < 0 ... decrease condition has been satisfied -AutoINC_stat times
          ! AutoINC_stat > 0 ... increase condition has been satisfied AutoINC_stat times
          if(to_be_decreased) then
            autoinc_stat = min(autoinc_stat,0) - 1
          else if(to_be_increased) then
            autoinc_stat = max(autoinc_stat,0) + 1
          else
            autoinc_stat = 0
          end if
 
          if( autoinc_stat <= -painc%NRtimes_s ) then
            timeinc0 = painc%ainc_Rs*timeinc0
            if(myrank == 0) write(*,'(2(A,E10.3))') 'time increment is decreased from ', time_inc_base, ' to ', timeinc0
          else if( autoinc_stat >= painc%NRtimes_l ) then
            timeinc0 = dmin1(painc%ainc_Rl*timeinc0,stepinfo%maxdt)
            if(myrank == 0) write(*,'(2(A,E10.3))') 'time increment is increased from ', time_inc_base, ' to ', timeinc0
          end if
        end if
      end if
    end if
 
    if( timeinc0 < stepinfo%mindt ) then
      if(myrank == 0) then
        write(*,'(2(A,E10.3))') 'Error: current time increment ',timeinc0,' is smaller than lower bound',stepinfo%mindt
        call fstr_timeinc_printstatus_final(.false.)
      endif
      stop
    end if
 
    endtime = stepinfo%starttime + stepinfo%elapsetime
    time_inc_base = dmin1(timeinc0,endtime-current_time)
    if( stepinfo%timepoint_id > 0 ) then
      remain = get_remain_to_next_timepoints(current_time,stepinfo%starttime,fstrparam%timepoints(stepinfo%timepoint_id))
      time_inc = dmin1(time_inc_base,remain)
    else
      time_inc = time_inc_base
    end if
 
  end subroutine
end module m_fstr_timeinc