GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/richards/annulus/problem.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 100 106 94.3%
Functions: 7 14 50.0%
Branches: 110 212 51.9%
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 #ifndef DUMUX_RICHARDS_ANNULUS_PROBLEM_HH
9 #define DUMUX_RICHARDS_ANNULUS_PROBLEM_HH
10
11 #include <iostream>
12 #include <cmath>
13
14 #include <dumux/common/boundarytypes.hh>
15 #include <dumux/common/properties.hh>
16 #include <dumux/common/parameters.hh>
17 #include <dumux/common/numeqvector.hh>
18
19 #include <dumux/common/integrate.hh>
20 #include <dumux/nonlinear/findscalarroot.hh>
21 #include <dumux/discretization/extrusion.hh>
22
23 #include <dumux/porousmediumflow/problem.hh>
24
25 namespace Dumux {
26
27 template<class TypeTag, class FluidSystem>
28 4 class RichardsAnnulusProblem : public PorousMediumFlowProblem<TypeTag>
29 {
30 using ParentType = PorousMediumFlowProblem<TypeTag>;
31
32 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
33 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
34 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
35 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
36
37 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
38 using BoundaryTypes = Dumux::BoundaryTypes<PrimaryVariables::size()>;
39 using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
40 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
41
42 public:
43 4 RichardsAnnulusProblem(std::shared_ptr<const GridGeometry> gridGeometry)
44
6/16
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 4 times.
✓ Branch 15 taken 4 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 12 : ParentType(gridGeometry)
45 {
46 // The default parameters correspond to the benchmark cases C1.1
47 // of Schnepf et al (2020) https://doi.org/10.3389/fpls.2020.00316
48
49 // material parameters
50
5/12
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 4 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 16 k_ = getParam<Scalar>(this->spatialParams().paramGroup() + ".Permeability");
51 static_assert(!FluidSystem::isCompressible(0), "Assumes incompressible fluid");
52 4 rho_ = FluidSystem::H2O::liquidDensity(0, 0);
53 static_assert(FluidSystem::viscosityIsConstant(0), "Assumes constant viscosity fluid");
54 4 mu_ = FluidSystem::H2O::liquidViscosity(0, 0);
55
56 // geometry
57
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 8 innerRadius_ = gridGeometry->bBoxMin()[0];
58
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 8 outerRadius_ = gridGeometry->bBoxMax()[0];
59
60 // boundary conditions (volumetric flow rate in m^2/s, pressure/psi in Pa)
61
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 const auto innerFlux = getParam<Scalar>("Problem.InnerFlux", 0.1); // positive means outflow (cm/day)
62 4 innerFlowRate_ = innerFlux*2*M_PI*innerRadius_*0.01/(24*60*60);
63
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 const auto outerFlux = getParam<Scalar>("Problem.OuterFlux", 0.0); // positive means outflow (cm/day)
64 4 outerFlowRate_ = outerFlux*2*M_PI*outerRadius_*0.01/(24*60*60);
65
66
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 const auto innerPressureHead = getParam<Scalar>("Problem.InnerPressureHeadInCm", -15000);
67 8 const auto innerPressure = headToPressure(innerPressureHead);
68 8 innerPsi_ = applyKirchhoffTransformation_(innerPressure);
69 4 source_ = (innerFlowRate_ + outerFlowRate_)/(M_PI*(innerRadius_*innerRadius_ - outerRadius_*outerRadius_));
70
71 // boundary condition types
72 8 enableOuterNeumann_ = getParam<bool>("Problem.EnableOuterNeumann", true);
73
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 enableInnerNeumann_ = getParam<bool>("Problem.EnableInnerNeumann", false);
74
75
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
 12 std::cout << "Inner flow rate: " << innerFlowRate_ << " m^2/s" << std::endl;
76
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
 12 std::cout << "Outer flow rate: " << outerFlowRate_ << " m^2/s" << std::endl;
77
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 8 std::cout << "Critical inner pressure: " << innerPressure << " Pa "
78
2/4
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 16 << "(-> head: " << innerPressureHead << " cm)" << std::endl;
79
80 // initial condition
81
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 const auto initialHead = getParam<Scalar>("Problem.InitialPressureHeadInCm", -100);
82 8 initialPressure_ = headToPressure(initialHead);
83 4 }
84
85 844 BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
86 {
87
1/2
✓ Branch 0 taken 844 times.
✗ Branch 1 not taken.
 844 BoundaryTypes values;
88
7/8
✓ Branch 0 taken 844 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 422 times.
✓ Branch 3 taken 422 times.
✓ Branch 4 taken 422 times.
✓ Branch 5 taken 422 times.
✓ Branch 6 taken 422 times.
✓ Branch 7 taken 422 times.
 844 if (enableOuterNeumann_ && globalPos[0] > this->gridGeometry().bBoxMax()[0] - 1e-10)
89 values.setAllNeumann();
90
5/8
✓ Branch 0 taken 409 times.
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 409 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 409 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 409 times.
 422 else if (enableInnerNeumann_ && globalPos[0] < this->gridGeometry().bBoxMin()[0] + 1e-10)
91 values.setAllNeumann();
92 else
93 values.setAllDirichlet();
94 844 return values;
95 }
96
97 // negative means extraction
98 NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
99 {
100
2/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 17400 times.
✓ Branch 3 taken 245400 times.
 262800 if (enableSource_)
101 34800 return { -rho_*source_ };
102 else
103 return { 0.0 };
104 }
105
106 // positive means outflow
107 NumEqVector neumannAtPos(const GlobalPosition& globalPos) const
108 {
109
6/12
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 422 times.
✓ Branch 7 taken 409 times.
✓ Branch 8 taken 422 times.
✓ Branch 9 taken 409 times.
✓ Branch 10 taken 422 times.
✓ Branch 11 taken 409 times.
 2493 if (globalPos[0] > this->gridGeometry().bBoxMax()[0] - 1e-10)
110 844 return { rho_*outerFlowRate_/(2*M_PI*outerRadius_) };
111 else
112 818 return { rho_*innerFlowRate_/(2*M_PI*innerRadius_) };
113 }
114
115 PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
116 13 { return exactSolution(globalPos); }
117
118 4808 PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
119 {
120
2/2
✓ Branch 0 taken 4206 times.
✓ Branch 1 taken 602 times.
 4808 if (useStationaryAnalyticInitialSolution_)
121 4206 return exactSolution(globalPos);
122 else
123 {
124 1204 return { initialPressure_ };
125 }
126 }
127
128 /*! \brief Analytical solution
129
130 Analytical solution of stationary incompressible
131 Richards equations with source term B:
132
133 ψ = -µ/k A/(2π) ln(r) + r^2/4 µ/k B + C
134
135 Three constants to fix with three boundary conditions:
136 - flux at inner boundary (r=r_i): q_i = 2πr_i k/µ ∂ψ/∂r
137 - flux at outer boundary (r=r_o): q_o = -2πr_o k/µ ∂ψ/∂r
138 - ψ = T(p) at inner boundary -> ψ = ψ_ri
139
140 Yields the solution in terms of q_i, q_o, ψ_ri:
141
142 ψ = ψ_ri - µ/k A/(2π) ln(r/r_i) + (r^2 - r_i^2)/4 µ/k B
143 where
144 A = - q_i + πr_i^2 B
145 B = (q_i + q_o)/(π(r_i^2 - r_o^2)) [source]
146 */
147 43819 PrimaryVariables exactSolution(const GlobalPosition& globalPos, const Scalar innerPsi) const
148 {
149 43819 const auto r = globalPos[0];
150 43819 const auto ri2 = innerRadius_*innerRadius_;
151 43819 const auto A = -innerFlowRate_ + M_PI*ri2*source_;
152 43819 const auto psi = innerPsi
153 43819 - A*mu_/(2*M_PI*k_)*std::log(r/innerRadius_)
154 43819 + source_*0.25*(r*r - ri2)*mu_/k_;
155
156 131457 return { applyInverseKirchhoffTransformation_(psi) };
157 }
158
159 PrimaryVariables exactSolution(const GlobalPosition& globalPos) const
160
1/4
✓ Branch 1 taken 8400 times.
✗ Branch 2 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 12619 { return exactSolution(globalPos, innerPsi_); }
161
162 ✗ Scalar nonwettingReferencePressure() const
163 ✗ { return 1.0e5; };
164
165 template<class SolutionVector, class GridVariables>
166 58 Scalar computeTime(const SolutionVector& p, const GridVariables& vars, const Scalar initialSaturation) const
167 {
168
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 const auto& gg = this->gridGeometry();
169
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 auto fvGeometry = localView(gg);
170
2/6
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 174 auto elemVolVars = localView(vars.curGridVolVars());
171 58 Scalar totalWaterVolume = 0.0;
172 58 Scalar refWaterVolume = 0.0;
173 using Extrusion = Extrusion_t<GridGeometry>;
174
4/8
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 19858 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 19800 times.
✗ Branch 10 not taken.
 39716 for (const auto& element : elements(gg.gridView()))
175 {
176
1/2
✓ Branch 1 taken 19800 times.
✗ Branch 2 not taken.
 19800 fvGeometry.bindElement(element);
177
1/2
✓ Branch 1 taken 19800 times.
✗ Branch 2 not taken.
 19800 elemVolVars.bindElement(element, fvGeometry, p);
178
4/4
✓ Branch 0 taken 39600 times.
✓ Branch 1 taken 19800 times.
✓ Branch 2 taken 39600 times.
✓ Branch 3 taken 19800 times.
 118800 for (const auto& scv : scvs(fvGeometry))
179 {
180 39600 const auto& volVars = elemVolVars[scv];
181 118800 totalWaterVolume += Extrusion::volume(fvGeometry, scv)*volVars.porosity()*volVars.saturation(0);
182 118800 refWaterVolume += Extrusion::volume(fvGeometry, scv)*volVars.porosity()*initialSaturation;
183 }
184 }
185
186
1/2
✓ Branch 0 taken 58 times.
✗ Branch 1 not taken.
 116 return (refWaterVolume - totalWaterVolume)/(innerFlowRate_ + outerFlowRate_);
187 }
188
189 template<class SolutionVector, class GridVariables>
190 50 void updateStorageDerivative(const SolutionVector& p, const SolutionVector& pOld,
191 const GridVariables& vars, const Scalar dt,
192 std::vector<Scalar>& storageDerivative) const
193 {
194
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 const auto& gg = this->gridGeometry();
195
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 auto fvGeometry = localView(gg);
196
2/6
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 150 auto elemVolVars = localView(vars.curGridVolVars());
197
3/6
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 50 times.
✗ Branch 7 not taken.
 150 auto oldElemVolVars = localView(vars.prevGridVolVars());
198
2/6
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 100 storageDerivative.assign(storageDerivative.size(), 0.0);
199
3/8
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 50 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 50 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 150 auto volumes = storageDerivative;
200 using Extrusion = Extrusion_t<GridGeometry>;
201
4/8
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 15050 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 15000 times.
✗ Branch 10 not taken.
 30100 for (const auto& element : elements(gg.gridView()))
202 {
203
1/2
✓ Branch 1 taken 15000 times.
✗ Branch 2 not taken.
 15000 fvGeometry.bindElement(element);
204
1/2
✓ Branch 1 taken 15000 times.
✗ Branch 2 not taken.
 15000 elemVolVars.bindElement(element, fvGeometry, p);
205
1/2
✓ Branch 1 taken 15000 times.
✗ Branch 2 not taken.
 15000 oldElemVolVars.bindElement(element, fvGeometry, pOld);
206
4/4
✓ Branch 0 taken 30000 times.
✓ Branch 1 taken 15000 times.
✓ Branch 2 taken 30000 times.
✓ Branch 3 taken 15000 times.
 90000 for (const auto& scv : scvs(fvGeometry))
207 {
208 30000 const auto& volVars = elemVolVars[scv];
209 30000 const auto& oldVolVars = oldElemVolVars[scv];
210 120000 storageDerivative[scv.dofIndex()] += Extrusion::volume(fvGeometry, scv)*(volVars.saturation(0) - oldVolVars.saturation(0))/dt;
211 90000 volumes[scv.dofIndex()] += Extrusion::volume(fvGeometry, scv);
212 }
213 }
214
215
4/4
✓ Branch 0 taken 15050 times.
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 15050 times.
✓ Branch 3 taken 50 times.
 30150 for (int i = 0; i < storageDerivative.size(); ++i)
216 45150 storageDerivative[i] /= volumes[i];
217 50 }
218
219 // return saturation provided pressure in Pa
220 const Scalar saturation(Scalar const pw) const
221 {
222
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 8 const auto& pcKrSwCurve = this->spatialParams().pcKrSwCurve();
223 4 const Scalar pc = this->nonwettingReferencePressure() - pw;
224
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 return pcKrSwCurve.sw(pc);
225 }
226
227 /*!
228 * \brief convert pressure head in cm to pressure in Pa
229 */
230 ✗ Scalar headToPressure(const Scalar head) const
231 {
232
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 10 return this->nonwettingReferencePressure() + head / 100.0 * rho_ * 9.81;
233 }
234
235 /*!
236 * \brief convert pressure in Pa to potential psi in Pa (Kirchhoff transformation)
237 */
238 Scalar pressureToPsi(const Scalar pw) const
239 {
240 104 return applyKirchhoffTransformation_(pw);
241 }
242
243 /*!
244 * \brief Disable the source term (for the instationary problem)
245 */
246 ✗ void disableSource()
247
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 { enableSource_ = false; }
248
249 /*!
250 * \brief Disable using the stationary initial solution (for the instationary problem)
251 */
252 ✗ void disableStationaryInitialSolution()
253
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 { useStationaryAnalyticInitialSolution_ = false; }
254
255 /*!
256 * \brief Enable using Neumann boundary conditions on the inner boundary
257 */
258 ✗ void enableInnerNeumannBC()
259
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 { enableInnerNeumann_ = true; }
260
261 private:
262
263 /*!
264 * \brief relative wetting phase permeability
265 */
266 739348195 Scalar kr_(Scalar pw) const
267 {
268 1478696390 const auto& pcKrSwCurve = this->spatialParams().pcKrSwCurve();
269 739348195 const Scalar pc = this->nonwettingReferencePressure() - pw;
270 739348195 const Scalar sw = pcKrSwCurve.sw(pc);
271 739348195 return pcKrSwCurve.krw(sw);
272 }
273
274 /*!
275 * \brief transformed variable p->psi
276 */
277 Scalar applyKirchhoffTransformation_(Scalar pw) const
278 {
279
6/12
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 50 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 50 times.
✗ Branch 17 not taken.
 56 return integrateScalarFunction(
280 370655719 [&](const double p){ return kr_(p); }, -1.5e6, pw, 1e-20
281 );
282 }
283
284 /*!
285 * \brief inverse transformation psi->p
286 */
287 Scalar applyInverseKirchhoffTransformation_(Scalar psi) const
288 {
289 43819 const auto residual = [&, psi](const auto& p){
290
3/4
✗ Branch 3 not taken.
✓ Branch 4 taken 43819 times.
✓ Branch 6 taken 1746651 times.
✓ Branch 7 taken 128898 times.
 1919368 return psi - applyKirchhoffTransformation_(p);
291 };
292
293 43819 return findScalarRootBrent(-1.5e6, 1.0e6, residual, 1e-14, 200000);
294 }
295
296 Scalar k_, mu_, rho_;
297 Scalar innerRadius_, outerRadius_, source_;
298 Scalar innerFlowRate_, outerFlowRate_, innerPsi_;
299 Scalar initialPressure_;
300 bool enableOuterNeumann_, enableInnerNeumann_;
301 bool useStationaryAnalyticInitialSolution_ = true;
302 bool enableSource_ = true;
303 };
304
305 } // end namespace Dumux
306
307 #endif
308