GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/experimental/timestepping/multistagetimestepper.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 53 55 96.4%
Functions: 41 41 100.0%
Branches: 38 83 45.8%
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-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Experimental
10 * \brief A time stepper performing a single time step of a transient simulation
11 */
12 #ifndef DUMUX_TIMESTEPPING_MULTISTAGE_TIMESTEPPER_HH
13 #define DUMUX_TIMESTEPPING_MULTISTAGE_TIMESTEPPER_HH
14
15 #include <memory>
16 #include <vector>
17 #include <cmath>
18 #include <iostream>
19
20 #include <dumux/io/format.hh>
21
22 #include <dumux/common/variablesbackend.hh>
23 #include <dumux/common/timeloop.hh>
24 #include <dumux/experimental/timestepping/multistagemethods.hh>
25
26 namespace Dumux::Experimental {
27
28 //! Data object for the parameters of a given stage
29 template<class Scalar>
30 2082 class MultiStageParams
31 {
32 struct Params {
33 Scalar alpha, betaDt, timeAtStage, dtFraction;
34 bool skipTemporal, skipSpatial;
35 };
36 public:
37 //! Extract params for stage i from method m
38 2082 MultiStageParams(const MultiStageMethod<Scalar>& m, std::size_t i, const Scalar t, const Scalar dt)
39
1/2
✓ Branch 1 taken 2082 times.
✗ Branch 2 not taken.
 2082 : size_(i+1)
40 {
41
1/2
✓ Branch 1 taken 2082 times.
✗ Branch 2 not taken.
 2082 params_.resize(size_);
42
2/2
✓ Branch 0 taken 4173 times.
✓ Branch 1 taken 2082 times.
 6255 for (std::size_t k = 0; k < size_; ++k)
43 {
44
1/2
✓ Branch 1 taken 4173 times.
✗ Branch 2 not taken.
 4173 auto& p = params_[k];
45
1/2
✓ Branch 1 taken 4173 times.
✗ Branch 2 not taken.
 4173 p.alpha = m.temporalWeight(i, k);
46
1/2
✓ Branch 1 taken 4173 times.
✗ Branch 2 not taken.
 4173 p.betaDt = m.spatialWeight(i, k)*dt;
47
1/2
✓ Branch 1 taken 4173 times.
✗ Branch 2 not taken.
 4173 p.timeAtStage = t + m.timeStepWeight(k)*dt;
48
1/2
✓ Branch 1 taken 4173 times.
✗ Branch 2 not taken.
 4173 p.dtFraction = m.timeStepWeight(k);
49
50 using std::abs;
51 4173 p.skipTemporal = (abs(p.alpha) < 1e-6);
52 4173 p.skipSpatial = (abs(p.betaDt) < 1e-6);
53 }
54 2082 }
55
56 5081853 std::size_t size () const
57
7/9
✓ Branch 1 taken 2158 times.
✓ Branch 2 taken 2246 times.
✓ Branch 3 taken 2226 times.
✓ Branch 4 taken 2154 times.
✓ Branch 5 taken 20 times.
✓ Branch 6 taken 54 times.
✓ Branch 7 taken 20 times.
✗ Branch 8 not taken.
✗ Branch 0 not taken.
 5081853 { return size_; }
58
59 //! weights of the temporal operator residual (\f$ \alpha_{ik} \f$)
60 5075280 Scalar temporalWeight (std::size_t k) const
61
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39936 times.
 10118672 { return params_[k].alpha; }
62
63 //! weights of the spatial operator residual (\f$ \beta_{ik} \f$)
64 5074994 Scalar spatialWeight (std::size_t k) const
65
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39936 times.
 10118386 { return params_[k].betaDt; }
66
67 //! the time at which we have to evaluate the operators
68 55 Scalar timeAtStage (std::size_t k) const
69 4159 { return params_[k].timeAtStage; }
70
71 //! the fraction of a time step corresponding to the k-th stage
72 15 Scalar timeStepFraction (std::size_t k) const
73 15 { return params_[k].dtFraction; }
74
75 //! If \f$ \alpha_{ik} = 0\f$
76 2338 Scalar skipTemporal (std::size_t k) const
77
3/4
✓ Branch 0 taken 2266 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 54 times.
✗ Branch 3 not taken.
 2338 { return params_[k].skipTemporal; }
78
79 //! If \f$ \beta_{ik} = 0\f$
80 2338 Scalar skipSpatial (std::size_t k) const
81
3/4
✓ Branch 0 taken 2034 times.
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 54 times.
 2338 { return params_[k].skipSpatial; }
82
83 private:
84 std::size_t size_;
85 std::vector<Params> params_;
86 };
87
88 /*!
89 * \brief Time stepping with a multi-stage method
90 * \note We limit ourselves to "diagonally" implicit multi-stage methods where solving
91 * a stage can only depend on the values of the same stage and stages before
92 * but not future stages (which would require solving larger linear systems)
93 */
94 template<class PDESolver, class Scalar = double>
95 class MultiStageTimeStepper
96 {
97 using Variables = typename PDESolver::Variables;
98 using StageParams = MultiStageParams<Scalar>;
99 using Backend = VariablesBackend<Variables>;
100
101 public:
102
103 /*!
104 * \brief The constructor
105 * \param pdeSolver Solver class for solving a PDE in each stage
106 * \param msMethod The multi-stage method which is to be used for time integration
107 * \param paramGroup A parameter group in which we look for parameters
108 */
109 16 MultiStageTimeStepper(std::shared_ptr<PDESolver> pdeSolver,
110 std::shared_ptr<const MultiStageMethod<Scalar>> msMethod,
111 const std::string& paramGroup = "")
112 16 : pdeSolver_(pdeSolver)
113 16 , msMethod_(msMethod)
114 {
115
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 initParams_(paramGroup);
116
117
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 if (pdeSolver_->verbosity() >= 1)
118 16 std::cout << "Initialize time stepper with method " << msMethod_->name()
119
7/12
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 14 times.
✓ Branch 7 taken 2 times.
✓ Branch 9 taken 16 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 16 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 16 times.
✗ Branch 16 not taken.
 78 << Fmt::format(" ({} stage{})", msMethod_->numStages(), (msMethod_->numStages() > 1 ? "s" : ""))
120 16 << std::endl;
121 16 }
122
123 /*!
124 * \brief Advance one time step of the given time loop
125 * \param vars The variables object at the current time level.
126 * \param t The current time level
127 * \param dt The time step size to be used
128 * \note We expect the time level in vars to correspond to the given time `t`
129 */
130 2077 void step(Variables& vars, const Scalar t, const Scalar dt)
131 {
132 // make sure there are no traces of previous stages
133 2077 pdeSolver_->assembler().clearStages();
134
135 // do time integration
136 2077 bool converged = step_(vars, t, dt);
137
138 // clear traces of previously registered stages
139 2077 pdeSolver_->assembler().clearStages();
140
141 // if the solver didn't converge we can't recover
142
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2077 times.
 2077 if (!converged)
143 ✗ DUNE_THROW(NumericalProblem, "Solver did not converge!");
144 2077 }
145
146 /*!
147 * \brief Advance one time step of the given time loop (adaptive time stepping on solver failure)
148 * \param vars The variables object at the current time level.
149 * \param timeLoop An instance of a time loop
150 * \note We expect the time level in vars to correspond to the given time `t`
151 */
152 void step(Variables& vars, TimeLoopBase<Scalar>& timeLoop)
153 {
154 // make sure there are no traces of previous stages
155 pdeSolver_->assembler().clearStages();
156
157 // try solving the non-linear system
158 for (std::size_t i = 0; i <= maxTimeStepDivisions_; ++i)
159 {
160 // try solving the non-linear system
161 bool converged = step_(vars, timeLoop.time(), timeLoop.timeStepSize());
162 if (converged)
163 {
164 // clear traces of previously registered stages
165 pdeSolver_->assembler().clearStages();
166 return;
167 }
168
169 else if (!converged && i < maxTimeStepDivisions_)
170 {
171 // set solution to previous solution & reset time step
172 Backend::update(vars, pdeSolver_->assembler().prevSol());
173 pdeSolver_->assembler().resetTimeStep(Backend::dofs(vars));
174
175 if (pdeSolver_->verbosity() >= 1)
176 {
177 const auto dt = timeLoop.timeStepSize();
178 std::cout << Fmt::format("Solver did not converge with dt = {} seconds. ", dt)
179 << Fmt::format("Retrying with time step of dt = {} seconds.\n", dt*retryTimeStepReductionFactor_);
180 }
181
182 // try again with dt = dt * retryTimeStepReductionFactor_
183 timeLoop.setTimeStepSize(timeLoop.timeStepSize() * retryTimeStepReductionFactor_);
184 }
185
186 else
187 {
188 pdeSolver_->assembler().clearStages();
189 DUNE_THROW(NumericalProblem,
190 Fmt::format("Solver didn't converge after {} time-step divisions; dt = {}.\n",
191 maxTimeStepDivisions_, timeLoop.timeStepSize()));
192 }
193 }
194
195 DUNE_THROW(Dune::InvalidStateException, "Unreachable");
196 }
197
198 private:
199 2077 bool step_(Variables& vars, const Scalar t, const Scalar dt)
200 {
201
2/2
✓ Branch 1 taken 2082 times.
✓ Branch 2 taken 2077 times.
 4159 for (auto stageIdx = 1UL; stageIdx <= msMethod_->numStages(); ++stageIdx)
202 {
203 // prepare the assembler for this stage
204
1/2
✓ Branch 2 taken 2082 times.
✗ Branch 3 not taken.
 2092 pdeSolver_->assembler().prepareStage(
205 vars,
206 2082 std::make_shared<StageParams>(*msMethod_, stageIdx, t, dt)
207 );
208
209 // assemble & solve
210 2082 bool converged = pdeSolver_->apply(vars);
211
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2082 times.
 2082 if (!converged)
212 ✗ return false;
213 }
214
215 2077 return true;
216 }
217
218 16 void initParams_(const std::string& group = "")
219 {
220 16 maxTimeStepDivisions_ = getParamFromGroup<std::size_t>(group, "TimeStepper.MaxTimeStepDivisions", 10);
221 16 retryTimeStepReductionFactor_ = getParamFromGroup<Scalar>(group, "TimeStepper.RetryTimeStepReductionFactor", 0.5);
222 16 }
223
224 std::shared_ptr<PDESolver> pdeSolver_;
225 std::shared_ptr<const MultiStageMethod<Scalar>> msMethod_;
226
227 Scalar maxTimeStepDivisions_;
228 Scalar retryTimeStepReductionFactor_;
229 };
230
231 } // end namespace Dumux::Experimental
232
233 #endif
234