!> FIXME : Add documentation module gs2_time implicit none private 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 contains !> 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 endif 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) ! TUNITS: DT(KY)=TUNITS(KY).CODE_DT ! 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. ! WUNITS: WUNITS(KY)=AKY(KY)*TUNITS(KY)/2 ! 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]. !CMRend if (wstar_units) then wunits = 1.0 where (aky /= 0.0) tunits = 2.0/aky elsewhere 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 else 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 reset=(code_dt>code_dt_cfl) 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 endif 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() ! Limit order in use if passed maximum_order is less ! than order set by available history order = min(order, get_option_with_default(maximum_order, order)) this_dt = get_option_with_default(target_dt, code_dt) ! We have one coefficient for each point of history we keep at ! this order allocate(coefficients(order)) select case (order) case (3) !Most common case - 3rd order coefficients ! see Tatsuno & Dorland, Physics of Plasmas 13, 092107 (2006) ! http://www.cscamm.umd.edu/publications/TatsunoDorland-PoP06_CS-06-11.pdf ! 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. ) / & (code_dt_prev1*code_dt_prev2) 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