GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/experimental/timestepping/multistagemethods.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 58 63 92.1%
Functions: 19 22 86.4%
Branches: 8 112 7.1%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 // SPDX-FileCopyrightInfo: Copyright © dune-pdelab developers, see LICENSE.md (permalink below)
7 // SPDX-License-Identifier: GPL-3.0-or-later
8 // The code is based on the implementation of time stepping method parameters
9 // in dune-pdelab (https://archive.softwareheritage.org/swh:1:cnt:9c3d72412f8e4d48d84a0090b2bb4362b5a0d843)
10 // licensed under GPL-3.0-or-later, see their LICENSE.md for a full list of copyright holders at
11 // https://archive.softwareheritage.org/swh:1:cnt:b11b484e74eefe20c74f7043309b1b02853df2eb.
12 // Modifications (different interface naming, comments, types) are licensed under GPL-3.0-or-later.
13 //
14 /*!
15 * \file
16 * \ingroup Experimental
17 * \brief Parameters for different multistage time stepping methods
18 * \note See e.g. https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods
19 */
20 #ifndef DUMUX_TIMESTEPPING_MULTISTAGE_METHODS_HH
21 #define DUMUX_TIMESTEPPING_MULTISTAGE_METHODS_HH
22
23 #include <cmath>
24 #include <string>
25 #include <array>
26
27 namespace Dumux::Experimental {
28
29 /*!
30 * \brief Abstract interface for one-step multi-stage method parameters in Shu/Osher form.
31 *
32 * This implementation is based on the Shu/Osher form from:
33 * Chi W. Shu and Stanley Osher. Efficient implementation of essentially
34 * non- oscillatory shock-capturing schemes. J. Comput. Phys., 77:439–471, 1988.
35 * https://doi.org/10.1016/0021-9991(88)90177-5.
36 * To this end Eq. (2.12) is extended for implicit schemes.
37 *
38 * We consider the general PDE form
39 *
40 * \f[
41 * \begin{equation}
42 * \frac{\partial M(x, t)}{\partial t} - R(x, t) = 0,
43 * \end{equation}
44 * \f]
45 *
46 * where \f$ M(x, t)\f$ is the temporal operator/residual in terms of the solution \f$ x \f$,
47 * and \f$ R(x, t)\f$ is the discrete spatial operator/residual.
48 * \f$ M(x)\f$ usually corresponds to the conserved quantity (e.g. mass), and \f$ R(x)\f$
49 * contains the rest of the residual. We can then construct \f$ m \f$-stage time discretization methods.
50 * Integrating from time \f$ t^n\f$ to \f$ t^{n+1}\f$ with time step size \f$ \Delta t^n\f$, we solve:
51 *
52 * \f[
53 * \begin{aligned}
54 * x^{(0)} &= u^n\\
55 * \sum_{k=0}^i \left[ \alpha_{ik} M\left(x^{(k)}, t^n + d_k\Delta t^n\right)
56 * + \beta_{ik}\Delta t^n R \left(x^{(k)}, t^n+d_k\Delta t^n \right)\right]
57 * &= 0 & \forall i \in \{1,\ldots,m\} \\
58 * x^{n+1} &= x^{(m)}
59 * \end{aligned}
60 * \f]
61 * where \f$ x^{(k)} \f$ denotes the intermediate solution at stage \f$ k \f$.
62 * Dependent on the number of stages \f$ m \f$, and the coefficients \f$ \alpha, \beta, d\f$,
63 * schemes with different properties and order of accuracy can be constructed.
64 *
65 * That the summation only goes up to \f$ i \f$ in stage \f$ i \f$ means that we
66 * restrict ourselves to diagonally-implicit Runge-Kutta schemes (DIRK)
67 * and explicit schemes.
68 *
69 * Note that when computing the Jacobian of the residual with respect
70 * to \f$ x^{(k)} \f$ at stage \f$ k \f$, only the terms containing the solution of
71 * the current stage \f$ k \f$ contribute to the derivatives.
72 */
73 template<class Scalar>
74 class MultiStageMethod
75 {
76 public:
77 virtual bool implicit () const = 0;
78
79 virtual std::size_t numStages () const = 0;
80
81 //! weights of the temporal operator residual (\f$ \alpha_{ik} \f$)
82 virtual Scalar temporalWeight (std::size_t i, std::size_t k) const = 0;
83
84 //! weights of the spatial operator residual (\f$ \beta_{ik} \f$)
85 virtual Scalar spatialWeight (std::size_t i, std::size_t k) const = 0;
86
87 //! time step weights for each stage (\f$ d_k \f$)
88 virtual Scalar timeStepWeight (std::size_t k) const = 0;
89
90 virtual std::string name () const = 0;
91
92 ✗ virtual ~MultiStageMethod() = default;
93 };
94
95 //! Multi-stage time stepping scheme implementations
96 namespace MultiStage {
97
98 /*!
99 * \brief A theta time stepping scheme
100 * theta=1.0 is an implicit Euler scheme,
101 * theta=0.0 an explicit Euler scheme,
102 * theta=0.5 is a Cranck-Nicholson scheme
103 */
104 template<class Scalar>
105 13 class Theta : public MultiStageMethod<Scalar>
106 {
107 public:
108 14 explicit Theta(const Scalar theta)
109 : paramAlpha_{{-1.0, 1.0}}
110 1 , paramBeta_{{1.0-theta, theta}}
111 15 , paramD_{{0.0, 1.0}}
112 {}
113
114 10289041 bool implicit () const final
115 20578082 { return paramBeta_[1] > 0.0; }
116
117 4178 std::size_t numStages () const final
118 4178 { return 1; }
119
120 4150 Scalar temporalWeight (std::size_t, std::size_t k) const final
121 8300 { return paramAlpha_[k]; }
122
123 4150 Scalar spatialWeight (std::size_t, std::size_t k) const final
124 8300 { return paramBeta_[k]; }
125
126 8300 Scalar timeStepWeight (std::size_t k) const final
127 16600 { return paramD_[k]; }
128
129 2 std::string name () const override
130 4 { return "theta scheme"; }
131
132 private:
133 std::array<Scalar, 2> paramAlpha_;
134 std::array<Scalar, 2> paramBeta_;
135 std::array<Scalar, 2> paramD_;
136 };
137
138 /*!
139 * \brief An explicit Euler time stepping scheme
140 */
141 template<class Scalar>
142 6 class ExplicitEuler final : public Theta<Scalar>
143 {
144 public:
145 6 ExplicitEuler() : Theta<Scalar>(0.0) {}
146
147 3 std::string name () const final
148
2/28
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
 8 { return "explicit Euler"; }
149 };
150
151 /*!
152 * \brief An implicit Euler time stepping scheme
153 */
154 template<class Scalar>
155 20 class ImplicitEuler final : public Theta<Scalar>
156 {
157 public:
158 20 ImplicitEuler() : Theta<Scalar>(1.0) {}
159
160 10 std::string name () const final
161
2/28
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
 22 { return "implicit Euler"; }
162 };
163
164 /*!
165 * \brief Classical explicit fourth order Runge-Kutta scheme
166 */
167 template<class Scalar>
168 1 class RungeKuttaExplicitFourthOrder final : public MultiStageMethod<Scalar>
169 {
170 public:
171 1 RungeKuttaExplicitFourthOrder()
172 : paramAlpha_{{{-1.0, 1.0, 0.0, 0.0, 0.0},
173 {-1.0, 0.0, 1.0, 0.0, 0.0},
174 {-1.0, 0.0, 0.0, 1.0, 0.0},
175 {-1.0, 0.0, 0.0, 0.0, 1.0}}}
176 , paramBeta_{{{0.5, 0.0, 0.0, 0.0, 0.0},
177 {0.0, 0.5, 0.0, 0.0, 0.0},
178 {0.0, 0.0, 1.0, 0.0, 0.0},
179 {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0, 0.0}}}
180 1 , paramD_{{0.0, 0.5, 0.5, 1.0, 1.0}}
181 1 {}
182
183 ✗ bool implicit () const final
184 ✗ { return false; }
185
186 7 std::size_t numStages () const final
187 7 { return 4; }
188
189 14 Scalar temporalWeight (std::size_t i, std::size_t k) const final
190 14 { return paramAlpha_[i-1][k]; }
191
192 14 Scalar spatialWeight (std::size_t i, std::size_t k) const final
193 42 { return paramBeta_[i-1][k]; }
194
195 28 Scalar timeStepWeight (std::size_t k) const final
196 56 { return paramD_[k]; }
197
198 1 std::string name () const final
199
2/28
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
 4 { return "explicit Runge-Kutta 4th order"; }
200
201 private:
202 std::array<std::array<Scalar, 5>, 4> paramAlpha_;
203 std::array<std::array<Scalar, 5>, 4> paramBeta_;
204 std::array<Scalar, 5> paramD_;
205 };
206
207 /*!
208 * \brief Third order DIRK scheme
209 * \note see Alexander (1977) https://doi.org/10.1137/0714068 (Theorem 5)
210 */
211 template<class Scalar>
212 1 class DIRKThirdOrderAlexander final : public MultiStageMethod<Scalar>
213 {
214 public:
215 1 DIRKThirdOrderAlexander()
216 1 {
217 1 constexpr Scalar alpha = []{
218 // Newton iteration for alpha
219 Scalar alpha = 0.4358665215; // initial guess
220 for (int i = 0; i < 10; ++i)
221 alpha = alpha - (alpha*(alpha*alpha-3.0*(alpha-0.5))-1.0/6.0)/(3.0*alpha*(alpha-2.0)+1.5);
222 return alpha;
223 }();
224
225 1 constexpr Scalar tau2 = (1.0+alpha)*0.5;
226 1 constexpr Scalar b1 = -(6.0*alpha*alpha -16.0*alpha + 1.0)*0.25;
227 1 constexpr Scalar b2 = (6*alpha*alpha - 20.0*alpha + 5.0)*0.25;
228
229 1 paramD_ = {{0.0, alpha, tau2, 1.0}};
230 1 paramAlpha_ = {{
231 {-1.0, 1.0, 0.0, 0.0},
232 {-1.0, 0.0, 1.0, 0.0},
233 {-1.0, 0.0, 0.0, 1.0}
234 }};
235 1 paramBeta_ = {{
236 {0.0, alpha, 0.0, 0.0},
237 {0.0, tau2-alpha, alpha, 0.0},
238 {0.0, b1, b2, alpha}
239 }};
240 }
241
242 ✗ bool implicit () const final
243 ✗ { return true; }
244
245 6 std::size_t numStages () const final
246 6 { return 3; }
247
248 9 Scalar temporalWeight (std::size_t i, std::size_t k) const final
249 27 { return paramAlpha_[i-1][k]; }
250
251 9 Scalar spatialWeight (std::size_t i, std::size_t k) const final
252 27 { return paramBeta_[i-1][k]; }
253
254 18 Scalar timeStepWeight (std::size_t k) const final
255 36 { return paramD_[k]; }
256
257 1 std::string name () const final
258
2/28
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
 4 { return "diagonally implicit Runge-Kutta 3rd order (Alexander)"; }
259
260 private:
261 std::array<std::array<Scalar, 4>, 3> paramAlpha_;
262 std::array<std::array<Scalar, 4>, 3> paramBeta_;
263 std::array<Scalar, 4> paramD_;
264 };
265
266 } // end namespace MultiStage
267 } // end namespace Dumux::Experimental
268
269 #endif
270