GCC Code Coverage Report

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