diff --git a/docs/References.bib b/docs/References.bib index 673bba495350..3362b0c7bab9 100644 --- a/docs/References.bib +++ b/docs/References.bib @@ -1134,6 +1134,18 @@ @book{Hairer1993 year = 1993 } +@Inbook{Hairer1996, + author="Hairer, Ernst and Wanner, Gerhard", + title="Solving Ordinary Differential Equations II: Stiff and + Differential-Algebraic Problems", + year="1996", + publisher="Springer Berlin Heidelberg", + address="Berlin, Heidelberg", + isbn="978-3-642-05221-7", + doi="10.1007/978-3-642-05221-7_1", + url="https://doi.org/10.1007/978-3-642-05221-7_1" +} + @article{Handmer2014qha, author = "Handmer, Casey J. and Szilagyi, B.", title = "{Spectral Characteristic Evolution: A New Algorithm for diff --git a/docs/Tutorials/Imex.md b/docs/Tutorials/Imex.md new file mode 100644 index 000000000000..ffe71ece601f --- /dev/null +++ b/docs/Tutorials/Imex.md @@ -0,0 +1,283 @@ +\cond NEVER +Distributed under the MIT License. +See LICENSE.txt for details. +\endcond +# IMEX {#tutorial_imex} + +\tableofcontents + +### Introduction to IMEX + +IMEX (implicit-explicit) integration is a method for stabilizing the numerical +evolution of stiff differential equations. While "stiff" is not a precisely +defined term, it generally refers to equations where the time step required to +achieve numerical stability is much shorter than the timescale of the evolution +of the solution, often by orders of magnitude. An example of a stiff equation +is + +\f{equation} + \label{eq:example-equation} + \frac{dy}{dt} = -k (y - f(t)) +\f} + +with $k$ very large. After a brief initial transient, the solution is very +well-approximated as $y = f(t) + O(k^{-1})$, but will be unstable if evolved +using an explicit time stepper with a step size larger than approximately +$k^{-1}$, even if $f$ is slowly varying. + +To see the source of the stability problem, consider evolving +$\eqref{eq:example-equation}$ using Euler's method. For a step from $t_0$ to +$t_1 = t_0 + \Delta t$, we find + +\f{equation} + \label{eq:example-Euler} + y_1 = y_0 + \Delta t \frac{dy}{dt}(y_0, t_0) + = y_0 - k \Delta t (y_0 - f(t_0)), +\f} + +giving the divergence from the known infinite-$k$ solution of + +\f{equation} + y_1 - f(t_1) = y_0 - f(t_1) - k \Delta t (y_0 - f(t_0)). +\f} + +If $k \Delta t$, is large, the last term will dominate and the deviation from +the approximate solution will increase by roughly that factor every step. + +To improve this, we can evolve using an implicit time stepper. Taking the same +step as above with the backwards Euler method gives + +\f{equation} + \label{eq:example-backwards-Euler} + y_1 = y_0 + \Delta t \frac{dy}{dt}(y_1, t_1) + = y_0 - k \Delta t (y_1 - f(t_1)), +\f} + +for a deviation of + +\f{equation} + y_1 - f(t_1) = \frac{y_0 - f(t_1)}{1 + k \Delta t}. +\f} + +This goes to zero as $k \Delta t$ becomes large, showing convergence to the +infinite-$k$ solution. + +The unconditional (linear) stability achievable by implicit integrators allows +them to solve problems not feasible with explicit methods, but it comes at a +significant price. While the explicit update $\eqref{eq:example-Euler}$ +directly gave $y_1$, the implicit update $\eqref{eq:example-backwards-Euler}$ +required the equation to be solved for the updated quantity. In this example +that was fairly easy, but, in a real application, finding an analytic solution +is likely to be difficult, if not impossible. Implicit integrators are +therefore usually used with a numerical solve of the update equation. Even for +simple systems, this root-finding procedure adds significantly to the +computational cost of the method, and for a system like a large PDE evolution +it is intractable. + +To mitigate this problem, we turn to the hybrid IMEX methods. The idea behind +IMEX is to split the derivative into two terms, one without stability problems +that will be treated explicitly, and another that will be treated implicitly. +Notationally, we will write + +\f{equation} + \frac{dy}{dt}(y, t) = E(y, t) + I(y, t). +\f} + +The simplest IMEX integrator combines the two forms of Euler's method from +above: + +\f{equation} + y_1 = y_0 + \Delta t (E(y_0, t_0) + I(y_1, t_1)). +\f} + +(This integrator is not implemented in SpECTRE.) In our example above, we can, +as an example, split $\eqref{eq:example-equation}$ as + +\f{align} + E(y, t) &= k f(t) & + I(y, t) &= -k y, +\f} + +which gives + +\f{equation} + y_1 = y_0 + k \Delta t (f(t_0) - y_1) + = \frac{y_0 + k \Delta t\, f(t_0)}{1 + k \Delta t}. +\f} + +In the limit where $k \Delta t$ goes to infinity, this gives $y_1 = f(t_0)$, +which is as small a deviation from the infinite-$k$ solution as possible when +$f$ is treated explicitly. + +Again, in this simple case there is no downside to using an IMEX integrator +over an explicit one, or an implicit method over IMEX. The gain appears in +real cases where the implicit-step equation cannot be solved exactly. + +#### Properties of IMEX time steppers + +Any IMEX time stepper can be used as either an implicit or an explicit time +stepper by choosing $E$ or $I$ to be zero. As such, many significant +properties of an IMEX time stepper are actually properties of one or the other +part. We won't discuss the standard properties of explicit time steppers here, +but most of them, such as error estimation, are used in SpECTRE ignoring the +implicit part. + +As the entire point of using implicit methods is to increase stability, the +most important properties of the implicit part of a time stepper are stability +properties. There are many stability classifications, but we will only discuss +the most important here, and those only qualitatively. More information can be +found in \cite Hairer1996 and \cite Kennedy2016. + +The basic desire for an implicit method is that as the equation being evolved +becomes increasingly stiff the evolution remains stable. (Linear) stability +for an infinitely large step size (or, equivalently, an infinitely stiff +equation) is known as *A-stability*. Adams methods above second order cannot +be A-stable, so IMEX work generally focuses on Runge-Kutta schemes. All IMEX +time steppers implemented in SpECTRE are A-stable. + +A stronger stability condition is known as *L-stability*. L-stability is +similar to A-stability, but instead of merely requiring that an analytically +decaying equation does not numerically diverge in the large-step-size limit, we +require that the solution reaches zero in a single step. Since IMEX +applications can take steps orders of magnitude larger than the decay timescale +of stiff terms, this property is often desirable. + +In addition to the stability properties, there is one interesting property of +the IMEX stepper as a whole. A wide class of time steppers, including all +explicit (global time-stepping) methods implementable in the current SpECTRE +interface, preserve (linear) conserved quantities. A linear conserved quantity +is a linear combination of the evolved variables that is analytically constant +under evolution, independent of the initial conditions. In physics, a common +example is the integral of the density, i.e., the total mass. Such quantities +will be numerically preserved during evolution to the level of roundoff error, +even if that is much smaller than the truncation error in the solution. + +Both the explicit and implicit portions of all IMEX time steppers that we will +consider are conservative in this manner. However, the combination of the two +parts into an IMEX scheme does not necessarily preserve this property. In +general, both the explicit and implicit parts of the split derivative must +individually conserve something for the full IMEX method to conserve it as +well. Some IMEX schemes, however, do preserve conserved quantities, +independent of how the derivative is split between the explicit and implicit +parts. Such methods are termed *stiffly accurate*. (This property is also +useful in deriving and analyzing implicit schemes, but that is not of great +importance to users of the methods.) + +### IMEX in SpECTRE + +SpECTRE supports IMEX integration, with restrictions to reduce the cost of the +solve of the implicit equation as much as possible. Additionally, not all +integration schemes extend to IMEX integration, so the list of available time +steppers is limited for IMEX evolutions. + +#### Mathematical restrictions + +For problems of interest to us, the evolved variables consist of several fields +coupled by the system derivative. The implicit solve is more expensive when +performed on more variables, so there is a method to restrict the solve to a +subset of the system if only some variables are affected by stiff terms. The +system can define one or more *implicit sectors*, which are subsets of the +evolved variables with implicit terms that can be solved independently from one +another. (If sectors depend on variables from other sectors, the evolution +system defines the order in which the implicit updates are applied.) + +As an example, for a fluid coupled to neutrinos, the neutrinos obey stiff +equations under some conditions. If we consider neutrino flavors to be +non-interacting, they can each be placed in their own implicit sector and +solved independently. The fluid equations will usually be non-stiff, so the +fluid variables will not be in any sector, even though they will still be used +as arguments for computing the stiff sources. + +Within each sector, restrictions are placed on the form of the implicit +derivative (the $I$ above). When evolving a hyperbolic PDE, in general, the +equation for an implicit step becomes an elliptic PDE. This is far too complex +and expensive to solve for every integration substep. We therefore restrict +the implicit derivative to only contain source terms, i.e., terms that depend +on the sector variables only in a pointwise manner, rather than through +derivatives or similar. (The terms may still depend on the derivatives of +non-sector variables.) This results in each point in the domain having an +independent set of algebraic equations to solve. + +There is a small complexity in handling non-autonomous (time-dependent) systems +of equations, such as those for evolutions with control systems. Some methods +are defined in a way such that the values of time for the explicit and implicit +parts of a substep appear to be inconsistent. SpECTRE always uses the time +derived from the explicit portion of the method, which is equivalent to +treating time as an additional evolved variable with a constant explicit +derivative of 1 and implicit derivative of 0. + +SpECTRE supports using an analytic solution of the implicit-step equation as +well as two modes for numerical root-finding: implicit and semi-implicit. The +analytic mode can be chosen in the definition of the implicit sector if the +form of the implicit source allows an analytic solution. The numerical methods +can be toggled at runtime. The fully implicit method does a numerical +root-find, while the semi-implicit one linearizes the equation and solves that. +Semi-implicit solves are faster, and, experimentally, they still stabilize the +evolution fairly well. Using a semi-implicit solve instead of a fully implicit +one does not affect conservation properties of the method. Both numerical +solvers require the jacobian of the implicit source to be coded, but an +analytic solution does not. + +#### Creation of SpECTRE IMEX executables + +The SpECTRE IMEX interface is defined by two protocols: \link +imex::protocols::ImexSystem \endlink and \link imex::protocols::ImplicitSector +\endlink. The normal evolution-system interface now only contains the explicit +portion of the derivative, and the implicit portion is defined by the contents +of the `implicit_sectors` typelist. See those protocols for details. + +In order to use an IMEX system, the implicit solver must be explicitly +instantiated for each sector. This is done in a separate source file in the +system directory. As an example, the instantiation file for the sectors in one +of the tests is + +\include tests/Unit/Helpers/Evolution/Imex/DoImplicitStepInstantiate.cpp + +In the executable itself, four sets of changes are necessary. First, the +time-stepper type, usually defined as a type alias near the start of the +metavariables, must be changed from `TimeStepper` to `ImexTimeStepper` +(IMEX-LTS is not currently supported), and an entry must be added to the +`factory_creation` list: + +``` +tmpl::pair +``` + +(The `TimeStepper` line can be removed, although doing so is not necessary.) + +Second, the IMEX actions must be added after the corresponding explicit +actions: + +* In the initialization phase, the `imex::Initialize` mutator must be + applied by adding it to the arguments of + `Initialization::Actions::InitializeItems`. It uses objects initialized by + `Initialization::TimeStepperHistory`, so must appear after that. + +* In the step actions, +``` +imex::Actions::RecordTimeStepperData +``` + must be added after each occurrence of + `Actions::RecordTimeStepperData`. (There may be only one.) + +* Again in the step actions, +``` +imex::Actions::DoImplicitStep +``` + must be added after each occurrence of `Actions::UpdateU`. (Again, + there may be only one.) + +Third, the `imex::ImplicitDenseOutput` dense output postprocessor must +be added to the argument of the `evolution::Actions::RunEventsAndDenseTriggers` +action. It must appear before any other preprocessors that use the evolved +variables. + +Finally, header includes must be added for all these things. The required +headers are +``` +#include "Evolution/Imex/Actions/DoImplicitStep.hpp" +#include "Evolution/Imex/Actions/RecordTimeStepperData.hpp" +#include "Evolution/Imex/ImplicitDenseOutput.hpp" +#include "Evolution/Imex/Initialize.hpp" +#include "Time/TimeSteppers/ImexTimeStepper.hpp" +``` diff --git a/docs/Tutorials/Tutorials.md b/docs/Tutorials/Tutorials.md index 7595d072adef..4cd7617bb7f8 100644 --- a/docs/Tutorials/Tutorials.md +++ b/docs/Tutorials/Tutorials.md @@ -20,3 +20,4 @@ See LICENSE.txt for details. checkpointing and restarting - \subpage tutorial_events_and_triggers - How to specify code to run in an input file +- \subpage tutorial_imex - How to add IMEX support to an evolution system diff --git a/src/Evolution/Imex/Protocols/ImexSystem.hpp b/src/Evolution/Imex/Protocols/ImexSystem.hpp index d8da425bb0e8..9bfb14006011 100644 --- a/src/Evolution/Imex/Protocols/ImexSystem.hpp +++ b/src/Evolution/Imex/Protocols/ImexSystem.hpp @@ -15,8 +15,10 @@ namespace imex::protocols { /// conforming to protocols::ImplicitSector, each of which describes /// an implicit solve to be performed during time steps. /// -/// For efficiency, the tensors in each sector are required to be -/// adjacent in the system's variables. +/// For efficiency, the tags in the `tensors` type alias of each +/// sector are required to be adjacent in the system's variables. +/// +/// \snippet DoImplicitStepSector.hpp ImexSystem struct ImexSystem { template struct test { diff --git a/src/Evolution/Imex/Protocols/ImplicitSector.hpp b/src/Evolution/Imex/Protocols/ImplicitSector.hpp index 8623d365b8b4..1cc780fe1403 100644 --- a/src/Evolution/Imex/Protocols/ImplicitSector.hpp +++ b/src/Evolution/Imex/Protocols/ImplicitSector.hpp @@ -87,10 +87,10 @@ namespace imex::protocols { /// All `Variables` in the DataBox, including the sources and source /// jacobian, will be initialized to zero with a single grid point. /// -/// \snippet Test_SolveImplicitSector.cpp ImplicitSector -/// \snippet Test_SolveImplicitSector.cpp initial_guess +/// \snippet DoImplicitStepSector.hpp simple_sector /// -/// Examples of definitions of the implicit source and jacobian: +/// Examples of definitions of a complicated implicit source and +/// jacobian: /// /// \snippet Test_SolveImplicitSector.cpp source /// \snippet Test_SolveImplicitSector.cpp jacobian diff --git a/tests/Unit/Evolution/Imex/Test_SolveImplicitSector.cpp b/tests/Unit/Evolution/Imex/Test_SolveImplicitSector.cpp index 04fb14f512ad..59f4be52bda4 100644 --- a/tests/Unit/Evolution/Imex/Test_SolveImplicitSector.cpp +++ b/tests/Unit/Evolution/Imex/Test_SolveImplicitSector.cpp @@ -135,7 +135,6 @@ struct Preparer { }; // End stuff only used for DataBox handling -// [initial_guess] struct AnalyticSolution { using return_tags = tmpl::list; using argument_tags = tmpl::list; @@ -159,7 +158,6 @@ struct AnalyticSolution { return {get(var1).size(), imex::GuessResult::ExactSolution}; } }; -// [initial_guess] struct InitialGuess { using return_tags = tmpl::list; @@ -261,7 +259,6 @@ struct Jacobian { } }; -// [ImplicitSector] template struct ImplicitSector : tt::ConformsTo { using tensors = tmpl::list; @@ -288,7 +285,6 @@ struct ImplicitSector : tt::ConformsTo { using solve_attempts = tmpl::list; }; -// [ImplicitSector] // ::tensors doesn't depend on the template parameter using sector_variables_tag = Tags::Variables::tensors>; diff --git a/tests/Unit/Helpers/Evolution/Imex/DoImplicitStepSector.hpp b/tests/Unit/Helpers/Evolution/Imex/DoImplicitStepSector.hpp index e8a9370a1c29..5df2c5235b46 100644 --- a/tests/Unit/Helpers/Evolution/Imex/DoImplicitStepSector.hpp +++ b/tests/Unit/Helpers/Evolution/Imex/DoImplicitStepSector.hpp @@ -25,6 +25,7 @@ // DoImplicitStepInstantiate.cpp). namespace do_implicit_step_helpers { template +// [simple_sector] struct Sector : tt::ConformsTo { using tensors = tmpl::list; @@ -61,6 +62,7 @@ struct Sector : tt::ConformsTo { using solve_attempts = tmpl::list; }; +// [simple_sector] struct Var1 : db::SimpleTag { using type = Scalar; @@ -70,10 +72,15 @@ struct Var2 : db::SimpleTag { using type = Scalar; }; +// [ImexSystem] struct System : tt::ConformsTo { using variables_tag = Tags::Variables>; + + // Explicit evolution stuff here... + using implicit_sectors = tmpl::list, Sector>; }; +// [ImexSystem] struct NonautonomousSector : tt::ConformsTo { using tensors = tmpl::list;