GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/nonlinear/newtonsolver.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 316 365 86.6%
Functions: 5033 9717 51.8%
Branches: 436 945 46.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 //
7 /*!
8 * \file
9 * \ingroup Nonlinear
10 * \brief Reference implementation of a Newton solver.
11 */
12 #ifndef DUMUX_NEWTON_SOLVER_HH
13 #define DUMUX_NEWTON_SOLVER_HH
14
15 #include <cmath>
16 #include <memory>
17 #include <iostream>
18 #include <type_traits>
19 #include <algorithm>
20 #include <numeric>
21
22 #include <dune/common/timer.hh>
23 #include <dune/common/exceptions.hh>
24 #include <dune/common/parallel/mpicommunication.hh>
25 #include <dune/common/parallel/mpihelper.hh>
26 #include <dune/common/std/type_traits.hh>
27 #include <dune/common/indices.hh>
28 #include <dune/common/hybridutilities.hh>
29
30 #include <dune/istl/bvector.hh>
31 #include <dune/istl/multitypeblockvector.hh>
32
33 #include <dumux/common/parameters.hh>
34 #include <dumux/common/exceptions.hh>
35 #include <dumux/common/typetraits/vector.hh>
36 #include <dumux/common/typetraits/isvalid.hh>
37 #include <dumux/common/timeloop.hh>
38 #include <dumux/common/pdesolver.hh>
39 #include <dumux/common/variablesbackend.hh>
40
41 #include <dumux/io/format.hh>
42
43 #include <dumux/linear/matrixconverter.hh>
44 #include <dumux/assembly/partialreassembler.hh>
45
46 #include "newtonconvergencewriter.hh"
47 #include "primaryvariableswitchadapter.hh"
48
49 namespace Dumux::Detail::Newton {
50
51 // Helper boolean that states if the assembler exports grid variables
52 template<class Assembler> using AssemblerGridVariablesType = typename Assembler::GridVariables;
53 template<class Assembler>
54 inline constexpr bool assemblerExportsGridVariables
55 = Dune::Std::is_detected_v<AssemblerGridVariablesType, Assembler>;
56
57 // helper struct to define the variables on which the privarswitch should operate
58 template<class Assembler, bool exportsGridVars = assemblerExportsGridVariables<Assembler>>
59 struct PriVarSwitchVariablesType { using Type = typename Assembler::GridVariables; };
60
61 // if assembler does not export them, use an empty class. These situations either mean
62 // that there is no privarswitch, or, it is handled by a derived implementation.
63 template<class Assembler>
64 struct PriVarSwitchVariablesType<Assembler, false>
65 { using Type = struct EmptyGridVariables {}; };
66
67 // Helper alias to deduce the variables types used in the privarswitch adapter
68 template<class Assembler>
69 using PriVarSwitchVariables
70 = typename PriVarSwitchVariablesType<Assembler, assemblerExportsGridVariables<Assembler>>::Type;
71
72 //! helper struct detecting if an assembler supports partial reassembly
73 struct supportsPartialReassembly
74 {
75 template<class Assembler>
76 auto operator()(Assembler&& a)
77 -> decltype(a.assembleJacobianAndResidual(std::declval<const typename Assembler::SolutionVector&>(),
78 std::declval<const PartialReassembler<Assembler>*>()))
79 {}
80 };
81
82 // helpers to implement max relative shift
83 template<class C> using dynamicIndexAccess = decltype(std::declval<C>()[0]);
84 template<class C> using staticIndexAccess = decltype(std::declval<C>()[Dune::Indices::_0]);
85 template<class C> static constexpr auto hasDynamicIndexAccess = Dune::Std::is_detected<dynamicIndexAccess, C>{};
86 template<class C> static constexpr auto hasStaticIndexAccess = Dune::Std::is_detected<staticIndexAccess, C>{};
87
88 template<class V, class Scalar, class Reduce, class Transform>
89 ✗ auto hybridInnerProduct(const V& v1, const V& v2, Scalar init, Reduce&& r, Transform&& t)
90 -> std::enable_if_t<hasDynamicIndexAccess<V>(), Scalar>
91 {
92 ✗ return std::inner_product(v1.begin(), v1.end(), v2.begin(), init, std::forward<Reduce>(r), std::forward<Transform>(t));
93 }
94
95 template<class V, class Scalar, class Reduce, class Transform>
96 auto hybridInnerProduct(const V& v1, const V& v2, Scalar init, Reduce&& r, Transform&& t)
97 -> std::enable_if_t<hasStaticIndexAccess<V>() && !hasDynamicIndexAccess<V>(), Scalar>
98 {
99 using namespace Dune::Hybrid;
100 41516 forEach(std::make_index_sequence<V::N()>{}, [&](auto i){
101 42237 init = r(init, hybridInnerProduct(v1[Dune::index_constant<i>{}], v2[Dune::index_constant<i>{}], init, std::forward<Reduce>(r), std::forward<Transform>(t)));
102 });
103 6679 return init;
104 }
105
106 // Maximum relative shift at a degree of freedom.
107 // For (primary variables) values below 1.0 we use
108 // an absolute shift.
109 template<class Scalar, class V>
110 ✗ auto maxRelativeShift(const V& v1, const V& v2)
111 -> std::enable_if_t<Dune::IsNumber<V>::value, Scalar>
112 {
113 using std::abs; using std::max;
114
20/36
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 55445702 times.
✓ Branch 7 taken 50422706 times.
✓ Branch 8 taken 55445702 times.
✓ Branch 9 taken 50422706 times.
✓ Branch 10 taken 62678270 times.
✓ Branch 11 taken 43190189 times.
✓ Branch 12 taken 31 times.
✓ Branch 13 taken 20 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 605554 times.
✓ Branch 17 taken 3921831 times.
✓ Branch 18 taken 605554 times.
✓ Branch 19 taken 3921831 times.
✓ Branch 20 taken 3644972 times.
✓ Branch 21 taken 882413 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 16416 times.
✓ Branch 27 taken 68704 times.
✓ Branch 28 taken 16416 times.
✓ Branch 29 taken 68704 times.
✓ Branch 30 taken 77824 times.
✓ Branch 31 taken 7296 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
 277029600 return abs(v1 - v2)/max<Scalar>(1.0, abs(v1 + v2)*0.5);
115 }
116
117 // Maximum relative shift for generic vector types.
118 // Recursively calls maxRelativeShift until Dune::IsNumber is true.
119 template<class Scalar, class V>
120 5958 auto maxRelativeShift(const V& v1, const V& v2)
121 -> std::enable_if_t<!Dune::IsNumber<V>::value, Scalar>
122 {
123 43235931 return hybridInnerProduct(v1, v2, Scalar(0.0),
124
18/24
✓ Branch 0 taken 46639358 times.
✓ Branch 1 taken 40830828 times.
✓ Branch 2 taken 3016136 times.
✓ Branch 3 taken 31125891 times.
✓ Branch 4 taken 13220837 times.
✓ Branch 5 taken 2154721 times.
✓ Branch 6 taken 6199126 times.
✓ Branch 7 taken 11162818 times.
✓ Branch 8 taken 412178 times.
✓ Branch 9 taken 4151882 times.
✓ Branch 10 taken 211935 times.
✓ Branch 11 taken 3520652 times.
✓ Branch 12 taken 9289 times.
✓ Branch 13 taken 329887 times.
✓ Branch 14 taken 672 times.
✓ Branch 15 taken 38240 times.
✓ Branch 16 taken 4751 times.
✓ Branch 17 taken 80369 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
 166412564 [](const auto& a, const auto& b){ using std::max; return max(a, b); },
125
14/14
✓ Branch 0 taken 55445702 times.
✓ Branch 1 taken 50422706 times.
✓ Branch 2 taken 605554 times.
✓ Branch 3 taken 7048461 times.
✓ Branch 4 taken 43700981 times.
✓ Branch 5 taken 105489 times.
✓ Branch 6 taken 1898723 times.
✓ Branch 7 taken 1990 times.
✓ Branch 8 taken 149632 times.
✓ Branch 9 taken 3592221 times.
✓ Branch 10 taken 672 times.
✓ Branch 11 taken 38240 times.
✓ Branch 13 taken 4751 times.
✓ Branch 14 taken 80369 times.
 163095491 [](const auto& a, const auto& b){ return maxRelativeShift<Scalar>(a, b); }
126 5958 );
127 }
128
129 template<class To, class From>
130 ✗ void assign(To& to, const From& from)
131 {
132 if constexpr (hasStaticIndexAccess<To>() && hasStaticIndexAccess<To>() && !hasDynamicIndexAccess<From>() && !hasDynamicIndexAccess<From>())
133 {
134 using namespace Dune::Hybrid;
135 4180 forEach(std::make_index_sequence<To::N()>{}, [&](auto i){
136 158 assign(to[Dune::index_constant<i>{}], from[Dune::index_constant<i>{}]);
137 });
138 }
139
140 else if constexpr (std::is_assignable<To&, From>::value)
141 1606188 to = from;
142
143 else if constexpr (hasDynamicIndexAccess<To>() && hasDynamicIndexAccess<From>())
144
2/4
✓ Branch 0 taken 1601395 times.
✓ Branch 1 taken 3283 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1604678 for (decltype(to.size()) i = 0; i < to.size(); ++i)
145 4804185 assign(to[i], from[i]);
146
147 else if constexpr (hasDynamicIndexAccess<To>() && Dune::IsNumber<From>::value)
148 {
149 assert(to.size() == 1);
150 assign(to[0], from);
151 }
152
153 else if constexpr (Dune::IsNumber<To>::value && hasDynamicIndexAccess<From>())
154 {
155 assert(from.size() == 1);
156 assign(to, from[0]);
157 }
158
159 else
160 DUNE_THROW(Dune::Exception, "Values are not assignable to each other!");
161 ✗ }
162
163 } // end namespace Dumux::Detail::Newton
164
165 namespace Dumux {
166
167 /*!
168 * \ingroup Nonlinear
169 * \brief An implementation of a Newton solver
170 * \tparam Assembler the assembler
171 * \tparam LinearSolver the linear solver
172 * \tparam Comm the communication object used to communicate with all processes
173 * \note If you want to specialize only some methods but are happy with the
174 * defaults of the reference solver, derive your solver from
175 * this class and simply overload the required methods.
176 */
177 template <class Assembler, class LinearSolver,
178 class Reassembler = PartialReassembler<Assembler>,
179 class Comm = Dune::Communication<Dune::MPIHelper::MPICommunicator> >
180 class NewtonSolver : public PDESolver<Assembler, LinearSolver>
181 {
182 using ParentType = PDESolver<Assembler, LinearSolver>;
183
184 protected:
185 using Backend = VariablesBackend<typename ParentType::Variables>;
186 using SolutionVector = typename Backend::DofVector;
187 using ResidualVector = typename Assembler::ResidualType;
188 using LinearAlgebraNativeBackend = VariablesBackend<ResidualVector>;
189 private:
190 using Scalar = typename Assembler::Scalar;
191 using JacobianMatrix = typename Assembler::JacobianMatrix;
192 using ConvergenceWriter = ConvergenceWriterInterface<SolutionVector, ResidualVector>;
193 using TimeLoop = TimeLoopBase<Scalar>;
194
195 // enable models with primary variable switch
196 // TODO: Always use ParentType::Variables once we require assemblers to export variables
197 static constexpr bool assemblerExportsVariables = Detail::PDESolver::assemblerExportsVariables<Assembler>;
198 using PriVarSwitchVariables
199 = std::conditional_t<assemblerExportsVariables,
200 typename ParentType::Variables,
201 Detail::Newton::PriVarSwitchVariables<Assembler>>;
202 using PrimaryVariableSwitchAdapter = Dumux::PrimaryVariableSwitchAdapter<PriVarSwitchVariables>;
203
204 public:
205 using typename ParentType::Variables;
206 using Communication = Comm;
207
208 495 NewtonSolver(std::shared_ptr<Assembler> assembler,
209 std::shared_ptr<LinearSolver> linearSolver,
210
3/8
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 311 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 324 const Communication& comm = Dune::MPIHelper::getCommunication(),
211 const std::string& paramGroup = "",
212 const std::string& paramGroupName = "Newton",
213 int verbosity = 2)
214 : ParentType(assembler, linearSolver)
215 , endIterMsgStream_(std::ostringstream::out)
216 , comm_(comm)
217 , paramGroup_(paramGroup)
218 , solverName_(paramGroupName)
219
10/28
✓ Branch 3 taken 479 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 479 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 479 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 479 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 479 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 479 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 479 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 479 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 452 times.
✓ Branch 26 taken 27 times.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
 1980 , priVarSwitchAdapter_(std::make_unique<PrimaryVariableSwitchAdapter>(paramGroup))
220 {
221
9/14
✓ Branch 0 taken 452 times.
✓ Branch 1 taken 27 times.
✓ Branch 2 taken 452 times.
✓ Branch 3 taken 27 times.
✓ Branch 5 taken 452 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 452 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 452 times.
✓ Branch 11 taken 27 times.
✓ Branch 12 taken 452 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 1458 verbosity_ = comm_.rank() == 0 ? getParamFromGroup<int>(paramGroup, solverName_ + ".Verbosity", verbosity) : 0;
222
223
1/2
✓ Branch 1 taken 479 times.
✗ Branch 2 not taken.
 495 initParams_(paramGroup);
224
225 // set the linear system (matrix & residual) in the assembler
226
2/4
✓ Branch 1 taken 479 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 474 times.
✗ Branch 5 not taken.
 985 this->assembler().setLinearSystem();
227
228 // set a different default for the linear solver residual reduction
229 // within the Newton the linear solver doesn't need to solve too exact
230
3/9
✓ Branch 1 taken 479 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 474 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 990 this->linearSolver().setResidualReduction(getParamFromGroup<Scalar>(paramGroup, "LinearSolver.ResidualReduction", 1e-6));
231
232 // initialize the partial reassembler
233
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 457 times.
 495 if (enablePartialReassembly_)
234
4/8
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 22 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 22 times.
 44 partialReassembler_ = std::make_unique<Reassembler>(this->assembler());
235 495 }
236
237 //! the communicator for parallel runs
238 const Communication& comm() const
239 { return comm_; }
240
241 /*!
242 * \brief Set the maximum acceptable difference of any primary variable
243 * between two iterations for declaring convergence.
244 *
245 * \param tolerance The maximum relative shift between two Newton
246 * iterations at which the scheme is considered finished
247 */
248 ✗ void setMaxRelativeShift(Scalar tolerance)
249 479 { shiftTolerance_ = tolerance; }
250
251 /*!
252 * \brief Set the maximum acceptable absolute residual for declaring convergence.
253 *
254 * \param tolerance The maximum absolute residual at which
255 * the scheme is considered finished
256 */
257 ✗ void setMaxAbsoluteResidual(Scalar tolerance)
258 479 { residualTolerance_ = tolerance; }
259
260 /*!
261 * \brief Set the maximum acceptable residual norm reduction.
262 *
263 * \param tolerance The maximum reduction of the residual norm
264 * at which the scheme is considered finished
265 */
266 ✗ void setResidualReduction(Scalar tolerance)
267 479 { reductionTolerance_ = tolerance; }
268
269 /*!
270 * \brief Set the number of iterations at which the Newton method
271 * should aim at.
272 *
273 * This is used to control the time-step size. The heuristic used
274 * is to scale the last time-step size by the deviation of the
275 * number of iterations used from the target steps.
276 *
277 * \param targetSteps Number of iterations which are considered "optimal"
278 */
279 ✗ void setTargetSteps(int targetSteps)
280 479 { targetSteps_ = targetSteps; }
281
282 /*!
283 * \brief Set the number of minimum iterations for the Newton
284 * method.
285 *
286 * \param minSteps Minimum number of iterations
287 */
288 ✗ void setMinSteps(int minSteps)
289 479 { minSteps_ = minSteps; }
290
291 /*!
292 * \brief Set the number of iterations after which the Newton
293 * method gives up.
294 *
295 * \param maxSteps Number of iterations after we give up
296 */
297 ✗ void setMaxSteps(int maxSteps)
298 479 { maxSteps_ = maxSteps; }
299
300 /*!
301 * \brief Run the Newton method to solve a non-linear system.
302 * Does time step control when the Newton fails to converge
303 * \param vars The variables object representing the current state of the
304 * numerical solution (primary and possibly secondary variables).
305 * \param timeLoop The time loop.
306 */
307 14285 void solve(Variables& vars, TimeLoop& timeLoop) override
308 {
309 if constexpr (!assemblerExportsVariables)
310 {
311
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 14260 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 14260 times.
 28540 if (this->assembler().isStationaryProblem())
312 ✗ DUNE_THROW(Dune::InvalidStateException, "Using time step control with stationary problem makes no sense!");
313 }
314
315 // try solving the non-linear system
316
1/2
✓ Branch 0 taken 14551 times.
✗ Branch 1 not taken.
 14574 for (std::size_t i = 0; i <= maxTimeStepDivisions_; ++i)
317 {
318 // linearize & solve
319 14574 const bool converged = solve_(vars);
320
321
2/2
✓ Branch 0 taken 304 times.
✓ Branch 1 taken 14262 times.
 14574 if (converged)
322 return;
323
324
1/2
✓ Branch 0 taken 289 times.
✗ Branch 1 not taken.
 289 else if (!converged && i < maxTimeStepDivisions_)
325 {
326 if constexpr (assemblerExportsVariables)
327 ✗ DUNE_THROW(Dune::NotImplemented, "Time step reset for new assembly methods");
328 else
329 {
330 // set solution to previous solution & reset time step
331
6/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 2 times.
 582 Backend::update(vars, this->assembler().prevSol());
332
3/4
✓ Branch 0 taken 272 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 270 times.
✗ Branch 3 not taken.
 574 this->assembler().resetTimeStep(Backend::dofs(vars));
333
334
2/2
✓ Branch 0 taken 287 times.
✓ Branch 1 taken 2 times.
 289 if (verbosity_ >= 1)
335 {
336 287 const auto dt = timeLoop.timeStepSize();
337 287 std::cout << Fmt::format("{} solver did not converge with dt = {} seconds. ", solverName_, dt)
338
5/14
✓ Branch 3 taken 287 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 287 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 287 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 287 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 287 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 1024 << Fmt::format("Retrying with time step of dt = {} seconds.\n", dt*retryTimeStepReductionFactor_);
339 }
340
341 // try again with dt = dt * retryTimeStepReductionFactor_
342 289 timeLoop.setTimeStepSize(timeLoop.timeStepSize() * retryTimeStepReductionFactor_);
343 }
344 }
345
346 else
347 {
348 ✗ DUNE_THROW(NumericalProblem,
349 Fmt::format("{} solver didn't converge after {} time-step divisions; dt = {}.\n",
350 solverName_, maxTimeStepDivisions_, timeLoop.timeStepSize()));
351 }
352 }
353 }
354
355 /*!
356 * \brief Run the Newton method to solve a non-linear system.
357 * The solver is responsible for all the strategic decisions.
358 * \param vars The variables object representing the current state of the
359 * numerical solution (primary and possibly secondary variables).
360 */
361 445 void solve(Variables& vars) override
362 {
363 445 const bool converged = solve_(vars);
364
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 417 times.
 445 if (!converged)
365
8/28
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 3 times.
✗ Branch 12 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 3 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 3 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 3 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 3 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 3 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✓ Branch 31 taken 3 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
 36 DUNE_THROW(NumericalProblem,
366 Fmt::format("{} solver didn't converge after {} iterations.\n", solverName_, numSteps_));
367 442 }
368
369 /*!
370 * \brief Run the Newton method to solve a non-linear system.
371 * The solver is responsible for all the strategic decisions.
372 * \param vars The variables object representing the current state of the
373 * numerical solution (primary and possibly secondary variables).
374 * \post If converged, the `Variables` will represent the solution. If convergence
375 * fails, they are in some intermediate, undefined state.
376 */
377 78 bool apply(Variables& vars) override
378 {
379 78 return solve_(vars);
380 }
381
382 /*!
383 * \brief Called before the Newton method is applied to an
384 * non-linear system of equations.
385 *
386 * \param initVars The variables representing the initial solution
387 */
388 15097 virtual void newtonBegin(Variables& initVars)
389 {
390 15097 numSteps_ = 0;
391
392 if constexpr (hasPriVarsSwitch<PriVarSwitchVariables>)
393 {
394 if constexpr (assemblerExportsVariables)
395 priVarSwitchAdapter_->initialize(Backend::dofs(initVars), initVars);
396 else // this assumes assembly with solution (i.e. Variables=SolutionVector)
397 12816 priVarSwitchAdapter_->initialize(initVars, this->assembler().gridVariables());
398 }
399
400
401
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 25 times.
 15097 const auto& initSol = Backend::dofs(initVars);
402
403 // write the initial residual if a convergence writer was set
404
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 15033 times.
 15097 if (convergenceWriter_)
405 {
406 62 this->assembler().assembleResidual(initVars);
407
408 // dummy vector, there is no delta before solving the linear system
409
1/4
✓ Branch 2 taken 31 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 124 ResidualVector delta = LinearAlgebraNativeBackend::zeros(Backend::size(initSol));
410
4/8
✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 31 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 31 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 31 times.
✗ Branch 11 not taken.
 124 convergenceWriter_->write(initSol, delta, this->assembler().residual());
411 }
412
413
2/2
✓ Branch 0 taken 1291 times.
✓ Branch 1 taken 13773 times.
 15097 if (enablePartialReassembly_)
414 {
415 2582 partialReassembler_->resetColors();
416 1291 resizeDistanceFromLastLinearization_(initSol, distanceFromLastLinearization_);
417 }
418 15097 }
419
420 /*!
421 * \brief Returns true if another iteration should be done.
422 *
423 * \param varsCurrentIter The variables of the current Newton iteration
424 * \param converged if the Newton method's convergence criterion was met in this step
425 */
426 77353 virtual bool newtonProceed(const Variables &varsCurrentIter, bool converged)
427 {
428
2/2
✓ Branch 0 taken 47094 times.
✓ Branch 1 taken 30117 times.
 77353 if (numSteps_ < minSteps_)
429 return true;
430
2/2
✓ Branch 0 taken 32322 times.
✓ Branch 1 taken 14772 times.
 47174 else if (converged)
431 return false; // we are below the desired tolerance
432
2/2
✓ Branch 0 taken 228 times.
✓ Branch 1 taken 32094 times.
 32369 else if (numSteps_ >= maxSteps_)
433 {
434 // We have exceeded the allowed number of steps. If the
435 // maximum relative shift was reduced by a factor of at least 4,
436 // we proceed even if we are above the maximum number of steps.
437
1/2
✓ Branch 0 taken 228 times.
✗ Branch 1 not taken.
 228 if (enableShiftCriterion_)
438 228 return shift_*4.0 < lastShift_;
439 else
440 ✗ return reduction_*4.0 < lastReduction_;
441 }
442
443 return true;
444 }
445
446 /*!
447 * \brief Indicates the beginning of a Newton iteration.
448 */
449 59412 virtual void newtonBeginStep(const Variables& vars)
450 {
451 62295 lastShift_ = shift_;
452
4/4
✓ Branch 0 taken 15051 times.
✓ Branch 1 taken 47135 times.
✓ Branch 2 taken 13 times.
✓ Branch 3 taken 25 times.
 62295 if (numSteps_ == 0)
453 {
454 15084 lastReduction_ = 1.0;
455 }
456 else
457 {
458 47211 lastReduction_ = reduction_;
459 }
460 59412 }
461
462 /*!
463 * \brief Assemble the linear system of equations \f$\mathbf{A}x - b = 0\f$.
464 *
465 * \param vars The current iteration's variables
466 */
467 62333 virtual void assembleLinearSystem(const Variables& vars)
468 {
469 62403 assembleLinearSystem_(this->assembler(), vars);
470
471
2/2
✓ Branch 0 taken 5089 times.
✓ Branch 1 taken 54144 times.
 59357 if (enablePartialReassembly_)
472 10178 partialReassembler_->report(comm_, endIterMsgStream_);
473 62333 }
474
475 /*!
476 * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
477 *
478 * Throws Dumux::NumericalProblem if the linear solver didn't
479 * converge.
480 *
481 * If the linear solver doesn't accept multitype matrices we copy the matrix
482 * into a 1x1 block BCRS matrix for solving.
483 *
484 * \param deltaU The difference between the current and the next solution
485 */
486 62333 void solveLinearSystem(ResidualVector& deltaU)
487 {
488 62333 bool converged = false;
489
490 try
491 {
492
2/2
✓ Branch 0 taken 15064 times.
✓ Branch 1 taken 47160 times.
 62333 if (numSteps_ == 0)
493
4/8
✓ Branch 1 taken 15046 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15046 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 15046 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 15046 times.
✗ Branch 11 not taken.
 60406 initialResidual_ = this->linearSolver().norm(this->assembler().residual());
494
495 // solve by calling the appropriate implementation depending on whether the linear solver
496 // is capable of handling MultiType matrices or not
497
2/2
✓ Branch 1 taken 62196 times.
✓ Branch 2 taken 28 times.
 62333 converged = solveLinearSystem_(deltaU);
498 }
499 28 catch (const Dune::Exception &e)
500 {
501
2/2
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 2 times.
 28 if (verbosity_ >= 1)
502
3/6
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 26 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 26 times.
✗ Branch 10 not taken.
 52 std::cout << solverName_ << ": Caught exception from the linear solver: \"" << e.what() << "\"\n";
503
504 28 converged = false;
505 }
506
507 // make sure all processes converged
508 62333 int convergedRemote = converged;
509
4/4
✓ Branch 0 taken 5932 times.
✓ Branch 1 taken 56292 times.
✓ Branch 2 taken 5932 times.
✓ Branch 3 taken 56292 times.
 124666 if (comm_.size() > 1)
510 5932 convergedRemote = comm_.min(converged);
511
512
2/2
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 62172 times.
 62333 if (!converged)
513 {
514
7/16
✓ Branch 2 taken 52 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 52 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 52 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 52 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 52 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 52 times.
✗ Branch 24 not taken.
✓ Branch 27 taken 52 times.
✗ Branch 28 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
 520 DUNE_THROW(NumericalProblem, "Linear solver did not converge");
515 ++numLinearSolverBreakdowns_;
516 }
517
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 62171 times.
 62281 else if (!convergedRemote)
518 {
519
7/16
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 1 times.
✗ Branch 24 not taken.
✓ Branch 27 taken 1 times.
✗ Branch 28 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
 10 DUNE_THROW(NumericalProblem, "Linear solver did not converge on a remote process");
520 ++numLinearSolverBreakdowns_; // we keep correct count for process 0
521 }
522 62280 }
523
524 /*!
525 * \brief Update the current solution with a delta vector.
526 *
527 * The error estimates required for the newtonConverged() and
528 * newtonProceed() methods should be updated inside this method.
529 *
530 * Different update strategies, such as line search and chopped
531 * updates can be implemented. The default behavior is just to
532 * subtract deltaU from uLastIter, i.e.
533 * \f[ u^{k+1} = u^k - \Delta u^k \f]
534 *
535 * \param vars The variables after the current iteration
536 * \param uLastIter The solution vector after the last iteration
537 * \param deltaU The delta as calculated from solving the linear
538 * system of equations. This parameter also stores
539 * the updated solution.
540 */
541 62280 void newtonUpdate(Variables& vars,
542 const SolutionVector& uLastIter,
543 const ResidualVector& deltaU)
544 {
545
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 62171 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 62280 if (enableShiftCriterion_ || enablePartialReassembly_)
546 62280 newtonUpdateShift_(uLastIter, deltaU);
547
548
2/2
✓ Branch 0 taken 5089 times.
✓ Branch 1 taken 57082 times.
 62280 if (enablePartialReassembly_) {
549 // Determine the threshold 'eps' that is used for the partial reassembly.
550 // Every entity where the primary variables exhibit a relative shift
551 // summed up since the last linearization above 'eps' will be colored
552 // red yielding a reassembly.
553 // The user can provide three parameters to influence the threshold:
554 // 'minEps' by 'Newton.ReassemblyMinThreshold' (1e-1*shiftTolerance_ default)
555 // 'maxEps' by 'Newton.ReassemblyMaxThreshold' (1e2*shiftTolerance_ default)
556 // 'omega' by 'Newton.ReassemblyShiftWeight' (1e-3 default)
557 // The threshold is calculated from the currently achieved maximum
558 // relative shift according to the formula
559 // eps = max( minEps, min(maxEps, omega*shift) ).
560 // Increasing/decreasing 'minEps' leads to less/more reassembly if
561 // 'omega*shift' is small, i.e., for the last Newton iterations.
562 // Increasing/decreasing 'maxEps' leads to less/more reassembly if
563 // 'omega*shift' is large, i.e., for the first Newton iterations.
564 // Increasing/decreasing 'omega' leads to more/less first and last
565 // iterations in this sense.
566 using std::max;
567 using std::min;
568
2/2
✓ Branch 0 taken 2953 times.
✓ Branch 1 taken 2136 times.
 5089 auto reassemblyThreshold = max(reassemblyMinThreshold_,
569 5089 min(reassemblyMaxThreshold_,
570
2/2
✓ Branch 0 taken 2797 times.
✓ Branch 1 taken 2292 times.
 5089 shift_*reassemblyShiftWeight_));
571
572 5089 updateDistanceFromLastLinearization_(uLastIter, deltaU);
573 10178 partialReassembler_->computeColors(this->assembler(),
574 5089 distanceFromLastLinearization_,
575 reassemblyThreshold);
576
577 // set the discrepancy of the red entities to zero
578
4/4
✓ Branch 0 taken 5877638 times.
✓ Branch 1 taken 5089 times.
✓ Branch 2 taken 5877638 times.
✓ Branch 3 taken 5089 times.
 5882727 for (unsigned int i = 0; i < distanceFromLastLinearization_.size(); i++)
579
6/6
✓ Branch 0 taken 2417136 times.
✓ Branch 1 taken 3460502 times.
✓ Branch 2 taken 2417136 times.
✓ Branch 3 taken 3460502 times.
✓ Branch 4 taken 2417136 times.
✓ Branch 5 taken 3460502 times.
 17632914 if (partialReassembler_->dofColor(i) == EntityColor::red)
580 4834272 distanceFromLastLinearization_[i] = 0;
581 }
582
583
2/2
✓ Branch 0 taken 1096 times.
✓ Branch 1 taken 61075 times.
 62280 if (useLineSearch_)
584 1096 lineSearchUpdate_(vars, uLastIter, deltaU);
585
586
2/2
✓ Branch 0 taken 931 times.
✓ Branch 1 taken 60144 times.
 61184 else if (useChop_)
587 931 choppedUpdate_(vars, uLastIter, deltaU);
588
589 else
590 {
591
1/4
✓ Branch 1 taken 54015 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 174547 auto uCurrentIter = uLastIter;
592 60253 Backend::axpy(-1.0, deltaU, uCurrentIter);
593
1/2
✓ Branch 1 taken 60099 times.
✗ Branch 2 not taken.
 60253 solutionChanged_(vars, uCurrentIter);
594
595
2/2
✓ Branch 0 taken 169 times.
✓ Branch 1 taken 59975 times.
 60253 if (enableResidualCriterion_)
596
1/2
✓ Branch 1 taken 169 times.
✗ Branch 2 not taken.
 185 computeResidualReduction_(vars);
597 }
598 62279 }
599
600 /*!
601 * \brief Indicates that one Newton iteration was finished.
602 *
603 * \param vars The variables after the current Newton iteration
604 * \param uLastIter The solution at the beginning of the current Newton iteration
605 */
606 62279 virtual void newtonEndStep(Variables &vars,
607 const SolutionVector &uLastIter)
608 {
609 if constexpr (hasPriVarsSwitch<PriVarSwitchVariables>)
610 {
611 if constexpr (assemblerExportsVariables)
612 priVarSwitchAdapter_->invoke(Backend::dofs(vars), vars);
613 else // this assumes assembly with solution (i.e. Variables=SolutionVector)
614 55696 priVarSwitchAdapter_->invoke(vars, this->assembler().gridVariables());
615 }
616
617 62256 ++numSteps_;
618
619
2/2
✓ Branch 0 taken 59177 times.
✓ Branch 1 taken 2970 times.
 62256 if (verbosity_ >= 1)
620 {
621
2/2
✓ Branch 0 taken 56495 times.
✓ Branch 1 taken 2682 times.
 59286 if (enableDynamicOutput_)
622 56604 std::cout << '\r'; // move cursor to beginning of line
623
624 59286 const auto width = Fmt::formatted_size("{}", maxSteps_);
625
2/6
✓ Branch 3 taken 59177 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 59177 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 59306 std::cout << Fmt::format("{} iteration {:{}} done", solverName_, numSteps_, width);
626
627
1/2
✓ Branch 0 taken 59177 times.
✗ Branch 1 not taken.
 59286 if (enableShiftCriterion_)
628
2/6
✓ Branch 3 taken 59177 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 59177 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 118572 std::cout << Fmt::format(", maximum relative shift = {:.4e}", shift_);
629
3/4
✓ Branch 0 taken 763 times.
✓ Branch 1 taken 58414 times.
✓ Branch 2 taken 763 times.
✗ Branch 3 not taken.
 59286 if (enableResidualCriterion_ && enableAbsoluteResidualCriterion_)
630
2/6
✓ Branch 3 taken 763 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 763 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1558 std::cout << Fmt::format(", residual = {:.4e}", residualNorm_);
631
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 58414 times.
 58507 else if (enableResidualCriterion_)
632 ✗ std::cout << Fmt::format(", residual reduction = {:.4e}", reduction_);
633
634
4/8
✓ Branch 2 taken 59177 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 59177 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4757 times.
✓ Branch 8 taken 54420 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 64043 std::cout << endIterMsgStream_.str() << "\n";
635 }
636
4/10
✓ Branch 1 taken 62147 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 62147 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 62147 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 62147 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 124512 endIterMsgStream_.str("");
637
638 // When the Newton iterations are done: ask the model to check whether it makes sense
639 // TODO: how do we realize this? -> do this here in the Newton solver
640 // model_().checkPlausibility();
641 62256 }
642
643 /*!
644 * \brief Called if the Newton method ended
645 * (not known yet if we failed or succeeded)
646 */
647 15020 virtual void newtonEnd() {}
648
649 /*!
650 * \brief Returns true if the error of the solution is below the
651 * tolerance.
652 */
653 78435 virtual bool newtonConverged() const
654 {
655 // in case the model has a priVar switch and some some primary variables
656 // actually switched their state in the last iteration, enforce another iteration
657
4/4
✓ Branch 0 taken 77153 times.
✓ Branch 1 taken 1140 times.
✓ Branch 2 taken 16306 times.
✓ Branch 3 taken 1140 times.
 95881 if (priVarSwitchAdapter_->switched())
658 return false;
659
660
3/4
✓ Branch 0 taken 77153 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 76110 times.
✓ Branch 3 taken 1043 times.
 77295 if (enableShiftCriterion_ && !enableResidualCriterion_)
661 {
662 76232 return shift_ <= shiftTolerance_;
663 }
664
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1043 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1063 else if (!enableShiftCriterion_ && enableResidualCriterion_)
665 {
666 ✗ if(enableAbsoluteResidualCriterion_)
667 ✗ return residualNorm_ <= residualTolerance_;
668 else
669 ✗ return reduction_ <= reductionTolerance_;
670 }
671
2/2
✓ Branch 0 taken 970 times.
✓ Branch 1 taken 73 times.
 1063 else if (satisfyResidualAndShiftCriterion_)
672 {
673
1/2
✓ Branch 0 taken 970 times.
✗ Branch 1 not taken.
 970 if(enableAbsoluteResidualCriterion_)
674 970 return shift_ <= shiftTolerance_
675
4/4
✓ Branch 0 taken 544 times.
✓ Branch 1 taken 426 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 540 times.
 1400 && residualNorm_ <= residualTolerance_;
676 else
677 ✗ return shift_ <= shiftTolerance_
678 ✗ && reduction_ <= reductionTolerance_;
679 }
680
2/4
✓ Branch 0 taken 73 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 73 times.
✗ Branch 3 not taken.
 93 else if(enableShiftCriterion_ && enableResidualCriterion_)
681 {
682
1/2
✓ Branch 0 taken 73 times.
✗ Branch 1 not taken.
 93 if(enableAbsoluteResidualCriterion_)
683 93 return shift_ <= shiftTolerance_
684
4/4
✓ Branch 0 taken 59 times.
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 45 times.
✓ Branch 3 taken 14 times.
 150 || residualNorm_ <= residualTolerance_;
685 else
686 ✗ return shift_ <= shiftTolerance_
687 ✗ || reduction_ <= reductionTolerance_;
688 }
689 else
690 {
691 ✗ return shift_ <= shiftTolerance_
692 ✗ || reduction_ <= reductionTolerance_
693 ✗ || residualNorm_ <= residualTolerance_;
694 }
695
696 return false;
697 }
698
699 /*!
700 * \brief Called if the Newton method broke down.
701 * This method is called _after_ newtonEnd()
702 */
703 256 virtual void newtonFail(Variables& u) {}
704
705 /*!
706 * \brief Called if the Newton method ended successfully
707 * This method is called _after_ newtonEnd()
708 */
709 14373 virtual void newtonSucceed() {}
710
711 /*!
712 * \brief output statistics / report
713 */
714 42 void report(std::ostream& sout = std::cout) const
715 {
716 sout << '\n'
717 42 << solverName_ << " statistics\n"
718 << "----------------------------------------------\n"
719 168 << "-- Total iterations: " << totalWastedIter_ + totalSucceededIter_ << '\n'
720 84 << "-- Total wasted iterations: " << totalWastedIter_ << '\n'
721 84 << "-- Total succeeded iterations: " << totalSucceededIter_ << '\n'
722 168 << "-- Average iterations per solve: " << std::setprecision(3) << double(totalSucceededIter_) / double(numConverged_) << '\n'
723 84 << "-- Number of linear solver breakdowns: " << numLinearSolverBreakdowns_ << '\n'
724 42 << std::endl;
725 42 }
726
727 /*!
728 * \brief reset the statistics
729 */
730 void resetReport()
731 {
732 totalWastedIter_ = 0;
733 totalSucceededIter_ = 0;
734 numConverged_ = 0;
735 numLinearSolverBreakdowns_ = 0;
736 }
737
738 /*!
739 * \brief Report the options and parameters this Newton is configured with
740 */
741 468 void reportParams(std::ostream& sout = std::cout) const
742 {
743 936 sout << "\n" << solverName_ << " solver configured with the following options and parameters:\n";
744 // options
745
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 441 times.
 479 if (useLineSearch_) sout << " -- " << solverName_ << ".UseLineSearch = true\n";
746
2/2
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 437 times.
 483 if (useChop_) sout << " -- " << solverName_ << ".EnableChop = true\n";
747
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 434 times.
 486 if (enablePartialReassembly_) sout << " -- " << solverName_ << ".EnablePartialReassembly = true\n";
748
2/2
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 439 times.
 485 if (enableAbsoluteResidualCriterion_) sout << " -- " << solverName_ << ".EnableAbsoluteResidualCriterion = true\n";
749
1/2
✓ Branch 0 taken 452 times.
✗ Branch 1 not taken.
 936 if (enableShiftCriterion_) sout << " -- " << solverName_ << ".EnableShiftCriterion = true (relative shift convergence criterion)\n";
750
2/2
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 439 times.
 485 if (enableResidualCriterion_) sout << " -- " << solverName_ << ".EnableResidualCriterion = true\n";
751
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 447 times.
 473 if (satisfyResidualAndShiftCriterion_) sout << " -- " << solverName_ << ".SatisfyResidualAndShiftCriterion = true\n";
752 // parameters
753
1/2
✓ Branch 0 taken 452 times.
✗ Branch 1 not taken.
 1404 if (enableShiftCriterion_) sout << " -- " << solverName_ << ".MaxRelativeShift = " << shiftTolerance_ << '\n';
754
2/2
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 439 times.
 502 if (enableAbsoluteResidualCriterion_) sout << " -- " << solverName_ << ".MaxAbsoluteResidual = " << residualTolerance_ << '\n';
755
2/2
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 439 times.
 502 if (enableResidualCriterion_) sout << " -- " << solverName_ << ".ResidualReduction = " << reductionTolerance_ << '\n';
756 1404 sout << " -- " << solverName_ << ".MinSteps = " << minSteps_ << '\n';
757 1404 sout << " -- " << solverName_ << ".MaxSteps = " << maxSteps_ << '\n';
758 1404 sout << " -- " << solverName_ << ".TargetSteps = " << targetSteps_ << '\n';
759
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 434 times.
 468 if (enablePartialReassembly_)
760 {
761 54 sout << " -- " << solverName_ << ".ReassemblyMinThreshold = " << reassemblyMinThreshold_ << '\n';
762 54 sout << " -- " << solverName_ << ".ReassemblyMaxThreshold = " << reassemblyMaxThreshold_ << '\n';
763 54 sout << " -- " << solverName_ << ".ReassemblyShiftWeight = " << reassemblyShiftWeight_ << '\n';
764 }
765 1404 sout << " -- " << solverName_ << ".RetryTimeStepReductionFactor = " << retryTimeStepReductionFactor_ << '\n';
766 1404 sout << " -- " << solverName_ << ".MaxTimeStepDivisions = " << maxTimeStepDivisions_ << '\n';
767 468 sout << std::endl;
768 468 }
769
770 /*!
771 * \brief Suggest a new time-step size based on the old time-step
772 * size.
773 *
774 * The default behavior is to suggest the old time-step size
775 * scaled by the ratio between the target iterations and the
776 * iterations required to actually solve the last time-step.
777 */
778 ✗ Scalar suggestTimeStepSize(Scalar oldTimeStep) const
779 {
780 // be aggressive reducing the time-step size but
781 // conservative when increasing it. the rationale is
782 // that we want to avoid failing in the next Newton
783 // iteration which would require another linearization
784 // of the problem.
785
2/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 66 times.
✓ Branch 3 taken 13549 times.
 13615 if (numSteps_ > targetSteps_) {
786 66 Scalar percent = Scalar(numSteps_ - targetSteps_)/targetSteps_;
787 66 return oldTimeStep/(1.0 + percent);
788 }
789
790 13549 Scalar percent = Scalar(targetSteps_ - numSteps_)/targetSteps_;
791 13549 return oldTimeStep*(1.0 + percent/1.2);
792 }
793
794 /*!
795 * \brief Specify the verbosity level
796 */
797 void setVerbosity(int val)
798 { verbosity_ = val; }
799
800 /*!
801 * \brief Return the verbosity level
802 */
803 ✗ int verbosity() const
804 ✗ { return verbosity_ ; }
805
806 /*!
807 * \brief Specify whether line search is enabled or not
808 */
809 void setUseLineSearch(bool val = true)
810 { useLineSearch_ = val; }
811
812 /*!
813 * \brief Return whether line search is enabled or not
814 */
815 bool useLineSearch() const
816 { return useLineSearch_; }
817
818 /*!
819 * \brief Returns the parameter group
820 */
821 const std::string& paramGroup() const
822 { return paramGroup_; }
823
824 /*!
825 * \brief Attach a convergence writer to write out intermediate results after each iteration
826 */
827 void attachConvergenceWriter(std::shared_ptr<ConvergenceWriter> convWriter)
828
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 { convergenceWriter_ = convWriter; }
829
830 /*!
831 * \brief Detach the convergence writer to stop the output
832 */
833 void detachConvergenceWriter()
834 { convergenceWriter_ = nullptr; }
835
836 /*!
837 * \brief Return the factor for reducing the time step after a Newton iteration has failed
838 */
839 Scalar retryTimeStepReductionFactor() const
840 { return retryTimeStepReductionFactor_; }
841
842 /*!
843 * \brief Set the factor for reducing the time step after a Newton iteration has failed
844 */
845 void setRetryTimeStepReductionFactor(const Scalar factor)
846 { retryTimeStepReductionFactor_ = factor; }
847
848 protected:
849
850 /*!
851 * \brief Update solution-dependent quantities like grid variables after the solution has changed.
852 * \todo TODO: In case we stop support for old-style grid variables / assemblers at one point,
853 * this would become obsolete as only the update call to the backend would remain.
854 */
855 60761 virtual void solutionChanged_(Variables& vars, const SolutionVector& uCurrentIter)
856 {
857 62406 Backend::update(vars, uCurrentIter);
858
859 if constexpr (!assemblerExportsVariables)
860 124522 this->assembler().updateGridVariables(Backend::dofs(vars));
861 60761 }
862
863 1344 void computeResidualReduction_(const Variables& vars)
864 {
865 // we assume that the assembler works on solution vectors
866 // if it doesn't export the variables type
867 if constexpr (!assemblerExportsVariables)
868
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 2688 this->assembler().assembleResidual(Backend::dofs(vars));
869 else
870
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 5 times.
 9 this->assembler().assembleResidual(vars);
871
872
8/8
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 5 times.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 5 times.
 5408 residualNorm_ = this->linearSolver().norm(this->assembler().residual());
873
874 1343 reduction_ = residualNorm_/initialResidual_;
875 1343 }
876
877 ✗ bool enableResidualCriterion() const
878 ✗ { return enableResidualCriterion_; }
879
880 //! optimal number of iterations we want to achieve
881 int targetSteps_;
882 //! minimum number of iterations we do
883 int minSteps_;
884 //! maximum number of iterations we do before giving up
885 int maxSteps_;
886 //! actual number of steps done so far
887 int numSteps_;
888
889 // residual criterion variables
890 Scalar reduction_;
891 Scalar residualNorm_;
892 Scalar lastReduction_;
893 Scalar initialResidual_;
894
895 // shift criterion variables
896 Scalar shift_;
897 Scalar lastShift_;
898
899 //! message stream to be displayed at the end of iterations
900 std::ostringstream endIterMsgStream_;
901
902
903 private:
904
905 /*!
906 * \brief Run the Newton method to solve a non-linear system.
907 * The solver is responsible for all the strategic decisions.
908 */
909 15097 bool solve_(Variables& vars)
910 {
911 try
912 {
913 // newtonBegin may manipulate the solution
914
1/2
✓ Branch 1 taken 15064 times.
✗ Branch 2 not taken.
 15097 newtonBegin(vars);
915
916 // the given solution is the initial guess
917
3/7
✓ Branch 1 taken 15046 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 66 times.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 30257 auto uLastIter = Backend::dofs(vars);
918
8/13
✓ Branch 1 taken 15046 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15046 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 12875 times.
✓ Branch 7 taken 15 times.
✓ Branch 8 taken 12875 times.
✓ Branch 9 taken 15 times.
✓ Branch 10 taken 66 times.
✓ Branch 11 taken 15 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 58246 ResidualVector deltaU = LinearAlgebraNativeBackend::zeros(Backend::size(Backend::dofs(vars)));
919
2/4
✓ Branch 1 taken 11842 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
 15120 Detail::Newton::assign(deltaU, Backend::dofs(vars));
920
921 // setup timers
922 15097 Dune::Timer assembleTimer(false);
923 15097 Dune::Timer solveTimer(false);
924 15097 Dune::Timer updateTimer(false);
925
926 // execute the method as long as the solver thinks
927 // that we should do another iteration
928 15097 bool converged = false;
929
3/4
✓ Branch 1 taken 77211 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 62224 times.
✓ Branch 4 taken 14987 times.
 77353 while (newtonProceed(vars, converged))
930 {
931 // notify the solver that we're about to start
932 // a new iteration
933
1/2
✓ Branch 1 taken 62224 times.
✗ Branch 2 not taken.
 62333 newtonBeginStep(vars);
934
935 // make the current solution to the old one
936
2/2
✓ Branch 0 taken 47160 times.
✓ Branch 1 taken 15064 times.
 62333 if (numSteps_ > 0)
937
2/4
✓ Branch 1 taken 47123 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 38 times.
✗ Branch 5 not taken.
 47284 uLastIter = Backend::dofs(vars);
938
939
4/4
✓ Branch 0 taken 59252 times.
✓ Branch 1 taken 2972 times.
✓ Branch 2 taken 56570 times.
✓ Branch 3 taken 2682 times.
 62333 if (verbosity_ >= 1 && enableDynamicOutput_)
940 56679 std::cout << "Assemble: r(x^k) = dS/dt + div F - q; M = grad r"
941
1/2
✓ Branch 1 taken 56570 times.
✗ Branch 2 not taken.
 56679 << std::flush;
942
943 ///////////////
944 // assemble
945 ///////////////
946
947 // linearize the problem at the current solution
948
1/2
✓ Branch 0 taken 62224 times.
✗ Branch 1 not taken.
 62333 assembleTimer.start();
949
1/2
✓ Branch 1 taken 62224 times.
✗ Branch 2 not taken.
 62333 assembleLinearSystem(vars);
950 62333 assembleTimer.stop();
951
952 ///////////////
953 // linear solve
954 ///////////////
955
956 // Clear the current line using an ansi escape
957 // sequence. for an explanation see
958 // http://en.wikipedia.org/wiki/ANSI_escape_code
959 62333 const char clearRemainingLine[] = { 0x1b, '[', 'K', 0 };
960
961
4/4
✓ Branch 0 taken 59252 times.
✓ Branch 1 taken 2972 times.
✓ Branch 2 taken 56570 times.
✓ Branch 3 taken 2682 times.
 62333 if (verbosity_ >= 1 && enableDynamicOutput_)
962 std::cout << "\rSolve: M deltax^k = r"
963
1/2
✓ Branch 3 taken 56570 times.
✗ Branch 4 not taken.
 170037 << clearRemainingLine << std::flush;
964
965 // solve the resulting linear equation system
966
1/2
✓ Branch 0 taken 62224 times.
✗ Branch 1 not taken.
 62333 solveTimer.start();
967
968 // set the delta vector to zero before solving the linear system!
969 62333 deltaU = 0;
970
971
2/2
✓ Branch 1 taken 62171 times.
✓ Branch 2 taken 53 times.
 62333 solveLinearSystem(deltaU);
972 62280 solveTimer.stop();
973
974 ///////////////
975 // update
976 ///////////////
977
4/4
✓ Branch 0 taken 59201 times.
✓ Branch 1 taken 2970 times.
✓ Branch 2 taken 56519 times.
✓ Branch 3 taken 2682 times.
 62280 if (verbosity_ >= 1 && enableDynamicOutput_)
978 std::cout << "\rUpdate: x^(k+1) = x^k - deltax^k"
979
1/2
✓ Branch 3 taken 56519 times.
✗ Branch 4 not taken.
 169884 << clearRemainingLine << std::flush;
980
981
1/2
✓ Branch 0 taken 62171 times.
✗ Branch 1 not taken.
 62280 updateTimer.start();
982 // update the current solution (i.e. uOld) with the delta
983 // (i.e. u). The result is stored in u
984
2/2
✓ Branch 1 taken 62170 times.
✓ Branch 2 taken 1 times.
 62280 newtonUpdate(vars, uLastIter, deltaU);
985 62279 updateTimer.stop();
986
987 // tell the solver that we're done with this iteration
988
2/2
✓ Branch 1 taken 62147 times.
✓ Branch 2 taken 23 times.
 62279 newtonEndStep(vars, uLastIter);
989
990 // if a convergence writer was specified compute residual and write output
991
2/2
✓ Branch 0 taken 104 times.
✓ Branch 1 taken 62043 times.
 62256 if (convergenceWriter_)
992 {
993
2/4
✓ Branch 1 taken 104 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 104 times.
✗ Branch 5 not taken.
 208 this->assembler().assembleResidual(vars);
994
4/10
✓ Branch 1 taken 104 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 104 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 104 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 104 times.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 416 convergenceWriter_->write(Backend::dofs(vars), deltaU, this->assembler().residual());
995 }
996
997 // detect if the method has converged
998
1/2
✓ Branch 1 taken 62147 times.
✗ Branch 2 not taken.
 62256 converged = newtonConverged();
999 }
1000
1001 // tell solver we are done
1002
1/2
✓ Branch 1 taken 14987 times.
✗ Branch 2 not taken.
 15020 newtonEnd();
1003
1004 // reset state if Newton failed
1005
3/4
✓ Branch 1 taken 14987 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 215 times.
✓ Branch 4 taken 14772 times.
 15020 if (!newtonConverged())
1006 {
1007 215 totalWastedIter_ += numSteps_;
1008
1/2
✓ Branch 1 taken 215 times.
✗ Branch 2 not taken.
 215 newtonFail(vars);
1009 return false;
1010 }
1011
1012 14805 totalSucceededIter_ += numSteps_;
1013 14805 numConverged_++;
1014
1015 // tell solver we converged successfully
1016
1/2
✓ Branch 1 taken 14772 times.
✗ Branch 2 not taken.
 14805 newtonSucceed();
1017
1018
2/2
✓ Branch 0 taken 14089 times.
✓ Branch 1 taken 683 times.
 14805 if (verbosity_ >= 1) {
1019 14122 const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
1020
3/8
✓ Branch 1 taken 14089 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14089 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 14089 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 28244 std::cout << Fmt::format("Assemble/solve/update time: {:.2g}({:.2f}%)/{:.2g}({:.2f}%)/{:.2g}({:.2f}%)\n",
1021 14122 assembleTimer.elapsed(), 100*assembleTimer.elapsed()/elapsedTot,
1022 14122 solveTimer.elapsed(), 100*solveTimer.elapsed()/elapsedTot,
1023 14122 updateTimer.elapsed(), 100*updateTimer.elapsed()/elapsedTot);
1024 }
1025 return true;
1026
1027 }
1028 154 catch (const NumericalProblem &e)
1029 {
1030
2/2
✓ Branch 0 taken 75 times.
✓ Branch 1 taken 2 times.
 77 if (verbosity_ >= 1)
1031
3/6
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 75 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 75 times.
✗ Branch 10 not taken.
 150 std::cout << solverName_ << ": Caught exception: \"" << e.what() << "\"\n";
1032
1033 77 totalWastedIter_ += numSteps_;
1034
1035
1/2
✓ Branch 1 taken 77 times.
✗ Branch 2 not taken.
 77 newtonFail(vars);
1036 return false;
1037 }
1038 }
1039
1040 //! assembleLinearSystem_ for assemblers that support partial reassembly
1041 template<class A>
1042 ✗ auto assembleLinearSystem_(const A& assembler, const Variables& vars)
1043 -> typename std::enable_if_t<decltype(isValid(Detail::Newton::supportsPartialReassembly())(assembler))::value, void>
1044 {
1045 166293 this->assembler().assembleJacobianAndResidual(vars, partialReassembler_.get());
1046 ✗ }
1047
1048 //! assembleLinearSystem_ for assemblers that don't support partial reassembly
1049 template<class A>
1050 ✗ auto assembleLinearSystem_(const A& assembler, const Variables& vars)
1051 -> typename std::enable_if_t<!decltype(isValid(Detail::Newton::supportsPartialReassembly())(assembler))::value, void>
1052 {
1053
1/2
✗ Branch 4 not taken.
✓ Branch 5 taken 53 times.
 13621 this->assembler().assembleJacobianAndResidual(vars);
1054 ✗ }
1055
1056 /*!
1057 * \brief Update the maximum relative shift of the solution compared to
1058 * the previous iteration. Overload for "normal" solution vectors.
1059 *
1060 * \param uLastIter The current iterative solution
1061 * \param deltaU The difference between the current and the next solution
1062 */
1063 62280 virtual void newtonUpdateShift_(const SolutionVector &uLastIter,
1064 const ResidualVector &deltaU)
1065 {
1066
1/4
✓ Branch 1 taken 55441 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 124509 auto uNew = uLastIter;
1067
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 20 times.
 62280 Backend::axpy(-1.0, deltaU, uNew);
1068
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 20 times.
 62280 shift_ = Detail::Newton::maxRelativeShift<Scalar>(uLastIter, uNew);
1069
1070
4/4
✓ Branch 0 taken 5928 times.
✓ Branch 1 taken 56243 times.
✓ Branch 2 taken 5928 times.
✓ Branch 3 taken 56243 times.
 124560 if (comm_.size() > 1)
1071
1/2
✓ Branch 1 taken 5912 times.
✗ Branch 2 not taken.
 5928 shift_ = comm_.max(shift_);
1072 62280 }
1073
1074 1096 virtual void lineSearchUpdate_(Variables &vars,
1075 const SolutionVector &uLastIter,
1076 const ResidualVector &deltaU)
1077 {
1078 1096 Scalar lambda = 1.0;
1079
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 1096 auto uCurrentIter = uLastIter;
1080
1081 3 while (true)
1082 {
1083 1168 Backend::axpy(-lambda, deltaU, uCurrentIter);
1084
1/2
✓ Branch 1 taken 1159 times.
✗ Branch 2 not taken.
 1168 solutionChanged_(vars, uCurrentIter);
1085
1086
3/3
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1163 times.
✓ Branch 2 taken 1 times.
 1168 computeResidualReduction_(vars);
1087
1088
4/4
✓ Branch 0 taken 87 times.
✓ Branch 1 taken 1080 times.
✓ Branch 2 taken 15 times.
✓ Branch 3 taken 72 times.
 1167 if (reduction_ < lastReduction_ || lambda <= lineSearchMinRelaxationFactor_)
1089 {
1090
4/11
✓ Branch 2 taken 1089 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1095 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1089 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 1096 endIterMsgStream_ << Fmt::format(", residual reduction {:.4e}->{:.4e}@lambda={:.4f}", lastReduction_, reduction_, lambda);
1091
1/2
✓ Branch 0 taken 495 times.
✗ Branch 1 not taken.
 1590 return;
1092 }
1093
1094 // try with a smaller update and reset solution
1095 72 lambda *= 0.5;
1096
1/2
✓ Branch 1 taken 69 times.
✗ Branch 2 not taken.
 72 uCurrentIter = uLastIter;
1097 }
1098 }
1099
1100 //! \note method must update the gridVariables, too!
1101 ✗ virtual void choppedUpdate_(Variables& vars,
1102 const SolutionVector& uLastIter,
1103 const ResidualVector& deltaU)
1104 {
1105 ✗ DUNE_THROW(Dune::NotImplemented,
1106 "Chopped " << solverName_ << " solver update strategy not implemented.");
1107 }
1108
1109 /*!
1110 * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
1111 *
1112 * Throws Dumux::NumericalProblem if the linear solver didn't converge.
1113 */
1114 62333 virtual bool solveLinearSystem_(ResidualVector& deltaU)
1115 {
1116
5/8
✗ Branch 0 not taken.
✓ Branch 1 taken 5964 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5964 times.
✓ Branch 4 taken 721 times.
✓ Branch 5 taken 5964 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 5964 times.
 137426 assert(this->checkSizesOfSubMatrices(this->assembler().jacobian()) && "Matrix blocks have wrong sizes!");
1117
1118 62368 return this->linearSolver().solve(
1119 124666 this->assembler().jacobian(),
1120 deltaU,
1121 62368 this->assembler().residual()
1122 249367 );
1123 }
1124
1125 //! initialize the parameters by reading from the parameter tree
1126 495 void initParams_(const std::string& group = "")
1127 {
1128
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 useLineSearch_ = getParamFromGroup<bool>(group, solverName_ + ".UseLineSearch", false);
1129
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 lineSearchMinRelaxationFactor_ = getParamFromGroup<Scalar>(group, solverName_ + ".LineSearchMinRelaxationFactor", 0.125);
1130
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 useChop_ = getParamFromGroup<bool>(group, solverName_ + ".EnableChop", false);
1131
3/4
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 468 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 11 times.
 495 if(useLineSearch_ && useChop_)
1132 ✗ DUNE_THROW(Dune::InvalidStateException, "Use either linesearch OR chop!");
1133
1134
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 enableAbsoluteResidualCriterion_ = getParamFromGroup<bool>(group, solverName_ + ".EnableAbsoluteResidualCriterion", false);
1135
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 enableShiftCriterion_ = getParamFromGroup<bool>(group, solverName_ + ".EnableShiftCriterion", true);
1136
5/10
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 466 times.
✓ Branch 7 taken 13 times.
✓ Branch 8 taken 479 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 1468 enableResidualCriterion_ = getParamFromGroup<bool>(group, solverName_ + ".EnableResidualCriterion", false) || enableAbsoluteResidualCriterion_;
1137
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 satisfyResidualAndShiftCriterion_ = getParamFromGroup<bool>(group, solverName_ + ".SatisfyResidualAndShiftCriterion", false);
1138
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 enableDynamicOutput_ = getParamFromGroup<bool>(group, solverName_ + ".EnableDynamicOutput", true);
1139
1140
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 479 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 495 if (!enableShiftCriterion_ && !enableResidualCriterion_)
1141 {
1142 ✗ DUNE_THROW(Dune::NotImplemented,
1143 "at least one of " << solverName_ << ".EnableShiftCriterion or "
1144 << solverName_ << ".EnableResidualCriterion has to be set to true");
1145 }
1146
1147
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 setMaxRelativeShift(getParamFromGroup<Scalar>(group, solverName_ + ".MaxRelativeShift", 1e-8));
1148
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 setMaxAbsoluteResidual(getParamFromGroup<Scalar>(group, solverName_ + ".MaxAbsoluteResidual", 1e-5));
1149
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 setResidualReduction(getParamFromGroup<Scalar>(group, solverName_ + ".ResidualReduction", 1e-5));
1150
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 setTargetSteps(getParamFromGroup<int>(group, solverName_ + ".TargetSteps", 10));
1151
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 479 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 setMinSteps(getParamFromGroup<int>(group, solverName_ + ".MinSteps", 2));
1152
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 479 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 setMaxSteps(getParamFromGroup<int>(group, solverName_ + ".MaxSteps", 18));
1153
1154
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 enablePartialReassembly_ = getParamFromGroup<bool>(group, solverName_ + ".EnablePartialReassembly", false);
1155
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 reassemblyMinThreshold_ = getParamFromGroup<Scalar>(group, solverName_ + ".ReassemblyMinThreshold", 1e-1*shiftTolerance_);
1156
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 reassemblyMaxThreshold_ = getParamFromGroup<Scalar>(group, solverName_ + ".ReassemblyMaxThreshold", 1e2*shiftTolerance_);
1157
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 reassemblyShiftWeight_ = getParamFromGroup<Scalar>(group, solverName_ + ".ReassemblyShiftWeight", 1e-3);
1158
1159
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 495 maxTimeStepDivisions_ = getParamFromGroup<std::size_t>(group, solverName_ + ".MaxTimeStepDivisions", 10);
1160
2/6
✓ Branch 2 taken 479 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 479 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 990 retryTimeStepReductionFactor_ = getParamFromGroup<Scalar>(group, solverName_ + ".RetryTimeStepReductionFactor", 0.5);
1161
1162 495 numSteps_ = 0;
1163
1164 // output a parameter report
1165
2/2
✓ Branch 0 taken 452 times.
✓ Branch 1 taken 27 times.
 495 if (verbosity_ >= 2)
1166 468 reportParams();
1167 495 }
1168
1169 template<class SolA, class SolB>
1170 5089 void updateDistanceFromLastLinearization_(const SolA& u, const SolB& uDelta)
1171 {
1172 if constexpr (Dune::IsNumber<SolA>::value)
1173 {
1174 ✗ auto nextPriVars = u;
1175 ✗ nextPriVars -= uDelta;
1176
1177 // add the current relative shift for this degree of freedom
1178 ✗ auto shift = Detail::Newton::maxRelativeShift<Scalar>(u, nextPriVars);
1179 ✗ distanceFromLastLinearization_[0] += shift;
1180 }
1181 else
1182 {
1183
2/2
✓ Branch 0 taken 5877638 times.
✓ Branch 1 taken 5089 times.
 5882727 for (std::size_t i = 0; i < u.size(); ++i)
1184 {
1185 5877638 const auto& currentPriVars(u[i]);
1186 5877638 auto nextPriVars(currentPriVars);
1187 11755276 nextPriVars -= uDelta[i];
1188
1189 // add the current relative shift for this degree of freedom
1190 5877638 auto shift = Detail::Newton::maxRelativeShift<Scalar>(currentPriVars, nextPriVars);
1191 11755276 distanceFromLastLinearization_[i] += shift;
1192 }
1193 }
1194 5089 }
1195
1196 template<class ...ArgsA, class...ArgsB>
1197 ✗ void updateDistanceFromLastLinearization_(const Dune::MultiTypeBlockVector<ArgsA...>& uLastIter,
1198 const Dune::MultiTypeBlockVector<ArgsB...>& deltaU)
1199 {
1200 ✗ DUNE_THROW(Dune::NotImplemented, "Reassembly for MultiTypeBlockVector");
1201 }
1202
1203 template<class Sol>
1204 ✗ void resizeDistanceFromLastLinearization_(const Sol& u, std::vector<Scalar>& dist)
1205 {
1206
0/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2582 dist.assign(Backend::size(u), 0.0);
1207 ✗ }
1208
1209 template<class ...Args>
1210 ✗ void resizeDistanceFromLastLinearization_(const Dune::MultiTypeBlockVector<Args...>& u,
1211 std::vector<Scalar>& dist)
1212 {
1213 ✗ DUNE_THROW(Dune::NotImplemented, "Reassembly for MultiTypeBlockVector");
1214 }
1215
1216 //! The communication object
1217 Communication comm_;
1218
1219 //! the verbosity level
1220 int verbosity_;
1221
1222 Scalar shiftTolerance_;
1223 Scalar reductionTolerance_;
1224 Scalar residualTolerance_;
1225
1226 // time step control
1227 std::size_t maxTimeStepDivisions_;
1228 Scalar retryTimeStepReductionFactor_;
1229
1230 // further parameters
1231 bool useLineSearch_;
1232 Scalar lineSearchMinRelaxationFactor_;
1233 bool useChop_;
1234 bool enableAbsoluteResidualCriterion_;
1235 bool enableShiftCriterion_;
1236 bool enableResidualCriterion_;
1237 bool satisfyResidualAndShiftCriterion_;
1238 bool enableDynamicOutput_;
1239
1240 //! the parameter group problem prefix for getting parameters from the parameter tree
1241 std::string paramGroup_;
1242 //! the parameter group for getting parameters from the parameter tree
1243 std::string solverName_;
1244
1245 // infrastructure for partial reassembly
1246 bool enablePartialReassembly_;
1247 std::unique_ptr<Reassembler> partialReassembler_;
1248 std::vector<Scalar> distanceFromLastLinearization_;
1249 Scalar reassemblyMinThreshold_;
1250 Scalar reassemblyMaxThreshold_;
1251 Scalar reassemblyShiftWeight_;
1252
1253 // statistics for the optional report
1254 std::size_t totalWastedIter_ = 0; //! Newton steps in solves that didn't converge
1255 std::size_t totalSucceededIter_ = 0; //! Newton steps in solves that converged
1256 std::size_t numConverged_ = 0; //! total number of converged solves
1257 std::size_t numLinearSolverBreakdowns_ = 0; //! total number of linear solves that failed
1258
1259 //! the class handling the primary variable switch
1260 std::unique_ptr<PrimaryVariableSwitchAdapter> priVarSwitchAdapter_;
1261
1262 //! convergence writer
1263 std::shared_ptr<ConvergenceWriter> convergenceWriter_ = nullptr;
1264 };
1265
1266 } // end namespace Dumux
1267
1268 #endif
1269