Afivo 0.3
All Classes Namespaces Functions Variables Pages
m_af_advance.f90
1!> Module with methods to perform time integration
3#include "cpp_macros.h"
4 use m_af_types
5
6 implicit none
7 private
8
9 integer, parameter, public :: af_num_integrators = 8
10
11 !< [time_integration_schemes]
12 !> Forward Euler method
13 integer, parameter, public :: af_forward_euler = 1
14 !> Heun's method (AKA modified Euler's method, explicit trapezoidal rule), CFL
15 !> coefficient of 1. See e.g. https://en.wikipedia.org/wiki/Heun%27s_method
16 integer, parameter, public :: af_heuns_method = 2
17 !> Midpoint method, see e.g. https://en.wikipedia.org/wiki/Midpoint_method
18 integer, parameter, public :: af_midpoint_method = 3
19 !> Optimal 3-stage third-order SSPRK method (Shu & Osher), CFL coefficient of
20 !> 1, see e.g. https://doi.org/10.1137/S0036142902419284
21 integer, parameter, public :: af_ssprk33_method = 4
22 !> Optimal 4-stage third-order SSPRK method (Ruuth & Spiteri), CFL coefficient
23 !> of 2, see e.g. https://doi.org/10.1137/S0036142902419284
24 integer, parameter, public :: af_ssprk43_method = 5
25 !> 1st order IMEX method, forward Euler and then implicit solve
26 integer, parameter, public :: af_imex_euler = 6
27 !> Trapezoidal IMEX method (2nd order)
28 integer, parameter, public :: af_imex_trapezoidal = 7
29 !> Classic 4th order Runge Kutta method, see e.g.
30 !> https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
31 integer, parameter, public :: af_rk4_method = 8
32 !< [time_integration_schemes]
33
34 character(len=af_nlen), public :: af_integrator_names(af_num_integrators) = &
35 [character(len=af_nlen) :: "forward_euler", "heuns_method", &
36 "midpoint_method", "ssprk33", "ssprk43", "imex_euler", &
37 "imex_trapezoidal", "rk4"]
38
39 !> How many steps the time integrators take
40 integer, parameter, public :: &
41 af_advance_num_steps(af_num_integrators) = [1, 2, 2, 3, 4, 1, 2, 4]
42
43 !> How many variable copies are required for the time integrators
44 integer, parameter :: req_copies(af_num_integrators) = af_advance_num_steps
45
46 !> Whether an implicit solver is required for the scheme
47 logical, parameter :: req_implicit(af_num_integrators) = &
48 [.false., .false., .false., .false., .false., .true., .true., .false.]
49
50 interface
51 !> Interface for a generic forward Euler scheme for time integration
52 !>
53 !> This method should advance the solution over a time dt. The method
54 !> assumes that several copies are stored for the variables to be
55 !> integrated. It should then operate on these different copies, which
56 !> correspond to temporal states. In this way, higher-order time
57 !> integration schemes can be constructed.
58 !>
59 !> The meaning of the temporal states is as follows. For an equation y' =
60 !> f(y), the method should perform:
61 !> y_out = sum(w_prev * y_prev) + dt * f(y_deriv).
62 !>
63 !> If the index of the variable `y` is `i`, then the index of `y_out` is
64 !> `i+s_out`, etc.
65 !>
66 !> In case of IMEX schemes, the time step dt_stiff (which is then not
67 !> always equal to dt) should be used for stiff terms. The equation to be
68 !> solved should be interpreted as:
69 !>
70 !> d/dt y = F0(y) + F1(y)
71 !>
72 !> where F0 is the non-stiff part and F1 is the stiff part
73 subroutine subr_feuler(tree, dt, dt_stiff, dt_lim, time, s_deriv, n_prev, &
74 s_prev, w_prev, s_out, i_step, n_steps)
75 import
76 type(af_t), intent(inout) :: tree
77 real(dp), intent(in) :: dt !< Time step for regular terms
78 real(dp), intent(in) :: dt_stiff !< Time step for stiff terms (IMEX)
79 real(dp), intent(inout) :: dt_lim !< Computed time step limit
80 real(dp), intent(in) :: time !< Current time
81 integer, intent(in) :: s_deriv !< State to compute derivatives from
82 integer, intent(in) :: n_prev !< Number of previous states
83 integer, intent(in) :: s_prev(n_prev) !< Previous states
84 real(dp), intent(in) :: w_prev(n_prev) !< Weights of previous states
85 integer, intent(in) :: s_out !< Output state
86 integer, intent(in) :: i_step !< Step of the integrator
87 integer, intent(in) :: n_steps !< Total number of steps
88 end subroutine subr_feuler
89
90 !> Interface for a implicit solver for integrating stiff terms over time
91 !>
92 !> This method should solve an equation of the form
93 !>
94 !> F(y_out, dt_stiff) = sum(w_prev * y_prev),
95 !>
96 !> where dt_stiff is the time step.
97 subroutine subr_implicit(tree, dt_stiff, time, n_prev, s_prev, w_prev, s_out)
98 import
99 type(af_t), intent(inout) :: tree
100 real(dp), intent(in) :: dt_stiff !< Time step for stiff terms
101 real(dp), intent(in) :: time !< Current time
102 integer, intent(in) :: n_prev !< Number of previous states
103 integer, intent(in) :: s_prev(n_prev) !< Previous states
104 real(dp), intent(in) :: w_prev(n_prev) !< Weights of previous states
105 integer, intent(in) :: s_out !< Output state
106 end subroutine subr_implicit
107 end interface
108
109 public :: af_advance
110
111contains
112
113 !> Advance solution over dt using time_integrator
114 !>
115 !> The user should supply a forward Euler method as documented in subr_feuler.
116 !> The indices of the cell-centered variables that will be operated on should
117 !> also be provided, so that higher-order schemes can be constructed
118 !> automatically from the forward Euler method.
119 !>
120 !> @todo check whether time in sub-steps is accurate with some test
121 subroutine af_advance(tree, dt, dt_lim, time, i_cc, time_integrator, &
122 forward_euler, implicit_solver)
123 type(af_t), intent(inout) :: tree
124 real(dp), intent(in) :: dt !< Current time step
125 real(dp), intent(out) :: dt_lim !< Time step limit
126 real(dp), intent(inout) :: time !< Current time
127 integer, intent(in) :: i_cc(:) !< Index of cell-centered variables
128 !> One of the pre-defined time integrators (e.g. af_heuns_method)
129 integer, intent(in) :: time_integrator
130 !> Forward Euler method provided by the user
131 procedure(subr_feuler) :: forward_euler
132 !> Implicit solver to be used with IMEX schemes
133 procedure(subr_implicit), optional :: implicit_solver
134 integer :: n_steps
135
136 real(dp), parameter :: third = 1/3.0_dp
137 real(dp), parameter :: sixth = 1/6.0_dp
138
139 if (time_integrator < 1 .or. time_integrator > af_num_integrators) &
140 error stop "Invalid time integrator"
141
142 if (any(tree%cc_num_copies(i_cc) < req_copies(time_integrator))) &
143 error stop "Not enough copies available"
144
145 if (req_implicit(time_integrator) .and. .not. present(implicit_solver)) &
146 error stop "implicit_solver required"
147
148 n_steps = af_advance_num_steps(time_integrator)
149 dt_lim = 1e100_dp
150
151 select case (time_integrator)
152 case (af_forward_euler)
153 call forward_euler(tree, dt, dt, dt_lim, time, 0, &
154 1, [0], [1.0_dp], 0, 1, n_steps)
155 case (af_midpoint_method)
156 call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time, 0, &
157 1, [0], [1.0_dp], 1, 1, n_steps)
158 call forward_euler(tree, dt, dt, dt_lim, time+0.5_dp*dt, 1, &
159 1, [0], [1.0_dp], 0, 2, n_steps)
160 case (af_heuns_method)
161 call forward_euler(tree, dt, dt, dt_lim, time, 0, &
162 1, [0], [1.0_dp], 1, 1, n_steps)
163 call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time+dt, 1, &
164 2, [0, 1], [0.5_dp, 0.5_dp], 0, 2, n_steps)
165 case (af_ssprk33_method)
166 call forward_euler(tree, dt, dt, dt_lim, time, 0, &
167 1, [0], [1.0_dp], 1, 1, n_steps)
168 call forward_euler(tree, 0.25_dp*dt, 0.25_dp*dt, dt_lim, time+dt, 1, &
169 2, [0, 1], [0.75_dp, 0.25_dp], 2, 2, n_steps)
170 call forward_euler(tree, 2*third*dt, 2*third*dt, dt_lim, time+0.5_dp*dt, 2, &
171 2, [0, 2], [third, 2*third], 0, 3, n_steps)
172 case (af_ssprk43_method)
173 call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time, 0, &
174 1, [0], [1.0_dp], 1, 1, n_steps)
175 call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time+0.5_dp*dt, 1, &
176 1, [1], [1.0_dp], 2, 2, n_steps)
177 call forward_euler(tree, sixth*dt, sixth*dt, dt_lim, time+dt, 2, &
178 2, [0, 2], [2*third, third], 3, 3, n_steps)
179 call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time+0.5_dp*dt, 3, &
180 1, [3], [1.0_dp], 0, 4, n_steps)
181 case (af_imex_euler)
182 call forward_euler(tree, dt, 0*dt, dt_lim, time, 0, &
183 1, [0], [1.0_dp], 0, 1, n_steps)
184 call implicit_solver(tree, dt, time, 1, [0], [1.0_dp], 0)
185 case (af_imex_trapezoidal)
186 ! Compute y* = y_n + dt*F0(y_n) + 0.5*dt * (F1(y_n) + F1(y*))
187 call forward_euler(tree, dt, 0.5_dp*dt, dt_lim, time, 0, &
188 1, [0], [1.0_dp], 1, 1, n_steps)
189 call implicit_solver(tree, 0.5_dp*dt, time, 1, [1], [1.0_dp], 1)
190
191 ! Compute y_n+1 = y_n + 0.5*dt * (F(y_n) + F(y*))
192 call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time, 0, &
193 1, [0], [1.0_dp], 0, 1, n_steps)
194 call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time, 1, &
195 1, [0], [1.0_dp], 0, 2, n_steps)
196 case (af_rk4_method)
197 ! This looks different than the standard formulation in most textbooks.
198 ! The idea is to construct the states needed for the derivatives, and
199 ! then take a linear combination. Note the negative coefficient used in
200 ! the last step.
201 call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time, 0, &
202 1, [0], [1.0_dp], 1, 1, n_steps)
203 call forward_euler(tree, 0.5_dp*dt, 0.5_dp*dt, dt_lim, time+0.5_dp*dt, 1, &
204 1, [0], [1.0_dp], 2, 2, n_steps)
205 call forward_euler(tree, dt, dt, dt_lim, time+0.5_dp*dt, 2, &
206 1, [0], [1.0_dp], 3, 3, n_steps)
207 call forward_euler(tree, sixth*dt, sixth*dt, dt_lim, time+dt, 3, &
208 4, [0, 1, 2, 3], [-third, third, 2*third, third], 0, 4, n_steps)
209 case default
210 error stop "Unknown time integrator"
211 end select
212
213 time = time + dt
214 end subroutine af_advance
215
216end module m_af_advance
Interface for a generic forward Euler scheme for time integration.
Interface for a implicit solver for integrating stiff terms over time.
Module with methods to perform time integration.
Definition: m_af_advance.f90:2
This module contains the basic types and constants that are used in the NDIM-dimensional version of A...
Definition: m_af_types.f90:3
Type which stores all the boxes and levels, as well as some information about the number of boxes,...
Definition: m_af_types.f90:326