GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 117 156 75.0%
Functions: 56 77 72.7%
Branches: 265 676 39.2%

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