| 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 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 39 times.
✗ Branch 2 not taken.
|
39 | 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 39 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 39 times.
✗ Branch 4 not taken.
|
39 | 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 | 931 | void choppedUpdate_(Variables &varsCurrentIter, | |
| 64 | const SolutionVector &uLastIter, | ||
| 65 | const ResidualVector &deltaU) final | ||
| 66 | { | ||
| 67 |
1/4✓ Branch 1 taken 931 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
1862 | auto uCurrentIter = uLastIter; |
| 68 | 931 | Backend::axpy(-1.0, deltaU, uCurrentIter); | |
| 69 | |||
| 70 | // do not clamp anything after 5 iterations | ||
| 71 |
2/2✓ Branch 0 taken 753 times.
✓ Branch 1 taken 178 times.
|
931 | if (this->numSteps_ <= 4) |
| 72 | { | ||
| 73 | // clamp saturation change to at most 20% per iteration | ||
| 74 |
2/4✓ Branch 1 taken 753 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 753 times.
✗ Branch 5 not taken.
|
1506 | const auto& gridGeometry = this->assembler().gridGeometry(); |
| 75 |
1/2✓ Branch 1 taken 753 times.
✗ Branch 2 not taken.
|
885 | auto fvGeometry = localView(gridGeometry); |
| 76 |
14/18✓ Branch 1 taken 753 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 753 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 169257 times.
✓ Branch 7 taken 264 times.
✓ Branch 8 taken 27456 times.
✓ Branch 9 taken 132 times.
✓ Branch 10 taken 141444 times.
✓ Branch 11 taken 27456 times.
✓ Branch 12 taken 28512 times.
✓ Branch 13 taken 132 times.
✓ Branch 14 taken 55968 times.
✓ Branch 15 taken 132 times.
✓ Branch 17 taken 28512 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 28512 times.
✗ Branch 21 not taken.
|
452919 | for (const auto& element : elements(gridGeometry.gridView())) |
| 77 | { | ||
| 78 |
1/2✓ Branch 1 taken 197280 times.
✗ Branch 2 not taken.
|
197280 | fvGeometry.bindElement(element); |
| 79 | |||
| 80 |
4/4✓ Branch 0 taken 454320 times.
✓ Branch 1 taken 197280 times.
✓ Branch 2 taken 342720 times.
✓ Branch 3 taken 85680 times.
|
1080000 | for (auto&& scv : scvs(fvGeometry)) |
| 81 | { | ||
| 82 |
1/2✓ Branch 1 taken 26288 times.
✗ Branch 2 not taken.
|
454320 | auto dofIdxGlobal = scv.dofIndex(); |
| 83 | |||
| 84 | // calculate the old wetting phase saturation | ||
| 85 |
3/6✓ Branch 1 taken 26288 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26288 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26288 times.
✗ Branch 8 not taken.
|
1362960 | const auto& spatialParams = this->assembler().problem().spatialParams(); |
| 86 |
1/2✓ Branch 1 taken 26288 times.
✗ Branch 2 not taken.
|
454320 | const auto elemSol = elementSolution(element, uCurrentIter, gridGeometry); |
| 87 | |||
| 88 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
454320 | const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol); |
| 89 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
454320 | const Scalar pcMin = fluidMatrixInteraction.pc(1.0); |
| 90 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 454320 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
454320 | const Scalar pw = uLastIter[dofIdxGlobal][pressureIdx]; |
| 91 | using std::max; | ||
| 92 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 454320 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 454320 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 454320 times.
|
1362960 | const Scalar pn = max(this->assembler().problem().nonwettingReferencePressure(), pw + pcMin); |
| 93 | 454320 | const Scalar pcOld = pn - pw; | |
| 94 |
3/5✗ Branch 0 not taken.
✓ Branch 1 taken 454320 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 376978 times.
✓ Branch 4 taken 77342 times.
|
454320 | 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.
|
454320 | const Scalar pwMin = pn - fluidMatrixInteraction.pc(SwOld - 0.2); |
| 98 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
454320 | const Scalar pwMax = pn - fluidMatrixInteraction.pc(SwOld + 0.2); |
| 99 | |||
| 100 | // clamp the result | ||
| 101 | using std::clamp; | ||
| 102 |
3/6✓ Branch 0 taken 454320 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 454320 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 454320 times.
✗ Branch 5 not taken.
|
1817280 | uCurrentIter[dofIdxGlobal][pressureIdx] = clamp(uCurrentIter[dofIdxGlobal][pressureIdx], pwMin, pwMax); |
| 103 | } | ||
| 104 | } | ||
| 105 | } | ||
| 106 | |||
| 107 | // update the variables | ||
| 108 |
1/2✓ Branch 1 taken 931 times.
✗ Branch 2 not taken.
|
931 | this->solutionChanged_(varsCurrentIter, uCurrentIter); |
| 109 | |||
| 110 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 931 times.
|
931 | if (this->enableResidualCriterion()) |
| 111 | ✗ | this->computeResidualReduction_(varsCurrentIter); | |
| 112 | 931 | } | |
| 113 | }; | ||
| 114 | |||
| 115 | } // end namespace Dumux | ||
| 116 | |||
| 117 | #endif | ||
| 118 |