GCC Code Coverage Report

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