| 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 | * \brief Invasion state class for the two-phase PNM. | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_PNM_2P_INVASIONSTATE_HH | ||
| 13 | #define DUMUX_PNM_2P_INVASIONSTATE_HH | ||
| 14 | |||
| 15 | #include <vector> | ||
| 16 | #include <type_traits> | ||
| 17 | #include <dune/common/std/type_traits.hh> | ||
| 18 | #include <dumux/common/parameters.hh> | ||
| 19 | #include <dumux/porenetwork/common/labels.hh> | ||
| 20 | |||
| 21 | namespace Dumux::PoreNetwork { | ||
| 22 | |||
| 23 | /*! | ||
| 24 | * \ingroup PNMTwoPModel | ||
| 25 | * \brief This class updates the invasion state for the two-phase PNM. | ||
| 26 | */ | ||
| 27 | template<class P> | ||
| 28 | class TwoPInvasionState | ||
| 29 | { | ||
| 30 | using Problem = P; | ||
| 31 | |||
| 32 | template <class T> | ||
| 33 | using GlobalCapillaryPressureDetector = decltype(std::declval<T>().globalCapillaryPressure()); | ||
| 34 | |||
| 35 | template<class T> | ||
| 36 | static constexpr bool hasGlobalCapillaryPressure() | ||
| 37 | { return Dune::Std::is_detected<GlobalCapillaryPressureDetector, T>::value; } | ||
| 38 | |||
| 39 | enum class EventType {invasion, snapOff, none}; | ||
| 40 | |||
| 41 | public: | ||
| 42 | |||
| 43 | 8 | TwoPInvasionState(const Problem& problem) : problem_(problem) | |
| 44 | { | ||
| 45 | // initialize the invasion state | ||
| 46 |
1/2✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
12 | invadedCurrentIteration_.resize(problem.gridGeometry().gridView().size(0)); |
| 47 |
1/2✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
12 | invadedPreviousTimeStep_.resize(problem.gridGeometry().gridView().size(0)); |
| 48 | |||
| 49 |
4/4✓ Branch 3 taken 282 times.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 282 times.
✓ Branch 6 taken 4 times.
|
576 | for (auto&& element : elements(problem.gridGeometry().gridView())) |
| 50 | { | ||
| 51 | 846 | const auto eIdx = problem.gridGeometry().elementMapper().index(element); | |
| 52 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 282 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 282 times.
|
564 | invadedCurrentIteration_[eIdx] = problem.initialInvasionState(element); |
| 53 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 282 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 282 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 282 times.
|
1128 | invadedPreviousTimeStep_[eIdx] = invadedCurrentIteration_[eIdx]; |
| 54 | } | ||
| 55 | |||
| 56 | 12 | numThroatsInvaded_ = std::count(invadedCurrentIteration_.begin(), invadedCurrentIteration_.end(), true); | |
| 57 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
8 | verbose_ = getParamFromGroup<bool>(problem.paramGroup(), "InvasionState.Verbosity", true); |
| 58 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | restrictToGlobalCapillaryPressure_ = getParamFromGroup<bool>(problem.paramGroup(), "InvasionState.RestrictInvasionToGlobalCapillaryPressure", false); |
| 59 | |||
| 60 | if constexpr (hasGlobalCapillaryPressure<Problem>()) | ||
| 61 | { | ||
| 62 | if (restrictToGlobalCapillaryPressure_) | ||
| 63 | std::cout << "\n *** Invasion behavior is restricted by a global capillary pressure defined in the problem! *** \n" << std::endl; | ||
| 64 | else | ||
| 65 | std::cout << "\n *** WARNING: global capillary pressure defined in the problem but InvasionState.RestrictInvasionToGlobalCapillaryPressure is set to false.\n" | ||
| 66 | << " Invasion behavior will NOT be restricted! ***\n" << std::endl; | ||
| 67 | } | ||
| 68 | 4 | } | |
| 69 | |||
| 70 | //! Return whether a given throat is invaded or not. | ||
| 71 | template<class Element> | ||
| 72 | 959820 | bool invaded(const Element& element) const | |
| 73 | { | ||
| 74 | 2879460 | const auto eIdx = problem_.gridGeometry().elementMapper().index(element); | |
| 75 | 1919640 | return invadedCurrentIteration_[eIdx]; | |
| 76 | } | ||
| 77 | |||
| 78 | //! Return the number of currently invaded throats | ||
| 79 | std::size_t numThroatsInvaded() const | ||
| 80 | { return numThroatsInvaded_; } | ||
| 81 | |||
| 82 | //! Update the invasion state of all throats. This is done after each Newton step by a call from the Newton solver. | ||
| 83 | template<class SolutionVector, class GridVolumeVariables, class GridFluxVariablesCache> | ||
| 84 | 1619 | bool update(const SolutionVector& sol, const GridVolumeVariables& gridVolVars, GridFluxVariablesCache& gridFluxVarsCache) | |
| 85 | { | ||
| 86 | 1619 | hasChangedInCurrentIteration_ = false; | |
| 87 | 3238 | auto fvGeometry = localView(problem_.gridGeometry()); | |
| 88 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
1619 | auto elemVolVars = localView(gridVolVars); |
| 89 | 1619 | auto elemFluxVarsCache = localView(gridFluxVarsCache); | |
| 90 |
5/6✓ Branch 3 taken 124713 times.
✓ Branch 4 taken 1619 times.
✓ Branch 5 taken 124713 times.
✓ Branch 6 taken 1619 times.
✓ Branch 8 taken 124713 times.
✗ Branch 9 not taken.
|
129570 | for (auto&& element : elements(problem_.gridGeometry().gridView())) |
| 91 | { | ||
| 92 |
1/2✓ Branch 1 taken 124713 times.
✗ Branch 2 not taken.
|
124713 | fvGeometry.bindElement(element); |
| 93 |
1/2✓ Branch 1 taken 124713 times.
✗ Branch 2 not taken.
|
124713 | elemVolVars.bind(element, fvGeometry, sol); |
| 94 |
1/2✓ Branch 1 taken 124713 times.
✗ Branch 2 not taken.
|
124713 | elemFluxVarsCache.bind(element, fvGeometry, elemVolVars); |
| 95 | |||
| 96 |
2/2✓ Branch 0 taken 124713 times.
✓ Branch 1 taken 124713 times.
|
249426 | for (auto&& scvf : scvfs(fvGeometry)) |
| 97 | { | ||
| 98 | // checks if invasion or snap-off occurred after Newton iteration step | ||
| 99 |
3/4✓ Branch 1 taken 124713 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 417 times.
✓ Branch 4 taken 124296 times.
|
124713 | if (const auto invasionResult = invasionSwitch_(element, elemVolVars, elemFluxVarsCache[scvf]); invasionResult) |
| 100 | { | ||
| 101 | 417 | hasChangedInCurrentIteration_ = true; | |
| 102 | if constexpr (GridFluxVariablesCache::cachingEnabled) | ||
| 103 | { | ||
| 104 | const auto eIdx = problem_.gridGeometry().elementMapper().index(element); | ||
| 105 | gridFluxVarsCache.cache(eIdx, scvf.index()).update(problem_, element, fvGeometry, elemVolVars, scvf, invadedCurrentIteration_[eIdx]); | ||
| 106 | } | ||
| 107 | } | ||
| 108 | } | ||
| 109 | } | ||
| 110 | 4857 | numThroatsInvaded_ = std::count(invadedCurrentIteration_.begin(), invadedCurrentIteration_.end(), true); | |
| 111 |
1/2✓ Branch 0 taken 1619 times.
✗ Branch 1 not taken.
|
3238 | return hasChangedInCurrentIteration_; |
| 112 | } | ||
| 113 | |||
| 114 | //! Restore the old invasion state after a Newton iteration has failed. | ||
| 115 | void reset() | ||
| 116 | { | ||
| 117 | 36 | hasChangedInCurrentIteration_ = false; | |
| 118 | 36 | invadedCurrentIteration_ = invadedPreviousTimeStep_; | |
| 119 | } | ||
| 120 | |||
| 121 | //! Return whether an invasion or snap-off occurred anywhere. Can be used, e.g., for output file writing control. | ||
| 122 | ✗ | bool hasChanged() const | |
| 123 | ✗ | { return hasChangedComparedToPreviousTimestep_; } | |
| 124 | |||
| 125 | //! Return whether an invasion or snap-off occurred anywhere during the current Newton iteration. | ||
| 126 | bool hasChangedInCurrentIteration() const | ||
| 127 | { return hasChangedInCurrentIteration_; } | ||
| 128 | |||
| 129 | //! This is called after the Newton method has successfully finished one time step. | ||
| 130 | 432 | void advance() | |
| 131 | { | ||
| 132 | 432 | hasChangedComparedToPreviousTimestep_ = (invadedPreviousTimeStep_ != invadedCurrentIteration_); | |
| 133 | 432 | invadedPreviousTimeStep_ = invadedCurrentIteration_; | |
| 134 | 432 | } | |
| 135 | |||
| 136 | template<class SolutionVector, class GridVolumeVariables, class GridFluxVariablesCache> | ||
| 137 | 451 | void checkIfCapillaryPressureIsCloseToEntryPressure(const SolutionVector& sol, | |
| 138 | const GridVolumeVariables& gridVolVars, | ||
| 139 | const GridFluxVariablesCache& gridFluxVarsCache) const | ||
| 140 | { | ||
| 141 | using Scalar = typename SolutionVector::block_type::value_type; | ||
| 142 |
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 accuracyCriterion = getParamFromGroup<Scalar>(problem_.paramGroup(), "InvasionState.AccuracyCriterion", -1.0); |
| 143 | |||
| 144 |
1/2✓ Branch 0 taken 451 times.
✗ Branch 1 not taken.
|
451 | if (accuracyCriterion < 0.0) |
| 145 | 451 | return; | |
| 146 | |||
| 147 | ✗ | auto fvGeometry = localView(problem_.gridGeometry()); | |
| 148 | ✗ | auto elemVolVars = localView(gridVolVars); | |
| 149 | ✗ | auto elemFluxVarsCache = localView(gridFluxVarsCache); | |
| 150 | ✗ | for (auto&& element : elements(problem_.gridGeometry().gridView())) | |
| 151 | { | ||
| 152 | // Only consider throats which have been invaded during the current time step | ||
| 153 | ✗ | const auto eIdx = problem_.gridGeometry().elementMapper().index(element); | |
| 154 | ✗ | if (!invadedCurrentIteration_[eIdx] || invadedPreviousTimeStep_[eIdx] == invadedCurrentIteration_[eIdx]) | |
| 155 | continue; | ||
| 156 | |||
| 157 | ✗ | fvGeometry.bindElement(element); | |
| 158 | ✗ | elemVolVars.bind(element, fvGeometry, sol); | |
| 159 | ✗ | elemFluxVarsCache.bind(element, fvGeometry, elemVolVars); | |
| 160 | |||
| 161 | ✗ | for (auto&& scvf : scvfs(fvGeometry)) | |
| 162 | { | ||
| 163 | // checks if pc is close enough to the entry pressure value | ||
| 164 | ✗ | const auto& fluxVarsCache = elemFluxVarsCache[scvf]; | |
| 165 | |||
| 166 | using std::max; | ||
| 167 | ✗ | const Scalar pc = max(elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure()); | |
| 168 | |||
| 169 | ✗ | if (pc < accuracyCriterion * fluxVarsCache.pcEntry()) | |
| 170 | ✗ | DUNE_THROW(NumericalProblem, "At element " << eIdx << ": pc " << pc << " too far away form pcEntry " << fluxVarsCache.pcEntry()); | |
| 171 | } | ||
| 172 | } | ||
| 173 | } | ||
| 174 | |||
| 175 | private: | ||
| 176 | |||
| 177 | //! The switch for determining the invasion state of a pore throat. Called at the end of each Newton step. | ||
| 178 | template<class Element, class ElementVolumeVariables, class FluxVariablesCache> | ||
| 179 | 124713 | auto invasionSwitch_(const Element& element, | |
| 180 | const ElementVolumeVariables& elemVolVars, | ||
| 181 | const FluxVariablesCache& fluxVarsCache) | ||
| 182 | |||
| 183 | { | ||
| 184 | using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type; | ||
| 185 | 124713 | const auto& gridGeometry = problem_.gridGeometry(); | |
| 186 | 124713 | const auto& spatialParams = problem_.spatialParams(); | |
| 187 | 249426 | const auto eIdx = gridGeometry.elementMapper().index(element); | |
| 188 |
4/4✓ Branch 0 taken 4 times.
✓ Branch 1 taken 124709 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 124709 times.
|
249426 | bool invadedBeforeSwitch = invadedCurrentIteration_[eIdx]; |
| 189 | 124713 | bool invadedAfterSwitch = invadedBeforeSwitch; | |
| 190 | |||
| 191 | // Result type, containing the local scv index of the pore from which the invasion/snap-off occurred | ||
| 192 | // Evaluates to 'false' if no invasion/snap-off occurred | ||
| 193 | 417 | struct Result | |
| 194 | { | ||
| 195 | std::uint8_t localScvIdxWithCriticalPc; | ||
| 196 | Scalar criticalPc; | ||
| 197 | EventType event = EventType::none; | ||
| 198 | |||
| 199 | ✗ | operator bool() const | |
| 200 | ✗ | { return event != EventType::none; } | |
| 201 | }; | ||
| 202 | |||
| 203 | // Block non-wetting phase flux out of the outlet | ||
| 204 |
7/14✓ Branch 0 taken 4 times.
✓ Branch 1 taken 124709 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.
✓ Branch 12 taken 4 times.
✗ Branch 13 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 4 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
124713 | static const auto blockNonwettingPhase = getParamFromGroup<std::vector<int>>(problem_.paramGroup(), "InvasionState.BlockNonwettingPhaseAtThroatLabel", std::vector<int>{}); |
| 205 |
2/8✗ Branch 0 not taken.
✓ Branch 1 taken 124713 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 124713 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
249426 | if (!blockNonwettingPhase.empty() && std::find(blockNonwettingPhase.begin(), blockNonwettingPhase.end(), gridGeometry.throatLabel(eIdx)) != blockNonwettingPhase.end()) |
| 206 | { | ||
| 207 | ✗ | invadedCurrentIteration_[eIdx] = false; | |
| 208 | ✗ | return Result{}; // nothing happened | |
| 209 | } | ||
| 210 | |||
| 211 | //Determine whether throat gets invaded or snap-off occurs | ||
| 212 | 249426 | const std::array<Scalar, 2> pc = { elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure() }; | |
| 213 | 374139 | const auto pcMax = std::max_element(pc.begin(), pc.end()); | |
| 214 | 124713 | const Scalar pcEntry = fluxVarsCache.pcEntry(); | |
| 215 | 124713 | const Scalar pcSnapoff = fluxVarsCache.pcSnapoff(); | |
| 216 | |||
| 217 | // check if there is a user-specified global capillary pressure which needs to be obeyed | ||
| 218 | 124713 | if (maybeRestrictToGlobalCapillaryPressure_(pcEntry)) | |
| 219 | { | ||
| 220 | if (*pcMax > pcEntry) | ||
| 221 | { | ||
| 222 | std::cout << "Throat " << eIdx << " would have been invaded by pc of " << *pcMax << "but a global capillary pressure restricion was set in the problem."; | ||
| 223 | std::cout << ". pcEntry: " << spatialParams.pcEntry(element, elemVolVars) << std::endl; | ||
| 224 | } | ||
| 225 | |||
| 226 | invadedCurrentIteration_[eIdx] = false; | ||
| 227 | return Result{}; //nothing happened | ||
| 228 | } | ||
| 229 | |||
| 230 |
2/2✓ Branch 0 taken 122998 times.
✓ Branch 1 taken 1715 times.
|
124713 | if (*pcMax > pcEntry) |
| 231 | invadedAfterSwitch = true; | ||
| 232 |
2/2✓ Branch 0 taken 98117 times.
✓ Branch 1 taken 24881 times.
|
122998 | else if (*pcMax <= pcSnapoff) |
| 233 | 98117 | invadedAfterSwitch = false; | |
| 234 | |||
| 235 |
4/4✓ Branch 0 taken 12238 times.
✓ Branch 1 taken 112475 times.
✓ Branch 2 taken 12238 times.
✓ Branch 3 taken 112475 times.
|
249426 | invadedCurrentIteration_[eIdx] = invadedAfterSwitch; |
| 236 | |||
| 237 |
2/2✓ Branch 0 taken 124296 times.
✓ Branch 1 taken 417 times.
|
124713 | if (invadedBeforeSwitch == invadedAfterSwitch) |
| 238 | 124296 | return Result{}; // nothing happened | |
| 239 | else | ||
| 240 | { | ||
| 241 | 834 | Result result; | |
| 242 |
4/4✓ Branch 0 taken 126 times.
✓ Branch 1 taken 291 times.
✓ Branch 2 taken 126 times.
✓ Branch 3 taken 291 times.
|
834 | result.localScvIdxWithCriticalPc = std::distance(pc.begin(), pcMax); |
| 243 | 417 | result.criticalPc = *pcMax; | |
| 244 |
2/2✓ Branch 0 taken 126 times.
✓ Branch 1 taken 291 times.
|
417 | result.event = !invadedBeforeSwitch && invadedAfterSwitch ? EventType::invasion : EventType::snapOff; |
| 245 | |||
| 246 |
1/2✓ Branch 0 taken 417 times.
✗ Branch 1 not taken.
|
417 | if (verbose_) |
| 247 | { | ||
| 248 | 417 | const auto wPhaseIdx = spatialParams.template wettingPhase<typename ElementVolumeVariables::VolumeVariables::FluidSystem>(element, elemVolVars); | |
| 249 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 417 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 417 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 417 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 417 times.
|
1668 | const std::array sw = { elemVolVars[0].saturation(wPhaseIdx), elemVolVars[1].saturation(wPhaseIdx) }; |
| 250 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 417 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 417 times.
|
834 | const auto vIdx = gridGeometry.gridView().indexSet().subIndex(element, result.localScvIdxWithCriticalPc, 1); |
| 251 |
2/2✓ Branch 0 taken 291 times.
✓ Branch 1 taken 126 times.
|
417 | if (result.event == EventType::invasion) |
| 252 | { | ||
| 253 | 873 | std::cout << "Throat " << eIdx << " was invaded from pore " << vIdx << " :"; | |
| 254 | 582 | std::cout << " pc: " << *pcMax; | |
| 255 | 582 | std::cout << ", pcEntry: " << spatialParams.pcEntry(element, elemVolVars); | |
| 256 | 873 | std::cout << ", sw: " << sw[result.localScvIdxWithCriticalPc] << std::endl; | |
| 257 | } | ||
| 258 | else | ||
| 259 | { | ||
| 260 | 378 | std::cout << "Snap-off occurred at throat " << eIdx << " from pore " << vIdx << " :"; | |
| 261 | 252 | std::cout << " pc: " << *pcMax; | |
| 262 | 252 | std::cout << ", pcSnapoff: " << spatialParams.pcSnapoff(element, elemVolVars); | |
| 263 | 378 | std::cout << ", sw: " << sw[result.localScvIdxWithCriticalPc] << std::endl; | |
| 264 | } | ||
| 265 | } | ||
| 266 | |||
| 267 | 417 | return result; | |
| 268 | } | ||
| 269 | } | ||
| 270 | |||
| 271 | //! If the user has specified a global capillary pressure, check if it is lower than the given entry capillary pressure. | ||
| 272 | //! This may be needed to exactly reproduce pc-S curves given by static network models. | ||
| 273 | template<class Scalar> | ||
| 274 | ✗ | bool maybeRestrictToGlobalCapillaryPressure_(const Scalar pcEntry) const | |
| 275 | { | ||
| 276 | if constexpr (hasGlobalCapillaryPressure<Problem>()) | ||
| 277 | return restrictToGlobalCapillaryPressure_ && (pcEntry > problem_.globalCapillaryPressure()); | ||
| 278 | else | ||
| 279 | ✗ | return false; | |
| 280 | } | ||
| 281 | |||
| 282 | std::vector<bool> invadedCurrentIteration_; | ||
| 283 | std::vector<bool> invadedPreviousTimeStep_; | ||
| 284 | bool hasChangedInCurrentIteration_ = false; | ||
| 285 | bool hasChangedComparedToPreviousTimestep_ = false; | ||
| 286 | std::size_t numThroatsInvaded_; | ||
| 287 | bool verbose_; | ||
| 288 | bool restrictToGlobalCapillaryPressure_; | ||
| 289 | |||
| 290 | const Problem& problem_; | ||
| 291 | }; | ||
| 292 | |||
| 293 | } // end namespace Dumux::PoreNetwork | ||
| 294 | |||
| 295 | #endif | ||
| 296 |