gs2_time.f90 Source File


Source Code

Source Code

!> FIXME : Add documentation
module gs2_time

  implicit none


  public :: init_gs2_time, finish_gs2_time
  public :: user_dt, code_dt, update_time, code_dt_old
  public :: time_history_available, dt_not_set
  public :: code_dt_prev1, code_dt_prev2
  public :: user_time, code_time, check_time_step_too_large
  public :: save_dt_min, save_dt_max, save_dt, save_dt_cfl, write_dt
  public :: init_tstart, init_delt
  public :: code_dt_cfl, code_dt_min, user2code, code_dt_max
  public :: get_adams_bashforth_coefficients
  public :: wunits, woutunits, tunits

  real, parameter :: dt_not_set = -1.0
  real :: user_dt, code_dt 
  real :: code_dt_prev1 = dt_not_set
  real :: code_dt_prev2 = dt_not_set
  real :: user_dt_cfl = 1e8
  real :: code_dt_cfl = 1e8

  real :: code_dt_min, user_dt_min
  real :: code_dt_max, user_dt_max

  ! added May 18, 2009 to take care of problems
  ! in exb_shear calculation after change in time step size
  real :: code_dt_old = 0.

! GGH points out that this initialization is not necessary (we think)
  real :: user_time = dt_not_set
  real :: code_time = dt_not_set

  real, dimension (:), allocatable :: wunits, woutunits, tunits


  !> Initialise the gs2_time module
  !> Note this doesn't actually initialise everything (for example
  !> it doesn't call init_tstart) due to module dependencies.
  subroutine init_gs2_time
    use run_parameters, only: delt_option_switch, delt_option_auto, delt
    use gs2_save, only: init_dt
    implicit none
    real :: delt_prev1, delt_prev2, delt0, delt1, delt2, delt_max, user_delt_max
    integer :: istatus

    user_delt_max = delt
    delt_prev1 = dt_not_set
    delt_prev2 = dt_not_set

    if (delt_option_switch == delt_option_auto) then
       call init_dt (delt0, delt1, delt2, delt_max, istatus, dt_not_set)
       if (delt0 /= dt_not_set) delt = delt0
       if (delt1 /= dt_not_set) delt_prev1 = delt1
       if (delt2 /= dt_not_set) delt_prev2 = delt2
       if (delt_max /= dt_not_set) user_delt_max = delt_max

    call init_delt (delt, delt_prev1, delt_prev2)
    call user2code (user_delt_max, code_dt_max)

    ! omega_* normalization of time:
    call adjust_time_norm
  end subroutine init_gs2_time

  !> Finalise the gs2_time module
  subroutine finish_gs2_time
    implicit none
    ! Note we don't (currently) reset module variables like user_time

    if(allocated(wunits)) deallocate (wunits)
    if(allocated(woutunits)) deallocate (woutunits)
    if(allocated(tunits)) deallocate (tunits)

  end subroutine finish_gs2_time

  !> FIXME : Add documentation
  subroutine adjust_time_norm
    use file_utils, only: error_unit
    use mp, only: proc0
    use kt_grids, only: aky, naky
    use run_parameters, only: wstar_units
    implicit none

    if(.not. allocated(wunits)) allocate (wunits(naky))
    if(.not. allocated(woutunits)) allocate (woutunits(naky))
    if(.not. allocated(tunits)) allocate (tunits(naky))

    !CMR: Sep 2010
    ! Attempt to understand time normalisation variables, which are arrays(naky)
    !            This is a generally very useful variable to store ky dependent
    !            timestep in the code time normalisation.
    !            Used to multiply ky independent source terms on RHS of GKE.
    !            Auxiliary variable.  Used to save compute operations when
    !            evaluating source terms on RHS of GKE that are proportional to ky.
    !            !! The Mysterious factor 1/2 Explained !!
    !            The factor of 1/2 arises because those source terms were first
    !            specified using the normalisation Tref=mref vtref^2
    ! [R Numata et al, "AstroGK: Astrophysical gyrokinetics code", JCP, 2010].
    if (wstar_units) then
       wunits = 1.0
       where (aky /= 0.0)
          tunits = 2.0/aky
          tunits = 0.0
       end where
       if (any(tunits == 0.0) .and. proc0) then
          write (error_unit(), *) &
               "WARNING: wstar_units=.true. and aky=0.0: garbage results"
          print *, &
               "WARNING: wstar_units=.true. and aky=0.0: garbage results"
       end if
       tunits = 1.0
       wunits = aky/2.0
    end if
    woutunits = 1.0/tunits
  end subroutine adjust_time_norm

  !> FIXME : Add documentation
  subroutine init_tstart (tstart)
    implicit none
    real, intent (in) :: tstart

    user_time = tstart
    code_time = tstart

  end subroutine init_tstart

  !> Set the full timestep history.
  !! Inputs set the current, previous and last but one timesteps
  subroutine init_delt (delt, delt1, delt2)
    implicit none
    real, intent (in) :: delt, delt1, delt2
! delt* are user input, from the run_parameters module.
! In a perfect world, we could have a gs2_time namelist.
    user_dt = delt
    code_dt = delt
    code_dt_prev1 = delt1
    code_dt_prev2 = delt2
  end subroutine init_delt

  !> Check if the current time step is too big
  subroutine check_time_step_too_large(reset)
    implicit none
    logical, intent(out) :: reset
  end subroutine check_time_step_too_large

  !> FIXME : Add documentation  
  subroutine update_time
    implicit none
