| 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 | // | ||
| 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 |
1/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2082 times.
✗ Branch 5 not taken.
|
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/4✓ Branch 1 taken 2082 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 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 | 8346 | p.skipSpatial = (abs(p.betaDt) < 1e-6); | |
| 53 | } | ||
| 54 | 2082 | } | |
| 55 | |||
| 56 | ✗ | std::size_t size () const | |
| 57 | ✗ | { return size_; } | |
| 58 | |||
| 59 | //! weights of the temporal operator residual (\f$ \alpha_{ik} \f$) | ||
| 60 | Scalar temporalWeight (std::size_t k) const | ||
| 61 |
0/3✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
16144 | { return params_[k].alpha; } |
| 62 | |||
| 63 | //! weights of the spatial operator residual (\f$ \beta_{ik} \f$) | ||
| 64 | Scalar spatialWeight (std::size_t k) const | ||
| 65 |
0/5✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
15858 | { return params_[k].betaDt; } |
| 66 | |||
| 67 | //! the time at which we have to evaluate the operators | ||
| 68 | Scalar timeAtStage (std::size_t k) const | ||
| 69 | 6316 | { return params_[k].timeAtStage; } | |
| 70 | |||
| 71 | //! the fraction of a time step corresponding to the k-th stage | ||
| 72 | Scalar timeStepFraction (std::size_t k) const | ||
| 73 | 15 | { return params_[k].dtFraction; } | |
| 74 | |||
| 75 | //! If \f$ \alpha_{ik} = 0\f$ | ||
| 76 | Scalar skipTemporal (std::size_t k) const | ||
| 77 |
6/8✓ Branch 0 taken 2266 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 2266 times.
✓ Branch 3 taken 18 times.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 54 times.
✗ Branch 7 not taken.
|
4676 | { return params_[k].skipTemporal; } |
| 78 | |||
| 79 | //! If \f$ \beta_{ik} = 0\f$ | ||
| 80 | 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 | : pdeSolver_(pdeSolver) | ||
| 113 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
16 | , msMethod_(msMethod) |
| 114 | { | ||
| 115 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
16 | initParams_(paramGroup); |
| 116 | |||
| 117 |
2/4✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
|
32 | if (pdeSolver_->verbosity() >= 1) |
| 118 | 32 | std::cout << "Initialize time stepper with method " << msMethod_->name() | |
| 119 |
13/28✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 16 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 16 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 16 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 14 times.
✓ Branch 17 taken 2 times.
✓ Branch 19 taken 16 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 16 times.
✗ Branch 23 not taken.
✓ Branch 26 taken 16 times.
✗ Branch 27 not taken.
✓ Branch 29 taken 16 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 2 times.
✓ Branch 32 taken 14 times.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
|
94 | << Fmt::format(" ({} stage{})", msMethod_->numStages(), (msMethod_->numStages() > 1 ? "s" : "")) |
| 120 |
2/4✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 16 times.
|
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 | 6231 | 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 | 6231 | 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 2 taken 2082 times.
✓ Branch 3 taken 2077 times.
|
4159 | for (auto stageIdx = 1UL; stageIdx <= msMethod_->numStages(); ++stageIdx) |
| 202 | { | ||
| 203 | // prepare the assembler for this stage | ||
| 204 |
4/13✓ Branch 3 taken 2082 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 2082 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 10 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 10 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
4174 | pdeSolver_->assembler().prepareStage( |
| 205 | vars, | ||
| 206 | 4164 | std::make_shared<StageParams>(*msMethod_, stageIdx, t, dt) | |
| 207 | ); | ||
| 208 | |||
| 209 | // assemble & solve | ||
| 210 | 4164 | 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 | ✗ | void initParams_(const std::string& group = "") | |
| 219 | { | ||
| 220 | ✗ | maxTimeStepDivisions_ = getParamFromGroup<std::size_t>(group, "TimeStepper.MaxTimeStepDivisions", 10); | |
| 221 | ✗ | retryTimeStepReductionFactor_ = getParamFromGroup<Scalar>(group, "TimeStepper.RetryTimeStepReductionFactor", 0.5); | |
| 222 | ✗ | } | |
| 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 |