| 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 / relPerm <-> saturation relation | ||
| 11 | * using a linear relation smoothed at the upper and lower bounds for kr. | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_MATERIAL_FLUIDMATRIX_TWOP_SMOOTHED_LINEAR_LAW_HH | ||
| 14 | #define DUMUX_MATERIAL_FLUIDMATRIX_TWOP_SMOOTHED_LINEAR_LAW_HH | ||
| 15 | |||
| 16 | |||
| 17 | #include <cmath> | ||
| 18 | #include <algorithm> | ||
| 19 | |||
| 20 | #include <dumux/common/parameters.hh> | ||
| 21 | #include <dumux/common/spline.hh> | ||
| 22 | #include <dumux/material/fluidmatrixinteractions/2p/efftoabsdefaultpolicy.hh> | ||
| 23 | #include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh> | ||
| 24 | |||
| 25 | namespace Dumux::FluidMatrix { | ||
| 26 | |||
| 27 | /*! | ||
| 28 | * \ingroup Fluidmatrixinteractions | ||
| 29 | * \brief Implements a linear saturation-capillary pressure relation | ||
| 30 | * | ||
| 31 | * The entry pressure is reached at \f$\mathrm{\overline{S}_w = 1}\f$, the maximum | ||
| 32 | * capillary pressure is observed at \f$\mathrm{\overline{S}_w = 0}\f$. | ||
| 33 | * | ||
| 34 | * The relative permeabilities are 0 or 1 outside of the range of effective saturation. | ||
| 35 | * However, the transition between the linearly changing and the constant part is not smooth but with a kink. | ||
| 36 | * The Newton scheme does not like that. Therefore a smooth transition is accomplished by interpolating these | ||
| 37 | * regions with a spline. | ||
| 38 | * | ||
| 39 | * An example of the regularization of the relative permeability is shown below: | ||
| 40 | * \image html regularizedLinearKr.png | ||
| 41 | */ | ||
| 42 | template<class ScalarType, class EffToAbsPolicy = TwoPEffToAbsDefaultPolicy> | ||
| 43 | class SmoothedLinearLaw : public Adapter<SmoothedLinearLaw<ScalarType, EffToAbsPolicy>, PcKrSw> | ||
| 44 | { | ||
| 45 | |||
| 46 | public: | ||
| 47 | using Scalar = ScalarType; | ||
| 48 | using EffToAbsParams = typename EffToAbsPolicy::template Params<Scalar>; | ||
| 49 | using EffToAbs = EffToAbsPolicy; | ||
| 50 | |||
| 51 | /*! | ||
| 52 | * \brief Return whether this law is regularized | ||
| 53 | */ | ||
| 54 | static constexpr bool isRegularized() | ||
| 55 | { return false; } | ||
| 56 | |||
| 57 | /*! | ||
| 58 | * \brief Return the number of fluid phases | ||
| 59 | */ | ||
| 60 | static constexpr int numFluidPhases() | ||
| 61 | { return 2; } | ||
| 62 | |||
| 63 | /*! | ||
| 64 | * \brief The parameter type | ||
| 65 | * \tparam Scalar The scalar type | ||
| 66 | */ | ||
| 67 | struct Params | ||
| 68 | { | ||
| 69 | 2 | Params(Scalar pe, Scalar pcMax, Scalar krLowS, Scalar krHighS) | |
| 70 | 2 | : pe_(pe), pcMax_(pcMax), krLowS_(krLowS), krHighS_(krHighS) | |
| 71 | {} | ||
| 72 | |||
| 73 | ✗ | Scalar pe() const { return pe_; } | |
| 74 | void setPe(Scalar pe) { pe_ = pe; } | ||
| 75 | |||
| 76 | ✗ | Scalar pcMax() const { return pcMax_; } | |
| 77 | void setPcMax(Scalar pcMax) { pcMax_ = pcMax; } | ||
| 78 | |||
| 79 | ✗ | Scalar krLowS() const { return krLowS_; } | |
| 80 | void setKrLowS(Scalar krLowS) { krLowS_ = krLowS; } | ||
| 81 | |||
| 82 | ✗ | Scalar krHighS() const { return krHighS_; } | |
| 83 | void setKrHighS(Scalar krHighS) { krHighS_ = krHighS; } | ||
| 84 | |||
| 85 | bool operator== (const Params& p) const | ||
| 86 | { | ||
| 87 | return Dune::FloatCmp::eq(pe(), p.pe(), 1e-6) | ||
| 88 | && Dune::FloatCmp::eq(pcMax(), p.pcMax(), 1e-6) | ||
| 89 | && Dune::FloatCmp::eq(krLowS(), p.krLowS(), 1e-6) | ||
| 90 | && Dune::FloatCmp::eq(krHighS(), p.krHighS(), 1e-6); | ||
| 91 | } | ||
| 92 | |||
| 93 | private: | ||
| 94 | Scalar pe_, pcMax_, krLowS_, krHighS_; | ||
| 95 | }; | ||
| 96 | |||
| 97 | /*! | ||
| 98 | * \brief Deleted default constructor (so we are never in an undefined state) | ||
| 99 | * \note store owning pointers to laws instead if you need default-constructible objects | ||
| 100 | */ | ||
| 101 | SmoothedLinearLaw() = delete; | ||
| 102 | |||
| 103 | /*! | ||
| 104 | * \brief Construct from a subgroup from the global parameter tree | ||
| 105 | * \note This will give you nice error messages if a mandatory parameter is missing | ||
| 106 | */ | ||
| 107 | 2 | explicit SmoothedLinearLaw(const std::string& paramGroup) | |
| 108 | 2 | : SmoothedLinearLaw(makeParams(paramGroup), EffToAbs::template makeParams<Scalar>(paramGroup)) | |
| 109 | 2 | {} | |
| 110 | |||
| 111 | /*! | ||
| 112 | * \brief Construct from parameter structs | ||
| 113 | * \note More efficient constructor but you need to ensure all parameters are initialized | ||
| 114 | */ | ||
| 115 | 2 | SmoothedLinearLaw(const Params& params, | |
| 116 | const EffToAbsParams& effToAbsParams = {}) | ||
| 117 | : params_(params) | ||
| 118 | , effToAbsParams_(effToAbsParams) | ||
| 119 | 2 | , splineM_((1.0 - ((1.0 - params_.krHighS()) + params_.krLowS())/2.0 ) | |
| 120 | 2 | / (1.0 - (1.0 - params_.krHighS()) - params_.krLowS())) | |
| 121 | , splineLowS_(0.0, params_.krLowS(), // x1, x2 | ||
| 122 | 2 | 0.0, params_.krLowS()/2.0, // y1, y2 | |
| 123 | 0.0, splineM_) // m1, m2 | ||
| 124 | , splineHighS_(params_.krHighS(), 1.0, // x1, x2 | ||
| 125 | 2 | 1.0 - (1.0 - params_.krHighS())/2.0, 1.0, // y1, y2 | |
| 126 | 2 | splineM_, 0.0) // m1, m2 | |
| 127 | 2 | {} | |
| 128 | |||
| 129 | /*! | ||
| 130 | * \brief Construct from a subgroup from the global parameter tree | ||
| 131 | * \note This will give you nice error messages if a mandatory parameter is missing | ||
| 132 | */ | ||
| 133 | 2 | static Params makeParams(const std::string& paramGroup) | |
| 134 | { | ||
| 135 | 2 | const auto pe = getParamFromGroup<Scalar>(paramGroup, "SmoothedLinearLawPe"); | |
| 136 | 2 | const auto pcMax = getParamFromGroup<Scalar>(paramGroup, "SmoothedLinearLawPcMax"); | |
| 137 | 2 | const auto krLowS = getParamFromGroup<Scalar>(paramGroup, "SmoothedLinearLawKrLowS"); | |
| 138 | 2 | const auto krHighS = getParamFromGroup<Scalar>(paramGroup, "SmoothedLinearLawKrHighS"); | |
| 139 | 4 | return {pe, pcMax, krLowS, krHighS}; | |
| 140 | } | ||
| 141 | |||
| 142 | /*! | ||
| 143 | * \brief The capillary pressure-saturation curve | ||
| 144 | */ | ||
| 145 | Scalar pc(Scalar swe) const | ||
| 146 | { | ||
| 147 | 586991 | return (1.0 - swe)*(params_.pcMax() - params_.pe()) + params_.pe(); | |
| 148 | } | ||
| 149 | |||
| 150 | /*! | ||
| 151 | * \brief The inverse saturation-capillary pressure curve | ||
| 152 | */ | ||
| 153 | Scalar swe(Scalar pc) const | ||
| 154 | { | ||
| 155 | return 1.0 - (pc - params_.pe())/(params_.pcMax() - params_.pe()); | ||
| 156 | } | ||
| 157 | |||
| 158 | /*! | ||
| 159 | * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure | ||
| 160 | */ | ||
| 161 | Scalar endPointPc() const | ||
| 162 | { return params_.pe(); } | ||
| 163 | |||
| 164 | /*! | ||
| 165 | * \brief The partial derivative of the capillary pressure w.r.t. the effective saturation | ||
| 166 | */ | ||
| 167 | Scalar dpc_dswe(Scalar swe) const | ||
| 168 | { | ||
| 169 | return params_.pe() - params_.pcMax(); | ||
| 170 | } | ||
| 171 | |||
| 172 | /*! | ||
| 173 | * \brief The partial derivative of the effective saturation w.r.t. the capillary pressure | ||
| 174 | */ | ||
| 175 | Scalar dswe_dpc(Scalar pc) const | ||
| 176 | { | ||
| 177 | return 1.0/(params_.pe() - params_.pcMax()); | ||
| 178 | } | ||
| 179 | |||
| 180 | /*! | ||
| 181 | * \brief The relative permeability for the wetting phase. | ||
| 182 | */ | ||
| 183 | Scalar krw(Scalar swe) const | ||
| 184 | { | ||
| 185 | ✗ | return relperm_(swe); | |
| 186 | } | ||
| 187 | |||
| 188 | /*! | ||
| 189 | * \brief The relative permeability for the non-wetting phase. | ||
| 190 | */ | ||
| 191 | Scalar krn(Scalar swe) const | ||
| 192 | { | ||
| 193 | ✗ | Scalar sne = 1.0 - swe; | |
| 194 | ✗ | return relperm_(sne); | |
| 195 | } | ||
| 196 | |||
| 197 | /*! | ||
| 198 | * \brief Equality comparison with another instance | ||
| 199 | */ | ||
| 200 | bool operator== (const SmoothedLinearLaw<Scalar, EffToAbs>& o) const | ||
| 201 | { | ||
| 202 | return params_ == o.params_ | ||
| 203 | && effToAbsParams_ == o.effToAbsParams_; | ||
| 204 | } | ||
| 205 | |||
| 206 | const EffToAbsParams& effToAbsParams() const | ||
| 207 | { return effToAbsParams_; } | ||
| 208 | |||
| 209 | private: | ||
| 210 | |||
| 211 | 1170096 | Scalar relperm_(Scalar S) const | |
| 212 | { | ||
| 213 | 1170096 | const Scalar lowS = params_.krLowS(); | |
| 214 | 1170096 | const Scalar highS = params_.krHighS(); | |
| 215 | |||
| 216 | // check whether the saturation is unpyhsical | ||
| 217 |
2/2✓ Branch 0 taken 677178 times.
✓ Branch 1 taken 492918 times.
|
1170096 | if (S >= 1.0) |
| 218 | return 1.0; | ||
| 219 |
2/2✓ Branch 0 taken 184276 times.
✓ Branch 1 taken 492902 times.
|
677178 | else if (S <= 0.0) |
| 220 | return 0; | ||
| 221 | // check whether the permeability needs to be regularized | ||
| 222 |
2/2✓ Branch 0 taken 87422 times.
✓ Branch 1 taken 96854 times.
|
184276 | else if (S < lowS) |
| 223 | 87422 | return splineLowS_.eval(S); | |
| 224 |
2/2✓ Branch 0 taken 87406 times.
✓ Branch 1 taken 9448 times.
|
96854 | else if (S > highS) |
| 225 | 87406 | return splineHighS_.eval(S); | |
| 226 | |||
| 227 | // straight line for S \in [lowS, highS] | ||
| 228 | 9448 | return lowS/2.0 + splineM_*(S - lowS); | |
| 229 | } | ||
| 230 | |||
| 231 | Params params_; | ||
| 232 | EffToAbsParams effToAbsParams_; | ||
| 233 | Scalar splineM_; | ||
| 234 | Dumux::Spline<Scalar> splineLowS_; | ||
| 235 | Dumux::Spline<Scalar> splineHighS_; | ||
| 236 | }; | ||
| 237 | } // end namespace Dumux::FluidMatrix | ||
| 238 | |||
| 239 | #endif | ||
| 240 |