| 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 | * \brief Implementation of van Genuchten's capillary pressure-saturation relation for three phases. | ||
| 11 | */ | ||
| 12 | #ifndef PARKER_VANGENUCHTEN_3P_HH | ||
| 13 | #define PARKER_VANGENUCHTEN_3P_HH | ||
| 14 | |||
| 15 | #include <algorithm> | ||
| 16 | #include <dumux/common/optionalscalar.hh> | ||
| 17 | #include <dumux/common/parameters.hh> | ||
| 18 | #include <dumux/common/spline.hh> | ||
| 19 | #include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh> | ||
| 20 | #include <dumux/material/fluidmatrixinteractions/2p/noregularization.hh> | ||
| 21 | |||
| 22 | namespace Dumux::FluidMatrix { | ||
| 23 | |||
| 24 | struct ParkerVanGenuchten3PEffToAbsPolicy | ||
| 25 | { | ||
| 26 | /*! | ||
| 27 | * \brief The parameter type | ||
| 28 | * \tparam Scalar The scalar type | ||
| 29 | * \note The efftoabs policy need two parameters: \f$\mathrm{S_{w,r}}, \mathrm{S_{n,r}}\f$. | ||
| 30 | * For the respective formulas check out the description of the free function. | ||
| 31 | */ | ||
| 32 | template<class Scalar> | ||
| 33 | struct Params | ||
| 34 | { | ||
| 35 | 17 | Params(const Scalar swr = 0.0, const Scalar snr = 0.0, const Scalar sgr = 0.0) | |
| 36 | 17 | : swr_(swr), snr_(snr), sgr_(sgr) | |
| 37 | {} | ||
| 38 | |||
| 39 | /*! | ||
| 40 | * \brief Return the residual wetting saturation. | ||
| 41 | */ | ||
| 42 | ✗ | Scalar swr() const | |
| 43 | ✗ | { return swr_; } | |
| 44 | |||
| 45 | /*! | ||
| 46 | * \brief Set the residual wetting saturation. | ||
| 47 | */ | ||
| 48 | ✗ | void setSwr(Scalar v) | |
| 49 | 17 | { swr_ = v; } | |
| 50 | |||
| 51 | /*! | ||
| 52 | * \brief Return the residual nonwetting saturation. | ||
| 53 | */ | ||
| 54 | Scalar snr() const | ||
| 55 | { return snr_; } | ||
| 56 | |||
| 57 | /*! | ||
| 58 | * \brief Set the residual nonwetting saturation. | ||
| 59 | */ | ||
| 60 | ✗ | void setSnr(Scalar v) | |
| 61 | 17 | { snr_ = v; } | |
| 62 | /*! | ||
| 63 | |||
| 64 | * \brief Return the residual gas phase saturation. | ||
| 65 | */ | ||
| 66 | ✗ | Scalar sgr() const | |
| 67 | ✗ | { return sgr_; } | |
| 68 | |||
| 69 | /*! | ||
| 70 | * \brief Set the residual gas phase saturation. | ||
| 71 | */ | ||
| 72 | ✗ | void setSgr(Scalar v) | |
| 73 | 17 | { sgr_ = v; } | |
| 74 | |||
| 75 | bool operator== (const Params& p) const | ||
| 76 | { | ||
| 77 | return Dune::FloatCmp::eq(swr(), p.swr(), 1e-6) | ||
| 78 | && Dune::FloatCmp::eq(snr(), p.snr(), 1e-6) | ||
| 79 | && Dune::FloatCmp::eq(sgr(), p.sgr(), 1e-6); | ||
| 80 | } | ||
| 81 | private: | ||
| 82 | Scalar swr_; | ||
| 83 | Scalar snr_; | ||
| 84 | Scalar sgr_; | ||
| 85 | }; | ||
| 86 | |||
| 87 | /*! | ||
| 88 | * \brief Construct from a subgroup from the global parameter tree | ||
| 89 | * \note This will give you nice error messages if a mandatory parameter is missing | ||
| 90 | */ | ||
| 91 | template<class Scalar> | ||
| 92 | 17 | static Params<Scalar> makeParams(const std::string& paramGroup) | |
| 93 | { | ||
| 94 | 34 | Params<Scalar> params; | |
| 95 | 17 | params.setSwr(getParamFromGroup<Scalar>(paramGroup, "Swr", 0.0)); | |
| 96 | 17 | params.setSnr(getParamFromGroup<Scalar>(paramGroup, "Snr", 0.0)); | |
| 97 | 17 | params.setSgr(getParamFromGroup<Scalar>(paramGroup, "Sgr", 0.0)); | |
| 98 | 17 | return params; | |
| 99 | } | ||
| 100 | |||
| 101 | /*! | ||
| 102 | * \brief Convert an absolute wetting saturation to an effective one. | ||
| 103 | * | ||
| 104 | * \param sw Absolute saturation of the wetting phase \f$\mathrm{[{S}_w]}\f$. | ||
| 105 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 106 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 107 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 108 | * \return Effective saturation of the wetting phase. | ||
| 109 | */ | ||
| 110 | template<class Scalar> | ||
| 111 | static Scalar swToSwe(const Scalar sw, const Params<Scalar>& params) | ||
| 112 | { | ||
| 113 | 18992205 | return (sw - params.swr())/(1.0 - params.swr()); // TODO other residual saturations? | |
| 114 | } | ||
| 115 | |||
| 116 | /*! | ||
| 117 | * \brief Convert an effective wetting saturation to an absolute one. | ||
| 118 | * | ||
| 119 | * \param swe Effective saturation of the nonwetting phase \f$\mathrm{[\overline{S}_n]}\f$. | ||
| 120 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 121 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 122 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 123 | * \return Absolute saturation of the nonwetting phase. | ||
| 124 | */ | ||
| 125 | template<class Scalar> | ||
| 126 | static Scalar sweToSw(const Scalar swe, const Params<Scalar>& params) | ||
| 127 | { | ||
| 128 | 34 | return swe*(1.0 - params.swr()) + params.swr(); // TODO other residual saturations? | |
| 129 | } | ||
| 130 | |||
| 131 | /*! | ||
| 132 | * \brief Derivative of the effective saturation w.r.t. the absolute saturation. | ||
| 133 | * | ||
| 134 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 135 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 136 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 137 | * \return Derivative of the effective saturation w.r.t. the absolute saturation. | ||
| 138 | */ | ||
| 139 | template<class Scalar> | ||
| 140 | static Scalar dswe_dsw(const Params<Scalar>& params) | ||
| 141 | { | ||
| 142 | 68 | return 1.0/(1.0 - params.swr()); // TODO other residual saturations? | |
| 143 | } | ||
| 144 | |||
| 145 | /*! | ||
| 146 | * \brief Derivative of the absolute saturation w.r.t. the effective saturation. | ||
| 147 | * | ||
| 148 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 149 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 150 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 151 | * \return Derivative of the absolute saturation w.r.t. the effective saturation. | ||
| 152 | */ | ||
| 153 | template<class Scalar> | ||
| 154 | static Scalar dsw_dswe(const Params<Scalar>& params) | ||
| 155 | { | ||
| 156 | 17 | return 1.0 - params.swr(); // TODO other residual saturations? | |
| 157 | } | ||
| 158 | |||
| 159 | /*! | ||
| 160 | * \brief Convert an absolute nonwetting saturation to an effective one. | ||
| 161 | * | ||
| 162 | * \param sn Absolute saturation of the nonwetting phase \f$\mathrm{[{S}_n]}\f$. | ||
| 163 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 164 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 165 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 166 | * \return Effective saturation of the nonwetting phase. | ||
| 167 | */ | ||
| 168 | template<class Scalar> | ||
| 169 | ✗ | static Scalar snToSne(const Scalar sn, const Params<Scalar>& params) | |
| 170 | { | ||
| 171 | ✗ | return sn; // sne equals sn // TODO other residual saturations? | |
| 172 | } | ||
| 173 | |||
| 174 | /*! | ||
| 175 | * \brief Convert an absolute total liquid saturation to an effective one. | ||
| 176 | * | ||
| 177 | * \param st Absolute saturation of the total liquid phase (sw+sn) \f$\mathrm{[{S}_n]}\f$. | ||
| 178 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 179 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 180 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 181 | * \return Effective saturation of the nonwetting phase. | ||
| 182 | */ | ||
| 183 | template<class Scalar> | ||
| 184 | static Scalar stToSte(const Scalar st, const Params<Scalar>& params) | ||
| 185 | { | ||
| 186 | 18994016 | return (st-params.swr()) / (1-params.swr()); // TODO other residual saturations? | |
| 187 | } | ||
| 188 | |||
| 189 | /*! | ||
| 190 | * \brief Convert an effective wetting saturation to an absolute one. | ||
| 191 | * | ||
| 192 | * \param ste Effective total liquid (wetting + nonwetting) saturation | ||
| 193 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 194 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 195 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 196 | * \return Absolute saturation of the nonwetting phase. | ||
| 197 | */ | ||
| 198 | template<class Scalar> | ||
| 199 | static Scalar steToSt(const Scalar ste, const Params<Scalar>& params) | ||
| 200 | { | ||
| 201 | return ste*(1.0 - params.swr()) + params.swr(); // TODO other residual saturations? | ||
| 202 | } | ||
| 203 | |||
| 204 | /*! | ||
| 205 | * \brief Derivative of the effective saturation w.r.t. the absolute saturation. | ||
| 206 | * | ||
| 207 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 208 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 209 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 210 | * \return Derivative of the effective saturation w.r.t. the absolute saturation. | ||
| 211 | */ | ||
| 212 | template<class Scalar> | ||
| 213 | static Scalar dste_dst(const Params<Scalar>& params) | ||
| 214 | { | ||
| 215 | 34 | return 1.0/(1.0 - params.swr() /*- params.snr() - params.sgr()*/); // TODO other residual saturations? | |
| 216 | } | ||
| 217 | |||
| 218 | /*! | ||
| 219 | * \brief Derivative of the absolute saturation w.r.t. the effective saturation. | ||
| 220 | * | ||
| 221 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 222 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 223 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 224 | * \return Derivative of the absolute saturation w.r.t. the effective saturation. | ||
| 225 | */ | ||
| 226 | template<class Scalar> | ||
| 227 | static Scalar dst_dste(const Params<Scalar>& params) | ||
| 228 | { | ||
| 229 | 17 | return 1.0 - params.swr(); // TODO other residual saturations? | |
| 230 | } | ||
| 231 | }; | ||
| 232 | |||
| 233 | /*! | ||
| 234 | * \ingroup Fluidmatrixinteractions | ||
| 235 | * \brief Implementation of Parker/vanGenuchten's capillary pressure <-> | ||
| 236 | * saturation relation for three phases. This class bundles the "raw" curves | ||
| 237 | * as static members and doesn't concern itself converting | ||
| 238 | * absolute to effective saturations and vince versa. | ||
| 239 | */ | ||
| 240 | class ParkerVanGenuchten3P | ||
| 241 | { | ||
| 242 | |||
| 243 | public: | ||
| 244 | /*! | ||
| 245 | * \brief The parameter type | ||
| 246 | * \tparam Scalar The scalar type | ||
| 247 | * \note The Parker/vanGenuchten laws are parameterized with four parameters: \f$\mathrm{n, m, \alpha, l}\f$. | ||
| 248 | * | ||
| 249 | * - \f$\mathrm{\alpha}\f$ shape parameter \f$\mathrm{[1/Pa]}\f$ | ||
| 250 | * - \f$\mathrm{n}\f$ shape parameter \f$\mathrm{[-]}\f$ | ||
| 251 | * - \f$\mathrm{swr}\f$ wetting phase residual saturation \f$\mathrm{[-]}\f$ | ||
| 252 | * - \f$\mathrm{swr}\f$ nonwetting phase residual saturation \f$\mathrm{[-]}\f$ | ||
| 253 | * - \f$\mathrm{betaNw}\f$ scaling parameter \f$\mathrm{[-]}\f$ | ||
| 254 | * - \f$\mathrm{betaGn}\f$ scaling parameter \f$\mathrm{[-]}\f$ | ||
| 255 | * - \f$\mathrm{betaGw}\f$ scaling parameter \f$\mathrm{[-]}\f$ | ||
| 256 | * - \f$\mathrm{regardSnr}\f$ determines whether snr is considered for krn or not | ||
| 257 | */ | ||
| 258 | template<class Scalar> | ||
| 259 | struct Params | ||
| 260 | { | ||
| 261 | 17 | Params(Scalar alpha, Scalar n, Scalar swr = 0.0, Scalar snr = 0.0, | |
| 262 | Scalar betaNw = 1.0, Scalar betaGn = 1.0, Scalar betaGw = 1.0, bool regardSnr = false) | ||
| 263 | 17 | : alpha_(alpha), n_(n), m_(1.0 - 1.0/n), swr_(swr), snr_(snr) | |
| 264 | 34 | , betaNw_(betaNw), betaGn_(betaGn), betaGw_(betaGw), regardSnr_(regardSnr) | |
| 265 | {} | ||
| 266 | |||
| 267 | ✗ | Scalar alpha() const { return alpha_; } | |
| 268 | void setAlpha(Scalar alpha) { alpha_ = alpha; } | ||
| 269 | |||
| 270 | ✗ | Scalar m() const { return m_; } | |
| 271 | void setM(Scalar m) { m_ = m; n_ = 1.0/(1.0 - m); } | ||
| 272 | |||
| 273 | ✗ | Scalar n() const{ return n_; } | |
| 274 | void setN(Scalar n){ n_ = n; m_ = 1.0 - 1.0/n; } | ||
| 275 | |||
| 276 | ✗ | Scalar swr() const { return swr_; } | |
| 277 | void setSwr(Scalar swr) { swr_ = swr; } | ||
| 278 | |||
| 279 | ✗ | Scalar snr() const { return snr_; } | |
| 280 | void setSnr(Scalar swr) { snr_ = swr; } | ||
| 281 | |||
| 282 | ✗ | Scalar betaNw() const { return betaNw_; } | |
| 283 | void setBetaNw(Scalar betaNw) { betaNw_ = betaNw; } | ||
| 284 | |||
| 285 | ✗ | Scalar betaGn() const { return betaGn_; } | |
| 286 | void setBetaGn(Scalar betaGn) { betaGn_ = betaGn; } | ||
| 287 | |||
| 288 | ✗ | Scalar betaGw() const { return betaGw_; } | |
| 289 | void setBetaGw(Scalar betaGw) { betaGw_ = betaGw; } | ||
| 290 | |||
| 291 | ✗ | bool regardSnrForKrn() const { return regardSnr_; } | |
| 292 | void setRegardSnrForKrn(bool v) {regardSnr_ = v; } | ||
| 293 | |||
| 294 | bool operator== (const Params& p) const | ||
| 295 | { | ||
| 296 | return Dune::FloatCmp::eq(alpha_, p.alpha_, 1e-6) | ||
| 297 | && Dune::FloatCmp::eq(n_, p.n_, 1e-6) | ||
| 298 | && Dune::FloatCmp::eq(m_, p.m_, 1e-6) | ||
| 299 | && Dune::FloatCmp::eq(swr_, p.swr_, 1e-6) | ||
| 300 | && Dune::FloatCmp::eq(snr_, p.snr_, 1e-6) | ||
| 301 | && Dune::FloatCmp::eq(betaNw_, p.betaNw_, 1e-6) | ||
| 302 | && Dune::FloatCmp::eq(betaGn_, p.betaGn_, 1e-6) | ||
| 303 | && Dune::FloatCmp::eq(betaGw_, p.betaGw_, 1e-6) | ||
| 304 | && regardSnr_ == p.regardSnr_; | ||
| 305 | } | ||
| 306 | |||
| 307 | private: | ||
| 308 | Scalar alpha_, n_, m_, swr_, snr_, betaNw_, betaGn_, betaGw_; | ||
| 309 | bool regardSnr_; | ||
| 310 | }; | ||
| 311 | |||
| 312 | /*! | ||
| 313 | * \brief Construct from a subgroup from the global parameter tree | ||
| 314 | * \note This will give you nice error messages if a mandatory parameter is missing | ||
| 315 | */ | ||
| 316 | template<class Scalar = double> | ||
| 317 | 17 | static Params<Scalar> makeParams(const std::string& paramGroup) | |
| 318 | { | ||
| 319 | 17 | const auto n = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenN"); | |
| 320 | 17 | const auto alpha = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenAlpha"); | |
| 321 | 17 | const auto swr = getParamFromGroup<Scalar>(paramGroup, "Swr", 0.0); | |
| 322 | 17 | const auto snr = getParamFromGroup<Scalar>(paramGroup, "Snr", 0.0); | |
| 323 | 17 | const auto betaNw = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenBetaNw", 1.0); | |
| 324 | 17 | const auto betaGn = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenBetaGn", 1.0); | |
| 325 | 17 | const auto betaGw = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenBetaGw", 1.0); | |
| 326 | 17 | const auto regardSnr = getParamFromGroup<bool>(paramGroup, "ParkerVanGenuchtenRegardSnrForKrn", false); | |
| 327 | return Params<Scalar>(alpha, n, swr, snr, | ||
| 328 | 34 | betaNw, betaGn, betaGw, regardSnr ); | |
| 329 | } | ||
| 330 | |||
| 331 | /*! | ||
| 332 | * \brief The capillary pressure-saturation curve for the gas and wetting phase | ||
| 333 | * \param params Array of parameters | ||
| 334 | * \param swe Effective wetting phase saturation | ||
| 335 | */ | ||
| 336 | template<class Scalar> | ||
| 337 | 4256623 | static Scalar pcgw(Scalar swe, const Params<Scalar>& params) | |
| 338 | { | ||
| 339 |
2/4✓ Branch 0 taken 4256623 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4256623 times.
|
4256623 | assert(0 <= swe && swe <= 1); |
| 340 | 4256623 | return pc_(swe, params); | |
| 341 | } | ||
| 342 | |||
| 343 | /*! | ||
| 344 | * \brief The capillary pressure-saturation curve for the non-wettigng and wetting phase | ||
| 345 | * \param params Array of parameters | ||
| 346 | * \param swe Effective wetting phase saturation | ||
| 347 | */ | ||
| 348 | template<class Scalar> | ||
| 349 | 4248444 | static Scalar pcnw(Scalar swe, const Params<Scalar>& params) | |
| 350 | { | ||
| 351 |
2/4✓ Branch 0 taken 4248444 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4248444 times.
|
4248444 | assert(0 <= swe && swe <= 1); |
| 352 | 4248444 | return pc_(swe, params)/params.betaNw(); | |
| 353 | } | ||
| 354 | |||
| 355 | /*! | ||
| 356 | * \brief The capillary pressure-saturation curve for the gas and nonwetting phase | ||
| 357 | * \param params Array of parameters | ||
| 358 | * \param ste Effective total liquid (wetting + nonwetting) saturation | ||
| 359 | */ | ||
| 360 | template<class Scalar> | ||
| 361 | 1332557 | static Scalar pcgn(const Scalar ste, const Params<Scalar>& params) | |
| 362 | { | ||
| 363 |
2/4✓ Branch 0 taken 1332557 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1332557 times.
|
1332557 | assert(0 <= ste && ste <= 1); |
| 364 | 1332557 | return pc_(ste, params)/params.betaGn(); | |
| 365 | } | ||
| 366 | |||
| 367 | /*! | ||
| 368 | * \brief This function ensures a continuous transition from 2 to 3 phases and vice versa | ||
| 369 | * \param params Array of parameters | ||
| 370 | * \param sne Non-wetting liquid saturation | ||
| 371 | */ | ||
| 372 | template<class Scalar> | ||
| 373 | static Scalar pcAlpha(Scalar sne, const Params<Scalar>& params) | ||
| 374 | { | ||
| 375 | /* regularization */ | ||
| 376 | 9494993 | if (sne <= 0.001) | |
| 377 | 4607687 | sne = 0.0; | |
| 378 |
1/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 9494993 times.
|
9494993 | if (sne >= 1.0) |
| 379 | ✗ | sne = 1.0; | |
| 380 | |||
| 381 |
2/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4773277 times.
✓ Branch 5 taken 4721716 times.
|
9494993 | if (sne > params.snr()) |
| 382 | return 1.0; | ||
| 383 | else | ||
| 384 | { | ||
| 385 |
1/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4773277 times.
✗ Branch 5 not taken.
|
4773277 | if (params.snr() >= 0.001) |
| 386 | 4773277 | return sne/params.snr(); | |
| 387 | else | ||
| 388 | return 0.0; | ||
| 389 | }; | ||
| 390 | } | ||
| 391 | |||
| 392 | /*! | ||
| 393 | * \brief Returns the partial derivative of the capillary | ||
| 394 | * pressure to the effective saturation. | ||
| 395 | * \param swe Effective wetting phase saturation | ||
| 396 | * \param params Array of parameters | ||
| 397 | */ | ||
| 398 | template<class Scalar> | ||
| 399 | 34 | static Scalar dpcgw_dswe(const Scalar swe, const Params<Scalar>& params) | |
| 400 | { | ||
| 401 | using std::pow; | ||
| 402 | 34 | const Scalar powSeRegu = pow(swe, -1/params.m()); | |
| 403 | 34 | return - 1.0/params.alpha() * pow(powSeRegu - 1, 1.0/params.n() - 1)/params.n() | |
| 404 | 34 | * powSeRegu/swe/params.m()/params.betaGw(); | |
| 405 | } | ||
| 406 | |||
| 407 | /*! | ||
| 408 | * \brief Returns the partial derivative of the capillary | ||
| 409 | * pressure to the effective saturation. | ||
| 410 | * \param swe Effective wetting phase saturation | ||
| 411 | * \param params Array of parameters | ||
| 412 | */ | ||
| 413 | template<class Scalar> | ||
| 414 | 34 | static Scalar dpcnw_dswe(const Scalar swe, const Params<Scalar>& params) | |
| 415 | { | ||
| 416 | using std::pow; | ||
| 417 | 34 | const Scalar powSeRegu = pow(swe, -1/params.m()); | |
| 418 | 34 | return - 1.0/params.alpha() * pow(powSeRegu - 1, 1.0/params.n() - 1)/params.n() | |
| 419 | 34 | * powSeRegu/swe/params.m()/params.betaNw(); | |
| 420 | } | ||
| 421 | |||
| 422 | /*! | ||
| 423 | * \brief Returns the partial derivative of the capillary | ||
| 424 | * pressure to the effective saturation. | ||
| 425 | * \param ste Effective total liquid (wetting + nonwetting) saturation | ||
| 426 | * \param params Array of parameters | ||
| 427 | */ | ||
| 428 | template<class Scalar> | ||
| 429 | 34 | static Scalar dpcgn_dste(const Scalar ste, const Params<Scalar>& params) | |
| 430 | { | ||
| 431 | using std::pow; | ||
| 432 | 34 | const Scalar powSeRegu = pow(ste, -1/params.m()); | |
| 433 | 34 | return - 1.0/params.alpha() * pow(powSeRegu - 1, 1.0/params.n() - 1)/params.n() | |
| 434 | 34 | * powSeRegu/ste/params.m()/params.betaGn(); | |
| 435 | } | ||
| 436 | |||
| 437 | /*! | ||
| 438 | * \brief The relative permeability for the wetting phase of | ||
| 439 | * the medium implied by van Genuchten's | ||
| 440 | * parameterization. | ||
| 441 | * | ||
| 442 | * The permeability of water in a 3p system equals the standard 2p description. | ||
| 443 | * (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models" | ||
| 444 | * MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.) \cite delshad1989 <BR> | ||
| 445 | * | ||
| 446 | * \param swe Effective wetting phase saturation | ||
| 447 | * \param params Array of parameters. | ||
| 448 | */ | ||
| 449 | template<class Scalar> | ||
| 450 | 5582914 | static Scalar krw(const Scalar swe, const Params<Scalar>& params) | |
| 451 | { | ||
| 452 | using std::pow; | ||
| 453 | using std::sqrt; | ||
| 454 | 5582914 | const Scalar r = 1.0 - pow(1 - pow(swe, 1/params.m()), params.m()); | |
| 455 | 5582914 | return sqrt(swe)*r*r; | |
| 456 | } | ||
| 457 | |||
| 458 | /*! | ||
| 459 | * \brief The relative permeability for the nonwetting phase | ||
| 460 | * after the Model of Parker et al. (1987). | ||
| 461 | * | ||
| 462 | * See model 7 in "Comparison of the Three-Phase Oil Relative Permeability Models" | ||
| 463 | * MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83 \cite delshad1989 <BR> | ||
| 464 | * or more comprehensive in | ||
| 465 | * "Estimation of primary drainage three-phase relative permeability for organic | ||
| 466 | * liquid transport in the vadose zone", Leonardo I. Oliveira, Avery H. Demond, | ||
| 467 | * Journal of Contaminant Hydrology 66 (2003), 261-285 \cite oliveira2003 <BR> | ||
| 468 | * | ||
| 469 | * | ||
| 470 | * \param params Array of parameters. | ||
| 471 | * \param swe Effective wetting phase saturation | ||
| 472 | * \param sn Absolute nonwetting liquid saturation | ||
| 473 | * \param ste Effective total liquid (wetting + nonwetting) saturation | ||
| 474 | */ | ||
| 475 | template<class Scalar> | ||
| 476 | 5445755 | static Scalar krn(const Scalar swe, const Scalar sn, const Scalar ste, const Params<Scalar>& params) | |
| 477 | { | ||
| 478 | Scalar krn; | ||
| 479 | using std::pow; | ||
| 480 | 5445755 | krn = pow(1 - pow(swe, 1/params.m()), params.m()); | |
| 481 | 5445755 | krn -= pow(1 - pow(ste, 1/params.m()), params.m()); | |
| 482 | 5445755 | krn *= krn; | |
| 483 | |||
| 484 | using std::clamp; | ||
| 485 | using std::sqrt; | ||
| 486 |
2/2✓ Branch 0 taken 1543969 times.
✓ Branch 1 taken 3901786 times.
|
5445755 | if (params.regardSnrForKrn()) |
| 487 | { | ||
| 488 | // regard Snr in the permeability of the n-phase, see Helmig1997 | ||
| 489 |
2/2✓ Branch 0 taken 787493 times.
✓ Branch 1 taken 756476 times.
|
1543969 | const Scalar resIncluded = clamp(sn - params.snr()/ (1-params.swr()), 0.0, 1.0); |
| 490 | 1543969 | krn *= sqrt(resIncluded); | |
| 491 | } | ||
| 492 | else | ||
| 493 | 3901786 | krn *= sqrt(sn / (1 - params.swr())); // Hint: (ste - swe) = sn / (1-Swr) | |
| 494 | |||
| 495 | 5445755 | return krn; | |
| 496 | } | ||
| 497 | |||
| 498 | /*! | ||
| 499 | * \brief The relative permeability for the nonwetting phase | ||
| 500 | * of the medium implied by van Genuchten's | ||
| 501 | * parameterization. | ||
| 502 | * | ||
| 503 | * The permeability of gas in a 3p system equals the standard 2p description. | ||
| 504 | * (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models" | ||
| 505 | * MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.) \cite delshad1989 <BR> | ||
| 506 | * | ||
| 507 | * \param params Array of parameters. | ||
| 508 | * \param ste Effective total liquid (wetting + nonwetting) saturation | ||
| 509 | */ | ||
| 510 | template<class Scalar> | ||
| 511 | 1433733 | static Scalar krg(const Scalar ste, const Params<Scalar>& params) | |
| 512 | { | ||
| 513 |
2/4✓ Branch 0 taken 1433733 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1433733 times.
|
1433733 | assert(0 <= ste && ste <= 1); |
| 514 | using std::cbrt; | ||
| 515 | using std::pow; | ||
| 516 | 1433733 | return cbrt(1 - ste) * pow(1 - pow(ste, 1/params.m()), 2*params.m()); | |
| 517 | } | ||
| 518 | |||
| 519 | /*! | ||
| 520 | * \brief The derivative of the relative permeability for the | ||
| 521 | * gas phase in regard to the total liquid saturation of | ||
| 522 | * the medium as implied by the van Genuchten | ||
| 523 | * parameterization. | ||
| 524 | * | ||
| 525 | * \param ste The mobile total liquid saturation. | ||
| 526 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 527 | */ | ||
| 528 | template<class Scalar> | ||
| 529 | 17 | static Scalar dkrg_dste(const Scalar ste, const Params<Scalar>& params) | |
| 530 | { | ||
| 531 |
2/4✓ Branch 0 taken 17 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 17 times.
|
17 | assert(0 < ste && ste <= 1); |
| 532 | |||
| 533 | using std::pow; | ||
| 534 | 17 | const Scalar x = pow(ste, 1.0/params.m()); | |
| 535 | return | ||
| 536 | 17 | -pow(1.0 - x, 2*params.m()) | |
| 537 | 17 | *pow(1.0 - ste, -2.0/3) | |
| 538 | 17 | *(1.0/3 + 2*x/ste); | |
| 539 | } | ||
| 540 | |||
| 541 | /*! | ||
| 542 | * \brief The relative permeability for a phase. | ||
| 543 | * \param params Array of parameters. | ||
| 544 | * \param phaseIdx Indicator, The saturation of all phases. | ||
| 545 | * \param swe Effective wetting phase saturation | ||
| 546 | * \param sne Effective nonwetting saturation | ||
| 547 | */ | ||
| 548 | template<class Scalar> | ||
| 549 | static Scalar kr(const int phaseIdx, const Scalar swe, const Scalar sne, const Params<Scalar>& params) | ||
| 550 | { | ||
| 551 | switch (phaseIdx) | ||
| 552 | { | ||
| 553 | case 0: | ||
| 554 | return krw(params, swe, sne); | ||
| 555 | case 1: | ||
| 556 | return krn(params, swe, sne); | ||
| 557 | case 2: | ||
| 558 | return krg(params, swe, sne); | ||
| 559 | } | ||
| 560 | DUNE_THROW(Dune::InvalidStateException, | ||
| 561 | "Invalid phase index "); | ||
| 562 | } | ||
| 563 | |||
| 564 | private: | ||
| 565 | |||
| 566 | /*! | ||
| 567 | * \brief The standard van Genuchten two-phase pc-S relation either with respect to | ||
| 568 | * the effective wetting phase saturation Swe or the effective total liquid saturation Ste. | ||
| 569 | * \param se Effective wetting phase or total liquid saturation | ||
| 570 | * \param params Array of parameters. | ||
| 571 | */ | ||
| 572 | template<class Scalar> | ||
| 573 | 9837624 | const static Scalar pc_(const Scalar se, const Params<Scalar>& params) | |
| 574 | { | ||
| 575 | using std::pow; | ||
| 576 | 9837624 | return pow(pow(se, -1/params.m()) - 1, 1/params.n())/params.alpha(); | |
| 577 | } | ||
| 578 | |||
| 579 | }; | ||
| 580 | |||
| 581 | /*! | ||
| 582 | * \ingroup Fluidmatrixinteractions | ||
| 583 | * \brief A regularization for the ParkerVanGenuchten3PRegularization material law | ||
| 584 | * \note Regularization can be turned of by setting the threshold parameters | ||
| 585 | * out of range (runtime) or by replacing | ||
| 586 | * this class by NoRegularization (compile time). | ||
| 587 | */ | ||
| 588 | template <class Scalar> | ||
| 589 | 68 | class ParkerVanGenuchten3PRegularization | |
| 590 | { | ||
| 591 | using BaseLawParams = typename ParkerVanGenuchten3P::Params<Scalar>; | ||
| 592 | |||
| 593 | public: | ||
| 594 | //! Regularization parameters | ||
| 595 | template<class S> | ||
| 596 | struct Params | ||
| 597 | { | ||
| 598 | /*! | ||
| 599 | * \brief Set the threshold saturation below which the capillary pressure is regularized. | ||
| 600 | * | ||
| 601 | * Most problems are very sensitive to this value (e.g. making it smaller might | ||
| 602 | * result in very high capillary pressures) | ||
| 603 | */ | ||
| 604 | void setPcLowSwe(Scalar pcLowSwe) | ||
| 605 | { pcLowSwe_ = pcLowSwe; } | ||
| 606 | |||
| 607 | /*! | ||
| 608 | * \brief Threshold saturation below which the capillary pressure is regularized. | ||
| 609 | */ | ||
| 610 | Scalar pcLowSwe() const | ||
| 611 | { return pcLowSwe_; } | ||
| 612 | |||
| 613 | /*! | ||
| 614 | * \brief Set the threshold saturation above which the capillary pressure is regularized. | ||
| 615 | */ | ||
| 616 | void setPcHighSwe(Scalar pcHighSwe) | ||
| 617 | { pcHighSwe_ = pcHighSwe; } | ||
| 618 | |||
| 619 | /*! | ||
| 620 | * \brief Threshold saturation above which the capillary pressure is regularized. | ||
| 621 | * | ||
| 622 | * Most problems are very sensitive to this value (e.g. making it smaller might | ||
| 623 | * result in negative capillary pressures). | ||
| 624 | */ | ||
| 625 | Scalar pcHighSwe() const | ||
| 626 | { return pcHighSwe_; } | ||
| 627 | |||
| 628 | /*! | ||
| 629 | * \brief Set the threshold saturation below which the relative | ||
| 630 | * permeability of the nonwetting phase gets regularized. | ||
| 631 | */ | ||
| 632 | void setKrnLowSwe(Scalar krnLowSwe) | ||
| 633 | { krnLowSwe_ = krnLowSwe; } | ||
| 634 | |||
| 635 | /*! | ||
| 636 | * \brief Threshold saturation below which the relative | ||
| 637 | * permeability of the nonwetting phase gets regularized. | ||
| 638 | */ | ||
| 639 | Scalar krnLowSwe() const | ||
| 640 | { return krnLowSwe_; } | ||
| 641 | |||
| 642 | /*! | ||
| 643 | * \brief Set the threshold saturation below which the relative | ||
| 644 | * permeability of the nonwetting phase gets regularized. | ||
| 645 | */ | ||
| 646 | void setKrgLowSte(Scalar krgLowSte) | ||
| 647 | { krgLowSte_ = krgLowSte; } | ||
| 648 | |||
| 649 | /*! | ||
| 650 | * \brief Threshold saturation below which the relative | ||
| 651 | * permeability of the nonwetting phase gets regularized. | ||
| 652 | */ | ||
| 653 | Scalar krgLowSte() const | ||
| 654 | { return krgLowSte_; } | ||
| 655 | |||
| 656 | /*! | ||
| 657 | * \brief Set the threshold saturation above which the relative | ||
| 658 | * permeability of the wetting phase gets regularized. | ||
| 659 | */ | ||
| 660 | void setKrwHighSwe(Scalar krwHighSwe) | ||
| 661 | { krwHighSwe_ = krwHighSwe; } | ||
| 662 | |||
| 663 | /*! | ||
| 664 | * \brief Threshold saturation above which the relative | ||
| 665 | * permeability of the wetting phase gets regularized. | ||
| 666 | */ | ||
| 667 | Scalar krwHighSwe() const | ||
| 668 | { return krwHighSwe_; } | ||
| 669 | |||
| 670 | |||
| 671 | /*! | ||
| 672 | * \brief Choose whether to use a constant value for regularization of the | ||
| 673 | * pc-S curves or not | ||
| 674 | * \param input True or false | ||
| 675 | */ | ||
| 676 | void setConstRegularization(const bool input) | ||
| 677 | { constRegularization_ = input; } | ||
| 678 | |||
| 679 | /*! | ||
| 680 | * \brief Returns whether to use a constant value for regularization of the | ||
| 681 | * pc-S curves or not | ||
| 682 | */ | ||
| 683 | bool constRegularization() const | ||
| 684 | { return constRegularization_; } | ||
| 685 | |||
| 686 | private: | ||
| 687 | S pcLowSwe_ = 0.01; | ||
| 688 | S pcHighSwe_ = 0.99; | ||
| 689 | S krnLowSwe_ = 0.1; | ||
| 690 | S krwHighSwe_ = 0.9; | ||
| 691 | S krgLowSte_ = 1e-3; | ||
| 692 | bool constRegularization_ = false; | ||
| 693 | }; | ||
| 694 | |||
| 695 | //! Initialize the spline | ||
| 696 | template<class MaterialLaw> | ||
| 697 | 17 | void init(const MaterialLaw* m, const std::string& paramGroup) | |
| 698 | { | ||
| 699 | 17 | pcLowSwe_ = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenPcLowSweThreshold", 0.01); | |
| 700 | 17 | pcHighSwe_ = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenPcHighSweThreshold", 0.99); | |
| 701 | 17 | krwHighSwe_ = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenKrwHighSweThreshold", 0.9); | |
| 702 | 17 | krnLowSwe_ = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenKrnLowSweThreshold", 0.1); | |
| 703 | 17 | krgLowSte_ = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenKrgLowSteThreshold", 1e-3); | |
| 704 | 17 | constRegularization_ = getParamFromGroup<bool>(paramGroup, "VanGenuchtenConstantRegularization", false); | |
| 705 | |||
| 706 | 17 | initPcParameters_(m, pcLowSwe_, pcHighSwe_); | |
| 707 | 17 | initKrParameters_(m, krnLowSwe_, krwHighSwe_); | |
| 708 | 17 | } | |
| 709 | |||
| 710 | template<class MaterialLaw, class BaseParams, class EffToAbsParams> | ||
| 711 | void init(const MaterialLaw* m, const BaseParams& bp, const EffToAbsParams& etap, const Params<Scalar>& p) | ||
| 712 | { | ||
| 713 | pcLowSwe_ = p.pcLowSwe(); | ||
| 714 | pcHighSwe_ = p.pcHighSwe(); | ||
| 715 | krwHighSwe_ = p.krwHighSwe(); | ||
| 716 | krnLowSwe_ = p.krnLowSwe(); | ||
| 717 | krgLowSte_ = p.krgLowSte(); | ||
| 718 | constRegularization_ = p.constRegularization(); | ||
| 719 | |||
| 720 | initPcParameters_(m, pcLowSwe_, pcHighSwe_); | ||
| 721 | initKrParameters_(m, krnLowSwe_, krwHighSwe_); | ||
| 722 | } | ||
| 723 | |||
| 724 | /*! | ||
| 725 | * \brief Equality comparison with another instance | ||
| 726 | */ | ||
| 727 | bool operator== (const ParkerVanGenuchten3PRegularization& o) const | ||
| 728 | { | ||
| 729 | return Dune::FloatCmp::eq(pcLowSwe_, o.pcLowSwe_, 1e-6) | ||
| 730 | && Dune::FloatCmp::eq(pcHighSwe_, o.pcHighSwe_, 1e-6) | ||
| 731 | && Dune::FloatCmp::eq(krwHighSwe_, o.krwHighSwe_, 1e-6) | ||
| 732 | && Dune::FloatCmp::eq(krnLowSwe_, o.krnLowSwe_, 1e-6) | ||
| 733 | && constRegularization_ == o.constRegularization_; | ||
| 734 | } | ||
| 735 | |||
| 736 | /*! | ||
| 737 | * \brief The regularized capillary pressure-saturation curve for the gas and wetting phase | ||
| 738 | * regularized part: | ||
| 739 | * - low saturation: extend the \f$\mathrm{p_{cgw}(S_{we})}\f$ curve with the slope at the regularization point (i.e. no kink). | ||
| 740 | * - high saturation: connect the high regularization point with | ||
| 741 | * with a spline and continue linearly for \f$\mathrm{S_{we} > 1}\f$ | ||
| 742 | * \param swe Effective wetting phase saturation | ||
| 743 | */ | ||
| 744 | 9523513 | OptionalScalar<Scalar> pcgw(Scalar swe) const | |
| 745 | { | ||
| 746 | // if specified, a constant value is used for regularization | ||
| 747 | using std::clamp; | ||
| 748 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9523513 times.
|
9523513 | if (constRegularization_) |
| 749 | ✗ | swe = clamp(swe, 0.0, 1.0); | |
| 750 | |||
| 751 | // make sure that the capillary pressure observes a derivative | ||
| 752 | // != 0 for 'illegal' saturations. This is favourable for the | ||
| 753 | // newton solver (if the derivative is calculated numerically) | ||
| 754 | // in order to get the saturation moving to the right | ||
| 755 | // direction if it temporarily is in an 'illegal' range. | ||
| 756 |
2/2✓ Branch 0 taken 4556958 times.
✓ Branch 1 taken 4966555 times.
|
9523513 | if (swe < pcLowSwe_) |
| 757 | 4556958 | return pcgwLowSwePcgwValue_ + pcgwDerivativeLowSw_*(swe - pcLowSwe_); | |
| 758 | |||
| 759 |
2/2✓ Branch 0 taken 146644 times.
✓ Branch 1 taken 4819911 times.
|
4966555 | else if (swe > 1.0) |
| 760 | 146644 | return pcgwDerivativeHighSweEnd_*(swe - 1.0); | |
| 761 | |||
| 762 |
2/2✓ Branch 0 taken 563339 times.
✓ Branch 1 taken 4256572 times.
|
4819911 | else if (swe > pcHighSwe_) |
| 763 | 563339 | return pcgwSpline_.eval(swe); | |
| 764 | |||
| 765 | else | ||
| 766 | 4256572 | return {}; // no regularization | |
| 767 | } | ||
| 768 | |||
| 769 | /*! | ||
| 770 | * \brief The regularized capillary pressure-saturation curve for the nonwetting and wetting phase | ||
| 771 | * regularized part: | ||
| 772 | * - low saturation: extend the \f$\mathrm{p_{cnw}(S_{we})}\f$ curve with the slope at the regularization point (i.e. no kink). | ||
| 773 | * - high saturation: connect the high regularization point with | ||
| 774 | * with a spline and continue linearly for \f$\mathrm{S_{we} > 1}\f$ | ||
| 775 | * \param swe Effective wetting phase saturation | ||
| 776 | */ | ||
| 777 | 9494993 | OptionalScalar<Scalar> pcnw(Scalar swe) const | |
| 778 | { | ||
| 779 | // if specified, a constant value is used for regularization | ||
| 780 | using std::clamp; | ||
| 781 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9494993 times.
|
9494993 | if (constRegularization_) |
| 782 | ✗ | swe = clamp(swe, 0.0, 1.0); | |
| 783 | |||
| 784 | // make sure that the capillary pressure observes a derivative | ||
| 785 | // != 0 for 'illegal' saturations. This is favourable for the | ||
| 786 | // newton solver (if the derivative is calculated numerically) | ||
| 787 | // in order to get the saturation moving to the right | ||
| 788 | // direction if it temporarily is in an 'illegal' range. | ||
| 789 |
2/2✓ Branch 0 taken 4538178 times.
✓ Branch 1 taken 4956815 times.
|
9494993 | if (swe < pcLowSwe_) |
| 790 | 4538178 | return pcnwLowSwePcnwValue_ + pcnwDerivativeLowSw_*(swe - pcLowSwe_); | |
| 791 | |||
| 792 |
2/2✓ Branch 0 taken 146644 times.
✓ Branch 1 taken 4810171 times.
|
4956815 | else if (swe > 1.0) |
| 793 | 146644 | return pcnwDerivativeHighSweEnd_*(swe - 1.0); | |
| 794 | |||
| 795 |
2/2✓ Branch 0 taken 561778 times.
✓ Branch 1 taken 4248393 times.
|
4810171 | else if (swe > pcHighSwe_) |
| 796 | 561778 | return pcnwSpline_.eval(swe); | |
| 797 | |||
| 798 | else | ||
| 799 | 4248393 | return {}; // no regularization | |
| 800 | } | ||
| 801 | |||
| 802 | /*! | ||
| 803 | * \brief The regularized capillary pressure-saturation curve for the gas and nonwetting phase | ||
| 804 | * regularized part: | ||
| 805 | * - low saturation: extend the \f$\mathrm{p_{cgn}(S_{teff})}\f$ curve with the slope at the regularization point (i.e. no kink). | ||
| 806 | * - high saturation: connect the high regularization point with | ||
| 807 | * with a spline and continue linearly for \f$\mathrm{S_{teff} > 1}\f$ | ||
| 808 | * \param ste Effective total liquid (sw + sn) saturation | ||
| 809 | */ | ||
| 810 | 9494993 | OptionalScalar<Scalar> pcgn(Scalar ste) const | |
| 811 | { | ||
| 812 | // if specified, a constant value is used for regularization | ||
| 813 | using std::clamp; | ||
| 814 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9494993 times.
|
9494993 | if (constRegularization_) |
| 815 | ✗ | ste = clamp(ste, 0.0, 1.0); | |
| 816 | |||
| 817 | |||
| 818 | // make sure that the capillary pressure observes a derivative | ||
| 819 | // != 0 for 'illegal' saturations. This is favourable for the | ||
| 820 | // newton solver (if the derivative is calculated numerically) | ||
| 821 | // in order to get the saturation moving to the right | ||
| 822 | // direction if it temporarily is in an 'illegal' range. | ||
| 823 | 9494993 | const Scalar pcLowSte = pcLowSwe_; | |
| 824 | 9494993 | const Scalar pcHighSte = pcHighSwe_; | |
| 825 |
2/2✓ Branch 0 taken 3598985 times.
✓ Branch 1 taken 5896008 times.
|
9494993 | if (ste < pcLowSte) |
| 826 | 3598985 | return pcgnLowStePcgnValue_ + pcgnDerivativeLowSt_*(ste - pcLowSte); | |
| 827 | |||
| 828 |
2/2✓ Branch 0 taken 708442 times.
✓ Branch 1 taken 5187566 times.
|
5896008 | else if (ste > 1.0) |
| 829 | 708442 | return pcgnDerivativeHighSteEnd_*(ste - 1.0); | |
| 830 | |||
| 831 |
2/2✓ Branch 0 taken 3855060 times.
✓ Branch 1 taken 1332506 times.
|
5187566 | else if (ste > pcHighSte) |
| 832 | 3855060 | return pcgnSpline_.eval(ste); | |
| 833 | |||
| 834 | else | ||
| 835 | 1332506 | return {}; // no regularization | |
| 836 | } | ||
| 837 | |||
| 838 | /*! | ||
| 839 | * \brief This function ensures a continuous transition from 2 to 3 phases and vice versa | ||
| 840 | * \param sne Effective nonwetting liquid saturation | ||
| 841 | */ | ||
| 842 | ✗ | OptionalScalar<Scalar> pcAlpha(Scalar sne) const | |
| 843 | { | ||
| 844 | // no regularization | ||
| 845 | 9494993 | return {}; | |
| 846 | } | ||
| 847 | |||
| 848 | /*! | ||
| 849 | * \brief The regularized relative permeability for the wetting phase | ||
| 850 | * \param swe Effective wetting phase saturation | ||
| 851 | */ | ||
| 852 | ✗ | OptionalScalar<Scalar> krw(const Scalar swe) const | |
| 853 | { | ||
| 854 | ✗ | if (swe < 0.0) | |
| 855 | ✗ | return 0.0; | |
| 856 | ✗ | else if (swe > 1.0 - std::numeric_limits<Scalar>::epsilon()) | |
| 857 | ✗ | return 1.0; | |
| 858 | else | ||
| 859 | ✗ | return {}; // no regularization | |
| 860 | } | ||
| 861 | |||
| 862 | |||
| 863 | /*! | ||
| 864 | * \brief The regularized relative permeability for the nonwetting phase | ||
| 865 | * \param swe Effective wetting phase saturation | ||
| 866 | * \param sn Nonwetting saturation | ||
| 867 | * \param ste Effective total (wetting + nonwetting) saturation | ||
| 868 | */ | ||
| 869 | 9496991 | OptionalScalar<Scalar> krn(Scalar swe, const Scalar sn, Scalar ste) const | |
| 870 | { | ||
| 871 | using std::clamp; | ||
| 872 |
2/2✓ Branch 0 taken 5744736 times.
✓ Branch 1 taken 3752255 times.
|
9496991 | swe = clamp(swe, 0.0, 1.0); |
| 873 |
2/2✓ Branch 0 taken 6630198 times.
✓ Branch 1 taken 2866793 times.
|
9496991 | ste = clamp(ste, 0.0, 1.0); |
| 874 | |||
| 875 |
2/2✓ Branch 0 taken 4051236 times.
✓ Branch 1 taken 5445755 times.
|
9496991 | if (ste - swe <= 0.0) |
| 876 | 4051236 | return 0.0; | |
| 877 | else | ||
| 878 | 5445755 | return ParkerVanGenuchten3P::krn(swe, sn, ste, *baseLawParamsPtr_); | |
| 879 | } | ||
| 880 | |||
| 881 | /*! | ||
| 882 | * \brief The regularized relative permeability for the gas phase | ||
| 883 | * \param ste Effective total (wetting + nonwetting) saturation | ||
| 884 | */ | ||
| 885 | 9496991 | OptionalScalar<Scalar> krg(const Scalar ste) const | |
| 886 | { | ||
| 887 | //return 0 if there is no gas | ||
| 888 |
2/2✓ Branch 0 taken 4563032 times.
✓ Branch 1 taken 4933959 times.
|
9496991 | if (ste > 1.0 - std::numeric_limits<Scalar>::epsilon()) |
| 889 | 4563032 | return 0.0; | |
| 890 | |||
| 891 | // use linear regularization for very high gas saturations | ||
| 892 | // to avoid a kink in the curve and to maintain a slope for | ||
| 893 | // the Newton solver | ||
| 894 |
2/2✓ Branch 0 taken 3500243 times.
✓ Branch 1 taken 1433716 times.
|
4933959 | if (ste <= krgLowSte_) |
| 895 | 3500243 | return krgLowStkrgValue_ + krgDerivativeLowSt_*(ste - krgLowSte_); | |
| 896 | else | ||
| 897 | { | ||
| 898 | // For very low gas saturations: | ||
| 899 | // We use a scaling factor that decreases the gas phase permeability quite fast a very low gas phase | ||
| 900 | // saturations, thus making that phase virtually immobile. | ||
| 901 | // This prevents numerical issues related to the degeneration of the gas phase mass balance for the 3p3c model | ||
| 902 | // at very low gas phase saturations. | ||
| 903 | |||
| 904 | // get the absolute gas phase saturation | ||
| 905 | 1433716 | const Scalar st = ste*(1 - swr_) + swr_; | |
| 906 | 1433716 | const Scalar sg = 1.0 - st; | |
| 907 | |||
| 908 | // do not regularize | ||
| 909 |
2/2✓ Branch 0 taken 1117928 times.
✓ Branch 1 taken 315788 times.
|
1433716 | if (sg > 0.1) |
| 910 | 1117928 | return {}; | |
| 911 | |||
| 912 | // return original curve scaled by factor | ||
| 913 | using std::max; | ||
| 914 |
2/2✓ Branch 0 taken 164768 times.
✓ Branch 1 taken 151020 times.
|
315788 | const Scalar scalFact = max(0.0, (sg - sgr_)/(0.1 - sgr_)); |
| 915 | |||
| 916 | 315788 | return ParkerVanGenuchten3P::krg(ste, *baseLawParamsPtr_) * scalFact; | |
| 917 | } | ||
| 918 | } | ||
| 919 | |||
| 920 | /*! | ||
| 921 | * \brief The relative permeability for a phase. | ||
| 922 | * \param phaseIdx Indicator, The saturation of all phases. | ||
| 923 | * \param swe Effective wetting phase saturation | ||
| 924 | * \param sne Effective nonwetting saturation | ||
| 925 | */ | ||
| 926 | OptionalScalar<Scalar> kr(const int phaseIdx, const Scalar swe, const Scalar sne) const | ||
| 927 | { | ||
| 928 | switch (phaseIdx) | ||
| 929 | { | ||
| 930 | case 0: | ||
| 931 | return krw(swe, sne); | ||
| 932 | case 1: | ||
| 933 | return krn(swe, sne); | ||
| 934 | case 2: | ||
| 935 | return krg(swe, sne); | ||
| 936 | } | ||
| 937 | DUNE_THROW(Dune::InvalidStateException, | ||
| 938 | "Invalid phase index "); | ||
| 939 | } | ||
| 940 | |||
| 941 | private: | ||
| 942 | template<class MaterialLaw> | ||
| 943 | 17 | void initPcParameters_(const MaterialLaw* m, const Scalar lowSwe, const Scalar highSwe) | |
| 944 | { | ||
| 945 | 51 | const auto lowSw = MaterialLaw::EffToAbs::sweToSw(lowSwe, m->effToAbsParams()); | |
| 946 | 51 | const auto highSw = MaterialLaw::EffToAbs::sweToSw(highSwe, m->effToAbsParams()); | |
| 947 | 51 | const auto dsw_dswe = MaterialLaw::EffToAbs::dsw_dswe(m->effToAbsParams()); | |
| 948 | 51 | const auto dst_dste = MaterialLaw::EffToAbs::dst_dste(m->effToAbsParams()); | |
| 949 | |||
| 950 | 34 | baseLawParamsPtr_ = &m->basicParams(); | |
| 951 | |||
| 952 | // pcgw | ||
| 953 | 34 | pcgwLowSwePcgwValue_ = m->template pcgw<false>(lowSw, 0.0); | |
| 954 | 34 | pcgwDerivativeLowSw_ = m->template dpcgw_dsw<false>(lowSw, 0.0)*dsw_dswe; | |
| 955 | 34 | pcgwHighSwePcgwValue_ = m->template pcgw<false>(highSw, 0.0); | |
| 956 | 34 | pcgwDerivativeHighSweThreshold_ = m->template dpcgw_dsw<false>(highSw, 0.0)*dsw_dswe; | |
| 957 | 34 | pcgwDerivativeHighSweEnd_ = 2.0*(0.0 - m->template pcgw<false>(highSw, 0.0))/(1.0 - highSwe); | |
| 958 | 17 | pcgwSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1 | |
| 959 | pcgwHighSwePcgwValue_, 0, // y0, y1 | ||
| 960 | pcgwDerivativeHighSweThreshold_, pcgwDerivativeHighSweEnd_); // m0, m1 | ||
| 961 | |||
| 962 | // pcnw | ||
| 963 | 34 | pcnwLowSwePcnwValue_ = m->template pcnw<false>(lowSw, 0.0); | |
| 964 | 34 | pcnwDerivativeLowSw_ = m->template dpcnw_dsw<false>(lowSw, 0.0)*dsw_dswe; | |
| 965 | 34 | pcnwHighSwePcnwValue_ = m->template pcnw<false>(highSw, 0.0); | |
| 966 | 34 | pcnwDerivativeHighSweThreshold_ = m->template dpcnw_dsw<false>(highSw, 0.0); | |
| 967 | 34 | pcnwDerivativeHighSweEnd_ = 2.0*(0.0 - m->template pcnw<false>(highSw, 0.0))/(1.0 - highSwe); | |
| 968 | 17 | pcnwSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1 | |
| 969 | pcnwHighSwePcnwValue_, 0, // y0, y1 | ||
| 970 | pcnwDerivativeHighSweThreshold_, pcnwDerivativeHighSweEnd_); // m0, m1 | ||
| 971 | |||
| 972 | // pcgn | ||
| 973 | 34 | pcgnLowStePcgnValue_ = m->template pcgn<false>(lowSw, 0.0); | |
| 974 | 34 | pcgnDerivativeLowSt_ = m->template dpcgn_dst<false>(lowSw, 0.0)*dst_dste; | |
| 975 | 34 | pcgnHighSwePcgnValue_ = m->template pcgn<false>(highSw, 0.0); | |
| 976 | 34 | pcgnDerivativeHighSteThreshold_ = m->template dpcgn_dst<false>(highSw, 0.0); | |
| 977 | 34 | pcgnDerivativeHighSteEnd_ = 2.0*(0.0 - m->template pcgn<false>(highSw, 0.0))/(1.0 - highSwe); | |
| 978 | 17 | pcgnSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1 | |
| 979 | pcgnHighSwePcgnValue_, 0, // y0, y1 | ||
| 980 | pcgnDerivativeHighSteThreshold_, pcgnDerivativeHighSteEnd_); // m0, m1 | ||
| 981 | |||
| 982 | 17 | } | |
| 983 | |||
| 984 | template<class MaterialLaw> | ||
| 985 | ✗ | void initKrParameters_(const MaterialLaw* m, const Scalar lowSwe, const Scalar highSwe) | |
| 986 | { | ||
| 987 | ✗ | krgLowStkrgValue_ = ParkerVanGenuchten3P::krg(krgLowSte_, *baseLawParamsPtr_); | |
| 988 | ✗ | krgDerivativeLowSt_ = ParkerVanGenuchten3P::dkrg_dste(krgLowSte_, *baseLawParamsPtr_); | |
| 989 | |||
| 990 | ✗ | swr_ = m->effToAbsParams().swr(); | |
| 991 | ✗ | sgr_ = m->effToAbsParams().sgr(); | |
| 992 | ✗ | } | |
| 993 | |||
| 994 | Scalar krgLowStkrgValue_; | ||
| 995 | Scalar krgDerivativeLowSt_; | ||
| 996 | |||
| 997 | Scalar pcLowSwe_, pcHighSwe_; | ||
| 998 | Scalar krwHighSwe_, krnLowSwe_, krgLowSte_; | ||
| 999 | |||
| 1000 | // pcgw | ||
| 1001 | Scalar pcgwLowSwePcgwValue_; | ||
| 1002 | Scalar pcgwHighSwePcgwValue_; | ||
| 1003 | Scalar pcgwDerivativeLowSw_; | ||
| 1004 | Scalar pcgwDerivativeHighSweThreshold_; | ||
| 1005 | Scalar pcgwDerivativeHighSweEnd_; | ||
| 1006 | |||
| 1007 | // pcgn | ||
| 1008 | Scalar pcgnLowStePcgnValue_; | ||
| 1009 | Scalar pcgnHighSwePcgnValue_; | ||
| 1010 | Scalar pcgnDerivativeLowSt_; | ||
| 1011 | Scalar pcgnDerivativeHighSteThreshold_; | ||
| 1012 | Scalar pcgnDerivativeHighSteEnd_; | ||
| 1013 | |||
| 1014 | // pcnw | ||
| 1015 | Scalar pcnwLowSwePcnwValue_; | ||
| 1016 | Scalar pcnwHighSwePcnwValue_; | ||
| 1017 | Scalar pcnwDerivativeLowSw_; | ||
| 1018 | Scalar pcnwDerivativeHighSweThreshold_; | ||
| 1019 | Scalar pcnwDerivativeHighSweEnd_; | ||
| 1020 | |||
| 1021 | Spline<Scalar> pcgwSpline_; | ||
| 1022 | Spline<Scalar> pcnwSpline_; | ||
| 1023 | Spline<Scalar> pcgnSpline_; | ||
| 1024 | Spline<Scalar> krwSpline_; | ||
| 1025 | Spline<Scalar> krnSpline_; | ||
| 1026 | |||
| 1027 | Scalar swr_, sgr_; | ||
| 1028 | |||
| 1029 | bool constRegularization_; | ||
| 1030 | |||
| 1031 | const BaseLawParams* baseLawParamsPtr_; | ||
| 1032 | }; | ||
| 1033 | |||
| 1034 | /*! | ||
| 1035 | * \ingroup Fluidmatrixinteractions | ||
| 1036 | * \brief Parker van Genuchten material law | ||
| 1037 | */ | ||
| 1038 | template<class ScalarType, | ||
| 1039 | class BaseLaw, | ||
| 1040 | class Regularization = NoRegularization, | ||
| 1041 | class EffToAbsPolicy = ParkerVanGenuchten3PEffToAbsPolicy> | ||
| 1042 | class ParkerVanGenuchtenMaterialLaw : public Adapter<ParkerVanGenuchtenMaterialLaw<ScalarType, BaseLaw, Regularization, EffToAbsPolicy>, ThreePhasePcKrSw> | ||
| 1043 | { | ||
| 1044 | public: | ||
| 1045 | |||
| 1046 | using Scalar = ScalarType; | ||
| 1047 | |||
| 1048 | using BasicParams = typename BaseLaw::template Params<Scalar>; | ||
| 1049 | using EffToAbsParams = typename EffToAbsPolicy::template Params<Scalar>; | ||
| 1050 | using RegularizationParams = typename Regularization::template Params<Scalar>; | ||
| 1051 | |||
| 1052 | using EffToAbs = EffToAbsPolicy; | ||
| 1053 | |||
| 1054 | /*! | ||
| 1055 | * \brief Return whether this law is regularized | ||
| 1056 | */ | ||
| 1057 | static constexpr bool isRegularized() | ||
| 1058 | { return !std::is_same<Regularization, NoRegularization>::value; } | ||
| 1059 | |||
| 1060 | /*! | ||
| 1061 | * \brief Deleted default constructor (so we are never in an undefined state) | ||
| 1062 | * \note store owning pointers to laws instead if you need default-constructible objects | ||
| 1063 | */ | ||
| 1064 | ParkerVanGenuchtenMaterialLaw() = delete; | ||
| 1065 | |||
| 1066 | /*! | ||
| 1067 | * \brief Construct from a subgroup from the global parameter tree | ||
| 1068 | * \note This will give you nice error messages if a mandatory parameter is missing | ||
| 1069 | */ | ||
| 1070 | 17 | explicit ParkerVanGenuchtenMaterialLaw(const std::string& paramGroup) | |
| 1071 | 34 | : basicParams_(makeBasicParams(paramGroup)) | |
| 1072 | 17 | , effToAbsParams_(makeEffToAbsParams(paramGroup)) | |
| 1073 | { | ||
| 1074 | if constexpr (isRegularized()) | ||
| 1075 | 17 | regularization_.init(this, paramGroup); | |
| 1076 | 17 | } | |
| 1077 | |||
| 1078 | /*! | ||
| 1079 | * \brief Construct from parameter structs | ||
| 1080 | * \note More efficient constructor but you need to ensure all parameters are initialized | ||
| 1081 | */ | ||
| 1082 | ParkerVanGenuchtenMaterialLaw(const BasicParams& baseParams, | ||
| 1083 | const EffToAbsParams& effToAbsParams = {}, | ||
| 1084 | const RegularizationParams& regParams = {}) | ||
| 1085 | : basicParams_(baseParams) | ||
| 1086 | , effToAbsParams_(effToAbsParams) | ||
| 1087 | { | ||
| 1088 | if constexpr (isRegularized()) | ||
| 1089 | regularization_.init(this, baseParams, effToAbsParams, regParams); | ||
| 1090 | } | ||
| 1091 | |||
| 1092 | /*! | ||
| 1093 | * \brief The capillary pressure-saturation curve for the gas and wetting phase | ||
| 1094 | */ | ||
| 1095 | template<bool enableRegularization = isRegularized()> | ||
| 1096 | ✗ | Scalar pcgw(const Scalar sw, const Scalar /*dummySn*/) const | |
| 1097 | { | ||
| 1098 | 102 | const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_); | |
| 1099 | if constexpr (enableRegularization) | ||
| 1100 | { | ||
| 1101 | ✗ | const auto regularized = regularization_.pcgw(swe); | |
| 1102 | ✗ | if (regularized) | |
| 1103 | return regularized.value(); | ||
| 1104 | } | ||
| 1105 | |||
| 1106 | 51 | return BaseLaw::pcgw(swe, basicParams_); | |
| 1107 | } | ||
| 1108 | |||
| 1109 | /*! | ||
| 1110 | * \brief The capillary pressure-saturation curve for the nonwetting and wetting phase | ||
| 1111 | */ | ||
| 1112 | template<bool enableRegularization = isRegularized()> | ||
| 1113 | ✗ | Scalar pcnw(const Scalar sw, const Scalar /*dummySn*/) const | |
| 1114 | { | ||
| 1115 | 102 | const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_); | |
| 1116 | if constexpr (enableRegularization) | ||
| 1117 | { | ||
| 1118 | ✗ | const auto regularized = regularization_.pcnw(swe); | |
| 1119 | ✗ | if (regularized) | |
| 1120 | return regularized.value(); | ||
| 1121 | } | ||
| 1122 | |||
| 1123 | 17 | return BaseLaw::pcnw(swe, basicParams_); | |
| 1124 | } | ||
| 1125 | |||
| 1126 | /*! | ||
| 1127 | * \brief The capillary pressure-saturation curve for the gas and nonwetting phase | ||
| 1128 | * \param sw Wetting saturation | ||
| 1129 | * \param sn Nonwetting saturation | ||
| 1130 | */ | ||
| 1131 | template<bool enableRegularization = isRegularized()> | ||
| 1132 | 9494993 | Scalar pcgn(const Scalar sw, const Scalar sn) const | |
| 1133 | { | ||
| 1134 | 18990088 | const auto swe = EffToAbs::swToSwe(sw + sn, effToAbsParams_); | |
| 1135 | if constexpr (enableRegularization) | ||
| 1136 | { | ||
| 1137 | 9494993 | const auto regularized = regularization_.pcgn(swe); | |
| 1138 |
2/2✓ Branch 0 taken 1332506 times.
✓ Branch 1 taken 8162487 times.
|
9494993 | if (regularized) |
| 1139 | return regularized.value(); | ||
| 1140 | } | ||
| 1141 | |||
| 1142 | 1332523 | return BaseLaw::pcgn(swe, basicParams_); | |
| 1143 | } | ||
| 1144 | |||
| 1145 | /*! | ||
| 1146 | * \brief This function ensures a continuous transition from 2 to 3 phases and vice versa | ||
| 1147 | */ | ||
| 1148 | template<bool enableRegularization = isRegularized()> | ||
| 1149 | ✗ | Scalar pcAlpha(const Scalar /*dummySw*/, const Scalar sn) const | |
| 1150 | { | ||
| 1151 | 9494993 | const auto sne = EffToAbs::snToSne(sn, effToAbsParams_); | |
| 1152 | if constexpr (enableRegularization) | ||
| 1153 | { | ||
| 1154 | 18989986 | const auto regularized = regularization_.pcAlpha(sne); | |
| 1155 |
2/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4607687 times.
✓ Branch 5 taken 4887306 times.
|
9494993 | if (regularized) |
| 1156 | return regularized.value(); | ||
| 1157 | } | ||
| 1158 |
4/8✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4607687 times.
✓ Branch 5 taken 4887306 times.
✓ Branch 6 taken 383 times.
✓ Branch 7 taken 3854565 times.
|
18989986 | return BaseLaw::pcAlpha(sne, basicParams_); |
| 1159 | } | ||
| 1160 | |||
| 1161 | /*! | ||
| 1162 | * \brief The partial derivative of the capillary pressure w.r.t. the saturation | ||
| 1163 | */ | ||
| 1164 | template<bool enableRegularization = isRegularized()> | ||
| 1165 | ✗ | Scalar dpcgw_dsw(const Scalar sw, const Scalar /*dummyS*/) const | |
| 1166 | { | ||
| 1167 | 68 | const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_); | |
| 1168 | if constexpr (enableRegularization) | ||
| 1169 | { | ||
| 1170 | const auto regularized = regularization_.dpcgw_dswe(swe); | ||
| 1171 | if (regularized) | ||
| 1172 | return regularized.value()*EffToAbs::dswe_dsw(effToAbsParams_); | ||
| 1173 | } | ||
| 1174 | |||
| 1175 | 17 | return BaseLaw::dpcgw_dswe(swe, basicParams_)*EffToAbs::dswe_dsw(effToAbsParams_); | |
| 1176 | } | ||
| 1177 | |||
| 1178 | /*! | ||
| 1179 | * \brief The partial derivative of the capillary pressure w.r.t. the saturation | ||
| 1180 | */ | ||
| 1181 | template<bool enableRegularization = isRegularized()> | ||
| 1182 | ✗ | Scalar dpcnw_dsw(const Scalar sw, const Scalar /*dummySw*/) const | |
| 1183 | { | ||
| 1184 | 68 | const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_); | |
| 1185 | if constexpr (enableRegularization) | ||
| 1186 | { | ||
| 1187 | const auto regularized = regularization_.dpcnw_dswe(swe); | ||
| 1188 | if (regularized) | ||
| 1189 | return regularized.value()*EffToAbs::dswe_dsw(effToAbsParams_); | ||
| 1190 | } | ||
| 1191 | |||
| 1192 | 17 | return BaseLaw::dpcnw_dswe(swe, basicParams_)*EffToAbs::dswe_dsw(effToAbsParams_); | |
| 1193 | } | ||
| 1194 | |||
| 1195 | /*! | ||
| 1196 | * \brief The partial derivative of the capillary pressure w.r.t. the saturation | ||
| 1197 | */ | ||
| 1198 | template<bool enableRegularization = isRegularized()> | ||
| 1199 | ✗ | Scalar dpcgn_dst(const Scalar st, const Scalar /*dummySw*/) const | |
| 1200 | { | ||
| 1201 | 68 | const auto ste = EffToAbs::stToSte(st, effToAbsParams_); | |
| 1202 | if constexpr (enableRegularization) | ||
| 1203 | { | ||
| 1204 | const auto regularized = regularization_.dpcgn_dste(ste); | ||
| 1205 | if (regularized) | ||
| 1206 | return regularized.value()*EffToAbs::dswte_dst(effToAbsParams_); | ||
| 1207 | } | ||
| 1208 | |||
| 1209 | 17 | return BaseLaw::dpcgn_dste(ste, basicParams_)*EffToAbs::dste_dst(effToAbsParams_); | |
| 1210 | } | ||
| 1211 | |||
| 1212 | /*! | ||
| 1213 | * \brief The relative permeability for the wetting phase | ||
| 1214 | * \param sw Wetting saturation | ||
| 1215 | * \param sn Nonwetting saturation | ||
| 1216 | */ | ||
| 1217 | template<bool enableRegularization = isRegularized()> | ||
| 1218 | ✗ | Scalar krw(const Scalar sw, const Scalar sn) const | |
| 1219 | { | ||
| 1220 | ✗ | const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_); | |
| 1221 | if constexpr (enableRegularization) | ||
| 1222 | { | ||
| 1223 | ✗ | const auto regularized = regularization_.krw(swe); | |
| 1224 | ✗ | if (regularized) | |
| 1225 | return regularized.value(); | ||
| 1226 | } | ||
| 1227 | |||
| 1228 | ✗ | return BaseLaw::krw(swe, basicParams_); | |
| 1229 | } | ||
| 1230 | |||
| 1231 | /*! | ||
| 1232 | * \brief The relative permeability for the nonwetting phase | ||
| 1233 | * \param sw Wetting saturation | ||
| 1234 | * \param sn Nonwetting saturation | ||
| 1235 | */ | ||
| 1236 | template<bool enableRegularization = isRegularized()> | ||
| 1237 | 9496991 | Scalar krn(const Scalar sw, const Scalar sn) const | |
| 1238 | { | ||
| 1239 | 18993982 | const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_); | |
| 1240 | 18993982 | const auto ste = EffToAbs::stToSte(sw + sn, effToAbsParams_); | |
| 1241 | if constexpr (enableRegularization) | ||
| 1242 | { | ||
| 1243 | 9496991 | const auto regularized = regularization_.krn(swe, sn, ste); | |
| 1244 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9496991 times.
|
9496991 | if (regularized) |
| 1245 | return regularized.value(); | ||
| 1246 | } | ||
| 1247 | |||
| 1248 | ✗ | return BaseLaw::krn(swe, sn, ste, basicParams_); | |
| 1249 | } | ||
| 1250 | |||
| 1251 | /*! | ||
| 1252 | * \brief The relative permeability for the nonwetting phase | ||
| 1253 | * \param sw Wetting saturation | ||
| 1254 | * \param sn Nonwetting saturation | ||
| 1255 | */ | ||
| 1256 | template<bool enableRegularization = isRegularized()> | ||
| 1257 | 9496991 | Scalar krg(const Scalar sw, const Scalar sn) const | |
| 1258 | { | ||
| 1259 | 18993982 | const auto ste = EffToAbs::stToSte(sw + sn, effToAbsParams_); | |
| 1260 | if constexpr (enableRegularization) | ||
| 1261 | { | ||
| 1262 | 9496991 | const auto regularized = regularization_.krg(ste); | |
| 1263 |
2/2✓ Branch 0 taken 1117928 times.
✓ Branch 1 taken 8379063 times.
|
9496991 | if (regularized) |
| 1264 | return regularized.value(); | ||
| 1265 | } | ||
| 1266 | |||
| 1267 | 1117928 | return BaseLaw::krg(ste, basicParams_); | |
| 1268 | } | ||
| 1269 | |||
| 1270 | /*! | ||
| 1271 | * \brief The relative permeability for the nonwetting phase | ||
| 1272 | * \param phaseIdx Indicator, The saturation of all phases. | ||
| 1273 | * \param sw Wetting saturation | ||
| 1274 | * \param sn Nonwetting saturation | ||
| 1275 | */ | ||
| 1276 | template<bool enableRegularization = isRegularized()> | ||
| 1277 | 28484973 | Scalar kr(const int phaseIdx, const Scalar sw, const Scalar sn) const | |
| 1278 | { | ||
| 1279 |
3/4✓ Branch 0 taken 9494991 times.
✓ Branch 1 taken 9494991 times.
✓ Branch 2 taken 9494991 times.
✗ Branch 3 not taken.
|
28484973 | switch (phaseIdx) |
| 1280 | { | ||
| 1281 | 9494991 | case 0: | |
| 1282 | 9494991 | return krw(sw, sn); | |
| 1283 | 9494991 | case 1: | |
| 1284 | 9494991 | return krn(sw, sn); | |
| 1285 | 9494991 | case 2: | |
| 1286 | 9494991 | return krg(sw, sn); | |
| 1287 | } | ||
| 1288 | ✗ | DUNE_THROW(Dune::InvalidStateException, | |
| 1289 | "Invalid phase index "); | ||
| 1290 | } | ||
| 1291 | |||
| 1292 | /*! | ||
| 1293 | * \brief The derivative of the relative permeability for the nonwetting phase w.r.t. saturation | ||
| 1294 | * \param st Total (wetting + nonwetting) saturation | ||
| 1295 | */ | ||
| 1296 | template<bool enableRegularization = isRegularized()> | ||
| 1297 | Scalar dkrg_dst(const Scalar st) const | ||
| 1298 | { | ||
| 1299 | const auto ste = EffToAbs::stToSte(st, effToAbsParams_); | ||
| 1300 | if constexpr (enableRegularization) | ||
| 1301 | { | ||
| 1302 | const auto regularized = regularization_.dkrg_dste(ste); | ||
| 1303 | if (regularized) | ||
| 1304 | return regularized.value()*EffToAbs::dste_dst(effToAbsParams_); | ||
| 1305 | } | ||
| 1306 | |||
| 1307 | return BaseLaw::dkrg_dste(ste, basicParams_)*EffToAbs::dste_dst(effToAbsParams_); | ||
| 1308 | } | ||
| 1309 | |||
| 1310 | /*! | ||
| 1311 | * \brief Equality comparison with another instance | ||
| 1312 | */ | ||
| 1313 | bool operator== (const ParkerVanGenuchtenMaterialLaw& o) const | ||
| 1314 | { | ||
| 1315 | return basicParams_ == o.basicParams_ | ||
| 1316 | && effToAbsParams_ == o.effToAbsParams_ | ||
| 1317 | && regularization_ == o.regularization_; | ||
| 1318 | } | ||
| 1319 | |||
| 1320 | /*! | ||
| 1321 | * \brief Create the base law's parameters using | ||
| 1322 | * input file parameters | ||
| 1323 | */ | ||
| 1324 | static BasicParams makeBasicParams(const std::string& paramGroup) | ||
| 1325 | { | ||
| 1326 | 17 | return BaseLaw::template makeParams<Scalar>(paramGroup); | |
| 1327 | } | ||
| 1328 | |||
| 1329 | /*! | ||
| 1330 | * \brief Return the base law's parameters | ||
| 1331 | */ | ||
| 1332 | const BasicParams& basicParams() const | ||
| 1333 | 17 | { return basicParams_; } | |
| 1334 | |||
| 1335 | /*! | ||
| 1336 | * \brief Create the parameters of the EffToAbs policy using | ||
| 1337 | * input file parameters | ||
| 1338 | */ | ||
| 1339 | static EffToAbsParams makeEffToAbsParams(const std::string& paramGroup) | ||
| 1340 | { | ||
| 1341 | 17 | return EffToAbs::template makeParams<Scalar>(paramGroup); | |
| 1342 | } | ||
| 1343 | |||
| 1344 | /*! | ||
| 1345 | * \brief Return the parameters of the EffToAbs policy | ||
| 1346 | */ | ||
| 1347 | const EffToAbsParams& effToAbsParams() const | ||
| 1348 | 17 | { return effToAbsParams_; } | |
| 1349 | |||
| 1350 | private: | ||
| 1351 | BasicParams basicParams_; | ||
| 1352 | EffToAbsParams effToAbsParams_; | ||
| 1353 | Regularization regularization_; | ||
| 1354 | }; | ||
| 1355 | |||
| 1356 | /*! | ||
| 1357 | * \ingroup Fluidmatrixinteractions | ||
| 1358 | * \brief A configuration for using the ParkerVanGenuchten material law without regularization | ||
| 1359 | */ | ||
| 1360 | template<class Scalar> | ||
| 1361 | using ParkerVanGenuchten3PNoReg = ParkerVanGenuchtenMaterialLaw<Scalar, ParkerVanGenuchten3P, NoRegularization, ParkerVanGenuchten3PEffToAbsPolicy>; | ||
| 1362 | |||
| 1363 | /*! | ||
| 1364 | * \ingroup Fluidmatrixinteractions | ||
| 1365 | * \brief A default configuration for using the ParkerVanGenuchten material law | ||
| 1366 | */ | ||
| 1367 | template<class Scalar> | ||
| 1368 | using ParkerVanGenuchten3PDefault = ParkerVanGenuchtenMaterialLaw<Scalar, ParkerVanGenuchten3P, ParkerVanGenuchten3PRegularization<Scalar>, ParkerVanGenuchten3PEffToAbsPolicy>; | ||
| 1369 | |||
| 1370 | } // end namespace Dumux::FluidMatrix | ||
| 1371 | |||
| 1372 | #endif | ||
| 1373 |