! for using a changing-timestep AB3 algorithm for the nonlinear term  
    code_dt_prev2 = code_dt_prev1
    code_dt_prev1 = code_dt

! MAB+CMR, 21/5/09: set code_dt_old to code_dt BEFORE any changes in timestep
    code_dt_old = code_dt
    code_time = code_time + code_dt
    user_time = user_time + user_dt

  end subroutine update_time

  !> FIXME : Add documentation
  subroutine save_dt_cfl (delt_cfl)
    implicit none
    real, intent (in) :: delt_cfl
    code_dt_cfl = delt_cfl
    user_dt_cfl = delt_cfl

  end subroutine save_dt_cfl

  !> Sets the minimum timestep size allowed based on passed value.
  subroutine save_dt_min (dt_min)
    implicit none
    real, intent (in) :: dt_min

    user_dt_min = dt_min
    code_dt_min = dt_min

  end subroutine save_dt_min

  !> Sets the maximum timestep size allowed based on passed value.
  subroutine save_dt_max (dt_max)
    implicit none
    real, intent (in) :: dt_max

    user_dt_max = dt_max
    code_dt_max = dt_max

  end subroutine save_dt_max

  !> FIXME : Add documentation
  subroutine save_dt(delt)
    implicit none
    real, intent (in) :: delt

    code_dt = delt
    user_dt = delt
  end subroutine save_dt

  !> FIXME : Add documentation  
  subroutine write_dt
    implicit none
    if (user_dt_cfl > 0. .and. user_dt_cfl < 1.e7) &
         write(*,*) user_dt_cfl,' : ',user_dt
  end subroutine write_dt

  !> FIXME : Add documentation  
  subroutine user2code (usertime, codetime)
    implicit none
    real, intent (in) :: usertime
    real, intent (out) :: codetime
    codetime = usertime
  end subroutine user2code

  !> Reports how many previous time steps are available
  !> by checking if the previous code_dt values have been
  !> set to anything. Note we assume that if a time step is
  !> set then all more recent values have also been set.
  function time_history_available()
    implicit none
    integer :: time_history_available
    if (code_dt_prev2 .ne. dt_not_set) then
       time_history_available =  3
    else if (code_dt_prev1 .ne. dt_not_set) then
       time_history_available =  2
    else !Note here we assume code_dt is always set
       time_history_available =  1
  end function time_history_available

  !> Returns the current set of Adams Bashforth coefficients.
  !> The order is determined by how much time history we have
  !> available up to a maximum of 3rd order. The coefficients
  !> are generalised for variable timestep.
  function get_adams_bashforth_coefficients(maximum_order, target_dt) result(coefficients)
    use optionals, only: get_option_with_default
    implicit none
    !> If present then sets the maxmimum method order to use when
    !> calculating the coefficients. Limited by how much time history
    !> we have available and no higher than 3.
    integer, intent(in), optional :: maximum_order
    !> If present sets the intended time step size we want to take. Otherwise
    !> assume we want to step by code_dt.
    real, intent(in), optional :: target_dt
    !> The result contains up to three coefficients. The first entry
    !> corresponds to the coefficient for the latest source term, the
    !> second entry corresponds to the coefficient for the previous source
    !> term and the last entry is for the oldest source term.
    real, dimension(:), allocatable :: coefficients
    real :: this_dt
    integer :: order

    order = time_history_available()
    if (present(maximum_order)) then
       ! Limit order in use if passed maximum_order is less
       ! than order set by available history
       order = min(order, maximum_order)
    end if

    this_dt = get_option_with_default(target_dt, code_dt)

    ! We have one coefficient for each point of history we keep at
    ! this order

    select case (order)
    case (3) !Most common case - 3rd order coefficients
       ! see Tatsuno & Dorland, Physics of Plasmas 13, 092107 (2006)
       ! if this_dt = code_dt_prev1 = code_dt_prev2, this reduces to usual AB3
       ! coefficients 23/12, -4/3 and 5/12.
       coefficients(1) = (((this_dt + &
            code_dt_prev1)**2/code_dt_prev1)*( &
            (this_dt+code_dt_prev1)/3. + code_dt_prev2/2. )  - &
            code_dt_prev1*(code_dt_prev1/3. + code_dt_prev2/2.) ) / &
            (this_dt*(code_dt_prev1 + code_dt_prev2))

       coefficients(2) = - this_dt * ( this_dt/3. + &
            (code_dt_prev1+code_dt_prev2)/2. ) / &

       coefficients(3) = this_dt * ( this_dt/3. + &
            code_dt_prev1/2. ) / ( &
            code_dt_prev2*(code_dt_prev1+code_dt_prev2) )

    case (2) !2nd order coefficients

       coefficients(1) = (0.5*this_dt + code_dt_prev1) / code_dt_prev1

       coefficients(2) = -0.5*this_dt / code_dt_prev1

    case (1) !1st order coefficients - this is just explicit Euler

       coefficients(1) = 1.0
    end select

  end function get_adams_bashforth_coefficients
end module gs2_time