GCC Code Coverage Report

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