| 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 the capillary pressure and | ||
| 11 | * relative permeability <-> saturation relations according to Brooks and Corey. | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_MATERIAL_FLUIDMATRIX_BROOKS_COREY_HH | ||
| 14 | #define DUMUX_MATERIAL_FLUIDMATRIX_BROOKS_COREY_HH | ||
| 15 | |||
| 16 | #include <cmath> | ||
| 17 | #include <algorithm> | ||
| 18 | |||
| 19 | #include <dumux/common/parameters.hh> | ||
| 20 | #include <dumux/common/optionalscalar.hh> | ||
| 21 | #include <dumux/material/fluidmatrixinteractions/2p/materiallaw.hh> | ||
| 22 | |||
| 23 | namespace Dumux::FluidMatrix { | ||
| 24 | |||
| 25 | /*! | ||
| 26 | * \ingroup Fluidmatrixinteractions | ||
| 27 | * | ||
| 28 | * \brief Implementation of the Brooks-Corey capillary pressure <-> | ||
| 29 | * saturation relation. This class bundles the "raw" curves | ||
| 30 | * as static members and doesn't concern itself converting | ||
| 31 | * absolute to effective saturations and vice versa. | ||
| 32 | * | ||
| 33 | * For general info: EffToAbsLaw | ||
| 34 | * | ||
| 35 | *\see BrooksCoreyParams | ||
| 36 | */ | ||
| 37 | class BrooksCorey | ||
| 38 | { | ||
| 39 | public: | ||
| 40 | /*! | ||
| 41 | * \brief The parameter type | ||
| 42 | * \tparam Scalar The scalar type | ||
| 43 | * \note The Brooks Corey laws are parameterized with two parameters: \f$\mathrm{p_{ce}, \lambda}\f$, | ||
| 44 | * the capillary entry pressure in \f$\mathrm{[Pa]}\f$] and a dimensionless shape parameter, respectively. | ||
| 45 | */ | ||
| 46 | template<class Scalar> | ||
| 47 | struct Params | ||
| 48 | { | ||
| 49 | 7 | Params(Scalar pcEntry, Scalar lambda) | |
| 50 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
7 | : pcEntry_(pcEntry), lambda_(lambda) |
| 51 | {} | ||
| 52 | |||
| 53 | ✗ | Scalar pcEntry() const{ return pcEntry_; } | |
| 54 | void setPcEntry(Scalar pcEntry){ pcEntry_ = pcEntry; } | ||
| 55 | |||
| 56 | ✗ | Scalar lambda() const { return lambda_; } | |
| 57 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | void setLambda(Scalar lambda) { lambda_ = lambda; } |
| 58 | |||
| 59 | 29504 | bool operator== (const Params& p) const | |
| 60 | { | ||
| 61 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 29504 times.
|
29504 | return Dune::FloatCmp::eq(pcEntry(), p.pcEntry(), 1e-6) |
| 62 |
4/6✓ Branch 0 taken 26340 times.
✓ Branch 1 taken 3164 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26340 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 26340 times.
|
29504 | && Dune::FloatCmp::eq(lambda(), p.lambda(), 1e-6); |
| 63 | } | ||
| 64 | |||
| 65 | private: | ||
| 66 | Scalar pcEntry_, lambda_; | ||
| 67 | }; | ||
| 68 | |||
| 69 | /*! | ||
| 70 | * \brief Construct from a subgroup from the global parameter tree | ||
| 71 | * \note This will give you nice error messages if a mandatory parameter is missing | ||
| 72 | */ | ||
| 73 | template<class Scalar = double> | ||
| 74 | 69 | static Params<Scalar> makeParams(const std::string& paramGroup) | |
| 75 | { | ||
| 76 | 69 | const auto pcEntry = getParamFromGroup<Scalar>(paramGroup, "BrooksCoreyPcEntry"); | |
| 77 | 69 | const auto lambda = getParamFromGroup<Scalar>(paramGroup, "BrooksCoreyLambda"); | |
| 78 | 69 | return {pcEntry, lambda}; | |
| 79 | } | ||
| 80 | |||
| 81 | /*! | ||
| 82 | * \brief The capillary pressure-saturation curve according to Brooks & Corey. | ||
| 83 | * | ||
| 84 | * The Brooks-Corey empirical capillary pressure <-> saturation | ||
| 85 | * function is given by | ||
| 86 | * | ||
| 87 | * \f$\mathrm{ p_c = p_{ce}\overline{S}_w^{-1/\lambda} | ||
| 88 | * }\f$ | ||
| 89 | * | ||
| 90 | * \param swe Effective saturation of the wetting phase \f$\mathrm{[\overline{S}_w]}\f$ | ||
| 91 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 92 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 93 | and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 94 | * \return Capillary pressure calculated by Brooks & Corey constitutive relation. | ||
| 95 | * | ||
| 96 | * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, | ||
| 97 | * by clamping the input. | ||
| 98 | */ | ||
| 99 | template<class Scalar> | ||
| 100 | 25110822 | static Scalar pc(Scalar swe, const Params<Scalar>& params) | |
| 101 | { | ||
| 102 | using std::pow; | ||
| 103 | using std::clamp; | ||
| 104 | |||
| 105 |
1/2✓ Branch 0 taken 25110822 times.
✗ Branch 1 not taken.
|
25110822 | swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0 |
| 106 | |||
| 107 | 25110822 | return params.pcEntry()*pow(swe, -1.0/params.lambda()); | |
| 108 | } | ||
| 109 | |||
| 110 | /*! | ||
| 111 | * \brief The saturation-capillary pressure curve according to Brooks & Corey. | ||
| 112 | * | ||
| 113 | * This is the inverse of the capillary pressure-saturation curve: | ||
| 114 | * \f$\mathrm{ \overline{S}_w = (\frac{p_c}{p_{ce}})^{-\lambda}}\f$ | ||
| 115 | * | ||
| 116 | * \param pc Capillary pressure \f$\mathrm{[p_c]}\f$ in \f$\mathrm{[Pa]}\f$. | ||
| 117 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 118 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container | ||
| 119 | * is constructed accordingly. Afterwards the values are set there, too. | ||
| 120 | * \return Effective wetting phase saturation calculated as inverse of BrooksCorey constitutive relation. | ||
| 121 | * | ||
| 122 | * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, | ||
| 123 | * by clamping the input. | ||
| 124 | */ | ||
| 125 | template<class Scalar> | ||
| 126 | static Scalar swe(Scalar pc, const Params<Scalar>& params) | ||
| 127 | { | ||
| 128 | using std::pow; | ||
| 129 | using std::max; | ||
| 130 | |||
| 131 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1335568 times.
|
1335568 | pc = max(pc, 0.0); // the equation below is undefined for negative pcs |
| 132 | |||
| 133 | 1335568 | return pow(pc/params.pcEntry(), -params.lambda()); | |
| 134 | } | ||
| 135 | |||
| 136 | /*! | ||
| 137 | * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure | ||
| 138 | * | ||
| 139 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 140 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container | ||
| 141 | * is constructed accordingly. Afterwards the values are set there, too. | ||
| 142 | */ | ||
| 143 | template<class Scalar> | ||
| 144 | static Scalar endPointPc(const Params<Scalar>& params) | ||
| 145 |
6/6✓ Branch 0 taken 3392238 times.
✓ Branch 1 taken 2189101 times.
✓ Branch 2 taken 360 times.
✓ Branch 3 taken 29505 times.
✓ Branch 4 taken 360 times.
✓ Branch 5 taken 29504 times.
|
5676360 | { return params.pcEntry(); } |
| 146 | |||
| 147 | /*! | ||
| 148 | * \brief The partial derivative of the capillary | ||
| 149 | * pressure w.r.t. the effective saturation according to Brooks & Corey. | ||
| 150 | * | ||
| 151 | * This is equivalent to | ||
| 152 | * \f$\mathrm{\frac{\partial p_c}{\partial \overline{S}_w} = | ||
| 153 | * -\frac{p_{ce}}{\lambda} \overline{S}_w^{-1/\lambda - 1} | ||
| 154 | * }\f$ | ||
| 155 | * | ||
| 156 | * \param swe Effective saturation of the wetting phase \f$\mathrm{[\overline{S}_w]}\f$ | ||
| 157 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 158 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container | ||
| 159 | * is constructed accordingly. Afterwards the values are set there, too. | ||
| 160 | * \return Partial derivative of \f$\mathrm{[p_c]}\f$ w.r.t. effective saturation according to Brooks & Corey. | ||
| 161 | * | ||
| 162 | * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, | ||
| 163 | * by clamping the input. | ||
| 164 | */ | ||
| 165 | template<class Scalar> | ||
| 166 | 314 | static Scalar dpc_dswe(Scalar swe, const Params<Scalar>& params) | |
| 167 | { | ||
| 168 | using std::pow; | ||
| 169 | using std::clamp; | ||
| 170 | |||
| 171 |
1/2✓ Branch 0 taken 314 times.
✗ Branch 1 not taken.
|
314 | swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0 |
| 172 | |||
| 173 | 314 | return -params.pcEntry()/params.lambda() * pow(swe, -1.0/params.lambda() - 1.0); | |
| 174 | } | ||
| 175 | |||
| 176 | /*! | ||
| 177 | * \brief The partial derivative of the effective | ||
| 178 | * saturation w.r.t. the capillary pressure according to Brooks & Corey. | ||
| 179 | * | ||
| 180 | * \param pc Capillary pressure \f$\mathrm{[p_c]}\f$ in \f$\mathrm{[Pa]}\f$. | ||
| 181 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 182 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container | ||
| 183 | * is constructed accordingly. Afterwards the values are set there, too. | ||
| 184 | * \return Partial derivative of effective saturation w.r.t. \f$\mathrm{[p_c]}\f$ according to Brooks & Corey. | ||
| 185 | * | ||
| 186 | * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, | ||
| 187 | * by clamping the input. | ||
| 188 | */ | ||
| 189 | template<class Scalar> | ||
| 190 | 158 | static Scalar dswe_dpc(Scalar pc, const Params<Scalar>& params) | |
| 191 | { | ||
| 192 | using std::pow; | ||
| 193 | using std::max; | ||
| 194 | |||
| 195 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 158 times.
|
158 | pc = max(pc, 0.0); // the equation below is undefined for negative pcs |
| 196 | |||
| 197 | 158 | return -params.lambda()/params.pcEntry() * pow(pc/params.pcEntry(), - params.lambda() - 1.0); | |
| 198 | } | ||
| 199 | |||
| 200 | /*! | ||
| 201 | * \brief The relative permeability for the wetting phase of | ||
| 202 | * the medium implied by the Brooks-Corey | ||
| 203 | * parameterization. | ||
| 204 | * | ||
| 205 | * \param swe The mobile saturation of the wetting phase. | ||
| 206 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 207 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 208 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 209 | * \return Relative permeability of the wetting phase calculated as implied by Brooks & Corey. | ||
| 210 | * | ||
| 211 | * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, | ||
| 212 | * by clamping the input. | ||
| 213 | */ | ||
| 214 | template<class Scalar> | ||
| 215 | 22145587 | static Scalar krw(Scalar swe, const Params<Scalar>& params) | |
| 216 | { | ||
| 217 | using std::pow; | ||
| 218 | using std::clamp; | ||
| 219 | |||
| 220 |
1/2✓ Branch 0 taken 22145587 times.
✗ Branch 1 not taken.
|
22145587 | swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0 |
| 221 | |||
| 222 | 22145587 | return pow(swe, 2.0/params.lambda() + 3.0); | |
| 223 | } | ||
| 224 | |||
| 225 | /*! | ||
| 226 | * \brief The derivative of the relative permeability for the | ||
| 227 | * wetting phase with regard to the wetting saturation of the | ||
| 228 | * medium implied by the Brooks-Corey parameterization. | ||
| 229 | * | ||
| 230 | * \param swe The mobile saturation of the wetting phase. | ||
| 231 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 232 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 233 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 234 | * \return Derivative of the relative permeability of the wetting phase w.r.t. effective wetting phase | ||
| 235 | * saturation calculated as implied by Brooks & Corey. | ||
| 236 | * | ||
| 237 | * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, | ||
| 238 | * by clamping the input. | ||
| 239 | */ | ||
| 240 | template<class Scalar> | ||
| 241 | 80 | static Scalar dkrw_dswe(Scalar swe, const Params<Scalar>& params) | |
| 242 | { | ||
| 243 | using std::pow; | ||
| 244 | using std::clamp; | ||
| 245 | |||
| 246 |
1/2✓ Branch 0 taken 80 times.
✗ Branch 1 not taken.
|
80 | swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0 |
| 247 | |||
| 248 | 80 | return (2.0/params.lambda() + 3.0)*pow(swe, 2.0/params.lambda() + 2.0); | |
| 249 | } | ||
| 250 | |||
| 251 | /*! | ||
| 252 | * \brief The relative permeability for the non-wetting phase of | ||
| 253 | * the medium as implied by the Brooks-Corey | ||
| 254 | * parameterization. | ||
| 255 | * | ||
| 256 | * \param swe The mobile saturation of the wetting phase. | ||
| 257 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 258 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container | ||
| 259 | * is constructed accordingly. Afterwards the values are set there, too. | ||
| 260 | * \return Relative permeability of the non-wetting phase calculated as implied by Brooks & Corey. | ||
| 261 | * | ||
| 262 | * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, | ||
| 263 | * by clamping the input. | ||
| 264 | */ | ||
| 265 | template<class Scalar> | ||
| 266 | 22145587 | static Scalar krn(Scalar swe, const Params<Scalar>& params) | |
| 267 | { | ||
| 268 | using std::pow; | ||
| 269 | using std::clamp; | ||
| 270 | |||
| 271 |
1/2✓ Branch 0 taken 22145587 times.
✗ Branch 1 not taken.
|
22145587 | swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0 |
| 272 | |||
| 273 | 22145587 | const Scalar exponent = 2.0/params.lambda() + 1.0; | |
| 274 | 22145587 | const Scalar sne = 1.0 - swe; | |
| 275 | 22145587 | return sne*sne*(1.0 - pow(swe, exponent)); | |
| 276 | } | ||
| 277 | |||
| 278 | /*! | ||
| 279 | * \brief The derivative of the relative permeability for the | ||
| 280 | * non-wetting phase in regard to the wetting saturation of | ||
| 281 | * the medium as implied by the Brooks-Corey | ||
| 282 | * parameterization. | ||
| 283 | * | ||
| 284 | * \param swe The mobile saturation of the wetting phase. | ||
| 285 | * \param params A container object that is populated with the appropriate coefficients for the respective law. | ||
| 286 | * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, | ||
| 287 | * and then the params container is constructed accordingly. Afterwards the values are set there, too. | ||
| 288 | * \return Derivative of the relative permeability of the non-wetting phase w.r.t. effective wetting phase | ||
| 289 | * saturation calculated as implied by Brooks & Corey. | ||
| 290 | * | ||
| 291 | * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, | ||
| 292 | * by clamping the input. | ||
| 293 | */ | ||
| 294 | template<class Scalar> | ||
| 295 | 80 | static Scalar dkrn_dswe(Scalar swe, const Params<Scalar>& params) | |
| 296 | { | ||
| 297 | using std::pow; | ||
| 298 | using std::clamp; | ||
| 299 | |||
| 300 |
1/2✓ Branch 0 taken 80 times.
✗ Branch 1 not taken.
|
80 | swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0 |
| 301 | |||
| 302 | 80 | const auto lambdaInv = 1.0/params.lambda(); | |
| 303 | 80 | const auto swePow = pow(swe, 2*lambdaInv); | |
| 304 | 80 | return 2.0*(swe - 1.0)*(1.0 + (0.5 + lambdaInv)*swePow - (1.5 + lambdaInv)*swePow*swe); | |
| 305 | } | ||
| 306 | }; | ||
| 307 | |||
| 308 | /*! | ||
| 309 | * \ingroup Fluidmatrixinteractions | ||
| 310 | * \brief A regularization for the BrooksCorey material law | ||
| 311 | * \note Regularization can be turned of by setting the threshold parameters | ||
| 312 | * out of range (runtime) or by replacing | ||
| 313 | * this class by NoRegularization (compile time). | ||
| 314 | */ | ||
| 315 | template <class Scalar> | ||
| 316 | class BrooksCoreyRegularization | ||
| 317 | { | ||
| 318 | public: | ||
| 319 | //! Regularization parameters | ||
| 320 | template<class S> | ||
| 321 | struct Params | ||
| 322 | { | ||
| 323 | 1 | Params(S pcLowSwe = 0.01) : pcLowSwe_(pcLowSwe) {} | |
| 324 | |||
| 325 | ✗ | S pcLowSwe() const { return pcLowSwe_; } | |
| 326 | 1 | void setPcLowSwe(S pcLowSwe) { pcLowSwe_ = pcLowSwe; } | |
| 327 | |||
| 328 | private: | ||
| 329 | S pcLowSwe_ = 0.01; | ||
| 330 | }; | ||
| 331 | |||
| 332 | template<class MaterialLaw> | ||
| 333 | 67 | void init(const MaterialLaw* m, const std::string& paramGroup) | |
| 334 | { | ||
| 335 | 67 | pcLowSwe_ = getParamFromGroup<Scalar>(paramGroup, "BrooksCoreyPcLowSweThreshold", 0.01); | |
| 336 | 67 | entryPressure_ = getParamFromGroup<Scalar>(paramGroup, "BrooksCoreyPcEntry"); | |
| 337 | |||
| 338 | 67 | initPcParameters_(m, pcLowSwe_); | |
| 339 | 67 | } | |
| 340 | |||
| 341 | template<class MaterialLaw, class BaseParams, class EffToAbsParams> | ||
| 342 | ✗ | void init(const MaterialLaw* m, const BaseParams& bp, const EffToAbsParams& etap, const Params<Scalar>& p) | |
| 343 | { | ||
| 344 | 7 | pcLowSwe_ = p.pcLowSwe(); | |
| 345 | 7 | entryPressure_ = bp.pcEntry(); | |
| 346 | |||
| 347 | 7 | initPcParameters_(m, pcLowSwe_); | |
| 348 | ✗ | } | |
| 349 | |||
| 350 | /*! | ||
| 351 | * \brief Equality comparison with another instance | ||
| 352 | */ | ||
| 353 | 26340 | bool operator== (const BrooksCoreyRegularization& o) const | |
| 354 | { | ||
| 355 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 26340 times.
|
26340 | return Dune::FloatCmp::eq(pcLowSwe_, o.pcLowSwe_, 1e-6) |
| 356 |
3/6✓ Branch 0 taken 26340 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 26340 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 26340 times.
|
26340 | && Dune::FloatCmp::eq(entryPressure_, o.entryPressure_, 1e-6); |
| 357 | } | ||
| 358 | |||
| 359 | /*! | ||
| 360 | * \brief The regularized capillary pressure-saturation curve | ||
| 361 | * regularized part: | ||
| 362 | * - low saturation: extend the \f$\mathrm{p_c(S_w)}\f$ curve with the slope at the regularization point (i.e. no kink). | ||
| 363 | * - high saturation: continue linearly | ||
| 364 | */ | ||
| 365 | OptionalScalar<Scalar> pc(const Scalar swe) const | ||
| 366 | { | ||
| 367 | // make sure that the capillary pressure observes a derivative | ||
| 368 | // != 0 for 'illegal' saturations. This is favourable for the | ||
| 369 | // newton solver (if the derivative is calculated numerically) | ||
| 370 | // in order to get the saturation moving to the right | ||
| 371 | // direction if it temporarily is in an 'illegal' range. | ||
| 372 |
2/2✓ Branch 0 taken 861485 times.
✓ Branch 1 taken 137318109 times.
|
138179594 | if (swe <= pcLowSwe_) |
| 373 | 861485 | return pcLowSwePcValue_ + pcDerivativeLowSw_*(swe - pcLowSwe_); | |
| 374 | |||
| 375 |
2/2✓ Branch 0 taken 112207466 times.
✓ Branch 1 taken 25110643 times.
|
137318109 | else if (swe >= 1.0) |
| 376 | 112207466 | return pcDerivativeHighSwEnd_*(swe - 1.0) + entryPressure_; | |
| 377 | |||
| 378 | else | ||
| 379 | return {}; // no regularization | ||
| 380 | } | ||
| 381 | |||
| 382 | /*! | ||
| 383 | * \brief The regularized partial derivative of the capillary pressure w.r.t. the saturation | ||
| 384 | */ | ||
| 385 | OptionalScalar<Scalar> dpc_dswe(const Scalar swe) const | ||
| 386 | { | ||
| 387 |
2/2✓ Branch 0 taken 22 times.
✓ Branch 1 taken 178 times.
|
200 | if (swe <= pcLowSwe_) |
| 388 | 22 | return pcDerivativeLowSw_; | |
| 389 | |||
| 390 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 158 times.
|
178 | else if (swe >= 1.0) |
| 391 | 20 | return pcDerivativeHighSwEnd_; | |
| 392 | |||
| 393 | else | ||
| 394 | return {}; // no regularization | ||
| 395 | } | ||
| 396 | |||
| 397 | /*! | ||
| 398 | * \brief The regularized saturation-capillary pressure curve | ||
| 399 | */ | ||
| 400 | OptionalScalar<Scalar> swe(const Scalar pc) const | ||
| 401 | { | ||
| 402 |
2/2✓ Branch 0 taken 842696 times.
✓ Branch 1 taken 1346804 times.
|
2189500 | if (pc <= entryPressure_) |
| 403 | 842696 | return 1.0 + (pc - entryPressure_)/pcDerivativeHighSwEnd_; | |
| 404 | |||
| 405 |
2/2✓ Branch 0 taken 11236 times.
✓ Branch 1 taken 1335568 times.
|
1346804 | else if (pc >= pcLowSwePcValue_) |
| 406 | 11236 | return (pc - pcLowSwePcValue_)/pcDerivativeLowSw_ + pcLowSwe_; | |
| 407 | |||
| 408 | else | ||
| 409 | return {}; // no regularization | ||
| 410 | } | ||
| 411 | |||
| 412 | /*! | ||
| 413 | * \brief The regularized partial derivative of the saturation to the capillary pressure | ||
| 414 | */ | ||
| 415 | OptionalScalar<Scalar> dswe_dpc(const Scalar pc) const | ||
| 416 | { | ||
| 417 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 180 times.
|
200 | if (pc <= entryPressure_) |
| 418 | 20 | return 1.0/pcDerivativeHighSwEnd_; | |
| 419 | |||
| 420 |
2/2✓ Branch 0 taken 22 times.
✓ Branch 1 taken 158 times.
|
180 | else if (pc >= pcLowSwePcValue_) |
| 421 | 22 | return 1.0/pcDerivativeLowSw_; | |
| 422 | |||
| 423 | else | ||
| 424 | return {}; // no regularization | ||
| 425 | } | ||
| 426 | |||
| 427 | /*! | ||
| 428 | * \brief The regularized relative permeability for the wetting phase | ||
| 429 | */ | ||
| 430 | ✗ | OptionalScalar<Scalar> krw(const Scalar swe) const | |
| 431 | { | ||
| 432 |
2/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 132163249 times.
✓ Branch 3 taken 431369 times.
|
132594618 | if (swe <= 0.0) |
| 433 | ✗ | return 0.0; | |
| 434 |
2/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 22145587 times.
✓ Branch 3 taken 110017662 times.
|
132163249 | else if (swe >= 1.0) |
| 435 | ✗ | return 1.0; | |
| 436 | else | ||
| 437 | 22145587 | return {}; // no regularization | |
| 438 | } | ||
| 439 | |||
| 440 | /*! | ||
| 441 | * \brief The regularized derivative of the relative permeability for the wetting phase w.r.t. saturation | ||
| 442 | */ | ||
| 443 | ✗ | OptionalScalar<Scalar> dkrw_dswe(const Scalar swe) const | |
| 444 | { | ||
| 445 |
2/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 90 times.
✓ Branch 3 taken 10 times.
|
100 | if (swe <= 0.0) |
| 446 | ✗ | return 0.0; | |
| 447 |
2/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 80 times.
✓ Branch 3 taken 10 times.
|
90 | else if (swe >= 1.0) |
| 448 | ✗ | return 0.0; | |
| 449 | else | ||
| 450 | 80 | return {}; // no regularization | |
| 451 | } | ||
| 452 | |||
| 453 | /*! | ||
| 454 | * \brief The regularized relative permeability for the non-wetting phase | ||
| 455 | */ | ||
| 456 | ✗ | OptionalScalar<Scalar> krn(const Scalar swe) const | |
| 457 | { | ||
| 458 |
2/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 132163249 times.
✓ Branch 3 taken 431369 times.
|
132594618 | if (swe <= 0.0) |
| 459 | ✗ | return 1.0; | |
| 460 |
2/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 22145587 times.
✓ Branch 3 taken 110017662 times.
|
132163249 | else if (swe >= 1.0) |
| 461 | ✗ | return 0.0; | |
| 462 | else | ||
| 463 | 22145587 | return {}; // no regularization | |
| 464 | } | ||
| 465 | |||
| 466 | /*! | ||
| 467 | * \brief The regularized derivative of the relative permeability for the non-wetting phase w.r.t. saturation | ||
| 468 | */ | ||
| 469 | ✗ | OptionalScalar<Scalar> dkrn_dswe(const Scalar swe) const | |
| 470 | { | ||
| 471 |
2/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 90 times.
✓ Branch 3 taken 10 times.
|
100 | if (swe <= 0.0) |
| 472 | ✗ | return 0.0; | |
| 473 |
2/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 80 times.
✓ Branch 3 taken 10 times.
|
90 | else if (swe >= 1.0) |
| 474 | ✗ | return 0.0; | |
| 475 | else | ||
| 476 | 80 | return {}; // no regularization | |
| 477 | } | ||
| 478 | |||
| 479 | private: | ||
| 480 | template<class MaterialLaw> | ||
| 481 | 78 | void initPcParameters_(const MaterialLaw* m, const Scalar lowSwe) | |
| 482 | { | ||
| 483 | 156 | const auto lowSw = MaterialLaw::EffToAbs::sweToSw(lowSwe, m->effToAbsParams()); | |
| 484 | 156 | const auto highSw = MaterialLaw::EffToAbs::sweToSw(1.0, m->effToAbsParams()); | |
| 485 | 156 | const auto dsw_dswe = MaterialLaw::EffToAbs::dsw_dswe(m->effToAbsParams()); | |
| 486 | 78 | pcDerivativeLowSw_ = m->template dpc_dsw<false>(lowSw)*dsw_dswe; | |
| 487 | 78 | pcDerivativeHighSwEnd_ = m->template dpc_dsw<false>(highSw)*dsw_dswe; | |
| 488 | 78 | pcLowSwePcValue_ = m->template pc<false>(lowSw); | |
| 489 | 78 | } | |
| 490 | |||
| 491 | Scalar pcLowSwe_; | ||
| 492 | Scalar pcLowSwePcValue_; | ||
| 493 | Scalar entryPressure_; | ||
| 494 | Scalar pcDerivativeLowSw_; | ||
| 495 | Scalar pcDerivativeHighSwEnd_; | ||
| 496 | }; | ||
| 497 | |||
| 498 | /*! | ||
| 499 | * \ingroup Fluidmatrixinteractions | ||
| 500 | * \brief A default configuration for using the Brooks Corey material law | ||
| 501 | */ | ||
| 502 | template<typename Scalar = double> | ||
| 503 | using BrooksCoreyDefault = TwoPMaterialLaw<Scalar, BrooksCorey, BrooksCoreyRegularization<Scalar>, TwoPEffToAbsDefaultPolicy>; | ||
| 504 | |||
| 505 | /*! | ||
| 506 | * \ingroup Fluidmatrixinteractions | ||
| 507 | * \brief A default configuration without regularization for using the Brooks Corey material law | ||
| 508 | */ | ||
| 509 | template<typename Scalar = double> | ||
| 510 | using BrooksCoreyNoReg = TwoPMaterialLaw<Scalar, BrooksCorey, NoRegularization, TwoPEffToAbsDefaultPolicy>; | ||
| 511 | |||
| 512 | } // end namespace Dumux::FluidMatrix | ||
| 513 | |||
| 514 | #endif | ||
| 515 |