| 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 Flux variables cache for the two-phase-flow PNM | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_PNM_2P_FLUXVARIABLESCACHE_HH | ||
| 13 | #define DUMUX_PNM_2P_FLUXVARIABLESCACHE_HH | ||
| 14 | |||
| 15 | #include <array> | ||
| 16 | #include <algorithm> | ||
| 17 | #include <dune/common/reservedvector.hh> | ||
| 18 | #include <dumux/porenetwork/common/throatproperties.hh> | ||
| 19 | |||
| 20 | namespace Dumux::PoreNetwork { | ||
| 21 | |||
| 22 | /*! | ||
| 23 | * \ingroup PNMTwoPModel | ||
| 24 | * \brief Flux variables cache for the two-phase-flow PNM | ||
| 25 | * Store data required for flux calculation. | ||
| 26 | */ | ||
| 27 | template<class AdvectionType, int maxNumCorners = 4> | ||
| 28 |
1/4✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 6 taken 31352 times.
✗ Branch 7 not taken.
|
474801 | class TwoPFluxVariablesCache |
| 29 | { | ||
| 30 | using Scalar = typename AdvectionType::Scalar; | ||
| 31 | static constexpr auto numPhases = 2; | ||
| 32 | using NumCornerVector = Dune::ReservedVector<Scalar, maxNumCorners>; | ||
| 33 | |||
| 34 | public: | ||
| 35 | //! whether the cache needs an update when the solution changes | ||
| 36 | static bool constexpr isSolDependent = true; | ||
| 37 | |||
| 38 | template<class Problem, class Element, class FVElementGeometry, | ||
| 39 | class ElementVolumeVariables, class SubControlVolumeFace> | ||
| 40 | 959820 | void update(const Problem& problem, | |
| 41 | const Element& element, | ||
| 42 | const FVElementGeometry& fvGeometry, | ||
| 43 | const ElementVolumeVariables& elemVolVars, | ||
| 44 | const SubControlVolumeFace& scvf, | ||
| 45 | bool invaded) | ||
| 46 | { | ||
| 47 | 1919640 | const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(element); | |
| 48 |
1/2✓ Branch 0 taken 959820 times.
✗ Branch 1 not taken.
|
959820 | throatCrossSectionShape_ = fvGeometry.gridGeometry().throatCrossSectionShape(eIdx); |
| 49 |
1/2✓ Branch 0 taken 959820 times.
✗ Branch 1 not taken.
|
959820 | throatShapeFactor_ = fvGeometry.gridGeometry().throatShapeFactor(eIdx); |
| 50 |
6/6✓ Branch 0 taken 253285 times.
✓ Branch 1 taken 706535 times.
✓ Branch 2 taken 253285 times.
✓ Branch 3 taken 706535 times.
✓ Branch 4 taken 253285 times.
✓ Branch 5 taken 706535 times.
|
2879460 | pc_ = std::max(elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure()); |
| 51 | 1919640 | pcEntry_ = problem.spatialParams().pcEntry(element, elemVolVars); | |
| 52 | 1919640 | pcSnapoff_ = problem.spatialParams().pcSnapoff(element, elemVolVars); | |
| 53 | 1919640 | throatInscribedRadius_ = problem.spatialParams().throatInscribedRadius(element, elemVolVars); | |
| 54 | 1919640 | throatLength_ = problem.spatialParams().throatLength(element, elemVolVars); | |
| 55 | 959820 | invaded_ = invaded; | |
| 56 | 959820 | poreToPoreDistance_ = element.geometry().volume(); | |
| 57 | |||
| 58 | // get the non-wetting phase index | ||
| 59 | using FluidSystem = typename ElementVolumeVariables::VolumeVariables::FluidSystem; | ||
| 60 | 959820 | const auto& spatialParams = problem.spatialParams(); | |
| 61 | 959820 | nPhaseIdx_ = 1 - spatialParams.template wettingPhase<FluidSystem>(element, elemVolVars); | |
| 62 | |||
| 63 | // take the average surface tension of both adjacent pores TODO: is this correct? | ||
| 64 | 1919640 | surfaceTension_ = 0.5*(elemVolVars[0].surfaceTension() + elemVolVars[1].surfaceTension()); | |
| 65 | |||
| 66 | 959820 | const auto& cornerHalfAngles = spatialParams.cornerHalfAngles(element); | |
| 67 | 1919640 | wettingLayerArea_.clear(); wettingLayerArea_.resize(cornerHalfAngles.size()); | |
| 68 | 959820 | const Scalar totalThroatCrossSectionalArea = spatialParams.throatCrossSectionalArea(element, elemVolVars); | |
| 69 | |||
| 70 |
2/2✓ Branch 0 taken 92457 times.
✓ Branch 1 taken 867363 times.
|
959820 | if (invaded) // two-phase flow |
| 71 | { | ||
| 72 | 92457 | const Scalar theta = spatialParams.contactAngle(element, elemVolVars); | |
| 73 |
2/2✓ Branch 0 taken 369828 times.
✓ Branch 1 taken 92457 times.
|
462285 | for (int i = 0; i< cornerHalfAngles.size(); ++i) |
| 74 | 739656 | wettingLayerArea_[i] = Throat::wettingLayerCrossSectionalArea(curvatureRadius(), theta, cornerHalfAngles[i]); | |
| 75 | |||
| 76 | // make sure the wetting phase area does not exceed the total cross-section area | ||
| 77 | 369828 | throatCrossSectionalArea_[wPhaseIdx()] = std::min( | |
| 78 |
2/2✓ Branch 0 taken 57 times.
✓ Branch 1 taken 92400 times.
|
369828 | std::accumulate(wettingLayerArea_.begin(), wettingLayerArea_.end(), 0.0), |
| 79 | totalThroatCrossSectionalArea | ||
| 80 | ); | ||
| 81 | 277371 | throatCrossSectionalArea_[nPhaseIdx()] = totalThroatCrossSectionalArea - throatCrossSectionalArea_[wPhaseIdx()]; | |
| 82 | } | ||
| 83 | else // single-phase flow | ||
| 84 | { | ||
| 85 |
2/2✓ Branch 0 taken 3469452 times.
✓ Branch 1 taken 867363 times.
|
4336815 | for (int i = 0; i< cornerHalfAngles.size(); ++i) |
| 86 | 6938904 | wettingLayerArea_[i] = 0.0; | |
| 87 | |||
| 88 | 867363 | throatCrossSectionalArea_[wPhaseIdx()] = totalThroatCrossSectionalArea; | |
| 89 | 1734726 | throatCrossSectionalArea_[nPhaseIdx()] = 0.0; | |
| 90 | } | ||
| 91 | |||
| 92 |
2/2✓ Branch 0 taken 1919640 times.
✓ Branch 1 taken 959820 times.
|
2879460 | for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) |
| 93 | { | ||
| 94 | 1919640 | singlePhaseCache_.fill(problem, element, fvGeometry, scvf, elemVolVars, *this, phaseIdx); | |
| 95 | 1919640 | nonWettingPhaseCache_.fill(problem, element, fvGeometry, scvf, elemVolVars, *this, phaseIdx); | |
| 96 | 1919640 | wettingLayerCache_.fill(problem, element, fvGeometry, scvf, elemVolVars, *this, phaseIdx); | |
| 97 | } | ||
| 98 | |||
| 99 |
2/2✓ Branch 0 taken 1919640 times.
✓ Branch 1 taken 959820 times.
|
2879460 | for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) |
| 100 | { | ||
| 101 | 1919640 | transmissibility_[phaseIdx] = AdvectionType::calculateTransmissibility( | |
| 102 | problem, element, fvGeometry, scvf, elemVolVars, *this, phaseIdx | ||
| 103 | ); | ||
| 104 | } | ||
| 105 | 959820 | } | |
| 106 | |||
| 107 | /*! | ||
| 108 | * \brief Returns the throats's cross-sectional shape. | ||
| 109 | */ | ||
| 110 | ✗ | Throat::Shape throatCrossSectionShape() const | |
| 111 | ✗ | { return throatCrossSectionShape_; } | |
| 112 | |||
| 113 | /*! | ||
| 114 | * \brief Returns the throats's shape factor. | ||
| 115 | */ | ||
| 116 | ✗ | Scalar throatShapeFactor() const | |
| 117 | ✗ | { return throatShapeFactor_; } | |
| 118 | |||
| 119 | /*! | ||
| 120 | * \brief Returns the throats's transmissibility. | ||
| 121 | */ | ||
| 122 | Scalar transmissibility(const int phaseIdx) const | ||
| 123 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1641413 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1641413 times.
|
3282826 | { return transmissibility_[phaseIdx]; } |
| 124 | |||
| 125 | /*! | ||
| 126 | * \brief Returns the throats's cross-sectional area for a given phaseIdx. | ||
| 127 | */ | ||
| 128 | Scalar throatCrossSectionalArea(const int phaseIdx) const | ||
| 129 |
4/4✓ Branch 0 taken 33903 times.
✓ Branch 1 taken 28801 times.
✓ Branch 2 taken 33903 times.
✓ Branch 3 taken 28801 times.
|
2255904 | { return throatCrossSectionalArea_[phaseIdx]; } |
| 130 | |||
| 131 | /*! | ||
| 132 | * \brief Returns the throats's total cross-sectional area. | ||
| 133 | */ | ||
| 134 | Scalar throatCrossSectionalArea() const | ||
| 135 | ✗ | { return throatCrossSectionalArea_[0] + throatCrossSectionalArea_[1]; } | |
| 136 | |||
| 137 | /*! | ||
| 138 | * \brief Returns the throats's length. | ||
| 139 | */ | ||
| 140 | ✗ | Scalar throatLength() const | |
| 141 | ✗ | { return throatLength_; } | |
| 142 | |||
| 143 | /*! | ||
| 144 | * \brief Returns the throats's inscribed radius. | ||
| 145 | */ | ||
| 146 | ✗ | Scalar throatInscribedRadius() const | |
| 147 | ✗ | { return throatInscribedRadius_; } | |
| 148 | |||
| 149 | /*! | ||
| 150 | * \brief Returns the throats's entry capillary pressure. | ||
| 151 | */ | ||
| 152 | ✗ | Scalar pcEntry() const | |
| 153 | ✗ | { return pcEntry_; } | |
| 154 | |||
| 155 | /*! | ||
| 156 | * \brief Returns the throats's snap-off capillary pressure. | ||
| 157 | */ | ||
| 158 | ✗ | Scalar pcSnapoff() const | |
| 159 | ✗ | { return pcSnapoff_; } | |
| 160 | |||
| 161 | /*! | ||
| 162 | * \brief Returns the capillary pressure within the throat. | ||
| 163 | */ | ||
| 164 | Scalar pc() const | ||
| 165 | { return pc_; } | ||
| 166 | |||
| 167 | /*! | ||
| 168 | * \brief Returns the surface tension within the throat. | ||
| 169 | */ | ||
| 170 | Scalar surfaceTension() const | ||
| 171 | { return surfaceTension_; } | ||
| 172 | |||
| 173 | /*! | ||
| 174 | * \brief Returns true if the throat is invaded by the nonwetting phase. | ||
| 175 | */ | ||
| 176 | ✗ | bool invaded() const | |
| 177 | ✗ | { return invaded_; } | |
| 178 | |||
| 179 | /*! | ||
| 180 | * \brief Returns the curvature radius within the throat. | ||
| 181 | */ | ||
| 182 | ✗ | Scalar curvatureRadius() const | |
| 183 | 462285 | { return surfaceTension_ / pc_;} | |
| 184 | |||
| 185 | /*! | ||
| 186 | * \brief Returns the cross-sectional area of a wetting layer within | ||
| 187 | * one of the throat's corners. | ||
| 188 | */ | ||
| 189 | Scalar wettingLayerCrossSectionalArea(const int cornerIdx) const | ||
| 190 | 739656 | { return wettingLayerArea_[cornerIdx]; } | |
| 191 | |||
| 192 | /*! | ||
| 193 | * \brief Returns the index of the wetting phase. | ||
| 194 | */ | ||
| 195 | ✗ | std::size_t wPhaseIdx() const | |
| 196 | 92457 | { return 1 - nPhaseIdx_; } | |
| 197 | |||
| 198 | /*! | ||
| 199 | * \brief Returns the index of the nonwetting phase. | ||
| 200 | */ | ||
| 201 | ✗ | std::size_t nPhaseIdx() const | |
| 202 | ✗ | { return nPhaseIdx_; } | |
| 203 | |||
| 204 | /*! | ||
| 205 | * \brief Returns the throats's cached flow variables for single-phase flow. | ||
| 206 | */ | ||
| 207 | const auto& singlePhaseFlowVariables() const | ||
| 208 | { return singlePhaseCache_; } | ||
| 209 | |||
| 210 | /*! | ||
| 211 | * \brief Returns the throats's cached flow variables for the nonwetting phase. | ||
| 212 | */ | ||
| 213 | const auto& nonWettingPhaseFlowVariables() const | ||
| 214 | { return nonWettingPhaseCache_; } | ||
| 215 | |||
| 216 | /*! | ||
| 217 | * \brief Returns the throats's cached flow variables for the wetting phase. | ||
| 218 | */ | ||
| 219 | const auto& wettingLayerFlowVariables() const | ||
| 220 | 369828 | { return wettingLayerCache_; } | |
| 221 | |||
| 222 | /*! | ||
| 223 | * \brief Returns the throats's pore-to-pore-center distance. | ||
| 224 | */ | ||
| 225 | ✗ | Scalar poreToPoreDistance() const | |
| 226 | ✗ | { return poreToPoreDistance_; } | |
| 227 | |||
| 228 | private: | ||
| 229 | Throat::Shape throatCrossSectionShape_; | ||
| 230 | Scalar throatShapeFactor_; | ||
| 231 | std::array<Scalar, numPhases> transmissibility_; | ||
| 232 | std::array<Scalar, numPhases> throatCrossSectionalArea_; | ||
| 233 | Scalar throatLength_; | ||
| 234 | Scalar throatInscribedRadius_; | ||
| 235 | Scalar pcEntry_; | ||
| 236 | Scalar pcSnapoff_; | ||
| 237 | Scalar pc_; | ||
| 238 | Scalar surfaceTension_; | ||
| 239 | bool invaded_; | ||
| 240 | NumCornerVector wettingLayerArea_; | ||
| 241 | std::size_t nPhaseIdx_; | ||
| 242 | Scalar poreToPoreDistance_; | ||
| 243 | |||
| 244 | typename AdvectionType::Transmissibility::SinglePhaseCache singlePhaseCache_; | ||
| 245 | typename AdvectionType::Transmissibility::NonWettingPhaseCache nonWettingPhaseCache_; | ||
| 246 | typename AdvectionType::Transmissibility::WettingLayerCache wettingLayerCache_; | ||
| 247 | }; | ||
| 248 | |||
| 249 | } // end Dumux::PoreNetwork | ||
| 250 | |||
| 251 | #endif | ||
| 252 |