FIXME : Add documentation
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
complex, | intent(inout), | dimension (-ntgrid:,:,g_lo%llim_proc:) | :: | g1 | ||
complex, | intent(inout), | dimension (-ntgrid:,:,g_lo%llim_proc:) | :: | g2 | ||
complex, | intent(inout), | dimension (-ntgrid:,:,g_lo%llim_proc:) | :: | g3 | ||
complex, | intent(in), | dimension (-ntgrid:,:,:) | :: | phi | ||
complex, | intent(in), | dimension (-ntgrid:,:,:) | :: | apar | ||
complex, | intent(in), | dimension (-ntgrid:,:,:) | :: | bpar | ||
integer, | intent(in) | :: | istep | |||
real, | intent(in) | :: | bd | |||
logical, | intent(in), | optional | :: | nl |
subroutine add_explicit (g1, g2, g3, phi, apar, bpar, istep, bd, nl)
use theta_grid, only: ntgrid
use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, is_idx
use dist_fn_arrays, only: g
use gs2_time, only: save_dt_cfl, code_dt, check_time_step_too_large
use run_parameters, only: reset, immediate_reset
use unit_tests, only: debug_message
use mp, only: nb_max_allreduce, max_allreduce
use array_utils, only: copy, zero_array
implicit none
complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in out) :: g1, g2, g3
complex, dimension (-ntgrid:,:,:), intent (in) :: phi, apar, bpar
integer, intent (in) :: istep
real, intent (in) :: bd
logical, intent (in), optional :: nl
integer, parameter :: verb = 4
integer :: iglo, ik
real :: dt_cfl
real :: error
real :: target_dt
call debug_message(verb, 'nonlinear_terms::add_explicit starting')
if (initializing) then
if (present(nl)) then
dt_cfl = dt_cfl_default_large
call save_dt_cfl (dt_cfl)
! Ensure the module level max_vel is consistent with
! the selected cfl limiting timestep.
max_vel_components = 1.0/dt_cfl
end if
return
endif
!
! @todo Currently not self-starting. Need to fix this.
!
call debug_message(verb, 'nonlinear_terms::add_explicit copying old terms')
! The following if statement ensures that the nonlinear term is calculated only once
! per timestep. It guards against calculating the nonlinear term in two cases:
! 1. the second call to `timeadv` in a usual timestep
! 2. both calls to `timeadv` after the CFL condition check has triggered a reset
! of this timestep (when `immediate_reset = .true.`). As the nonlinear term is
! independent of `delt`, it does not need to be calculated again.
!
! N.B. When `immediate_reset = .false.`, the timestep is completed even though the
! CFL condition is broken. `delt` and the response matrix are changed at the end of
! the timestep; the following timestep is "normal", and the nonlinear term is
! calculated.
if (istep /= istep_last) then
istep_last = istep
! Shuffle time history along -- this is a place pointers might help
! avoid some copying.
! Should we only need to do this if we have a nonlinear simulation?
! It's possible we might have other explicit sources, in which case we
! probably always need to do this shuffling, however in that case this
! subroutine probably best belongs in a different module.
call copy(g2, g3)
call copy(g1, g2)
call debug_message(verb, 'nonlinear_terms::add_explicit copied old terms')
! if running nonlinearly, then compute the nonlinear term at grid points
! and store it in g1
if (present(nl)) then
call debug_message(verb, 'nonlinear_terms::add_explicit calling add_nl')
! If we're not using various limits then we'd like to zero
! out the maximum velocity however, we later take 1/max_vel
! so let's just set it small.
max_vel_components = epsilon(0.0)
call add_nl (g, g1, phi, apar, bpar)
! takes g1 at grid points and returns 2*g1 at cell centers
call nl_center (g1, bd)
! Currently the module level max_vel_components contains the
! maximum x and y velocities found from the nonlinear term
! on this processor. This can be used to find the cfl
! limting time step.
! Reset the x and y components of max_vel_components if we
! don't want to consider the cfl limit
if (.not. use_cfl_limit) then
max_vel_components(1:2) = epsilon(0.0)
end if
! Calculate the effective maximum velocity based on the error
! based time step limit (really we estimate the maximum time step
! to keep our errors below our target and then convert to an
! effective velocity).
if (use_order_based_error .and. (istep > istep_error_start)) then
! Get the current order based error estimate
error = estimate_error(g1, g2, g3)
! Based on the current error estimate and the expected error scaling
! work out what time step we expect to give us the target error.
! Note that we don't scale by cfl here as the user has direct control
! over error_target already. Also note that by taking the square root
! we are assuming that the error looks like dt^2. We should probably
! look at what order schemes we're actually comparing to determine
! the appropriate dependency.
target_dt =code_dt*sqrt(error_target/(error))
! Now we have an error limited time step lets convert this
! to an effective maximum velocity
max_vel_components(3) = 1.0/target_dt
end if
! Now we can reduce our local limiting velocities to get
! a global value for the maximum velocity. We only need to
! reduce if at least one time limit check is active.
if (use_cfl_limit .or. use_order_based_error) then
! Estimate the global cfl limit based on max_vel calculated
! in `add_nl`.If immediate_reset is true, check the max
! allreduce immediately. If not, overlap comms and work,
! and check for reset at the end of the timestep.
if(immediate_reset)then
call max_allreduce(max_vel_components)
! This call can update the value of reset.
call set_cfl_time_step_and_check_if_too_large
else
! Do a nonblocking reduce and check later in
! [[gs2_reinit:check_time_step]].
call nb_max_allreduce(max_vel_components, cfl_req_hand)
end if
else
! If we don't have any checks active then we better make sure that the
! limiting time step is set. We will have set max_vel earlier to
! epsilon(0.0) so we can just use the usual routine to set this, but
! can skip communications. Really we probably just need to do this once
! but for now just put this here to be on the safe side. The cost should
! be low.
call set_cfl_time_step_and_check_if_too_large
end if
if(reset) then
return !Return if resetting
endif
else
call zero_array(g1)
end if
end if
if (zip) then
do iglo = g_lo%llim_proc, g_lo%ulim_proc
ik = ik_idx(g_lo,iglo)
if (ik /= 1) then
g (:,1,iglo) = 0.
g (:,2,iglo) = 0.
g1(:,1,iglo) = 0.
g1(:,2,iglo) = 0.
end if
end do
end if
istep_last = istep
end subroutine add_explicit