Line data Source code
1 0 : \cond NEVER
2 : Distributed under the MIT License.
3 : See LICENSE.txt for details.
4 : \endcond
5 : # IMEX {#tutorial_imex}
6 :
7 : \tableofcontents
8 :
9 : ### Introduction to IMEX
10 :
11 : IMEX (implicit-explicit) integration is a method for stabilizing the numerical
12 : evolution of stiff differential equations. While "stiff" is not a precisely
13 : defined term, it generally refers to equations where the time step required to
14 : achieve numerical stability is much shorter than the timescale of the evolution
15 : of the solution, often by orders of magnitude. An example of a stiff equation
16 : is
17 :
18 : \f{equation}
19 : \label{eq:example-equation}
20 : \frac{dy}{dt} = -k (y - f(t))
21 : \f}
22 :
23 : with $k$ very large. After a brief initial transient, the solution is very
24 : well-approximated as $y = f(t) + O(k^{-1})$, but will be unstable if evolved
25 : using an explicit time stepper with a step size larger than approximately
26 : $k^{-1}$, even if $f$ is slowly varying.
27 :
28 : To see the source of the stability problem, consider evolving
29 : $\eqref{eq:example-equation}$ using Euler's method. For a step from $t_0$ to
30 : $t_1 = t_0 + \Delta t$, we find
31 :
32 : \f{equation}
33 : \label{eq:example-Euler}
34 : y_1 = y_0 + \Delta t \frac{dy}{dt}(y_0, t_0)
35 : = y_0 - k \Delta t (y_0 - f(t_0)),
36 : \f}
37 :
38 : giving the divergence from the known infinite-$k$ solution of
39 :
40 : \f{equation}
41 : y_1 - f(t_1) = y_0 - f(t_1) - k \Delta t (y_0 - f(t_0)).
42 : \f}
43 :
44 : If $k \Delta t$, is large, the last term will dominate and the deviation from
45 : the approximate solution will increase by roughly that factor every step.
46 :
47 : To improve this, we can evolve using an implicit time stepper. Taking the same
48 : step as above with the backwards Euler method gives
49 :
50 : \f{equation}
51 : \label{eq:example-backwards-Euler}
52 : y_1 = y_0 + \Delta t \frac{dy}{dt}(y_1, t_1)
53 : = y_0 - k \Delta t (y_1 - f(t_1)),
54 : \f}
55 :
56 : for a deviation of
57 :
58 : \f{equation}
59 : y_1 - f(t_1) = \frac{y_0 - f(t_1)}{1 + k \Delta t}.
60 : \f}
61 :
62 : This goes to zero as $k \Delta t$ becomes large, showing convergence to the
63 : infinite-$k$ solution.
64 :
65 : The unconditional (linear) stability achievable by implicit integrators allows
66 : them to solve problems not feasible with explicit methods, but it comes at a
67 : significant price. While the explicit update $\eqref{eq:example-Euler}$
68 : directly gave $y_1$, the implicit update $\eqref{eq:example-backwards-Euler}$
69 : required the equation to be solved for the updated quantity. In this example
70 : that was fairly easy, but, in a real application, finding an analytic solution
71 : is likely to be difficult, if not impossible. Implicit integrators are
72 : therefore usually used with a numerical solve of the update equation. Even for
73 : simple systems, this root-finding procedure adds significantly to the
74 : computational cost of the method, and for a system like a large PDE evolution
75 : it is intractable.
76 :
77 : To mitigate this problem, we turn to the hybrid IMEX methods. The idea behind
78 : IMEX is to split the derivative into two terms, one without stability problems
79 : that will be treated explicitly, and another that will be treated implicitly.
80 : Notationally, we will write
81 :
82 : \f{equation}
83 : \frac{dy}{dt}(y, t) = E(y, t) + I(y, t).
84 : \f}
85 :
86 : The simplest IMEX integrator combines the two forms of Euler's method from
87 : above:
88 :
89 : \f{equation}
90 : y_1 = y_0 + \Delta t (E(y_0, t_0) + I(y_1, t_1)).
91 : \f}
92 :
93 : (This integrator is not implemented in SpECTRE.) In our example above, we can,
94 : as an example, split $\eqref{eq:example-equation}$ as
95 :
96 : \f{align}
97 : E(y, t) &= k f(t) &
98 : I(y, t) &= -k y,
99 : \f}
100 :
101 : which gives
102 :
103 : \f{equation}
104 : y_1 = y_0 + k \Delta t (f(t_0) - y_1)
105 : = \frac{y_0 + k \Delta t\, f(t_0)}{1 + k \Delta t}.
106 : \f}
107 :
108 : In the limit where $k \Delta t$ goes to infinity, this gives $y_1 = f(t_0)$,
109 : which is as small a deviation from the infinite-$k$ solution as possible when
110 : $f$ is treated explicitly.
111 :
112 : Again, in this simple case there is no downside to using an IMEX integrator
113 : over an explicit one, or an implicit method over IMEX. The gain appears in
114 : real cases where the implicit-step equation cannot be solved exactly.
115 :
116 : #### Properties of IMEX time steppers
117 :
118 : Any IMEX time stepper can be used as either an implicit or an explicit time
119 : stepper by choosing $E$ or $I$ to be zero. As such, many significant
120 : properties of an IMEX time stepper are actually properties of one or the other
121 : part. We won't discuss the standard properties of explicit time steppers here,
122 : but most of them, such as error estimation, are used in SpECTRE ignoring the
123 : implicit part.
124 :
125 : As the entire point of using implicit methods is to increase stability, the
126 : most important properties of the implicit part of a time stepper are stability
127 : properties. There are many stability classifications, but we will only discuss
128 : the most important here, and those only qualitatively. More information can be
129 : found in \cite Hairer1996 and \cite Kennedy2016.
130 :
131 : The basic desire for an implicit method is that as the equation being evolved
132 : becomes increasingly stiff the evolution remains stable. (Linear) stability
133 : for an infinitely large step size (or, equivalently, an infinitely stiff
134 : equation) is known as *A-stability*. Adams methods above second order cannot
135 : be A-stable, so IMEX work generally focuses on Runge-Kutta schemes. All IMEX
136 : time steppers implemented in SpECTRE are A-stable.
137 :
138 : A stronger stability condition is known as *L-stability*. L-stability is
139 : similar to A-stability, but instead of merely requiring that an analytically
140 : decaying equation does not numerically diverge in the large-step-size limit, we
141 : require that the solution reaches zero in a single step. Since IMEX
142 : applications can take steps orders of magnitude larger than the decay timescale
143 : of stiff terms, this property is often desirable.
144 :
145 : In addition to the stability properties, there is one interesting property of
146 : the IMEX stepper as a whole. A wide class of time steppers, including all
147 : explicit (global time-stepping) methods implementable in the current SpECTRE
148 : interface, preserve (linear) conserved quantities. A linear conserved quantity
149 : is a linear combination of the evolved variables that is analytically constant
150 : under evolution, independent of the initial conditions. In physics, a common
151 : example is the integral of the density, i.e., the total mass. Such quantities
152 : will be numerically preserved during evolution to the level of roundoff error,
153 : even if that is much smaller than the truncation error in the solution.
154 :
155 : Both the explicit and implicit portions of all IMEX time steppers that we will
156 : consider are conservative in this manner. However, the combination of the two
157 : parts into an IMEX scheme does not necessarily preserve this property. In
158 : general, both the explicit and implicit parts of the split derivative must
159 : individually conserve something for the full IMEX method to conserve it as
160 : well. Some IMEX schemes, however, do preserve conserved quantities,
161 : independent of how the derivative is split between the explicit and implicit
162 : parts. Such methods are termed *stiffly accurate*. (This property is also
163 : useful in deriving and analyzing implicit schemes, but that is not of great
164 : importance to users of the methods.)
165 :
166 : ### IMEX in SpECTRE
167 :
168 : SpECTRE supports IMEX integration, with restrictions to reduce the cost of the
169 : solve of the implicit equation as much as possible. Additionally, not all
170 : integration schemes extend to IMEX integration, so the list of available time
171 : steppers is limited for IMEX evolutions.
172 :
173 : #### Mathematical restrictions
174 :
175 : For problems of interest to us, the evolved variables consist of several fields
176 : coupled by the system derivative. The implicit solve is more expensive when
177 : performed on more variables, so there is a method to restrict the solve to a
178 : subset of the system if only some variables are affected by stiff terms. The
179 : system can define one or more *implicit sectors*, which are subsets of the
180 : evolved variables with implicit terms that can be solved independently from one
181 : another. (If sectors depend on variables from other sectors, the evolution
182 : system defines the order in which the implicit updates are applied.)
183 :
184 : As an example, for a fluid coupled to neutrinos, the neutrinos obey stiff
185 : equations under some conditions. If we consider neutrino flavors to be
186 : non-interacting, they can each be placed in their own implicit sector and
187 : solved independently. The fluid equations will usually be non-stiff, so the
188 : fluid variables will not be in any sector, even though they will still be used
189 : as arguments for computing the stiff sources.
190 :
191 : Within each sector, restrictions are placed on the form of the implicit
192 : derivative (the $I$ above). When evolving a hyperbolic PDE, in general, the
193 : equation for an implicit step becomes an elliptic PDE. This is far too complex
194 : and expensive to solve for every integration substep. We therefore restrict
195 : the implicit derivative to only contain source terms, i.e., terms that depend
196 : on the sector variables only in a pointwise manner, rather than through
197 : derivatives or similar. (The terms may still depend on the derivatives of
198 : non-sector variables.) This results in each point in the domain having an
199 : independent set of algebraic equations to solve.
200 :
201 : There is a small complexity in handling non-autonomous (time-dependent) systems
202 : of equations, such as those for evolutions with control systems. Some methods
203 : are defined in a way such that the values of time for the explicit and implicit
204 : parts of a substep appear to be inconsistent. SpECTRE always uses the time
205 : derived from the explicit portion of the method, which is equivalent to
206 : treating time as an additional evolved variable with a constant explicit
207 : derivative of 1 and implicit derivative of 0.
208 :
209 : SpECTRE supports using an analytic solution of the implicit-step equation as
210 : well as two modes for numerical root-finding: implicit and semi-implicit. The
211 : analytic mode can be chosen in the definition of the implicit sector if the
212 : form of the implicit source allows an analytic solution. The numerical methods
213 : can be toggled at runtime. The fully implicit method does a numerical
214 : root-find, while the semi-implicit one linearizes the equation and solves that.
215 : Semi-implicit solves are faster, and, experimentally, they still stabilize the
216 : evolution fairly well. Using a semi-implicit solve instead of a fully implicit
217 : one does not affect conservation properties of the method. Both numerical
218 : solvers require the jacobian of the implicit source to be coded, but an
219 : analytic solution does not.
220 :
221 : #### Creation of SpECTRE IMEX executables
222 :
223 : The SpECTRE IMEX interface is defined by two protocols: \link
224 : imex::protocols::ImexSystem \endlink and \link imex::protocols::ImplicitSector
225 : \endlink. The normal evolution-system interface now only contains the explicit
226 : portion of the derivative, and the implicit portion is defined by the contents
227 : of the `implicit_sectors` typelist. See those protocols for details.
228 :
229 : In order to use an IMEX system, the implicit solver must be explicitly
230 : instantiated for each sector. This is done in a separate source file in the
231 : system directory. As an example, the instantiation file for the sectors in one
232 : of the tests is
233 :
234 : \include tests/Unit/Helpers/Evolution/Imex/DoImplicitStepInstantiate.cpp
235 :
236 : In the executable itself, four sets of changes are necessary. First, the
237 : time-stepper type, usually defined as a type alias near the start of the
238 : metavariables, must be changed from `TimeStepper` to `ImexTimeStepper`
239 : (IMEX-LTS is not currently supported), and an entry must be added to the
240 : `factory_creation` list:
241 :
242 : ```
243 : tmpl::pair<ImexTimeStepper, TimeSteppers::imex_time_steppers>
244 : ```
245 :
246 : (The `TimeStepper` line can be removed, although doing so is not necessary.)
247 :
248 : Second, the IMEX actions must be added after the corresponding explicit
249 : actions:
250 :
251 : * In the initialization phase, the `imex::Initialize<system>` mutator must be
252 : applied by adding it to the arguments of
253 : `Initialization::Actions::InitializeItems`. It uses objects initialized by
254 : `Initialization::TimeStepperHistory`, so must appear after that.
255 :
256 : * In the step actions,
257 : ```
258 : imex::Actions::RecordTimeStepperData<system>
259 : ```
260 : must be added after each occurrence of
261 : `Actions::RecordTimeStepperData<system>`. (There may be only one.)
262 :
263 : * Again in the step actions,
264 : ```
265 : imex::Actions::DoImplicitStep<system>
266 : ```
267 : must be added after each occurrence of `Actions::UpdateU<system>`. (Again,
268 : there may be only one.)
269 :
270 : Third, the `imex::ImplicitDenseOutput<system>` dense output postprocessor must
271 : be added to the argument of the `evolution::Actions::RunEventsAndDenseTriggers`
272 : action. It must appear before any other preprocessors that use the evolved
273 : variables.
274 :
275 : Finally, header includes must be added for all these things. The required
276 : headers are
277 : ```
278 : #include "Evolution/Imex/Actions/DoImplicitStep.hpp"
279 : #include "Evolution/Imex/Actions/RecordTimeStepperData.hpp"
280 : #include "Evolution/Imex/ImplicitDenseOutput.hpp"
281 : #include "Evolution/Imex/Initialize.hpp"
282 : #include "Time/TimeSteppers/ImexTimeStepper.hpp"
283 : ```
|
