m_af_advance.f90

#include "cpp_macros.h"
 
module m_af_advance
  use m_af_types
 
  implicit none
  private
 
  integer, parameter, public :: af_num_integrators  = 8
  integer, parameter, public :: af_forward_euler    = 1
  integer, parameter, public :: af_heuns_method     = 2
  integer, parameter, public :: af_midpoint_method  = 3
  integer, parameter, public :: af_ssprk33_method   = 4
  integer, parameter, public :: af_ssprk43_method   = 5
  integer, parameter, public :: af_imex_euler       = 6
  integer, parameter, public :: af_imex_trapezoidal = 7
  integer, parameter, public :: af_rk4_method       = 8
 
  character(len=af_nlen), public :: af_integrator_names(af_num_integrators) = &
       [character(len=af_nlen) :: "forward_euler", "heuns_method", &
       "midpoint_method", "ssprk33", "ssprk43", "imex_euler", &
       "imex_trapezoidal", "rk4"]
 
  integer, parameter, public :: &
       af_advance_num_steps(af_num_integrators) = [1, 2, 2, 3, 4, 1, 2, 4]
 
  integer, parameter :: req_copies(af_num_integrators) = af_advance_num_steps
 
  logical, parameter :: req_implicit(af_num_integrators) = &
       [.false., .false., .false., .false., .false., .true., .true., .false.]
 
  interface
 
     subroutine subr_feuler(tree, dt, dt_stiff, dt_lim, time, s_deriv, n_prev, &
          s_prev, w_prev, s_out, i_step, n_steps)
       import
       type(af_t), intent(inout) :: tree
       real(dp), intent(in)      :: dt
       real(dp), intent(in)      :: dt_stiff
       real(dp), intent(inout)   :: dt_lim
       real(dp), intent(in)      :: time
       integer, intent(in)       :: s_deriv
       integer, intent(in)       :: n_prev
       integer, intent(in)       :: s_prev(n_prev)
       real(dp), intent(in)      :: w_prev(n_prev)
       integer, intent(in)       :: s_out
       integer, intent(in)       :: i_step
       integer, intent(in)       :: n_steps
     end subroutine subr_feuler
 
     subroutine subr_implicit(tree, dt_stiff, time, n_prev, s_prev, w_prev, s_out)
       import
       type(af_t), intent(inout) :: tree
       real(dp), intent(in)      :: dt_stiff
       real(dp), intent(in)      :: time
       integer, intent(in)       :: n_prev
       integer, intent(in)       :: s_prev(n_prev)
       real(dp), intent(in)      :: w_prev(n_prev)
       integer, intent(in)       :: s_out
     end subroutine subr_implicit
  end interface
 
  public :: af_advance
 
contains
 
  subroutine af_advance(tree, dt, dt_lim, time, i_cc, time_integrator, &
       forward_euler, implicit_solver)
    type(af_t), intent(inout) :: tree
    real(dp), intent(in)      :: dt
    real(dp), intent(out)     :: dt_lim
    real(dp), intent(inout)   :: time
    integer, intent(in)       :: i_cc(:)
    integer, intent(in)       :: time_integrator
    procedure(subr_feuler)    :: forward_euler
    procedure(subr_implicit), optional  :: implicit_solver
    integer                   :: n_steps
 
    real(dp), parameter :: third = 1/3.0_dp
    real(dp), parameter :: sixth = 1/6.0_dp
 
    if (time_integrator < 1 .or. time_integrator > af_num_integrators) &
         error stop "Invalid time integrator"
 
    if (any(tree%cc_num_copies(i_cc) < req_copies(time_integrator))) &
         error stop "Not enough copies available"
 
    if (req_implicit(time_integrator) .and. .not. present(implicit_solver)) &
         error stop "implicit_solver required"
 
    n_steps = af_advance_num_steps(time_integrator)
    dt_lim = 1e100_dp
 
    select case (time_integrator)
    case (af_forward_euler)
       call forward_euler(tree, dt, dt, dt_lim, time, 0, &
            1, [0], [1.0_dp], 0, 1, n_steps)
    case (af_midpoint_method)
       call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time, 0, &
            1, [0], [1.0_dp], 1, 1, n_steps)
       call forward_euler(tree, dt, dt, dt_lim, time+0.5_dp*dt, 1, &
            1, [0], [1.0_dp], 0, 2, n_steps)
    case (af_heuns_method)
       call forward_euler(tree, dt, dt, dt_lim, time, 0, &
            1, [0], [1.0_dp], 1, 1, n_steps)
       call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time+dt, 1, &
            2, [0, 1], [0.5_dp, 0.5_dp], 0, 2, n_steps)
    case (af_ssprk33_method)
       call forward_euler(tree, dt, dt, dt_lim, time, 0, &
            1, [0], [1.0_dp], 1, 1, n_steps)
       call forward_euler(tree, 0.25_dp*dt, 0.25_dp*dt, dt_lim, time+dt, 1, &
            2, [0, 1], [0.75_dp, 0.25_dp], 2, 2, n_steps)
       call forward_euler(tree, 2*third*dt, 2*third*dt, dt_lim, time+0.5_dp*dt, 2, &
            2, [0, 2], [third, 2*third], 0, 3, n_steps)
    case (af_ssprk43_method)
       call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time, 0, &
            1, [0], [1.0_dp], 1, 1, n_steps)
       call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time+0.5_dp*dt, 1, &
            1, [1], [1.0_dp], 2, 2, n_steps)
       call forward_euler(tree, sixth*dt, sixth*dt, dt_lim, time+dt, 2, &
            2, [0, 2], [2*third, third], 3, 3, n_steps)
       call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time+0.5_dp*dt, 3, &
            1, [3], [1.0_dp], 0, 4, n_steps)
    case (af_imex_euler)
       call forward_euler(tree, dt, 0*dt, dt_lim, time, 0, &
            1, [0], [1.0_dp], 0, 1, n_steps)
       call implicit_solver(tree, dt, time, 1, [0], [1.0_dp], 0)
    case (af_imex_trapezoidal)
       ! Compute y* = y_n + dt*F0(y_n) + 0.5*dt * (F1(y_n) + F1(y*))
       call forward_euler(tree, dt, 0.5_dp*dt, dt_lim, time, 0, &
            1, [0], [1.0_dp], 1, 1, n_steps)
       call implicit_solver(tree, 0.5_dp*dt, time, 1, [1], [1.0_dp], 1)
 
       ! Compute y_n+1 = y_n + 0.5*dt * (F(y_n) + F(y*))
       call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time, 0, &
            1, [0], [1.0_dp], 0, 1, n_steps)
       call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time, 1, &
            1, [0], [1.0_dp], 0, 2, n_steps)
    case (af_rk4_method)
       ! This looks different than the standard formulation in most textbooks.
       ! The idea is to construct the states needed for the derivatives, and
       ! then take a linear combination. Note the negative coefficient used in
       ! the last step.
       call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time, 0, &
            1, [0], [1.0_dp], 1, 1, n_steps)
       call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time+0.5_dp*dt, 1, &
            1, [0], [1.0_dp], 2, 2, n_steps)
       call forward_euler(tree, dt, dt, dt_lim, time+0.5_dp*dt, 2, &
            1, [0], [1.0_dp], 3, 3, n_steps)
       call forward_euler(tree, sixth*dt, sixth*dt, dt_lim, time+dt, 3, &
            4, [0, 1, 2, 3], [-third, third, 2*third, third], 0, 4, n_steps)
    case default
       error stop "Unknown time integrator"
    end select
 
    time = time + dt
  end subroutine af_advance
 
end module m_af_advance