add_explicit Subroutine

private subroutine add_explicit(g1, g2, g3, phi, apar, bpar, istep, bd, nl)

FIXME : Add documentation

Arguments

Type IntentOptional 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

Contents

Source Code


Source Code

  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