| 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 A spline approximation wrapper for 2p material laws | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_MATERIAL_FLUIDMATRIX_TWOP_SPLINE_MATERIAL_LAW_HH | ||
| 13 | #define DUMUX_MATERIAL_FLUIDMATRIX_TWOP_SPLINE_MATERIAL_LAW_HH | ||
| 14 | |||
| 15 | #include <memory> // unique_ptr | ||
| 16 | |||
| 17 | #include <dumux/common/parameters.hh> | ||
| 18 | #include <dumux/common/monotonecubicspline.hh> | ||
| 19 | #include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh> | ||
| 20 | #include <dumux/material/fluidmatrixinteractions/2p/materiallaw.hh> | ||
| 21 | |||
| 22 | namespace Dumux::FluidMatrix { | ||
| 23 | |||
| 24 | /*! | ||
| 25 | * \ingroup Fluidmatrixinteractions | ||
| 26 | * \brief A spline approximation wrapper for 2p material laws | ||
| 27 | * \tparam TwoPMaterialLaw the type of material law to be wrapped | ||
| 28 | * \tparam approximatePcSwInverse if this is set true, the | ||
| 29 | * spline approximates sw(pc) and evaluating pc(sw) needs spline inversion. | ||
| 30 | * if this is false, the spline approximates pc(sw) and evaluating | ||
| 31 | * sw(pc) needs spline inversion. Spline inversion is rather expensive | ||
| 32 | * since it has to be done numerically. | ||
| 33 | */ | ||
| 34 | template<class TwoPMaterialLaw, bool approximatePcSwInverse = false> | ||
| 35 | class SplineTwoPMaterialLaw | ||
| 36 | : public TwoPMaterialLaw | ||
| 37 | , public Adapter<SplineTwoPMaterialLaw<TwoPMaterialLaw, approximatePcSwInverse>, PcKrSw> | ||
| 38 | { | ||
| 39 | public: | ||
| 40 | using Scalar = typename TwoPMaterialLaw::Scalar; | ||
| 41 | |||
| 42 | using BasicParams = typename TwoPMaterialLaw::BasicParams; | ||
| 43 | using EffToAbsParams = typename TwoPMaterialLaw::BasicParams; | ||
| 44 | using RegularizationParams = typename TwoPMaterialLaw::RegularizationParams; | ||
| 45 | |||
| 46 | /*! | ||
| 47 | * \brief Return the number of fluid phases | ||
| 48 | */ | ||
| 49 | static constexpr int numFluidPhases() | ||
| 50 | { return 2; } | ||
| 51 | |||
| 52 | /*! | ||
| 53 | * \brief We are always regularized in the sense that we replace | ||
| 54 | * the original curve by a cubic spline | ||
| 55 | */ | ||
| 56 | static constexpr bool isRegularized() | ||
| 57 | { return true; } | ||
| 58 | |||
| 59 | /*! | ||
| 60 | * \brief Deleted default constructor (so we are never in an undefined state) | ||
| 61 | * \note store owning pointers to laws instead if you need default-constructible objects | ||
| 62 | */ | ||
| 63 | SplineTwoPMaterialLaw() = delete; | ||
| 64 | |||
| 65 | /*! | ||
| 66 | * \brief Construct from a subgroup from the global parameter tree | ||
| 67 | * \note This will give you nice error messages if a mandatory parameter is missing | ||
| 68 | */ | ||
| 69 | 1 | explicit SplineTwoPMaterialLaw(const std::string& paramGroup) | |
| 70 |
3/10✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
1 | : TwoPMaterialLaw(paramGroup) |
| 71 | { | ||
| 72 | 1 | const std::array<Scalar, 2> defaultInterval{{ 0.01, 1.0 }}; | |
| 73 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | const auto sweInterval = getParamFromGroup<std::array<Scalar, 2>>(paramGroup, "SplineSweInterval", defaultInterval); |
| 74 | 3 | swInterval_ = {{ TwoPMaterialLaw::EffToAbs::sweToSw(sweInterval[0], this->effToAbsParams()), | |
| 75 | 3 | TwoPMaterialLaw::EffToAbs::sweToSw(sweInterval[1], this->effToAbsParams()) }}; | |
| 76 | 2 | swIntervalPc_ = {{ TwoPMaterialLaw::pc(swInterval_[1]), | |
| 77 | 2 | TwoPMaterialLaw::pc(swInterval_[0]) }}; | |
| 78 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | numSwSamples_ = getParamFromGroup<std::size_t>(paramGroup, "SplineNumSwSamples", 30); |
| 79 | |||
| 80 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
|
1 | pcSpline_ = makeSweSpline_( |
| 81 | 100 | [&](const Scalar s){ return TwoPMaterialLaw::pc(s); }, | |
| 82 | approximatePcSwInverse | ||
| 83 | ); | ||
| 84 | |||
| 85 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
|
1 | krwSpline_ = makeSweSpline_( |
| 86 | 100 | [&](const Scalar s){ return TwoPMaterialLaw::krw(s); } | |
| 87 | ); | ||
| 88 | |||
| 89 |
3/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1 | krnSpline_ = makeSweSpline_( |
| 90 | 100 | [&](const Scalar s){ return TwoPMaterialLaw::krn(s); } | |
| 91 | ); | ||
| 92 | 1 | } | |
| 93 | |||
| 94 | /*! | ||
| 95 | * \brief Construct from parameter structs | ||
| 96 | * \note More efficient constructor but you need to ensure all parameters are initialized | ||
| 97 | */ | ||
| 98 | SplineTwoPMaterialLaw(const std::array<Scalar, 2>& sweInterval, | ||
| 99 | std::size_t numSwSamples, | ||
| 100 | TwoPMaterialLaw&& twoP) | ||
| 101 | : TwoPMaterialLaw(std::move(twoP)) | ||
| 102 | , numSwSamples_(numSwSamples) | ||
| 103 | { | ||
| 104 | swInterval_ = {{ TwoPMaterialLaw::EffToAbs::sweToSw(sweInterval[0], this->effToAbsParams()), | ||
| 105 | TwoPMaterialLaw::EffToAbs::sweToSw(sweInterval[1], this->effToAbsParams()) }}; | ||
| 106 | swIntervalPc_ = {{ TwoPMaterialLaw::pc(swInterval_[1]), | ||
| 107 | TwoPMaterialLaw::pc(swInterval_[0]) }}; | ||
| 108 | } | ||
| 109 | |||
| 110 | /*! | ||
| 111 | * \brief The capillary pressure-saturation curve | ||
| 112 | */ | ||
| 113 | 3000000 | Scalar pc(const Scalar sw) const | |
| 114 | { | ||
| 115 |
8/8✓ Branch 0 taken 2999997 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 2999997 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 2999994 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 2999994 times.
✓ Branch 7 taken 3 times.
|
6000000 | if (sw > swInterval_[0] && sw < swInterval_[1]) |
| 116 | { | ||
| 117 | if constexpr (approximatePcSwInverse) | ||
| 118 | return pcSpline_->evalInverse(sw); | ||
| 119 | else | ||
| 120 | 5999988 | return pcSpline_->eval(sw); | |
| 121 | } | ||
| 122 | |||
| 123 | 12 | return TwoPMaterialLaw::pc(sw); | |
| 124 | } | ||
| 125 | |||
| 126 | /*! | ||
| 127 | * \brief The partial derivative of the capillary pressure w.r.t. the saturation | ||
| 128 | */ | ||
| 129 | 1000000 | Scalar dpc_dsw(const Scalar sw) const | |
| 130 | { | ||
| 131 |
8/8✓ Branch 0 taken 999999 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 999999 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 999998 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 999998 times.
✓ Branch 7 taken 1 times.
|
2000000 | if (sw > swInterval_[0] && sw < swInterval_[1]) |
| 132 | { | ||
| 133 | if constexpr (approximatePcSwInverse) | ||
| 134 | return 1.0/pcSpline_->evalDerivative(pcSpline_->evalInverse(sw)); | ||
| 135 | else | ||
| 136 | 1999996 | return pcSpline_->evalDerivative(sw); | |
| 137 | } | ||
| 138 | |||
| 139 | 2 | return TwoPMaterialLaw::dpc_dsw(sw); | |
| 140 | } | ||
| 141 | |||
| 142 | /*! | ||
| 143 | * \brief The saturation-capillary pressure curve | ||
| 144 | */ | ||
| 145 | 1000000 | Scalar sw(const Scalar pc) const | |
| 146 | { | ||
| 147 |
8/8✓ Branch 0 taken 999999 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 999999 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 999998 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 999998 times.
✓ Branch 7 taken 1 times.
|
2000000 | if (pc > swIntervalPc_[0] && pc < swIntervalPc_[1]) |
| 148 | { | ||
| 149 | if constexpr (approximatePcSwInverse) | ||
| 150 | return pcSpline_->eval(pc); | ||
| 151 | else | ||
| 152 | 1999996 | return pcSpline_->evalInverse(pc); | |
| 153 | } | ||
| 154 | |||
| 155 | 4 | return TwoPMaterialLaw::sw(pc); | |
| 156 | } | ||
| 157 | |||
| 158 | /*! | ||
| 159 | * \brief The partial derivative of the saturation to the capillary pressure | ||
| 160 | */ | ||
| 161 | 1000000 | Scalar dsw_dpc(const Scalar pc) const | |
| 162 | { | ||
| 163 |
8/8✓ Branch 0 taken 999999 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 999999 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 999998 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 999998 times.
✓ Branch 7 taken 1 times.
|
2000000 | if (pc > swIntervalPc_[0] && pc < swIntervalPc_[1]) |
| 164 | { | ||
| 165 | if constexpr (approximatePcSwInverse) | ||
| 166 | return pcSpline_->evalDerivative(pc); | ||
| 167 | else | ||
| 168 | 2999994 | return 1.0/pcSpline_->evalDerivative(pcSpline_->evalInverse(pc)); | |
| 169 | } | ||
| 170 | |||
| 171 | 4 | return TwoPMaterialLaw::dsw_dpc(pc); | |
| 172 | } | ||
| 173 | |||
| 174 | /*! | ||
| 175 | * \brief The relative permeability for the wetting phase | ||
| 176 | */ | ||
| 177 | 1000000 | Scalar krw(const Scalar sw) const | |
| 178 | { | ||
| 179 |
8/8✓ Branch 0 taken 999999 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 999999 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 999998 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 999998 times.
✓ Branch 7 taken 1 times.
|
2000000 | if (sw > swInterval_[0] && sw < swInterval_[1]) |
| 180 | 1999996 | return krwSpline_->eval(sw); | |
| 181 | |||
| 182 | 4 | return TwoPMaterialLaw::krw(sw); | |
| 183 | } | ||
| 184 | |||
| 185 | /*! | ||
| 186 | * \brief The derivative of the relative permeability for the wetting phase w.r.t. saturation | ||
| 187 | */ | ||
| 188 | 1000000 | Scalar dkrw_dsw(const Scalar sw) const | |
| 189 | { | ||
| 190 |
12/12✓ Branch 0 taken 999999 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 999999 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 999999 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 999998 times.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 999998 times.
✓ Branch 9 taken 1 times.
✓ Branch 10 taken 999998 times.
✓ Branch 11 taken 1 times.
|
3000000 | if (sw > swInterval_[0] && sw < swInterval_[1]) |
| 191 | 1999996 | return krwSpline_->evalDerivative(sw); | |
| 192 | |||
| 193 | 2 | return TwoPMaterialLaw::dkrw_dsw(sw); | |
| 194 | } | ||
| 195 | |||
| 196 | /*! | ||
| 197 | * \brief The relative permeability for the non-wetting phase | ||
| 198 | */ | ||
| 199 | 1000000 | Scalar krn(const Scalar sw) const | |
| 200 | { | ||
| 201 |
8/8✓ Branch 0 taken 999999 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 999999 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 999998 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 999998 times.
✓ Branch 7 taken 1 times.
|
2000000 | if (sw > swInterval_[0] && sw < swInterval_[1]) |
| 202 | 1999996 | return krnSpline_->eval(sw); | |
| 203 | |||
| 204 | 4 | return TwoPMaterialLaw::krn(sw); | |
| 205 | } | ||
| 206 | |||
| 207 | /*! | ||
| 208 | * \brief The derivative of the relative permeability for the non-wetting phase w.r.t. saturation | ||
| 209 | */ | ||
| 210 | 1000000 | Scalar dkrn_dsw(const Scalar sw) const | |
| 211 | { | ||
| 212 |
8/8✓ Branch 0 taken 999999 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 999999 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 999998 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 999998 times.
✓ Branch 7 taken 1 times.
|
2000000 | if (sw > swInterval_[0] && sw < swInterval_[1]) |
| 213 | 1999996 | return krnSpline_->evalDerivative(sw); | |
| 214 | |||
| 215 | 2 | return TwoPMaterialLaw::dkrn_dsw(sw); | |
| 216 | } | ||
| 217 | |||
| 218 | private: | ||
| 219 | template<class Function> | ||
| 220 | std::unique_ptr<MonotoneCubicSpline<Scalar>> | ||
| 221 | 3 | makeSweSpline_(const Function& f, bool invert = false) const | |
| 222 | { | ||
| 223 | 12 | const auto sw = linspace(swInterval_[0], swInterval_[1], numSwSamples_); | |
| 224 | |||
| 225 | 9 | auto values = sw; | |
| 226 | 12 | std::transform(sw.begin(), sw.end(), values.begin(), f); | |
| 227 | |||
| 228 | 3 | if (invert) | |
| 229 | ✗ | return std::make_unique<MonotoneCubicSpline<Scalar>>(values, sw); | |
| 230 | else | ||
| 231 | 3 | return std::make_unique<MonotoneCubicSpline<Scalar>>(sw, values); | |
| 232 | } | ||
| 233 | |||
| 234 | std::unique_ptr<MonotoneCubicSpline<Scalar>> pcSpline_; | ||
| 235 | std::unique_ptr<MonotoneCubicSpline<Scalar>> krwSpline_; | ||
| 236 | std::unique_ptr<MonotoneCubicSpline<Scalar>> krnSpline_; | ||
| 237 | |||
| 238 | std::array<Scalar, 2> swInterval_; | ||
| 239 | std::array<Scalar, 2> swIntervalPc_; | ||
| 240 | std::size_t numSwSamples_; | ||
| 241 | }; | ||
| 242 | |||
| 243 | } // end namespace Dumux::FluidMatrix | ||
| 244 | |||
| 245 | #endif | ||
| 246 |