Time integration

The m_af_advance module can be used to perform time integration.

Methods

Currently, the following explicit methods are included:

  !> Forward Euler method
  integer, parameter, public :: af_forward_euler    = 1
  !> Heun's method (AKA modified Euler's method, explicit trapezoidal rule), CFL
  !> coefficient of 1. See e.g. https://en.wikipedia.org/wiki/Heun%27s_method
  integer, parameter, public :: af_heuns_method     = 2
  !> Midpoint method, see e.g. https://en.wikipedia.org/wiki/Midpoint_method
  integer, parameter, public :: af_midpoint_method  = 3
  !> Optimal 3-stage third-order SSPRK method (Shu & Osher), CFL coefficient of
  !> 1, see e.g. https://doi.org/10.1137/S0036142902419284
  integer, parameter, public :: af_ssprk33_method   = 4
  !> Optimal 4-stage third-order SSPRK method (Ruuth & Spiteri), CFL coefficient
  !> of 2, see e.g. https://doi.org/10.1137/S0036142902419284
  integer, parameter, public :: af_ssprk43_method   = 5
  !> 1st order IMEX method, forward Euler and then implicit solve
  integer, parameter, public :: af_imex_euler       = 6
  !> Trapezoidal IMEX method (2nd order)
  integer, parameter, public :: af_imex_trapezoidal = 7
  !> Classic 4th order Runge Kutta method, see e.g.
  !> https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
  integer, parameter, public :: af_rk4_method       = 8

New integrators can be added relatively easily by modifying the af_advance routine.

How to use

Time integration can be performed with a call to m_af_advance::af_advance.

User-supplied forward Euler method

To use the built-in time integration, a forward_euler routine has to be provided, see m_af_advance::subr_feuler for details. This routine will then be used to construct the various time integration schemes.

For higher-order schemes, the idea is that multiple copies will be stored for variables that will be integrated, which correspond to different temporal states. The number of copies can be specified in m_af_core::af_add_cc_variable.

The user-provided forward Euler method then takes three extra arguments that indicate which temporal states to operate on. For a simple equation of the form

y' = f(y)

the method should provide a solution

y_out = y_prev + dt * f(y_deriv)

More details can be found at m_af_advance::subr_feuler, and an example from euler_gas_dynamics.f90 is shown below

  subroutine forward_euler(tree, dt, dt_stiff, dt_lim, time, s_deriv, n_prev, &
       s_prev, w_prev, s_out, i_step, n_steps)
    type(af_t), intent(inout) :: tree
    real(dp), intent(in)      :: dt
    real(dp), intent(in)      :: dt_stiff       !< Time step for stiff terms
    real(dp), intent(inout)   :: dt_lim
    real(dp), intent(in)      :: time
    integer, intent(in)       :: s_deriv
    integer, intent(in)       :: n_prev         !< Number of previous states
    integer, intent(in)       :: s_prev(n_prev) !< Previous states
    real(dp), intent(in)      :: w_prev(n_prev) !< Weights of previous states
    integer, intent(in)       :: s_out
    integer, intent(in)       :: i_step, n_steps
    real(dp)                  :: dummy_dt(0)
 
    call flux_generic_tree(tree, n_vars, variables, s_deriv, fluxes, dt_lim, &
         max_wavespeed, get_fluxes, flux_dummy_modify, flux_dummy_line_modify, &
         to_primitive, to_conservative, af_limiter_vanleer_t)
    call flux_update_densities(tree, dt, n_vars, variables, n_vars, variables, fluxes, &
         s_deriv, n_prev, s_prev, w_prev, s_out, flux_dummy_source, 0, dummy_dt)
  end subroutine forward_euler