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 |