GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/nonlinear/newtonsolver.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 352 382 92.1%
Functions: 5161 7756 66.5%
Branches: 379 872 43.5%
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 64464900 auto hybridInnerProduct(const V& v1, const V& v2, Scalar init, Reduce&& r, Transform&& t)
91 -> std::enable_if_t<hasDynamicIndexAccess<V>(), Scalar>
92 {
93 64464900 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 117025252 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 234050504 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 70575141 auto maxRelativeShift(const V& v1, const V& v2)
122 -> std::enable_if_t<!Dune::IsNumber<V>::value, Scalar>
123 {
124 70633685 return hybridInnerProduct(v1, v2, Scalar(0.0),
125
11/12
✓ Branch 0 taken 68582906 times.
✓ Branch 1 taken 43459983 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.
 248081963 [](const auto& a, const auto& b){ using std::max; return max(a, b); },
126
12/12
✓ Branch 0 taken 68582919 times.
✓ Branch 1 taken 43459970 times.
✓ Branch 2 taken 7117664 times.
✓ Branch 3 taken 50075110 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.
 175509119 [](const auto& a, const auto& b){ return maxRelativeShift<Scalar>(a, b); }
127 70570130 );
128 }
129
130 template<class To, class From>
131 1621809 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 1611541 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 9884 }
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 14465 void solve(Variables& vars, TimeLoop& timeLoop) override
310 {
311 if constexpr (!assemblerExportsVariables)
312 {
313
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14436 times.
 14450 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 14705 times.
✗ Branch 1 not taken.
 14732 for (std::size_t i = 0; i <= maxTimeStepDivisions_; ++i)
319 {
320 // linearize & solve
321 14732 const bool converged = solve_(vars);
322
323
2/2
✓ Branch 0 taken 282 times.
✓ Branch 1 taken 14438 times.
 14732 if (converged)
324 return;
325
326
1/2
✓ Branch 0 taken 267 times.
✗ Branch 1 not taken.
 267 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.
 267 Backend::update(vars, this->assembler().prevSol());
334
1/2
✓ Branch 0 taken 248 times.
✗ Branch 1 not taken.
 267 this->assembler().resetTimeStep(Backend::dofs(vars));
335
336
2/2
✓ Branch 0 taken 265 times.
✓ Branch 1 taken 2 times.
 267 if (verbosity_ >= 1)
337 {
338 265 const auto dt = timeLoop.timeStepSize();
339 265 std::cout << Fmt::format("{} solver did not converge with dt = {} seconds. ", solverName_, dt)
340
1/4
✓ Branch 2 taken 265 times.
✗ Branch 3 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 795 << Fmt::format("Retrying with time step of dt = {} seconds.\n", dt*retryTimeStepReductionFactor_);
341 }
342
343 // try again with dt = dt * retryTimeStepReductionFactor_
344 267 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.
 267 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 15606 virtual void newtonBegin(Variables& initVars)
391 {
392
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 25 times.
 15606 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 15606 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 15529 times.
 15606 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 14267 times.
 15606 if (enablePartialReassembly_)
416 {
417 1293 partialReassembler_->resetColors();
418 2586 resizeDistanceFromLastLinearization_(initSol, distanceFromLastLinearization_);
419 }
420 15606 }
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 79309 virtual bool newtonProceed(const Variables &varsCurrentIter, bool converged)
429 {
430
2/2
✓ Branch 0 taken 47768 times.
✓ Branch 1 taken 31349 times.
 79309 if (numSteps_ < minSteps_)
431 return true;
432
10/10
✓ Branch 0 taken 32487 times.
✓ Branch 1 taken 15285 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.
 47874 else if (converged)
433 return false; // we are below the desired tolerance
434
6/10
✓ Branch 0 taken 233 times.
✓ Branch 1 taken 32254 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.
 32544 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 230 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.
 230 if (enableShiftCriterion_)
440 230 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 63708 virtual void newtonBeginStep(const Variables& vars)
452 {
453 63708 lastShift_ = shift_;
454
4/4
✓ Branch 0 taken 15540 times.
✓ Branch 1 taken 48034 times.
✓ Branch 2 taken 15 times.
✓ Branch 3 taken 27 times.
 63708 if (numSteps_ == 0)
455 {
456 15580 lastReduction_ = 1.0;
457 }
458 else
459 {
460 48128 lastReduction_ = reduction_;
461 }
462 59909 }
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 63762 virtual void assembleLinearSystem(const Variables& vars)
470 {
471 63783 assembleLinearSystem_(this->assembler(), vars);
472
473
2/2
✓ Branch 0 taken 5103 times.
✓ Branch 1 taken 54627 times.
 59854 if (enablePartialReassembly_)
474 5103 partialReassembler_->report(comm_, endIterMsgStream_);
475 63762 }
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 63783 void solveLinearSystem(ResidualVector& deltaU)
489 {
490 63783 bool converged = false;
491
492 try
493 {
494
2/2
✓ Branch 0 taken 15560 times.
✓ Branch 1 taken 48077 times.
 63783 if (numSteps_ == 0)
495
1/2
✓ Branch 1 taken 15537 times.
✗ Branch 2 not taken.
 15606 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 63612 times.
✓ Branch 2 taken 4 times.
 63825 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 from the linear solver: \"" << e.what() << "\"\n";
505
506 4 converged = false;
507 }
508
509 // make sure all processes converged
510 63783 int convergedRemote = converged;
511
2/2
✓ Branch 0 taken 5976 times.
✓ Branch 1 taken 57640 times.
 63741 if (comm_.size() > 1)
512 5976 convergedRemote = comm_.min(converged);
513
514
2/2
✓ Branch 0 taken 29 times.
✓ Branch 1 taken 63608 times.
 63783 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 63586 times.
 63712 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 63753 }
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 63753 void newtonUpdate(Variables& vars,
544 const SolutionVector& uLastIter,
545 const ResidualVector& deltaU)
546 {
547
2/2
✓ Branch 0 taken 1772 times.
✓ Branch 1 taken 61835 times.
 63753 if (useLineSearch_)
548 1793 lineSearchUpdate_(vars, uLastIter, deltaU);
549
550
2/2
✓ Branch 0 taken 975 times.
✓ Branch 1 taken 60860 times.
 61960 else if (useChop_)
551 975 choppedUpdate_(vars, uLastIter, deltaU);
552
553 else
554 {
555 60985 auto uCurrentIter = uLastIter;
556 60985 Backend::axpy(-1.0, deltaU, uCurrentIter);
557
1/2
✓ Branch 1 taken 60815 times.
✗ Branch 2 not taken.
 60985 solutionChanged_(vars, uCurrentIter);
558
559
2/2
✓ Branch 0 taken 171 times.
✓ Branch 1 taken 60689 times.
 60985 if (enableResidualCriterion_)
560
1/2
✓ Branch 1 taken 171 times.
✗ Branch 2 not taken.
 189 computeResidualReduction_(vars);
561 60940 }
562
563
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 63606 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 63752 if (enableShiftCriterion_ || enablePartialReassembly_)
564 63794 newtonComputeShift_(Backend::dofs(vars), uLastIter);
565
566
2/2
✓ Branch 0 taken 5103 times.
✓ Branch 1 taken 58503 times.
 63752 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 63752 }
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 63752 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 63729 ++numSteps_;
622
623
2/2
✓ Branch 0 taken 60591 times.
✓ Branch 1 taken 2992 times.
 63729 if (verbosity_ >= 1)
624 {
625
2/2
✓ Branch 0 taken 57625 times.
✓ Branch 1 taken 2966 times.
 60737 if (enableDynamicOutput_)
626 57771 std::cout << '\r'; // move cursor to beginning of line
627
628 60737 const auto width = Fmt::formatted_size("{}", maxSteps_);
629 121474 std::cout << Fmt::format("{} iteration {:{}} done", solverName_, numSteps_, width);
630
631
1/2
✓ Branch 0 taken 60591 times.
✗ Branch 1 not taken.
 60737 if (enableShiftCriterion_)
632 121474 std::cout << Fmt::format(", maximum relative shift = {:.4e}", shift_);
633
3/4
✓ Branch 0 taken 59272 times.
✓ Branch 1 taken 1319 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 59272 times.
 60737 if (enableResidualCriterion_ || enableAbsoluteResidualCriterion_)
634 2674 std::cout << Fmt::format(", residual = {:.4e}, residual reduction = {:.4e}", residualNorm_, reduction_);
635
636
2/4
✓ Branch 2 taken 60591 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 60591 times.
✗ Branch 6 not taken.
 121474 std::cout << endIterMsgStream_.str() << "\n";
637 }
638 127458 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 63729 }
644
645 /*!
646 * \brief Called if the Newton method ended
647 * (not known yet if we failed or succeeded)
648 */
649 15542 virtual void newtonEnd() {}
650
651 /*!
652 * \brief Returns true if the error of the solution is below the
653 * tolerance.
654 */
655 80442 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 79102 times.
✓ Branch 1 taken 1148 times.
 80442 if (priVarSwitchAdapter_->switched())
660 return false;
661
662
3/4
✓ Branch 0 taken 79102 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 77280 times.
✓ Branch 3 taken 1822 times.
 79294 if (enableShiftCriterion_ && !enableResidualCriterion_)
663 {
664 77448 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 */
705 235 virtual void newtonFail(Variables& u) {}
706
707 /*!
708 * \brief Called if the Newton method ended successfully
709 * This method is called _after_ newtonEnd()
710 */
711 14893 virtual void newtonSucceed() {}
712
713 /*!
714 * \brief output statistics / report
715 */
716 45 void report(std::ostream& sout = std::cout) const
717 {
718 sout << '\n'
719 45 << solverName_ << " statistics\n"
720 << "----------------------------------------------\n"
721 90 << "-- Total iterations: " << totalWastedIter_ + totalSucceededIter_ << '\n'
722 45 << "-- Total wasted iterations: " << totalWastedIter_ << '\n'
723 45 << "-- Total succeeded iterations: " << totalSucceededIter_ << '\n'
724 45 << "-- Average iterations per solve: " << std::setprecision(3) << double(totalSucceededIter_) / double(numConverged_) << '\n'
725 45 << "-- Number of linear solver breakdowns: " << numLinearSolverBreakdowns_ << '\n'
726 45 << std::endl;
727 45 }
728
729 /*!
730 * \brief reset the statistics
731 */
732 void resetReport()
733 {
734 totalWastedIter_ = 0;
735 totalSucceededIter_ = 0;
736 numConverged_ = 0;
737 numLinearSolverBreakdowns_ = 0;
738 }
739
740 /*!
741 * \brief Report the options and parameters this Newton is configured with
742 */
743 521 void reportParams(std::ostream& sout = std::cout) const
744 {
745 521 sout << "\n" << solverName_ << " solver configured with the following options and parameters:\n";
746 // options
747
2/2
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 479 times.
 534 if (useLineSearch_) sout << " -- " << solverName_ << ".UseLineSearch = true\n";
748
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 476 times.
 537 if (useChop_) sout << " -- " << solverName_ << ".EnableChop = true\n";
749
2/2
✓ Branch 0 taken 19 times.
✓ Branch 1 taken 473 times.
 540 if (enablePartialReassembly_) sout << " -- " << solverName_ << ".EnablePartialReassembly = true\n";
750
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 476 times.
 543 if (enableAbsoluteResidualCriterion_) sout << " -- " << solverName_ << ".EnableAbsoluteResidualCriterion = true\n";
751
1/2
✓ Branch 0 taken 492 times.
✗ Branch 1 not taken.
 1042 if (enableShiftCriterion_) sout << " -- " << solverName_ << ".EnableShiftCriterion = true (relative shift convergence criterion)\n";
752
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 476 times.
 543 if (enableResidualCriterion_) sout << " -- " << solverName_ << ".EnableResidualCriterion = true\n";
753
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 486 times.
 527 if (satisfyResidualAndShiftCriterion_) sout << " -- " << solverName_ << ".SatisfyResidualAndShiftCriterion = true\n";
754 // parameters
755
1/2
✓ Branch 0 taken 492 times.
✗ Branch 1 not taken.
 1042 if (enableShiftCriterion_) sout << " -- " << solverName_ << ".MaxRelativeShift = " << shiftTolerance_ << '\n';
756
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 476 times.
 543 if (enableAbsoluteResidualCriterion_) sout << " -- " << solverName_ << ".MaxAbsoluteResidual = " << residualTolerance_ << '\n';
757
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 476 times.
 543 if (enableResidualCriterion_) sout << " -- " << solverName_ << ".ResidualReduction = " << reductionTolerance_ << '\n';
758 521 sout << " -- " << solverName_ << ".MinSteps = " << minSteps_ << '\n';
759 521 sout << " -- " << solverName_ << ".MaxSteps = " << maxSteps_ << '\n';
760 521 sout << " -- " << solverName_ << ".TargetSteps = " << targetSteps_ << '\n';
761
2/2
✓ Branch 0 taken 19 times.
✓ Branch 1 taken 473 times.
 521 if (enablePartialReassembly_)
762 {
763 19 sout << " -- " << solverName_ << ".ReassemblyMinThreshold = " << reassemblyMinThreshold_ << '\n';
764 19 sout << " -- " << solverName_ << ".ReassemblyMaxThreshold = " << reassemblyMaxThreshold_ << '\n';
765 19 sout << " -- " << solverName_ << ".ReassemblyShiftWeight = " << reassemblyShiftWeight_ << '\n';
766 }
767 521 sout << " -- " << solverName_ << ".RetryTimeStepReductionFactor = " << retryTimeStepReductionFactor_ << '\n';
768 521 sout << " -- " << solverName_ << ".MaxTimeStepDivisions = " << maxTimeStepDivisions_ << '\n';
769 521 sout << std::endl;
770 521 }
771
772 /*!
773 * \brief Suggest a new time-step size based on the old time-step
774 * size.
775 *
776 * The default behavior is to suggest the old time-step size
777 * scaled by the ratio between the target iterations and the
778 * iterations required to actually solve the last time-step.
779 */
780 14016 Scalar suggestTimeStepSize(Scalar oldTimeStep) const
781 {
782 // be aggressive reducing the time-step size but
783 // conservative when increasing it. the rationale is
784 // that we want to avoid failing in the next Newton
785 // iteration which would require another linearization
786 // of the problem.
787
3/4
✓ Branch 0 taken 66 times.
✓ Branch 1 taken 13940 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
 14016 if (numSteps_ > targetSteps_) {
788 66 Scalar percent = Scalar(numSteps_ - targetSteps_)/targetSteps_;
789 66 return oldTimeStep/(1.0 + percent);
790 }
791
792 13950 Scalar percent = Scalar(targetSteps_ - numSteps_)/targetSteps_;
793 13950 return oldTimeStep*(1.0 + percent/1.2);
794 }
795
796 /*!
797 * \brief Specify the verbosity level
798 */
799 void setVerbosity(int val)
800 { verbosity_ = val; }
801
802 /*!
803 * \brief Return the verbosity level
804 */
805 56 int verbosity() const
806
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_ ; }
807
808 /*!
809 * \brief Specify whether line search is enabled or not
810 */
811 5 void setUseLineSearch(bool val = true)
812 5 { useLineSearch_ = val; }
813
814 /*!
815 * \brief Return whether line search is enabled or not
816 */
817 bool useLineSearch() const
818 { return useLineSearch_; }
819
820 /*!
821 * \brief Returns the parameter group
822 */
823
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
824
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_; }
825
826 /*!
827 * \brief Attach a convergence writer to write out intermediate results after each iteration
828 */
829 4 void attachConvergenceWriter(std::shared_ptr<ConvergenceWriter> convWriter)
830
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 { convergenceWriter_ = convWriter; }
831
832 /*!
833 * \brief Detach the convergence writer to stop the output
834 */
835 void detachConvergenceWriter()
836 { convergenceWriter_ = nullptr; }
837
838 /*!
839 * \brief Return the factor for reducing the time step after a Newton iteration has failed
840 */
841 Scalar retryTimeStepReductionFactor() const
842 { return retryTimeStepReductionFactor_; }
843
844 /*!
845 * \brief Set the factor for reducing the time step after a Newton iteration has failed
846 */
847 void setRetryTimeStepReductionFactor(const Scalar factor)
848 { retryTimeStepReductionFactor_ = factor; }
849
850 protected:
851
852 /*!
853 * \brief Update solution-dependent quantities like grid variables after the solution has changed.
854 * \todo TODO: In case we stop support for old-style grid variables / assemblers at one point,
855 * this would become obsolete as only the update call to the backend would remain.
856 */
857
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.
 63850 virtual void solutionChanged_(Variables& vars, const SolutionVector& uCurrentIter)
858 {
859
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.
 63850 Backend::update(vars, uCurrentIter);
860
861 if constexpr (!assemblerExportsVariables)
862 63738 this->assembler().updateGridVariables(Backend::dofs(vars));
863 62859 }
864
865 2058 void computeResidualReduction_(const Variables& vars)
866 {
867 // we assume that the assembler works on solution vectors
868 // if it doesn't export the variables type
869 if constexpr (!assemblerExportsVariables)
870
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 2049 this->assembler().assembleResidual(Backend::dofs(vars));
871 else
872
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 5 times.
 9 this->assembler().assembleResidual(vars);
873
874
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 5 times.
 2057 residualNorm_ = this->linearSolver().norm(this->assembler().residual());
875
876 2057 reduction_ = residualNorm_/initialResidual_;
877 2048 }
878
879 975 bool enableResidualCriterion() const
880
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 975 times.
 975 { return enableResidualCriterion_; }
881
882 //! optimal number of iterations we want to achieve
883 int targetSteps_;
884 //! minimum number of iterations we do
885 int minSteps_;
886 //! maximum number of iterations we do before giving up
887 int maxSteps_;
888 //! actual number of steps done so far
889 int numSteps_;
890
891 // residual criterion variables
892 Scalar reduction_;
893 Scalar residualNorm_;
894 Scalar lastReduction_;
895 Scalar initialResidual_;
896
897 // shift criterion variables
898 Scalar shift_;
899 Scalar lastShift_;
900
901 //! message stream to be displayed at the end of iterations
902 std::ostringstream endIterMsgStream_;
903
904
905 private:
906
907 /*!
908 * \brief Run the Newton method to solve a non-linear system.
909 * The solver is responsible for all the strategic decisions.
910 */
911 15606 bool solve_(Variables& vars)
912 {
913 try
914 {
915 // newtonBegin may manipulate the solution
916
1/2
✓ Branch 1 taken 15560 times.
✗ Branch 2 not taken.
 15606 newtonBegin(vars);
917
918 // the given solution is the initial guess
919
1/2
✓ Branch 1 taken 15542 times.
✗ Branch 2 not taken.
 15606 auto uLastIter = Backend::dofs(vars);
920
1/2
✓ Branch 1 taken 15542 times.
✗ Branch 2 not taken.
 15606 ResidualVector deltaU = LinearAlgebraNativeBackend::zeros(Backend::size(Backend::dofs(vars)));
921
1/2
✓ Branch 1 taken 12329 times.
✗ Branch 2 not taken.
 25515 Detail::Newton::assign(deltaU, Backend::dofs(vars));
922
923 // setup timers
924 15606 Dune::Timer assembleTimer(false);
925 15606 Dune::Timer solveTimer(false);
926 15606 Dune::Timer updateTimer(false);
927
928 // execute the method as long as the solver thinks
929 // that we should do another iteration
930 15606 bool converged = false;
931
4/5
✓ Branch 1 taken 79133 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 63616 times.
✓ Branch 4 taken 15501 times.
✓ Branch 0 taken 10 times.
 79335 while (newtonProceed(vars, converged))
932 {
933 // notify the solver that we're about to start
934 // a new iteration
935
2/3
✓ Branch 1 taken 63632 times.
✗ Branch 2 not taken.
✓ Branch 0 taken 5 times.
 63783 newtonBeginStep(vars);
936
937 // make the current solution to the old one
938
2/2
✓ Branch 0 taken 48077 times.
✓ Branch 1 taken 15560 times.
 63783 if (numSteps_ > 0)
939
1/2
✓ Branch 1 taken 48040 times.
✗ Branch 2 not taken.
 48177 uLastIter = Backend::dofs(vars);
940
941
4/4
✓ Branch 0 taken 60643 times.
✓ Branch 1 taken 2994 times.
✓ Branch 2 taken 57677 times.
✓ Branch 3 taken 2966 times.
 63783 if (verbosity_ >= 1 && enableDynamicOutput_)
942
1/2
✓ Branch 1 taken 57677 times.
✗ Branch 2 not taken.
 57823 std::cout << "Assemble: r(x^k) = dS/dt + div F - q; M = grad r"
943
1/2
✓ Branch 0 taken 63637 times.
✗ Branch 1 not taken.
 63783 << std::flush;
944
945 ///////////////
946 // assemble
947 ///////////////
948
949 // linearize the problem at the current solution
950 63783 assembleTimer.start();
951
1/2
✓ Branch 1 taken 63637 times.
✗ Branch 2 not taken.
 63783 assembleLinearSystem(vars);
952 63783 assembleTimer.stop();
953
954 ///////////////
955 // linear solve
956 ///////////////
957
958 // Clear the current line using an ansi escape
959 // sequence. for an explanation see
960 // http://en.wikipedia.org/wiki/ANSI_escape_code
961 63783 const char clearRemainingLine[] = { 0x1b, '[', 'K', 0 };
962
963
4/4
✓ Branch 0 taken 60643 times.
✓ Branch 1 taken 2994 times.
✓ Branch 2 taken 57677 times.
✓ Branch 3 taken 2966 times.
 63783 if (verbosity_ >= 1 && enableDynamicOutput_)
964 std::cout << "\rSolve: M deltax^k = r"
965
3/6
✓ Branch 1 taken 57677 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57677 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 63637 times.
✗ Branch 7 not taken.
 63783 << clearRemainingLine << std::flush;
966
967 // solve the resulting linear equation system
968
2/2
✓ Branch 0 taken 46370500 times.
✓ Branch 1 taken 57731 times.
 47228826 solveTimer.start();
969
970 // set the delta vector to zero before solving the linear system!
971 63783 deltaU = 0;
972
973
2/2
✓ Branch 1 taken 63607 times.
✓ Branch 2 taken 30 times.
 63783 solveLinearSystem(deltaU);
974 63753 solveTimer.stop();
975
976 ///////////////
977 // update
978 ///////////////
979
4/4
✓ Branch 0 taken 60615 times.
✓ Branch 1 taken 2992 times.
✓ Branch 2 taken 57649 times.
✓ Branch 3 taken 2966 times.
 63753 if (verbosity_ >= 1 && enableDynamicOutput_)
980 std::cout << "\rUpdate: x^(k+1) = x^k - deltax^k"
981
3/6
✓ Branch 1 taken 57649 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57649 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 63607 times.
✗ Branch 7 not taken.
 63753 << clearRemainingLine << std::flush;
982
983 63753 updateTimer.start();
984 // update the current solution (i.e. uOld) with the delta
985 // (i.e. u). The result is stored in u
986
2/2
✓ Branch 1 taken 63606 times.
✓ Branch 2 taken 1 times.
 63753 newtonUpdate(vars, uLastIter, deltaU);
987 63752 updateTimer.stop();
988
989 // tell the solver that we're done with this iteration
990
2/2
✓ Branch 1 taken 63583 times.
✓ Branch 2 taken 23 times.
 63752 newtonEndStep(vars, uLastIter);
991
992 // if a convergence writer was specified compute residual and write output
993
2/2
✓ Branch 0 taken 104 times.
✓ Branch 1 taken 63479 times.
 63729 if (convergenceWriter_)
994 {
995
1/2
✓ Branch 1 taken 104 times.
✗ Branch 2 not taken.
 104 this->assembler().assembleResidual(vars);
996
1/2
✓ Branch 1 taken 104 times.
✗ Branch 2 not taken.
 104 convergenceWriter_->write(Backend::dofs(vars), deltaU, this->assembler().residual());
997 }
998
999 // detect if the method has converged
1000
1/2
✓ Branch 1 taken 63562 times.
✗ Branch 2 not taken.
 63729 converged = newtonConverged();
1001 }
1002
1003 // tell solver we are done
1004
1/2
✓ Branch 1 taken 15501 times.
✗ Branch 2 not taken.
 15542 newtonEnd();
1005
1006 // reset state if Newton failed
1007
3/5
✓ Branch 1 taken 15506 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 217 times.
✓ Branch 4 taken 15284 times.
✗ Branch 0 not taken.
 15552 if (!newtonConverged())
1008 {
1009 217 totalWastedIter_ += numSteps_;
1010
1/2
✓ Branch 1 taken 217 times.
✗ Branch 2 not taken.
 217 newtonFail(vars);
1011 ✗ return false;
1012 }
1013
1014 15335 totalSucceededIter_ += numSteps_;
1015 15335 numConverged_++;
1016
1017 // tell solver we converged successfully
1018
1/2
✓ Branch 1 taken 15284 times.
✗ Branch 2 not taken.
 15325 newtonSucceed();
1019
1020
2/2
✓ Branch 0 taken 14602 times.
✓ Branch 1 taken 687 times.
 15335 if (verbosity_ >= 1) {
1021 14648 const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
1022 14648 std::cout << Fmt::format("Assemble/solve/update time: {:.2g}({:.2f}%)/{:.2g}({:.2f}%)/{:.2g}({:.2f}%)\n",
1023
1/2
✓ Branch 2 taken 14602 times.
✗ Branch 3 not taken.
 14648 assembleTimer.elapsed(), 100*assembleTimer.elapsed()/elapsedTot,
1024 14648 solveTimer.elapsed(), 100*solveTimer.elapsed()/elapsedTot,
1025
1/2
✓ Branch 3 taken 29 times.
✗ Branch 4 not taken.
 29346 updateTimer.elapsed(), 100*updateTimer.elapsed()/elapsedTot);
1026 }
1027 return true;
1028
1029 28718 }
1030
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
 108 catch (const NumericalProblem &e)
1031 {
1032
2/2
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 2 times.
 54 if (verbosity_ >= 1)
1033
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";
1034
1035 54 totalWastedIter_ += numSteps_;
1036
1037
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 newtonFail(vars);
1038 return false;
1039 }
1040 }
1041
1042 //! assembleLinearSystem_ for assemblers that support partial reassembly
1043 template<class A>
1044 55903 auto assembleLinearSystem_(const A& assembler, const Variables& vars)
1045 -> typename std::enable_if_t<decltype(isValid(Detail::Newton::supportsPartialReassembly())(assembler))::value, void>
1046 {
1047 55903 this->assembler().assembleJacobianAndResidual(vars, partialReassembler_.get());
1048 }
1049
1050 //! assembleLinearSystem_ for assemblers that don't support partial reassembly
1051 template<class A>
1052 7734 auto assembleLinearSystem_(const A& assembler, const Variables& vars)
1053 -> typename std::enable_if_t<!decltype(isValid(Detail::Newton::supportsPartialReassembly())(assembler))::value, void>
1054 {
1055
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);
1056 }
1057
1058 /*!
1059 * \brief Update the maximum relative shift of the solution compared to
1060 * the previous iteration. Overload for "normal" solution vectors.
1061 *
1062 * \param uLastIter The current iterative solution
1063 * \param deltaU The difference between the current and the next solution
1064 */
1065 [[deprecated("Use computeShift_(u1, u2) instead")]]
1066 ✗ virtual void newtonUpdateShift_(const SolutionVector &uLastIter,
1067 const ResidualVector &deltaU)
1068 {
1069 ✗ auto uNew = uLastIter;
1070 ✗ Backend::axpy(-1.0, deltaU, uNew);
1071 ✗ newtonComputeShift_(uLastIter, uNew);
1072 ✗ }
1073
1074 /*!
1075 * \brief Update the maximum relative shift of one solution
1076 * compared to another.
1077 */
1078
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 20 times.
 63731 virtual void newtonComputeShift_(const SolutionVector &u1,
1079 const SolutionVector &u2)
1080 {
1081 63731 shift_ = Detail::Newton::maxRelativeShift<Scalar>(u1, u2);
1082
2/10
✓ Branch 0 taken 5972 times.
✓ Branch 1 taken 57613 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.
 63710 if (comm_.size() > 1)
1083 5972 shift_ = comm_.max(shift_);
1084 63731 }
1085
1086 /*!
1087 * \brief Use a line search update based on simple backtracking
1088 * \note method must update the gridVariables if the solution changes
1089 *
1090 * \param vars The variables to be updated
1091 * \param uLastIter The solution vector after the last iteration
1092 * \param deltaU The computed difference between the current and the next solution (full update)
1093 */
1094 1751 virtual void lineSearchUpdate_(Variables &vars,
1095 const SolutionVector &uLastIter,
1096 const ResidualVector &deltaU)
1097 {
1098 1751 Scalar lambda = 1.0;
1099 1842 auto uCurrentIter = uLastIter;
1100
1101 3 while (true)
1102 {
1103 3109 Backend::axpy(-lambda, deltaU, uCurrentIter);
1104
1/2
✓ Branch 1 taken 1860 times.
✗ Branch 2 not taken.
 1869 solutionChanged_(vars, uCurrentIter);
1105
1106
2/2
✓ Branch 1 taken 1859 times.
✓ Branch 2 taken 1 times.
 1869 computeResidualReduction_(vars);
1107
1108
4/4
✓ Branch 0 taken 143 times.
✓ Branch 1 taken 1725 times.
✓ Branch 2 taken 25 times.
✓ Branch 3 taken 118 times.
 1868 if (reduction_ < lastReduction_ || lambda <= lineSearchMinRelaxationFactor_)
1109 {
1110
1/4
✓ Branch 1 taken 1744 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 3501 endIterMsgStream_ << Fmt::format(", residual reduction {:.4e}->{:.4e}@lambda={:.4f}", lastReduction_, reduction_, lambda);
1111
1/2
✓ Branch 0 taken 596 times.
✗ Branch 1 not taken.
 1750 return;
1112 }
1113
1114 // try with a smaller update and reset solution
1115
1/2
✓ Branch 1 taken 91 times.
✗ Branch 2 not taken.
 118 lambda *= 0.5;
1116
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 1267 uCurrentIter = uLastIter;
1117 }
1118 1744 }
1119
1120 /*!
1121 * \brief Use a custom chopped update strategy (do not use the full update)
1122 * \note method must update the gridVariables if the solution changes
1123 *
1124 * \param vars The variables to be updated
1125 * \param uLastIter The solution vector after the last iteration
1126 * \param deltaU The computed difference between the current and the next solution (full update)
1127 */
1128 ✗ virtual void choppedUpdate_(Variables& vars,
1129 const SolutionVector& uLastIter,
1130 const ResidualVector& deltaU)
1131 {
1132 ✗ DUNE_THROW(Dune::NotImplemented,
1133 "Chopped " << solverName_ << " solver update strategy not implemented.");
1134 }
1135
1136 /*!
1137 * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
1138 *
1139 * Throws Dumux::NumericalProblem if the linear solver didn't converge.
1140 */
1141 63762 virtual bool solveLinearSystem_(ResidualVector& deltaU)
1142 {
1143
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.
 70504 assert(this->checkSizesOfSubMatrices(this->assembler().jacobian()) && "Matrix blocks have wrong sizes!");
1144
1145
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.
 63762 return this->linearSolver().solve(
1146
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.
 63762 this->assembler().jacobian(),
1147 deltaU,
1148
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.
 63762 this->assembler().residual()
1149 63737 );
1150 }
1151
1152 //! initialize the parameters by reading from the parameter tree
1153 549 void initParams_(const std::string& group = "")
1154 {
1155 1098 useLineSearch_ = getParamFromGroup<bool>(group, solverName_ + ".UseLineSearch", false);
1156 1098 lineSearchMinRelaxationFactor_ = getParamFromGroup<Scalar>(group, solverName_ + ".LineSearchMinRelaxationFactor", 0.125);
1157 1098 useChop_ = getParamFromGroup<bool>(group, solverName_ + ".EnableChop", false);
1158
3/4
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 507 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13 times.
 549 if(useLineSearch_ && useChop_)
1159 ✗ DUNE_THROW(Dune::InvalidStateException, "Use either linesearch OR chop!");
1160
1161 1098 enableAbsoluteResidualCriterion_ = getParamFromGroup<bool>(group, solverName_ + ".EnableAbsoluteResidualCriterion", false);
1162 1098 enableShiftCriterion_ = getParamFromGroup<bool>(group, solverName_ + ".EnableShiftCriterion", true);
1163
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_;
1164 1098 satisfyResidualAndShiftCriterion_ = getParamFromGroup<bool>(group, solverName_ + ".SatisfyResidualAndShiftCriterion", false);
1165 1098 enableDynamicOutput_ = getParamFromGroup<bool>(group, solverName_ + ".EnableDynamicOutput", true);
1166
1167
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 520 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 549 if (!enableShiftCriterion_ && !enableResidualCriterion_)
1168 {
1169 ✗ DUNE_THROW(Dune::NotImplemented,
1170 "at least one of " << solverName_ << ".EnableShiftCriterion or "
1171 << solverName_ << ".EnableResidualCriterion has to be set to true");
1172 }
1173
1174 1098 setMaxRelativeShift(getParamFromGroup<Scalar>(group, solverName_ + ".MaxRelativeShift", 1e-8));
1175 1098 setMaxAbsoluteResidual(getParamFromGroup<Scalar>(group, solverName_ + ".MaxAbsoluteResidual", 1e-5));
1176 1098 setResidualReduction(getParamFromGroup<Scalar>(group, solverName_ + ".ResidualReduction", 1e-5));
1177 1098 setTargetSteps(getParamFromGroup<int>(group, solverName_ + ".TargetSteps", 10));
1178 1098 setMinSteps(getParamFromGroup<int>(group, solverName_ + ".MinSteps", 2));
1179 1098 setMaxSteps(getParamFromGroup<int>(group, solverName_ + ".MaxSteps", 18));
1180
1181 1098 enablePartialReassembly_ = getParamFromGroup<bool>(group, solverName_ + ".EnablePartialReassembly", false);
1182 1098 reassemblyMinThreshold_ = getParamFromGroup<Scalar>(group, solverName_ + ".ReassemblyMinThreshold", 1e-1*shiftTolerance_);
1183 1098 reassemblyMaxThreshold_ = getParamFromGroup<Scalar>(group, solverName_ + ".ReassemblyMaxThreshold", 1e2*shiftTolerance_);
1184 1098 reassemblyShiftWeight_ = getParamFromGroup<Scalar>(group, solverName_ + ".ReassemblyShiftWeight", 1e-3);
1185
1186
1/2
✓ Branch 2 taken 520 times.
✗ Branch 3 not taken.
 549 maxTimeStepDivisions_ = getParamFromGroup<std::size_t>(group, solverName_ + ".MaxTimeStepDivisions", 10);
1187 1098 retryTimeStepReductionFactor_ = getParamFromGroup<Scalar>(group, solverName_ + ".RetryTimeStepReductionFactor", 0.5);
1188
1189 549 numSteps_ = 0;
1190
1191 // output a parameter report
1192
2/2
✓ Branch 0 taken 492 times.
✓ Branch 1 taken 28 times.
 549 if (verbosity_ >= 2)
1193 521 reportParams();
1194 549 }
1195
1196 template<class SolA, class SolB>
1197 5103 void updateDistanceFromLastLinearization_(const SolA& u, const SolB& uDelta)
1198 {
1199 if constexpr (Dune::IsNumber<SolA>::value)
1200 {
1201 ✗ auto nextPriVars = u;
1202 ✗ nextPriVars -= uDelta;
1203
1204 // add the current relative shift for this degree of freedom
1205 ✗ auto shift = Detail::Newton::maxRelativeShift<Scalar>(u, nextPriVars);
1206 ✗ distanceFromLastLinearization_[0] += shift;
1207 }
1208 else
1209 {
1210
2/2
✓ Branch 0 taken 5898178 times.
✓ Branch 1 taken 5103 times.
 5903281 for (std::size_t i = 0; i < u.size(); ++i)
1211 {
1212
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 5898178 const auto& currentPriVars(u[i]);
1213 5898178 auto nextPriVars(currentPriVars);
1214
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 5898178 nextPriVars -= uDelta[i];
1215
1216 // add the current relative shift for this degree of freedom
1217 5898178 auto shift = Detail::Newton::maxRelativeShift<Scalar>(currentPriVars, nextPriVars);
1218 5898178 distanceFromLastLinearization_[i] += shift;
1219 }
1220 }
1221 5103 }
1222
1223 template<class ...ArgsA, class...ArgsB>
1224 ✗ void updateDistanceFromLastLinearization_(const Dune::MultiTypeBlockVector<ArgsA...>& uLastIter,
1225 const Dune::MultiTypeBlockVector<ArgsB...>& deltaU)
1226 {
1227 ✗ DUNE_THROW(Dune::NotImplemented, "Reassembly for MultiTypeBlockVector");
1228 }
1229
1230 template<class Sol>
1231 1293 void resizeDistanceFromLastLinearization_(const Sol& u, std::vector<Scalar>& dist)
1232 {
1233 1293 dist.assign(Backend::size(u), 0.0);
1234 1293 }
1235
1236 template<class ...Args>
1237 ✗ void resizeDistanceFromLastLinearization_(const Dune::MultiTypeBlockVector<Args...>& u,
1238 std::vector<Scalar>& dist)
1239 {
1240 ✗ DUNE_THROW(Dune::NotImplemented, "Reassembly for MultiTypeBlockVector");
1241 }
1242
1243 //! The communication object
1244 Communication comm_;
1245
1246 //! the verbosity level
1247 int verbosity_;
1248
1249 Scalar shiftTolerance_;
1250 Scalar reductionTolerance_;
1251 Scalar residualTolerance_;
1252
1253 // time step control
1254 std::size_t maxTimeStepDivisions_;
1255 Scalar retryTimeStepReductionFactor_;
1256
1257 // further parameters
1258 bool useLineSearch_;
1259 Scalar lineSearchMinRelaxationFactor_;
1260 bool useChop_;
1261 bool enableAbsoluteResidualCriterion_;
1262 bool enableShiftCriterion_;
1263 bool enableResidualCriterion_;
1264 bool satisfyResidualAndShiftCriterion_;
1265 bool enableDynamicOutput_;
1266
1267 //! the parameter group problem prefix for getting parameters from the parameter tree
1268 std::string paramGroup_;
1269 //! the parameter group for getting parameters from the parameter tree
1270 std::string solverName_;
1271
1272 // infrastructure for partial reassembly
1273 bool enablePartialReassembly_;
1274 std::unique_ptr<Reassembler> partialReassembler_;
1275 std::vector<Scalar> distanceFromLastLinearization_;
1276 Scalar reassemblyMinThreshold_;
1277 Scalar reassemblyMaxThreshold_;
1278 Scalar reassemblyShiftWeight_;
1279
1280 // statistics for the optional report
1281 std::size_t totalWastedIter_ = 0; //! Newton steps in solves that didn't converge
1282 std::size_t totalSucceededIter_ = 0; //! Newton steps in solves that converged
1283 std::size_t numConverged_ = 0; //! total number of converged solves
1284 std::size_t numLinearSolverBreakdowns_ = 0; //! total number of linear solves that failed
1285
1286 //! the class handling the primary variable switch
1287 std::unique_ptr<PrimaryVariableSwitchAdapter> priVarSwitchAdapter_;
1288
1289 //! convergence writer
1290 std::shared_ptr<ConvergenceWriter> convergenceWriter_ = nullptr;
1291 };
1292
1293 } // end namespace Dumux
1294
1295 #endif
1296