| 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 PNMTwoPModel | ||
| 10 | * \copydoc Dumux::PoreNetwork::TwoPNewtonConsistencyChecks | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_PNM_NEWTON_CONSISTENCY_CHECKS | ||
| 13 | #define DUMUX_PNM_NEWTON_CONSISTENCY_CHECKS | ||
| 14 | |||
| 15 | #include <vector> | ||
| 16 | #include <iostream> | ||
| 17 | #include <dune/common/exceptions.hh> | ||
| 18 | #include <dumux/common/exceptions.hh> | ||
| 19 | #include <dumux/common/parameters.hh> | ||
| 20 | #include <dumux/discretization/localview.hh> | ||
| 21 | |||
| 22 | namespace Dumux::PoreNetwork { | ||
| 23 | |||
| 24 | #ifndef DOXYGEN | ||
| 25 | namespace Detail { | ||
| 26 | |||
| 27 | // Primary template | ||
| 28 | template <typename T, typename U = int> | ||
| 29 | struct SaturationIndex {}; | ||
| 30 | |||
| 31 | // Specialization for 2p | ||
| 32 | template <typename T> | ||
| 33 | struct SaturationIndex <T, decltype((void) T::saturationIdx, 0)> | ||
| 34 | { | ||
| 35 | static constexpr auto value = T::saturationIdx; | ||
| 36 | }; | ||
| 37 | |||
| 38 | // Specialization for 2pnc | ||
| 39 | template <typename T> | ||
| 40 | struct SaturationIndex <T, decltype((void) T::switchIdx, 0)> | ||
| 41 | { | ||
| 42 | static constexpr auto value = T::switchIdx; | ||
| 43 | }; | ||
| 44 | |||
| 45 | // Specialization for MPNC | ||
| 46 | template <typename T> | ||
| 47 | struct SaturationIndex <T, decltype((void) T::s0Idx, 0)> | ||
| 48 | { | ||
| 49 | static constexpr auto value = T::s0Idx; | ||
| 50 | }; | ||
| 51 | } // end namespace Detail | ||
| 52 | #endif | ||
| 53 | |||
| 54 | /*! | ||
| 55 | * \ingroup PNMTwoPModel | ||
| 56 | * \brief Consistency checks for the PNM two-phase model | ||
| 57 | */ | ||
| 58 | template<class GridVariables, class SolutionVector> | ||
| 59 | class TwoPNewtonConsistencyChecks | ||
| 60 | { | ||
| 61 | public: | ||
| 62 | |||
| 63 | /*! | ||
| 64 | * \brief Perform all checks. | ||
| 65 | */ | ||
| 66 | 451 | void performChecks(const GridVariables& gridVariables, const SolutionVector& uCurrentIter, const SolutionVector& prevSol) const | |
| 67 | { | ||
| 68 | 451 | checkRelativeSaturationShift(gridVariables, uCurrentIter, prevSol); | |
| 69 | 451 | checkIfValuesArePhysical(gridVariables, uCurrentIter); | |
| 70 | 902 | checkIfCapillaryPressureIsCloseToEntryPressure(gridVariables, uCurrentIter); | |
| 71 | 451 | } | |
| 72 | |||
| 73 | /*! | ||
| 74 | * \brief Checks if the relative shift of saturation between to consecutive time steps is below a given threshold. | ||
| 75 | */ | ||
| 76 | 451 | void checkRelativeSaturationShift(const GridVariables& gridVariables, const SolutionVector& uCurrentIter, const SolutionVector& prevSol) const | |
| 77 | { | ||
| 78 | using Scalar = typename SolutionVector::block_type::value_type; | ||
| 79 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 447 times.
|
451 | const auto& problem = gridVariables.curGridVolVars().problem(); |
| 80 | |||
| 81 |
5/8✓ Branch 0 taken 4 times.
✓ Branch 1 taken 447 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
|
451 | static const Scalar allowedSaturationChange = getParamFromGroup<Scalar>(problem.paramGroup(), "Newton.AllowedSaturationChange", -1.0); |
| 82 |
1/2✓ Branch 0 taken 451 times.
✗ Branch 1 not taken.
|
451 | if (allowedSaturationChange < 0.0) |
| 83 | 451 | return; | |
| 84 | |||
| 85 | ✗ | const auto& curGridVolVars = gridVariables.curGridVolVars(); | |
| 86 | ✗ | const auto& prevGridVolVars = gridVariables.prevGridVolVars(); | |
| 87 | |||
| 88 | ✗ | std::vector<bool> dofVisited(uCurrentIter.size(), false); | |
| 89 | |||
| 90 | ✗ | auto fvGeometry = localView(problem.gridGeometry()); | |
| 91 | ✗ | auto curElemVolVars = localView(curGridVolVars); | |
| 92 | ✗ | auto prevElemVolVars = localView(prevGridVolVars); | |
| 93 | |||
| 94 | ✗ | for (const auto& element : elements(problem.gridGeometry().gridView())) | |
| 95 | { | ||
| 96 | ✗ | fvGeometry.bindElement(element); | |
| 97 | ✗ | curElemVolVars.bindElement(element, fvGeometry, uCurrentIter); | |
| 98 | ✗ | prevElemVolVars.bindElement(element, fvGeometry, prevSol); | |
| 99 | |||
| 100 | ✗ | for (const auto& scv : scvs(fvGeometry)) | |
| 101 | { | ||
| 102 | ✗ | if (dofVisited[scv.dofIndex()]) | |
| 103 | continue; | ||
| 104 | ✗ | dofVisited[scv.dofIndex()] = true; | |
| 105 | |||
| 106 | ✗ | const Scalar satNew = curElemVolVars[scv].saturation(0); | |
| 107 | ✗ | const Scalar satOld = prevElemVolVars[scv].saturation(0); | |
| 108 | using std::abs; | ||
| 109 | ✗ | const Scalar deltaS = abs(satNew - satOld); | |
| 110 | ✗ | static const bool considerRelativeShift = getParamFromGroup<bool>(problem.paramGroup(), "Newton.SaturationChangeIsRelative", false); | |
| 111 | |||
| 112 | ✗ | if (considerRelativeShift) | |
| 113 | { | ||
| 114 | // satOld mgiht be (close to) zero, so have to take care of this | ||
| 115 | ✗ | if (satOld > 1e-3 && deltaS / satOld > allowedSaturationChange) | |
| 116 | ✗ | DUNE_THROW(NumericalProblem, "Saturation change too high at dof " << scv.dofIndex() << ", old sat. " << satOld << ", new sat. " << satNew << std::endl); | |
| 117 | } | ||
| 118 | ✗ | else if (deltaS > allowedSaturationChange) | |
| 119 | ✗ | DUNE_THROW(NumericalProblem, "Saturation change too high at dof " << scv.dofIndex() << ", old sat. " << satOld << ", new sat. " << satNew << std::endl); | |
| 120 | } | ||
| 121 | } | ||
| 122 | } | ||
| 123 | |||
| 124 | /*! | ||
| 125 | * \brief Checks if the saturation is between zero and one. | ||
| 126 | */ | ||
| 127 | 451 | void checkIfValuesArePhysical(const GridVariables& gridVariables, const SolutionVector& uCurrentIter) const | |
| 128 | { | ||
| 129 |
4/4✓ Branch 0 taken 4 times.
✓ Branch 1 taken 447 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 447 times.
|
902 | const auto& problem = gridVariables.curGridVolVars().problem(); |
| 130 | |||
| 131 |
5/8✓ Branch 0 taken 4 times.
✓ Branch 1 taken 447 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
|
451 | static const bool doPlausibilityCheck = getParamFromGroup<bool>(problem.paramGroup(), "Newton.PlausibilityCheck", false); |
| 132 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 451 times.
|
451 | if (!doPlausibilityCheck) |
| 133 | return; | ||
| 134 | |||
| 135 | ✗ | for (std::size_t i = 0; i < uCurrentIter.size(); ++i) | |
| 136 | { | ||
| 137 | ✗ | const auto& priVars = uCurrentIter[i]; | |
| 138 | using Indices = typename GridVariables::VolumeVariables::Indices; | ||
| 139 | static constexpr auto saturationOrMoleFractionIndex = Detail::SaturationIndex<Indices>::value; | ||
| 140 | ✗ | if (priVars[saturationOrMoleFractionIndex] < 0.0 || priVars[saturationOrMoleFractionIndex] > 1.0) | |
| 141 | { | ||
| 142 | ✗ | std::cout << "at dof " << i << ": saturation " << priVars[saturationOrMoleFractionIndex] << std::endl; | |
| 143 | ✗ | DUNE_THROW(NumericalProblem, "Saturation is below 0 or above 1"); | |
| 144 | } | ||
| 145 | } | ||
| 146 | } | ||
| 147 | |||
| 148 | /*! | ||
| 149 | * \brief Checks if the capillary pressure at pore from which a throat was invaded is sufficiently close to the throat's entry capillary pressure. | ||
| 150 | */ | ||
| 151 | ✗ | void checkIfCapillaryPressureIsCloseToEntryPressure(const GridVariables& gridVariables, const SolutionVector& uCurrentIter) const | |
| 152 | { | ||
| 153 | // this check is implemented in the invasion state itself | ||
| 154 | 902 | const auto& invasionState = gridVariables.gridFluxVarsCache().invasionState(); | |
| 155 | 1353 | invasionState.checkIfCapillaryPressureIsCloseToEntryPressure(uCurrentIter, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache()); | |
| 156 | ✗ | } | |
| 157 | }; | ||
| 158 | } // end namespace Dumux | ||
| 159 | #endif | ||
| 160 |