GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/material/fluidmatrixinteractions/2p/heatpipelaw.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 0 29 0.0%
Functions: 0 5 0.0%
Branches: 0 14 0.0%
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 <-> saturation relation
11 * for the heatpipe problem.
12 */
13 #ifndef DUMUX_MATERIAL_FLUIDMATRIX_TWOP_HEATPIPELAW_HH
14 #define DUMUX_MATERIAL_FLUIDMATRIX_TWOP_HEATPIPELAW_HH
15
16 #include <cmath>
17 #include <algorithm>
18
19 #include <dumux/common/parameters.hh>
20 #include <dumux/common/spline.hh>
21 #include <dumux/material/fluidmatrixinteractions/fluidmatrixinteraction.hh>
22 #include <dumux/material/fluidmatrixinteractions/2p/efftoabsdefaultpolicy.hh>
23
24 namespace Dumux::FluidMatrix {
25
26 /*!
27 * \ingroup Fluidmatrixinteractions
28 * \brief Implementation of the capillary pressure <-> saturation
29 * relation for the heatpipe problem.
30 */
31 template<class ScalarType, class EffToAbsPolicy = TwoPEffToAbsDefaultPolicy>
32 class HeatPipeLaw : Adapter<HeatPipeLaw<ScalarType, EffToAbsPolicy>, PcKrSw>
33 {
34
35 public:
36 using Scalar = ScalarType;
37 using EffToAbsParams = typename EffToAbsPolicy::template Params<Scalar>;
38 using EffToAbs = EffToAbsPolicy;
39
40 /*!
41 * \brief Return the number of fluid phases
42 */
43 static constexpr int numFluidPhases()
44 { return 2; }
45
46 /*!
47 * \brief Return whether this law is regularized
48 */
49 static constexpr bool isRegularized()
50 { return false; }
51
52 /*!
53 * \brief The parameter type
54 * \tparam Scalar The scalar type
55 */
56 struct Params
57 {
58 ✗ Params(Scalar gamma, Scalar p0)
59 ✗ : gamma_(gamma), p0_(p0)
60 {}
61
62 ✗ Scalar gamma() const { return gamma_; }
63 void setGamma(Scalar gamma) { gamma_ = gamma; }
64
65 ✗ Scalar p0() const { return p0_; }
66 void setP0(Scalar p0) { p0_ = p0; }
67
68 bool operator== (const Params& p) const
69 {
70 return Dune::FloatCmp::eq(gamma(), p.gamma(), 1e-6)
71 && Dune::FloatCmp::eq(p0(), p.p0(), 1e-6);
72 }
73
74 private:
75 Scalar gamma_, p0_;
76 };
77
78 /*!
79 * \brief Deleted default constructor (so we are never in an undefined state)
80 * \note store owning pointers to laws instead if you need default-constructible objects
81 */
82 HeatPipeLaw() = delete;
83
84 /*!
85 * \brief Construct from a subgroup from the global parameter tree
86 * \note This will give you nice error messages if a mandatory parameter is missing
87 */
88 explicit HeatPipeLaw(const std::string& paramGroup)
89 {
90 params_ = makeParams(paramGroup);
91 effToAbsParams_ = EffToAbs::template makeParams<Scalar>(paramGroup);
92 }
93
94 /*!
95 * \brief Construct from parameter structs
96 * \note More efficient constructor but you need to ensure all parameters are initialized
97 */
98 ✗ HeatPipeLaw(const Params& params,
99 const EffToAbsParams& effToAbsParams = {})
100 : params_(params)
101 , effToAbsParams_(effToAbsParams)
102 , spline_(eps_, 1.0, // x1, x2
103 eps_*eps_*eps_, 1.0, // y1, y2
104 ✗ 3.0*eps_*eps_, 0.0) // m1, m2
105 ✗ {}
106
107 /*!
108 * \brief Construct from a subgroup from the global parameter tree
109 * \note This will give you nice error messages if a mandatory parameter is missing
110 */
111 static Params makeParams(const std::string& paramGroup)
112 {
113 const auto gamma = getParamFromGroup<Scalar>(paramGroup, "HeatpipeLawGamma");
114 const auto p0 = getParamFromGroup<Scalar>(paramGroup, "HeatpipeLawP0");
115 return {gamma, p0};
116 }
117
118 /*!
119 * \brief The capillary pressure-saturation curve.
120 *
121 * \param sw Saturation of the wetting phase \f$\mathrm{S_w}\f$
122 */
123 ✗ Scalar pc(const Scalar sw) const
124 {
125 ✗ const Scalar swe = EffToAbs::swToSwe(sw, effToAbsParams_);
126 ✗ const Scalar sne = 1 - swe;
127 ✗ const Scalar p0Gamma = params_.p0()*params_.gamma();
128
129 // regularization
130 ✗ if (sne >= 1.0)
131 {
132 ✗ const Scalar y = p0Gamma*( (1.263*1.0 - 2.120)*1.0 + 1.417)*1.0;
133 ✗ const Scalar m = p0Gamma*((3*1.263*1.0 - 2*2.120)*1.0 + 1.417);
134 ✗ return (sne - 1)*m + y;
135 }
136 ✗ else if (sne <= 0.0)
137 ✗ return sne * p0Gamma*1.417;
138 else
139 ✗ return p0Gamma*((1.263*sne - 2.120)*sne + 1.417) * sne;
140 }
141
142 /*!
143 * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure
144 */
145 Scalar endPointPc() const
146 { return 0.0; }
147
148 /*!
149 * \brief The partial derivative of the capillary
150 * pressure w.r.t. the effective saturation.
151 *
152 * \param sw Saturation of the wetting phase \f$\mathrm{S_w}\f$
153 */
154 Scalar dpc_dsw(const Scalar sw) const
155 {
156 const Scalar swe = EffToAbs::swToSwe(sw, effToAbsParams_);
157 const Scalar sne = 1 - swe;
158 const Scalar p0Gamma = params_.p0()*params_.gamma();
159 if (sne > 1.0)
160 sne = 1.0;
161 else if (sne <= 0.0)
162 return -p0Gamma*1.417*EffToAbs::dswe_dsw(effToAbsParams_);
163 else
164 return - p0Gamma*((3*1.263*sne - 2*2.120)*sne + 1.417)*EffToAbs::dswe_dsw(effToAbsParams_);
165 }
166
167 /*!
168 * \brief The relative permeability for the wetting phase of
169 * the medium.
170 *
171 * \param sw Saturation of the wetting phase \f$\mathrm{S_w}\f$
172 */
173 Scalar krw(const Scalar sw) const
174 {
175 ✗ const Scalar swe = EffToAbs::swToSwe(sw, effToAbsParams_);
176 ✗ return kr_(swe);
177 }
178
179 /*!
180 * \brief The relative permeability for the non-wetting phase
181 * of the medium.
182 *
183 * \param sw Saturation of the wetting phase \f$\mathrm{S_w}\f$
184 */
185 Scalar krn(const Scalar sw) const
186 {
187 ✗ const Scalar swe = EffToAbs::swToSwe(sw, effToAbsParams_);
188 ✗ const Scalar sne = 1 - swe; // TODO does this make sense?
189 ✗ return kr_(sne);
190 }
191
192 /*!
193 * \brief Equality comparison with another instance
194 */
195 bool operator== (const HeatPipeLaw<Scalar, EffToAbs>& o) const
196 {
197 return params_ == o.params_
198 && effToAbsParams_ == o.effToAbsParams_;
199 }
200
201 const EffToAbsParams& effToAbsParams() const
202 { return effToAbsParams_; }
203
204 private:
205
206 ✗ Scalar kr_(Scalar s) const
207 {
208 ✗ if (s >= 1.0)
209 return 1;
210 ✗ else if (s <= 0.0)
211 return 0;
212 ✗ else if (s > eps_)
213 ✗ return spline_.eval(s); // regularize
214 else
215 ✗ return s*s*s;
216 }
217
218 Params params_;
219 EffToAbsParams effToAbsParams_;
220 static constexpr Scalar eps_ = 0.95;
221 Dumux::Spline<Scalar> spline_;
222 };
223 } // end namespace Dumux::FluidMatrix
224
225 #endif
226