GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/porousmediumflow/richards/newtonsolver.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 27 28 96.4%
Functions: 9 19 47.4%
Branches: 33 57 57.9%
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 RichardsModel
10 * \brief A Richards model Newton solver.
11 */
12
13 #ifndef DUMUX_RICHARDS_NEWTON_SOLVER_HH
14 #define DUMUX_RICHARDS_NEWTON_SOLVER_HH
15
16 #include <algorithm>
17
18 #include <dune/common/parallel/communication.hh>
19 #include <dune/common/parallel/mpihelper.hh>
20
21 #include <dumux/common/properties.hh>
22 #include <dumux/assembly/partialreassembler.hh>
23 #include <dumux/nonlinear/newtonsolver.hh>
24 #include <dumux/discretization/elementsolution.hh>
25
26 namespace Dumux {
27 /*!
28 * \ingroup RichardsModel
29 * \brief A Richards model specific Newton solver.
30 *
31 * This solver 'knows' what a 'physically meaningful' solution is
32 * and can thus do update smarter than the plain Newton solver.
33 *
34 */
35 template <class Assembler, class LinearSolver,
36 class Reassembler = PartialReassembler<Assembler>,
37 class Comm = Dune::Communication<Dune::MPIHelper::MPICommunicator> >
38
1/2
✓ Branch 1 taken 41 times.
✗ Branch 2 not taken.
 41 class RichardsNewtonSolver : public NewtonSolver<Assembler, LinearSolver, Reassembler, Comm>
39 {
40 using Scalar = typename Assembler::Scalar;
41 using ParentType = NewtonSolver<Assembler, LinearSolver, Reassembler, Comm>;
42 using Indices = typename Assembler::GridVariables::VolumeVariables::Indices;
43 enum { pressureIdx = Indices::pressureIdx };
44
45 using typename ParentType::Backend;
46 using typename ParentType::SolutionVector;
47 using typename ParentType::ResidualVector;
48
49 public:
50
2/4
✓ Branch 1 taken 41 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 41 times.
✗ Branch 4 not taken.
 41 using ParentType::ParentType;
51 using typename ParentType::Variables;
52
53 private:
54
55 /*!
56 * \brief Update the current solution of the Newton method
57 *
58 * \param varsCurrentIter The variables after the current Newton iteration \f$ u^{k+1} \f$
59 * \param uLastIter The solution after the last Newton iteration \f$ u^k \f$
60 * \param deltaU The vector of differences between the last
61 * iterative solution and the next one \f$ \Delta u^k \f$
62 */
63 975 void choppedUpdate_(Variables &varsCurrentIter,
64 const SolutionVector &uLastIter,
65 const ResidualVector &deltaU) final
66 {
67 975 auto uCurrentIter = uLastIter;
68 975 Backend::axpy(-1.0, deltaU, uCurrentIter);
69
70 // do not clamp anything after 5 iterations
71
2/2
✓ Branch 0 taken 791 times.
✓ Branch 1 taken 184 times.
 975 if (this->numSteps_ <= 4)
72 {
73 // clamp saturation change to at most 20% per iteration
74
1/2
✓ Branch 1 taken 791 times.
✗ Branch 2 not taken.
 791 const auto& gridGeometry = this->assembler().gridGeometry();
75 791 auto fvGeometry = localView(gridGeometry);
76
10/15
✓ Branch 1 taken 791 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 177331 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 162314 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 28574 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 13462 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 13392 times.
✓ Branch 16 taken 62 times.
✓ Branch 12 taken 15120 times.
✓ Branch 6 taken 14560 times.
✓ Branch 9 taken 12896 times.
 574163 for (const auto& element : elements(gridGeometry.gridView()))
77 {
78
1/2
✓ Branch 1 taken 85680 times.
✗ Branch 2 not taken.
 324688 fvGeometry.bindElement(element);
79
80
5/5
✓ Branch 1 taken 205184 times.
✓ Branch 2 taken 93216 times.
✓ Branch 3 taken 26288 times.
✓ Branch 4 taken 26288 times.
✓ Branch 0 taken 342720 times.
 667408 for (auto&& scv : scvs(fvGeometry))
81 {
82 462224 auto dofIdxGlobal = scv.dofIndex();
83
84 // calculate the old wetting phase saturation
85
1/2
✓ Branch 1 taken 26288 times.
✗ Branch 2 not taken.
 462224 const auto& spatialParams = this->assembler().problem().spatialParams();
86
2/2
✓ Branch 0 taken 574448 times.
✓ Branch 1 taken 65126 times.
 732790 const auto elemSol = elementSolution(element, uCurrentIter, gridGeometry);
87
88
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 462224 const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol);
89 462224 const Scalar pcMin = fluidMatrixInteraction.pc(1.0);
90
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 462224 times.
 462224 const Scalar pw = uLastIter[dofIdxGlobal][pressureIdx];
91 using std::max;
92
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 462224 times.
 462224 const Scalar pn = max(this->assembler().problem().nonwettingReferencePressure(), pw + pcMin);
93
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 462224 times.
✗ Branch 2 not taken.
 462224 const Scalar pcOld = pn - pw;
94
2/2
✓ Branch 0 taken 383809 times.
✓ Branch 1 taken 78415 times.
 462224 const Scalar SwOld = max(0.0, fluidMatrixInteraction.sw(pcOld));
95
96 // convert into minimum and maximum wetting phase pressures
97
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 462224 const Scalar pwMin = pn - fluidMatrixInteraction.pc(SwOld - 0.2);
98 462224 const Scalar pwMax = pn - fluidMatrixInteraction.pc(SwOld + 0.2);
99
100 // clamp the result
101 using std::clamp;
102
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 462224 times.
 924448 uCurrentIter[dofIdxGlobal][pressureIdx] = clamp(uCurrentIter[dofIdxGlobal][pressureIdx], pwMin, pwMax);
103 }
104 }
105 791 }
106
107 // update the variables
108
1/2
✓ Branch 1 taken 975 times.
✗ Branch 2 not taken.
 975 this->solutionChanged_(varsCurrentIter, uCurrentIter);
109
110
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 975 times.
 975 if (this->enableResidualCriterion())
111 ✗ this->computeResidualReduction_(varsCurrentIter);
112 975 }
113 };
114
115 } // end namespace Dumux
116
117 #endif
118