add_explicit Subroutine

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

FIXME : Add documentation


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


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
    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
    real, parameter :: zero = epsilon(0.0)
    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 = 1.e8
          call save_dt_cfl (dt_cfl)

          ! Ensure the module level max_vel is consistent with
          ! the selected cfl limiting timestep.
          max_vel = 1.0/dt_cfl
       end if

    ! @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.
       g3 = g2
       g2 = g1

       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(gryfx_zonal%on) then
             call add_nl_gryfx (g1)
             call add_nl (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 contains the maximum
          ! velocity found from the nonlinear term on this
          ! processor. This can be used to find the cfl limting time
          ! step.

          ! If we're not using the cfl limit then we'd like to zero
          ! out the maximum velocity however, we later take 1/max_vel
          ! so let's just set it small.
          if (.not. use_cfl_limit) then
             max_vel = 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 = max(max_vel, 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.
                call max_allreduce(max_vel)
                ! This call can update the value of reset.
                call set_cfl_time_step_and_check_if_too_large
                ! Do a nonblocking reduce and check later in
                ! [[gs2_reinit:check_time_step]].
                call nb_max_allreduce(max_vel, cfl_req_hand)
             end if
             ! 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

          g1 = 0.
       end if

#ifdef LOWFLOW
       ! do something

    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