Semi-implicit multistep ODE integrator \(
\begin{align}
v^{n+1} = \sum_{q=0}^{s-1} a_q v^{n-q} + \Delta t\left[\left(\sum_{q=0}^{s-1}b_q \hat E(t^{n}-q\Delta t, v^{n-q}) + \sum_{q=1}^{s} c_q \hat I( t^n - q\Delta t, v^{n-q})\right) + c_0\hat I(t^{n}+\Delta t, v^{n+1})\right]
\end{align}
\).
More...
|
| | ImExMultistep ()=default |
| | No memory allocation.
|
| |
| | ImExMultistep (ConvertsToMultistepTableau< value_type > tableau, const ContainerType ©able) |
| | Reserve memory for integration and construct Solver.
|
| |
| template<class ... Params> |
| void | construct (Params &&...ps) |
| | Perfect forward parameters to one of the constructors.
|
| |
| const ContainerType & | copyable () const |
| | Return an object of same size as the object used for construction.
|
| |
| template<class ExplicitRHS , class ImplicitRHS , class Solver > |
| void | init (const std::tuple< ExplicitRHS, ImplicitRHS, Solver > &ode, value_type t0, const ContainerType &u0, value_type dt) |
| | Initialize timestepper. Call before using the step function.
|
| |
| template<class ExplicitRHS , class ImplicitRHS , class Solver > |
| void | step (const std::tuple< ExplicitRHS, ImplicitRHS, Solver > &ode, value_type &t, ContainerType &u) |
| | Advance one timestep.
|
| |
template<class ContainerType>
struct dg::ImExMultistep< ContainerType >
Semi-implicit multistep ODE integrator \(
\begin{align}
v^{n+1} = \sum_{q=0}^{s-1} a_q v^{n-q} + \Delta t\left[\left(\sum_{q=0}^{s-1}b_q \hat E(t^{n}-q\Delta t, v^{n-q}) + \sum_{q=1}^{s} c_q \hat I( t^n - q\Delta t, v^{n-q})\right) + c_0\hat I(t^{n}+\Delta t, v^{n+1})\right]
\end{align}
\).
which discretizes
\[
\frac{\partial v}{\partial t} = \hat E(t,v) + \hat I(t,v)
\]
where \( \hat E \) contains the explicit and \( \hat I \) the implicit part of the equations. As an example, the coefficients for the 3-step, 3rd order "Karniadakis" scheme are
\[
a_0 = \frac{18}{11}\ a_1 = -\frac{9}{11}\ a_2 = \frac{2}{11} \\
b_0 = \frac{18}{11}\ b_1 = -\frac{18}{11}\ b_2 = \frac{6}{11} \\
c_0 = \frac{6}{11}\quad c_1 = c_2 = c_3 = 0 \\
\]
You can use your own coefficients defined as a dg::MultistepTableau or use one of the predefined coefficients in
We follow the naming convention as NAME-S-Q
- NAME is the author or name of the method
- S is the number of steps in the method
- Q is the global order of the method
| Name | Identifier | Description |
| ImEx-Euler-1-1 | dg::IMEX_EULER_1_1 | Explicit Euler combined with Implicit Euler |
| Euler | dg::IMEX_EULER_1_1 | For convenience |
| ImEx-Koto-2-2 | dg::IMEX_KOTO_2_2 | Koto T. Front. Math. China 2009, 4(1): 113-129 A stabilized 2nd order scheme with a large region of stability |
| ImEx-Adams-2-2 | dg::IMEX_ADAMS_2_2 | Hundsdorfer and Ruuth, Journal of Computational Physics 225 (2007) - Note
- (C=0.44)
|
| ImEx-BDF-X-X | dg::IMEX_BDF_X_X | The family of schems described in Hundsdorfer and Ruuth, Journal of Computational Physics 225 (2007)
The implicit part is a normal BDF scheme https://en.wikipedia.org/wiki/Backward_differentiation_formula while the explicit part equals the Minimal Projecting method by Alfeld, P., Math. Comput. 33.148 1195-1212 (1979) or extrapolated BDF in Hundsdorfer, W., Ruuth, S. J., & Spiteri, R. J. (2003). Monotonicity-preserving linear multistep methods. SIAM Journal on Numerical Analysis, 41(2), 605-623
- Note
- Possible values for X: 1 (C=1.00), 2 (C=0.63), 3 (C=0.39), 4 (C=0.22), 5 (C=0.09), 6
Note that X=3 is identical to the "Karniadakis" scheme
|
| Karniadakis | dg::IMEX_BDF_3_3 | The ImEx-BDF-3-3 scheme is identical to the widely used "Karniadakis" scheme Karniadakis, et al. J. Comput. Phys. 97 (1991) |
| ImEx-TVB-X-X | dg::IMEX_TVB_X_X | The family of schems described in < Hundsdorfer and Ruuth, Journal of Computational Physics 225 (2007)
The explicit part is a TVB scheme while the implicit part is optimized to maximize damping of high wavelength
- Note
- Possible values for X: 3 (C=0.54), 4 (C=0.46), 5 (C=0.38)
|
- Note
- the CFL coefficient C is given relative to the forward Euler method: \( \Delta t < C \Delta t_{FE}\).
- Attention
- The coefficient C is the one that ensures the TVD property of the scheme and is not directly related to the stability region of the scheme
The necessary Inversion in the implicit part is provided by the Implicit class.
The following code example demonstrates how to implement the method of manufactured solutions on a 2d partial differential equation with the dg library:
struct Explicit
{
Explicit( double nu): m_nu( nu) {}
void operator()( double t, const std::array<double,2>& T,
std::array<double,2>& Tp)
{
Tp[0] = m_nu*cos(t) - sin(t);
Tp[1] = m_nu*sin(t) + cos(t);
}
private:
double m_nu;
};
struct Implicit
{
Implicit( double nu): m_nu(nu) { }
void operator()( double t, const std::array<double,2>& T,
std::array<double,2>& Tp)
{
Tp[0] = -m_nu*T[0];
Tp[1] = -m_nu*T[1];
}
void operator()( double alpha, double t, std::array<double,2>& y, const
std::array<double,2>& rhs)
{
y[0] = rhs[0]/(1+alpha*m_nu);
y[1] = rhs[1]/(1+alpha*m_nu);
}
private:
double m_nu;
};
In the main function:
time = 0., y0 = init;
unsigned NT = NTs[u];
double dt = T/ (double)NT;
imex.init( std::tie(ex, im, im), time, y0, dt);
for( unsigned k=0; k<NT; k++)
imex.step( std::tie(ex, im, im), time, y0);
- Note
- In our experience the implicit treatment of diffusive or hyperdiffusive terms may significantly reduce the required number of time steps. This outweighs the increased computational cost of the additional matrix inversions. However, each PDE is different and general statements like this one should be treated with care.
- Note
- Uses only
blas1::axpby routines to integrate one step.
-
The difference between a multistep and a single step method like RungeKutta is that the multistep only takes one right-hand-side evaluation per step. This is advantageous if the right hand side is expensive to evaluate. Even though Runge Kutta methods can have a larger absolute timestep, if the effective timestep per rhs evaluation is compared, multistep methods generally win.
-
a disadvantage of multistep is that timestep adaption is not easily done.
- Template Parameters
-
| ContainerType | Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others
dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp) std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors) double (scalar) and other primitive types ...
If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them |
- See also
- See The dg dispatch system for a detailed explanation of our type dispatch system
1.12.0