| 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 Fluidmatrixinteractions | ||
| 10 | * \ingroup PoreNetworkModels | ||
| 11 | * \brief Pore-local pc-Sw curves for for platonic bodies | ||
| 12 | * (tetrahedron, cube, octahedron, dodecahedron, icosahedron). | ||
| 13 | */ | ||
| 14 | #ifndef DUMUX_PNM_2P_LOCAL_RULES_FOR_PLATONIC_BODY_HH | ||
| 15 | #define DUMUX_PNM_2P_LOCAL_RULES_FOR_PLATONIC_BODY_HH | ||
| 16 | |||
| 17 | #include <algorithm> | ||
| 18 | #include <cassert> | ||
| 19 | #include <cmath> | ||
| 20 | #include <optional> | ||
| 21 | #include <dune/common/math.hh> | ||
| 22 | #include <dumux/common/optionalscalar.hh> | ||
| 23 | #include <dumux/common/spline.hh> | ||
| 24 | #include <dumux/common/typetraits/typetraits.hh> | ||
| 25 | #include <dumux/porenetwork/common/poreproperties.hh> | ||
| 26 | #include "singleshapelocalrules.hh" | ||
| 27 | |||
| 28 | namespace Dumux::PoreNetwork::FluidMatrix { | ||
| 29 | |||
| 30 | /*! | ||
| 31 | * \brief The parameter type | ||
| 32 | * \tparam Scalar The scalar type | ||
| 33 | */ | ||
| 34 | template<class Scalar> | ||
| 35 | struct PlatonicBodyParams | ||
| 36 | { | ||
| 37 | PlatonicBodyParams() = default; | ||
| 38 | |||
| 39 | PlatonicBodyParams& setPoreInscribedRadius(Scalar r) { radius_ = r; return *this;} | ||
| 40 | PlatonicBodyParams& setPoreShape(Pore::Shape s) { shape_ = s; return *this;} | ||
| 41 | PlatonicBodyParams& setSurfaceTension(Scalar st) { surfaceTension_ = st; return *this;} | ||
| 42 | |||
| 43 | template<class SpatialParams, class Element, class SubControlVolume, class ElemSol> | ||
| 44 | 6042272 | PlatonicBodyParams(const SpatialParams& spatialParams, | |
| 45 | const Element& element, | ||
| 46 | const SubControlVolume& scv, | ||
| 47 | const ElemSol& elemSol) | ||
| 48 | 12084544 | : shape_(spatialParams.gridGeometry().poreGeometry()[scv.dofIndex()]) | |
| 49 | , radius_(spatialParams.poreInscribedRadius(element, scv, elemSol)) | ||
| 50 | 18126816 | , surfaceTension_(spatialParams.surfaceTension(element, scv, elemSol)) | |
| 51 | 6042272 | {} | |
| 52 | |||
| 53 | template<class SpatialParams, class Element, class SubControlVolume, class ElemSol> | ||
| 54 | void update(const SpatialParams& spatialParams, | ||
| 55 | const Element& element, | ||
| 56 | const SubControlVolume& scv, | ||
| 57 | const ElemSol& elemSol) | ||
| 58 | { | ||
| 59 | const auto& gridGeometry = spatialParams.gridGeometry(); | ||
| 60 | shape_ = gridGeometry.poreGeometry()[scv.dofIndex()]; | ||
| 61 | radius_ = spatialParams.poreInscribedRadius(element, scv, elemSol); | ||
| 62 | surfaceTension_ = spatialParams.surfaceTension(element, scv, elemSol); | ||
| 63 | } | ||
| 64 | |||
| 65 | ✗ | Pore::Shape poreShape() const { return shape_; } | |
| 66 | |||
| 67 | ✗ | Scalar poreInscribedRadius() const { return radius_; } | |
| 68 | |||
| 69 | ✗ | Scalar surfaceTension() const { return surfaceTension_; } | |
| 70 | |||
| 71 | bool operator== (const PlatonicBodyParams& p) const | ||
| 72 | { | ||
| 73 | return Dune::FloatCmp::eq(radius_, p.radius_, 1e-6) | ||
| 74 | && Dune::FloatCmp::eq(surfaceTension_, p.surfaceTension_, 1e-6) | ||
| 75 | && shape_ == p.shape_; | ||
| 76 | } | ||
| 77 | |||
| 78 | private: | ||
| 79 | Pore::Shape shape_; | ||
| 80 | Scalar radius_; | ||
| 81 | Scalar surfaceTension_; | ||
| 82 | }; | ||
| 83 | |||
| 84 | /*! | ||
| 85 | * \ingroup Fluidmatrixinteractions | ||
| 86 | * \ingroup PoreNetworkModels | ||
| 87 | * \brief Implementation of the simplified pore-local capillary pressure-saturation curve | ||
| 88 | * for platonic bodies (tetrahedron, cube, octahedron, dodecahedron, icosahedron). | ||
| 89 | * | ||
| 90 | * See Joekar-Niasar et al., 2010 and Sweijen et al., 2018. | ||
| 91 | */ | ||
| 92 | template<Pore::Shape shape> | ||
| 93 | struct TwoPLocalRulesPlatonicBody | ||
| 94 | { | ||
| 95 | template<class Scalar> | ||
| 96 | using Params = PlatonicBodyParams<Scalar>; | ||
| 97 | |||
| 98 | template<class SpatialParams, class Element, class SubControlVolume, class ElemSol> | ||
| 99 | static auto makeParams(const SpatialParams& spatialParams, | ||
| 100 | const Element& element, | ||
| 101 | const SubControlVolume& scv, | ||
| 102 | const ElemSol& elemSol) | ||
| 103 | { | ||
| 104 | using Scalar = std::decay_t<decltype(spatialParams.poreInscribedRadius(element, scv, elemSol))>; | ||
| 105 | return Params<Scalar>(spatialParams, element, scv, elemSol); | ||
| 106 | } | ||
| 107 | |||
| 108 | /*! | ||
| 109 | * \brief The capillary pressure-saturation curve | ||
| 110 | * | ||
| 111 | * \param sw Saturation of the wetting phase \f$\mathrm{S_{w,i}}\f$ at pore \f$i\f$ | ||
| 112 | * \param params The parameters container | ||
| 113 | */ | ||
| 114 | template<class Scalar> | ||
| 115 | 3156562 | static Scalar pc(Scalar sw, const Params<Scalar>& params) | |
| 116 | { | ||
| 117 | ✗ | assert(isPlatonicBody_(params.poreShape())); | |
| 118 | |||
| 119 | using std::clamp; | ||
| 120 |
1/2✓ Branch 0 taken 1619844 times.
✗ Branch 1 not taken.
|
3156562 | sw = clamp(sw, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0 |
| 121 | |||
| 122 | 3156562 | const Scalar poreRadius = params.poreInscribedRadius(); | |
| 123 | 3156562 | const Scalar sigma = params.surfaceTension(); | |
| 124 | // TODO incorporate contact angle!!! | ||
| 125 | |||
| 126 | using std::exp; | ||
| 127 | 3156562 | return 2.0*sigma / (poreRadius*(1.0 - exp(-expFactor_<Scalar>()*sw))); | |
| 128 | } | ||
| 129 | |||
| 130 | /*! | ||
| 131 | * \brief The wetting-phase saturation of a pore body | ||
| 132 | * | ||
| 133 | * \param pc The capillary pressure \f$\mathrm{p_{c,i}}\f$ at pore \f$i\f$ | ||
| 134 | * \param params The parameters container | ||
| 135 | */ | ||
| 136 | template<class Scalar> | ||
| 137 | 43131 | static Scalar sw(Scalar pc, const Params<Scalar>& params) | |
| 138 | { | ||
| 139 | ✗ | assert(isPlatonicBody_(params.poreShape())); | |
| 140 | |||
| 141 | using std::max; | ||
| 142 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42347 times.
|
43131 | pc = max(pc, 0.0); // the equation below is undefined for negative pcs |
| 143 | |||
| 144 | 43131 | const Scalar poreRadius = params.poreInscribedRadius(); | |
| 145 | 43131 | const Scalar sigma = params.surfaceTension(); | |
| 146 | |||
| 147 | using std::log; | ||
| 148 | 43131 | return -1.0/expFactor_<Scalar>()* log(1.0 - 2.0*sigma/(poreRadius*pc)); | |
| 149 | } | ||
| 150 | |||
| 151 | /*! | ||
| 152 | * \brief The partial derivative of the capillary | ||
| 153 | * pressure w.r.t. the wetting phase saturation. | ||
| 154 | * | ||
| 155 | * \param sw Saturation of the wetting phase \f$\mathrm{S_{w,i}}\f$ at pore \f$i\f$ | ||
| 156 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 157 | */ | ||
| 158 | template<class Scalar> | ||
| 159 | 38228 | static Scalar dpc_dsw(Scalar sw, const Params<Scalar>& params) | |
| 160 | { | ||
| 161 | ✗ | assert(isPlatonicBody_(params.poreShape())); | |
| 162 | |||
| 163 | using std::clamp; | ||
| 164 |
1/2✓ Branch 0 taken 19114 times.
✗ Branch 1 not taken.
|
38228 | sw = clamp(sw, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0 |
| 165 | |||
| 166 | 38228 | const Scalar sigma = params.surfaceTension(); | |
| 167 | 38228 | const Scalar poreRadius = params.poreInscribedRadius(); | |
| 168 | using std::exp; | ||
| 169 | 38228 | const Scalar e = exp(expFactor_<Scalar>()*sw); | |
| 170 | 38228 | return -(2.0*expFactor_<Scalar>()*sigma*e) / (poreRadius*(1.0-e)*(1.0-e)); | |
| 171 | } | ||
| 172 | |||
| 173 | /*! | ||
| 174 | * \brief The partial derivative of the wetting phase saturation | ||
| 175 | * w.r.t. the capillary pressure. | ||
| 176 | * | ||
| 177 | * \param pc The capillary pressure \f$\mathrm{p_{c,i}}\f$ at pore \f$i\f$ | ||
| 178 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 179 | */ | ||
| 180 | template<class Scalar> | ||
| 181 | 1568 | static Scalar dsw_dpc(Scalar pc, const Params<Scalar>& params) | |
| 182 | { | ||
| 183 | ✗ | assert(isPlatonicBody_(params.poreShape())); | |
| 184 | |||
| 185 | using std::max; | ||
| 186 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 784 times.
|
1568 | pc = max(pc, 0.0); // the equation below is undefined for negative pcs |
| 187 | |||
| 188 | 1568 | const Scalar sigma = params.surfaceTension(); | |
| 189 | 1568 | const Scalar poreRadius = params.poreInscribedRadius(); | |
| 190 | 1568 | return sigma / (expFactor_<Scalar>()*sigma*pc - 0.5*expFactor_<Scalar>()*poreRadius*pc*pc); | |
| 191 | } | ||
| 192 | |||
| 193 | private: | ||
| 194 | |||
| 195 | template<class Scalar> | ||
| 196 | static constexpr Scalar expFactor_() | ||
| 197 | { | ||
| 198 | if constexpr (shape == Pore::Shape::tetrahedron) | ||
| 199 | return 3.87; | ||
| 200 | else if constexpr (shape == Pore::Shape::cube) | ||
| 201 | return 6.83; | ||
| 202 | else if constexpr (shape == Pore::Shape::octahedron) | ||
| 203 | return 8.71; | ||
| 204 | else if constexpr (shape == Pore::Shape::dodecahedron) | ||
| 205 | return 22.87; | ||
| 206 | else if constexpr (shape == Pore::Shape::icosahedron) | ||
| 207 | return 24.11; | ||
| 208 | else | ||
| 209 | { | ||
| 210 | static_assert(AlwaysFalse<Scalar>::value, "Shape not supported"); | ||
| 211 | return 0; | ||
| 212 | } | ||
| 213 | } | ||
| 214 | |||
| 215 | static constexpr bool isPlatonicBody_(Pore::Shape s) | ||
| 216 | { | ||
| 217 | 3364178 | return s == Pore::Shape::tetrahedron || | |
| 218 | 1682089 | s == Pore::Shape::cube || | |
| 219 | 5898 | s == Pore::Shape::octahedron || | |
| 220 |
40/120✓ Branch 0 taken 98 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 98 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 98 times.
✓ Branch 7 taken 83126 times.
✓ Branch 8 taken 98 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 98 times.
✓ Branch 12 taken 98 times.
✓ Branch 13 taken 41563 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 98 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 16819 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 98 times.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 98 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 30 taken 98 times.
✗ Branch 31 not taken.
✓ Branch 32 taken 98 times.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✓ Branch 35 taken 98 times.
✓ Branch 36 taken 98 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 98 times.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✓ Branch 43 taken 98 times.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✓ Branch 48 taken 197 times.
✓ Branch 49 taken 42477 times.
✗ Branch 50 not taken.
✓ Branch 51 taken 197 times.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✓ Branch 54 taken 197 times.
✗ Branch 55 not taken.
✓ Branch 56 taken 197 times.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✓ Branch 59 taken 197 times.
✓ Branch 60 taken 197 times.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✓ Branch 63 taken 197 times.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✓ Branch 67 taken 197 times.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✓ Branch 72 taken 590 times.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✓ Branch 75 taken 590 times.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✓ Branch 78 taken 590 times.
✓ Branch 79 taken 1489738 times.
✓ Branch 80 taken 590 times.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✓ Branch 83 taken 590 times.
✓ Branch 84 taken 590 times.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✓ Branch 87 taken 590 times.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✓ Branch 91 taken 590 times.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✓ Branch 97 taken 392 times.
✗ Branch 98 not taken.
✗ Branch 99 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 102 not taken.
✓ Branch 103 taken 392 times.
✗ Branch 104 not taken.
✗ Branch 105 not taken.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ Branch 108 not taken.
✓ Branch 109 taken 788 times.
✗ Branch 110 not taken.
✗ Branch 111 not taken.
✗ Branch 112 not taken.
✗ Branch 113 not taken.
✗ Branch 114 not taken.
✓ Branch 115 taken 2960 times.
✗ Branch 116 not taken.
✗ Branch 117 not taken.
✗ Branch 118 not taken.
✗ Branch 119 not taken.
|
1682089 | s == Pore::Shape::dodecahedron || |
| 221 | s == Pore::Shape::icosahedron; | ||
| 222 | } | ||
| 223 | }; | ||
| 224 | |||
| 225 | /*! | ||
| 226 | * \ingroup Fluidmatrixinteractions | ||
| 227 | * \ingroup PoreNetworkModels | ||
| 228 | * \brief Two-phase rules for regularizing the pc-SW for platonic bodies. | ||
| 229 | */ | ||
| 230 | template<class Scalar, class BaseLaw> | ||
| 231 |
35/110✓ Branch 1 taken 41563 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 41563 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 41563 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 41563 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 41563 times.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✓ Branch 101 taken 1 times.
✗ Branch 102 not taken.
✓ Branch 104 taken 1 times.
✗ Branch 105 not taken.
✓ Branch 107 taken 1 times.
✗ Branch 108 not taken.
✓ Branch 110 taken 1 times.
✗ Branch 111 not taken.
✓ Branch 113 taken 1 times.
✗ Branch 114 not taken.
✓ Branch 116 taken 1 times.
✗ Branch 117 not taken.
✓ Branch 119 taken 1 times.
✗ Branch 120 not taken.
✓ Branch 122 taken 1 times.
✗ Branch 123 not taken.
✓ Branch 125 taken 1 times.
✗ Branch 126 not taken.
✓ Branch 128 taken 1 times.
✗ Branch 129 not taken.
✓ Branch 131 taken 1 times.
✗ Branch 132 not taken.
✓ Branch 134 taken 1 times.
✗ Branch 135 not taken.
✓ Branch 137 taken 1 times.
✗ Branch 138 not taken.
✓ Branch 140 taken 1 times.
✗ Branch 141 not taken.
✓ Branch 143 taken 1 times.
✗ Branch 144 not taken.
✓ Branch 146 taken 3062797 times.
✗ Branch 147 not taken.
✓ Branch 149 taken 3062797 times.
✗ Branch 150 not taken.
✓ Branch 152 taken 3062797 times.
✗ Branch 153 not taken.
✓ Branch 155 taken 3062797 times.
✗ Branch 156 not taken.
✓ Branch 158 taken 3062797 times.
✗ Branch 159 not taken.
✓ Branch 161 taken 1 times.
✗ Branch 162 not taken.
✓ Branch 164 taken 1 times.
✗ Branch 165 not taken.
✓ Branch 167 taken 1 times.
✗ Branch 168 not taken.
✓ Branch 170 taken 1 times.
✗ Branch 171 not taken.
✓ Branch 173 taken 1 times.
✗ Branch 174 not taken.
✓ Branch 176 taken 3 times.
✗ Branch 177 not taken.
✓ Branch 179 taken 3 times.
✗ Branch 180 not taken.
✓ Branch 182 taken 3 times.
✗ Branch 183 not taken.
✓ Branch 185 taken 3 times.
✗ Branch 186 not taken.
✓ Branch 188 taken 3 times.
✗ Branch 189 not taken.
|
15521835 | class TwoPLocalRulesPlatonicBodyRegularization |
| 232 | { | ||
| 233 | using ThisType = TwoPLocalRulesPlatonicBodyRegularization<Scalar, BaseLaw>; | ||
| 234 | using BaseLawParams = typename BaseLaw::template Params<Scalar>; | ||
| 235 | public: | ||
| 236 | |||
| 237 | /*! | ||
| 238 | * \brief The available options for regularizing the pc-SW curve at high wetting-phase saturations. | ||
| 239 | */ | ||
| 240 | enum class HighSwRegularizationMethod | ||
| 241 | { | ||
| 242 | linear, spline, powerLaw | ||
| 243 | }; | ||
| 244 | |||
| 245 | template<class S> | ||
| 246 |
5/10✓ Branch 1 taken 41564 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3062797 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
|
3104364 | struct Params |
| 247 | { | ||
| 248 | // export method type | ||
| 249 | using HighSwRegularizationMethod = ThisType::HighSwRegularizationMethod; | ||
| 250 | /*! | ||
| 251 | * \brief Set the threshold saturation below which the capillary pressure is regularized. | ||
| 252 | * | ||
| 253 | * Most problems are very sensitive to this value (e.g. making it smaller might | ||
| 254 | * result in very high capillary pressures) | ||
| 255 | */ | ||
| 256 | void setpcLowSw(S pcLowSw) | ||
| 257 | { pcLowSw_ = pcLowSw; } | ||
| 258 | |||
| 259 | /*! | ||
| 260 | * \brief Threshold saturation below which the capillary pressure is regularized. | ||
| 261 | */ | ||
| 262 | ✗ | S pcLowSw() const | |
| 263 | ✗ | { return pcLowSw_; } | |
| 264 | |||
| 265 | /*! | ||
| 266 | * \brief Set the threshold saturation above which the capillary pressure is regularized. | ||
| 267 | */ | ||
| 268 | void setpcHighSw(S pcHighSw) | ||
| 269 | { pcHighSw_ = pcHighSw; } | ||
| 270 | |||
| 271 | /*! | ||
| 272 | * \brief Threshold saturation above which the capillary pressure is regularized. | ||
| 273 | * | ||
| 274 | * Most problems are very sensitive to this value (e.g. making it smaller might | ||
| 275 | * result in negative capillary pressures). | ||
| 276 | */ | ||
| 277 | ✗ | S pcHighSw() const | |
| 278 | ✗ | { return pcHighSw_; } | |
| 279 | |||
| 280 | /*! | ||
| 281 | * \brief Set the regularization method for high saturations. | ||
| 282 | */ | ||
| 283 | ✗ | void setHighSwRegularizationMethod(HighSwRegularizationMethod method) | |
| 284 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
|
3 | { highSwRegularizationMethod_ = method; } |
| 285 | |||
| 286 | /*! | ||
| 287 | * \brief Return the regularization method for high saturations. | ||
| 288 | */ | ||
| 289 | ✗ | HighSwRegularizationMethod highSwRegularizationMethod() const | |
| 290 | ✗ | { return highSwRegularizationMethod_; } | |
| 291 | |||
| 292 | private: | ||
| 293 | S pcLowSw_ = 0.01; | ||
| 294 | S pcHighSw_ = 0.99; | ||
| 295 | HighSwRegularizationMethod highSwRegularizationMethod_ = HighSwRegularizationMethod::linear; | ||
| 296 | }; | ||
| 297 | |||
| 298 | |||
| 299 | template<class MaterialLaw> | ||
| 300 | void init(const MaterialLaw* m, const Params<Scalar>& p, const std::string& paramGroup = "") | ||
| 301 | { | ||
| 302 |
7/22✓ Branch 1 taken 41563 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 1 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 1 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 3062797 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 1 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 3 times.
✗ Branch 37 not taken.
|
3104367 | initPcParameters_(m, p, paramGroup); |
| 303 | } | ||
| 304 | |||
| 305 | 3072996 | OptionalScalar<Scalar> pc(const Scalar sw) const | |
| 306 | { | ||
| 307 | // make sure that the capilary pressure observes a | ||
| 308 | // derivative != 0 for 'illegal' saturations. This is | ||
| 309 | // required for example by newton solvers (if the | ||
| 310 | // derivative is calculated numerically) in order to get the | ||
| 311 | // saturation moving to the right direction if it | ||
| 312 | // temporarily is in an 'illegal' range. | ||
| 313 |
2/2✓ Branch 0 taken 17610 times.
✓ Branch 1 taken 1518888 times.
|
3072996 | if (sw < pcLowSw_) |
| 314 | 105660 | return pcLowSwPcValue_() + pcDerivativeLowSw_() * (sw - pcLowSw_); | |
| 315 | |||
| 316 |
2/2✓ Branch 0 taken 270465 times.
✓ Branch 1 taken 1248423 times.
|
3037776 | if (sw <= pcHighSw_) |
| 317 | 540930 | return {}; // standard | |
| 318 |
2/2✓ Branch 0 taken 221681 times.
✓ Branch 1 taken 1026742 times.
|
2496846 | else if (sw < 1.0) // regularized part below sw = 1.0 |
| 319 | { | ||
| 320 | using std::pow; | ||
| 321 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 221681 times.
|
443362 | if (highSwRegularizationMethod_ == HighSwRegularizationMethod::powerLaw) |
| 322 | ✗ | return pcHighSwPcValue_() * pow(((1.0-sw)/(1.0-pcHighSw_)), 1.0/3.0); | |
| 323 | |||
| 324 |
2/2✓ Branch 0 taken 131445 times.
✓ Branch 1 taken 90236 times.
|
443362 | else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::linear) |
| 325 | { | ||
| 326 | 525780 | const Scalar slope = -pcHighSwPcValue_() / (1.0 - pcHighSw_); | |
| 327 | 525780 | return pcHighSwPcValue_() + (sw - pcHighSw_) * slope; | |
| 328 | } | ||
| 329 | |||
| 330 |
1/2✓ Branch 0 taken 90236 times.
✗ Branch 1 not taken.
|
180472 | else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::spline) |
| 331 | 180472 | return pcSpline_().eval(sw); | |
| 332 | |||
| 333 | else | ||
| 334 | ✗ | DUNE_THROW(Dune::NotImplemented, "Regularization not method not implemented"); | |
| 335 | } | ||
| 336 | else // regularized part above sw = 1.0 | ||
| 337 | 2053484 | return pcDerivativeHighSwEnd_()*(sw - 1.0); | |
| 338 | } | ||
| 339 | |||
| 340 | /*! | ||
| 341 | * \brief The regularized saturation-capillary pressure curve | ||
| 342 | */ | ||
| 343 | 43163 | OptionalScalar<Scalar> sw(const Scalar pc) const | |
| 344 | { | ||
| 345 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 42355 times.
|
43163 | if (pc <= 0.0) |
| 346 | { | ||
| 347 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
16 | if (pcHighSw_ >= 1.0) |
| 348 | ✗ | return 1.0; // no regularization for high sw | |
| 349 | else | ||
| 350 | 16 | return pc/pcDerivativeHighSwEnd_() + 1.0; | |
| 351 | } | ||
| 352 | |||
| 353 | // low saturation | ||
| 354 | 86294 | if (pc > pcLowSwPcValue_()) | |
| 355 | 48 | return (pc - pcLowSwPcValue_())/pcDerivativeLowSw_() + pcLowSw_; | |
| 356 | |||
| 357 | // high saturation | ||
| 358 | 86262 | else if (pc <= pcHighSwPcValue_()) | |
| 359 | { | ||
| 360 | ✗ | if (highSwRegularizationMethod_ == HighSwRegularizationMethod::powerLaw) | |
| 361 | { | ||
| 362 | // invert power law | ||
| 363 | using Dune::power; | ||
| 364 | ✗ | return power(pc/pcHighSwPcValue_(), 3) * (pcHighSw_ - 1.0) + 1.0; | |
| 365 | } | ||
| 366 | |||
| 367 | ✗ | else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::linear) | |
| 368 | ✗ | return pc/pcDerivativeHighSwEnd_() + 1.0; | |
| 369 | |||
| 370 | ✗ | else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::spline) | |
| 371 | // invert spline | ||
| 372 | ✗ | return pcSpline_().intersectInterval(pcHighSw_, 1.0, 0.0, 0.0, 0.0, pc); | |
| 373 | } | ||
| 374 | |||
| 375 | // no regularization | ||
| 376 | 43131 | return {}; | |
| 377 | } | ||
| 378 | |||
| 379 | /*! | ||
| 380 | * \brief The regularized partial derivative of the capillary pressure w.r.t. the saturation | ||
| 381 | */ | ||
| 382 | 3200 | OptionalScalar<Scalar> dpc_dsw(const Scalar sw) const | |
| 383 | { | ||
| 384 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 1576 times.
|
3200 | if (sw <= pcLowSw_) |
| 385 | 96 | return pcDerivativeLowSw_(); | |
| 386 | |||
| 387 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1568 times.
|
3152 | else if (sw >= 1.0) |
| 388 | 16 | return pcDerivativeHighSwEnd_(); | |
| 389 | |||
| 390 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1568 times.
|
3136 | else if (sw > pcHighSw_) |
| 391 | { | ||
| 392 | ✗ | if (highSwRegularizationMethod_ == HighSwRegularizationMethod::powerLaw) | |
| 393 | { | ||
| 394 | using std::pow; | ||
| 395 | ✗ | return pcHighSwPcValue_()/3.0 * 1.0 /((pcHighSw_-1.0) * pow((sw-1.0)/(pcHighSw_-1.0), 2.0/3.0)); | |
| 396 | } | ||
| 397 | |||
| 398 | ✗ | else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::linear) | |
| 399 | ✗ | return pcDerivativeHighSwEnd_(); | |
| 400 | |||
| 401 | else | ||
| 402 | ✗ | return pcSpline_().evalDerivative(sw); | |
| 403 | } | ||
| 404 | |||
| 405 | else | ||
| 406 | 3136 | return {}; // no regularization | |
| 407 | } | ||
| 408 | |||
| 409 | /*! | ||
| 410 | * \brief The regularized partial derivative of the saturation to the capillary pressure | ||
| 411 | */ | ||
| 412 | 1600 | OptionalScalar<Scalar> dsw_dpc(const Scalar pc) const | |
| 413 | { | ||
| 414 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 792 times.
|
1600 | if (pc <= 0.0) |
| 415 | { | ||
| 416 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
16 | if (pcHighSw_ >= 1.0) |
| 417 | ✗ | return 0.0; | |
| 418 | else | ||
| 419 | 16 | return 1.0/pcDerivativeHighSwEnd_(); | |
| 420 | } | ||
| 421 | |||
| 422 | // derivative of the inverse of the function is one over derivative of the function | ||
| 423 | 3168 | else if (pc <= pcHighSwPcValue_()) | |
| 424 | { | ||
| 425 | ✗ | if (highSwRegularizationMethod_ == HighSwRegularizationMethod::powerLaw) | |
| 426 | { | ||
| 427 | using Dune::power; | ||
| 428 | ✗ | return (3.0*pcHighSw_ - 3.0) * power(pc, 2) / power(pcHighSwPcValue_(), 3); | |
| 429 | } | ||
| 430 | |||
| 431 | ✗ | else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::linear) | |
| 432 | ✗ | return 1.0/pcDerivativeHighSwEnd_(); | |
| 433 | |||
| 434 | else | ||
| 435 | ✗ | return 1.0/pcSpline_().evalDerivative(pcSpline_().intersectInterval(pcHighSw_, 1.0, 0.0, 0.0, 0.0, pc)); | |
| 436 | } | ||
| 437 | |||
| 438 | 3168 | else if (pc >= pcLowSwPcValue_()) | |
| 439 | 32 | return 1.0/pcDerivativeLowSw_(); | |
| 440 | |||
| 441 | else | ||
| 442 | 1568 | return {}; // no regularization | |
| 443 | } | ||
| 444 | |||
| 445 | template<class SpatialParams, class Element, class SubControlVolume, class ElemSol> | ||
| 446 | void updateParams(const SpatialParams& spatialParams, | ||
| 447 | const Element& element, | ||
| 448 | const SubControlVolume& scv, | ||
| 449 | const ElemSol& elemSol) | ||
| 450 | {} | ||
| 451 | |||
| 452 | private: | ||
| 453 | |||
| 454 | template<class MaterialLaw> | ||
| 455 | 6167171 | void initPcParameters_(const MaterialLaw* m, const Params<Scalar>& params, const std::string& paramGroup) | |
| 456 | { | ||
| 457 | 6167171 | highSwRegularizationMethod_ = params.highSwRegularizationMethod(); | |
| 458 | |||
| 459 | // maybe overwrite method using input | ||
| 460 |
7/14✓ Branch 0 taken 10 times.
✓ Branch 1 taken 3104357 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 10 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 10 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 10 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 10 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
|
6167228 | static const bool hasMethod = hasParamInGroup(paramGroup, "HighSwRegularizationMethod"); |
| 461 |
2/2✓ Branch 0 taken 1182963 times.
✓ Branch 1 taken 1921404 times.
|
6167171 | if (hasMethod) |
| 462 | { | ||
| 463 | // set method only once | ||
| 464 |
3/4✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1182960 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
|
2324370 | static const HighSwRegularizationMethod inputmethod = [&] |
| 465 | { | ||
| 466 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
6 | static const auto input = getParamFromGroup<std::string>(paramGroup, "HighSwRegularizationMethod"); |
| 467 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
3 | if (input == "Linear") |
| 468 | return HighSwRegularizationMethod::linear; | ||
| 469 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
|
3 | else if (input == "Spline") |
| 470 | return HighSwRegularizationMethod::spline; | ||
| 471 | ✗ | else if (input == "PowerLaw") | |
| 472 | return HighSwRegularizationMethod::powerLaw; | ||
| 473 | else | ||
| 474 | ✗ | DUNE_THROW(Dune::InvalidStateException, input << " is not a valid regularization method"); | |
| 475 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
|
13 | }(); |
| 476 | |||
| 477 | 2324363 | highSwRegularizationMethod_ = inputmethod; | |
| 478 | } | ||
| 479 | |||
| 480 |
4/6✓ Branch 0 taken 10 times.
✓ Branch 1 taken 3104357 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 10 times.
✗ Branch 7 not taken.
|
6167171 | static const bool splineZeroSlope = getParamFromGroup<bool>(paramGroup, "HighSwSplineZeroSlope", true); |
| 481 | 6167171 | highSwSplineZeroSlope_ = splineZeroSlope; | |
| 482 | |||
| 483 | // print name only once | ||
| 484 |
3/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 3104357 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
|
6167199 | [[maybe_unused]] static const bool printName = [&] |
| 485 | { | ||
| 486 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
10 | std::string name; |
| 487 | 10 | if (highSwRegularizationMethod_ == HighSwRegularizationMethod::linear) | |
| 488 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
7 | name = "Linear"; |
| 489 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
3 | else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::powerLaw) |
| 490 | ✗ | name = "PowerLaw"; | |
| 491 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
3 | else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::spline) |
| 492 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
3 | name = "Spline"; |
| 493 | else | ||
| 494 | ✗ | DUNE_THROW(Dune::NotImplemented, "Regularization not method not implemented"); | |
| 495 | |||
| 496 |
2/4✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
30 | std::cout << "\n*****\nUsing " << name << " as regularization method for high Sw\n*****" << std::endl; |
| 497 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
10 | return true; |
| 498 |
2/4✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
|
48 | }(); |
| 499 | |||
| 500 | using std::isnan; | ||
| 501 |
4/6✓ Branch 0 taken 10 times.
✓ Branch 1 taken 3104357 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 10 times.
✗ Branch 7 not taken.
|
6167171 | static const auto pcLowSwInput = getParamFromGroup<Scalar>(paramGroup, "RegularizationLowSw", params.pcLowSw()); |
| 502 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3104367 times.
|
6167171 | pcLowSw_ = !isnan(pcLowSwInput) ? pcLowSwInput : params.pcLowSw(); |
| 503 | |||
| 504 |
4/6✓ Branch 0 taken 10 times.
✓ Branch 1 taken 3104357 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 10 times.
✗ Branch 7 not taken.
|
6167171 | static const auto pcHighSwInput = getParamFromGroup<Scalar>(paramGroup, "RegularizationHighSw", std::numeric_limits<Scalar>::quiet_NaN()); |
| 505 |
1/2✓ Branch 0 taken 3104367 times.
✗ Branch 1 not taken.
|
6167171 | pcHighSw_ = !isnan(pcHighSwInput) ? pcHighSwInput : params.pcHighSw(); |
| 506 | |||
| 507 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 3104357 times.
|
6167171 | baseLawParamsPtr_ = &m->basicParams(); |
| 508 | |||
| 509 |
4/6✓ Branch 0 taken 10 times.
✓ Branch 1 taken 3104357 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 10 times.
✗ Branch 7 not taken.
|
6167171 | static const auto highSwFixedSlopeInput = getParamFromGroup<Scalar>(paramGroup, "RegularizationHighSwFixedSlope", std::numeric_limits<Scalar>::quiet_NaN()); |
| 510 | 6167171 | highSwFixedSlope_ = highSwFixedSlopeInput; | |
| 511 | |||
| 512 | // Note: we do not pre-calculate all end-point values (e.g., pcLowSwPcValue_ and pcDerivativeLowSw_) | ||
| 513 | // as done in, e.g., VanGenuchten. This is because this law will generally be instantiated only as | ||
| 514 | // a temporary object since all pore bodies generally have different parameters. | ||
| 515 | // We still cache above-mentioned values, but only if they are actually needed (see below). | ||
| 516 | 6167171 | } | |
| 517 | |||
| 518 | // the capillary pressure for the lower saturation threshold | ||
| 519 | Scalar pcLowSwPcValue_() const | ||
| 520 | { | ||
| 521 | // calculated value within first function call, used cached value otherwise | ||
| 522 |
27/40✓ Branch 0 taken 41563 times.
✓ Branch 1 taken 99 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 99 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 99 times.
✓ Branch 6 taken 16721 times.
✓ Branch 7 taken 99 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 99 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 99 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 99 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 99 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 1 times.
✓ Branch 24 taken 1 times.
✓ Branch 25 taken 8 times.
✓ Branch 26 taken 818 times.
✓ Branch 27 taken 8 times.
✓ Branch 28 taken 1 times.
✓ Branch 29 taken 8 times.
✓ Branch 30 taken 1 times.
✓ Branch 31 taken 8 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 396 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 396 times.
✗ Branch 36 not taken.
✓ Branch 37 taken 4 times.
✓ Branch 38 taken 4 times.
✓ Branch 39 taken 32 times.
|
60765 | if (!optionalPcLowSwPcValue_) |
| 523 | 59109 | optionalPcLowSwPcValue_ = BaseLaw::pc(pcLowSw_, *baseLawParamsPtr_); | |
| 524 |
20/20✓ Branch 0 taken 1 times.
✓ Branch 1 taken 41661 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 98 times.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 98 times.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 98 times.
✓ Branch 8 taken 1 times.
✓ Branch 9 taken 98 times.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 98 times.
✓ Branch 12 taken 1 times.
✓ Branch 13 taken 98 times.
✓ Branch 14 taken 1 times.
✓ Branch 15 taken 98 times.
✓ Branch 16 taken 4 times.
✓ Branch 17 taken 392 times.
✓ Branch 18 taken 4 times.
✓ Branch 19 taken 392 times.
|
59868 | return optionalPcLowSwPcValue_.value(); |
| 525 | } | ||
| 526 | |||
| 527 | // dpc_dsw for the lower saturation threshold | ||
| 528 | Scalar pcDerivativeLowSw_() const | ||
| 529 | { | ||
| 530 |
27/40✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 16721 times.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
✓ Branch 16 taken 817 times.
✓ Branch 17 taken 3 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 3 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 3 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 3 times.
✓ Branch 24 taken 1 times.
✓ Branch 25 taken 8 times.
✓ Branch 26 taken 1 times.
✓ Branch 27 taken 8 times.
✓ Branch 28 taken 1 times.
✓ Branch 29 taken 8 times.
✓ Branch 30 taken 1 times.
✓ Branch 31 taken 8 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 4 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 4 times.
✗ Branch 36 not taken.
✓ Branch 37 taken 12 times.
✓ Branch 38 taken 4 times.
✓ Branch 39 taken 32 times.
|
17650 | if (!optionalPcDerivativeLowSw_) |
| 531 | 17546 | optionalPcDerivativeLowSw_ = BaseLaw::dpc_dsw(pcLowSw_, *baseLawParamsPtr_); | |
| 532 | 17650 | return optionalPcDerivativeLowSw_.value(); | |
| 533 | } | ||
| 534 | |||
| 535 | // the capillary pressure for the upper saturation threshold | ||
| 536 | Scalar pcHighSwPcValue_() const | ||
| 537 | { | ||
| 538 |
13/70✓ Branch 0 taken 41563 times.
✓ Branch 1 taken 99 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 99 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 99 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 99 times.
✓ Branch 14 taken 999 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 98 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 98 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 98 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✓ Branch 29 taken 98 times.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✓ Branch 34 taken 130446 times.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✓ Branch 57 taken 396 times.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✓ Branch 61 taken 392 times.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
|
174584 | if (!optionalPcHighSwPcValue_) |
| 539 | 1289970 | optionalPcHighSwPcValue_ = BaseLaw::pc(pcHighSw_, *baseLawParamsPtr_); | |
| 540 |
11/30✗ Branch 0 not taken.
✓ Branch 1 taken 41662 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 99 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 99 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1098 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 98 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 98 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 98 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 98 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 130446 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 396 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 392 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
|
174584 | return optionalPcHighSwPcValue_.value(); |
| 541 | } | ||
| 542 | |||
| 543 | // dpc_dsw at Sw = 1.0 | ||
| 544 | 2053532 | Scalar pcDerivativeHighSwEnd_() const | |
| 545 | { | ||
| 546 | using std::isnan; | ||
| 547 |
1/2✓ Branch 0 taken 1026766 times.
✗ Branch 1 not taken.
|
2053532 | if (!isnan(highSwFixedSlope_)) |
| 548 | return highSwFixedSlope_; | ||
| 549 | else | ||
| 550 |
2/2✓ Branch 0 taken 1026726 times.
✓ Branch 1 taken 40 times.
|
4106984 | return (0.0 - pcHighSwPcValue_()) / (1.0 - pcHighSw_); |
| 551 | } | ||
| 552 | |||
| 553 | 90236 | const Spline<Scalar>& pcSpline_() const | |
| 554 | { | ||
| 555 |
0/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
180472 | if (!optionalPcSpline_) |
| 556 | { | ||
| 557 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
90236 | const auto slopes = highSwSplineZeroSlope_ ? std::array{0.0, 0.0} |
| 558 | ✗ | : std::array{BaseLaw::dpc_dsw(pcHighSw_, *baseLawParamsPtr_), pcDerivativeHighSwEnd_()}; | |
| 559 | |||
| 560 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
180472 | optionalPcSpline_ = Spline<Scalar>(pcHighSw_, 1.0, // x0, x1 |
| 561 | pcHighSwPcValue_(), 0, // y0, y1 | ||
| 562 |
0/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
270708 | slopes[0], slopes[1]); // m0, m1 |
| 563 | } | ||
| 564 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
90236 | return optionalPcSpline_.value(); |
| 565 | } | ||
| 566 | |||
| 567 | Scalar pcLowSw_, pcHighSw_; | ||
| 568 | mutable OptionalScalar<Scalar> optionalPcLowSwPcValue_, optionalPcDerivativeLowSw_; | ||
| 569 | mutable OptionalScalar<Scalar> optionalPcHighSwPcValue_; | ||
| 570 | HighSwRegularizationMethod highSwRegularizationMethod_; | ||
| 571 | const BaseLawParams* baseLawParamsPtr_; | ||
| 572 | mutable std::optional<Spline<Scalar>> optionalPcSpline_; | ||
| 573 | bool highSwSplineZeroSlope_; | ||
| 574 | Scalar highSwFixedSlope_; | ||
| 575 | }; | ||
| 576 | |||
| 577 | /*! | ||
| 578 | * \ingroup Fluidmatrixinteractions | ||
| 579 | * \brief A default configuration for using the VanGenuchten material law | ||
| 580 | */ | ||
| 581 | template<Pore::Shape shape, typename Scalar = double> | ||
| 582 | using TwoPLocalRulesPlatonicBodyDefault = SingleShapeTwoPLocalRules<Scalar, | ||
| 583 | TwoPLocalRulesPlatonicBody<shape>, | ||
| 584 | TwoPLocalRulesPlatonicBodyRegularization<Scalar, | ||
| 585 | TwoPLocalRulesPlatonicBody<shape>>>; | ||
| 586 | |||
| 587 | /*! | ||
| 588 | * \ingroup Fluidmatrixinteractions | ||
| 589 | * \brief A default configuration without regularization for using the VanGenuchten material law | ||
| 590 | */ | ||
| 591 | template<Pore::Shape shape, typename Scalar = double> | ||
| 592 | using TwoPLocalRulesPlatonicBodyNoReg = SingleShapeTwoPLocalRules<Scalar, TwoPLocalRulesPlatonicBody<shape>, Dumux::FluidMatrix::NoRegularization>; | ||
| 593 | |||
| 594 | } // end namespace Dumux::FluidMatrix | ||
| 595 | |||
| 596 | #endif | ||
| 597 |