GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/3p/convection/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 59 62 95.2%
Functions: 8 12 66.7%
Branches: 79 134 59.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 ThreePTests
10 * \brief Test for the ThreePModel in combination with the NI model for a convection problem.
11 */
12
13 #ifndef DUMUX_3PNI_CONVECTION_PROBLEM_HH
14 #define DUMUX_3PNI_CONVECTION_PROBLEM_HH
15
16 #include <algorithm>
17 #include <cmath>
18
19 #include <dumux/common/boundarytypes.hh>
20 #include <dumux/common/parameters.hh>
21 #include <dumux/common/properties.hh>
22 #include <dumux/common/numeqvector.hh>
23
24 #include <dumux/porousmediumflow/problem.hh>
25 #include <dumux/material/components/h2o.hh>
26
27 namespace Dumux {
28
29 /*!
30 * \ingroup ThreePTests
31 *
32 * \brief Test for the ThreePModel in combination with the NI model for a convection problem.
33 *
34 * The simulation domain is a tube where water with an elevated temperature is injected
35 * at a constant rate on the left hand side.
36 *
37 * Initially the domain is fully saturated with water at a constant temperature.
38 * On the left hand side water is injected at a constant rate and on the right
39 * hand side a Dirichlet boundary with constant pressure, saturation and
40 * temperature is applied.
41 *
42 * The results are compared to an analytical solution where a retarded front
43 * velocity is calculated as follows:
44 \f[
45 v_{Front}=\frac{q S_{water}}{\phi S_{total}}
46 \f]
47 */
48 template <class TypeTag>
49 class ThreePNIConvectionProblem : public PorousMediumFlowProblem<TypeTag>
50 {
51 using ParentType = PorousMediumFlowProblem<TypeTag>;
52
53 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
54 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
55 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
56 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
57 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
58 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
59 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
60 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
61
62 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
63 using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
64 using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
65 using VolumeVariables = typename GridVariables::GridVolumeVariables::VolumeVariables;
66
67 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
68 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
69 using IapwsH2O = Components::H2O<Scalar>;
70
71 // copy some indices for convenience
72 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
73 enum {
74 // index of the primary variables
75 pressureIdx = Indices::pressureIdx,
76 swIdx = Indices::swIdx,
77 snIdx = Indices::snIdx,
78 temperatureIdx = Indices::temperatureIdx,
79 wPhaseIdx = FluidSystem::wPhaseIdx,
80 conti0EqIdx = Indices::conti0EqIdx,
81 energyEqIdx = Indices::energyEqIdx
82 };
83
84 enum { dimWorld = GridView::dimensionworld };
85
86 using Element = typename GridView::template Codim<0>::Entity;
87 using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
88
89 public:
90 2 ThreePNIConvectionProblem(std::shared_ptr<const GridGeometry> gridGeometry)
91
8/24
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
 6 : ParentType(gridGeometry)
92 {
93 //initialize fluid system
94
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 FluidSystem::init();
95
96
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 2 name_ = getParam<std::string>("Problem.Name");
97
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 outputInterval_ = getParam<int>("Problem.OutputInterval");
98
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 darcyVelocity_ = getParam<Scalar>("Problem.DarcyVelocity");
99
100 2 temperatureHigh_ = 291.;
101 2 temperatureLow_ = 290.;
102 2 pressureHigh_ = 2e5;
103 2 pressureLow_ = 1e5;
104
105
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 5 temperatureExact_.resize(this->gridGeometry().numDofs());
106 2 }
107
108 //! Get exact temperature vector for output
109 const std::vector<Scalar>& getExactTemperature()
110 {
111
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return temperatureExact_;
112 }
113
114 //! Update the analytical temperature
115 86 void updateExactTemperature(const SolutionVector& curSol, Scalar time)
116 {
117 258 const auto someElement = *(elements(this->gridGeometry().gridView()).begin());
118
119 172 const auto someElemSol = elementSolution(someElement, curSol, this->gridGeometry());
120
1/2
✓ Branch 2 taken 86 times.
✗ Branch 3 not taken.
 86 const auto someInitSol = initialAtPos(someElement.geometry().center());
121
122
2/4
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
 344 const auto someFvGeometry = localView(this->gridGeometry()).bindElement(someElement);
123 86 const auto someScv = *(scvs(someFvGeometry).begin());
124
125 86 VolumeVariables volVars;
126
1/2
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
 86 volVars.update(someElemSol, *this, someElement, someScv);
127
128
2/4
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
 172 const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
129
1/2
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
 86 const auto densityW = volVars.density(wPhaseIdx);
130
3/6
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
 258 const auto heatCapacityW = IapwsH2O::liquidHeatCapacity(someInitSol[temperatureIdx], someInitSol[pressureIdx]);
131 86 const auto storageW = densityW*heatCapacityW*porosity;
132 86 const auto densityS = volVars.solidDensity();
133 86 const auto heatCapacityS = volVars.solidHeatCapacity();
134 86 const auto storageTotal = storageW + densityS*heatCapacityS*(1 - porosity);
135
2/4
✓ Branch 2 taken 86 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 86 times.
✗ Branch 6 not taken.
 172 std::cout << "storage: " << storageTotal << '\n';
136
137 using std::max;
138
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 86 times.
 86 time = max(time, 1e-10);
139 86 const Scalar retardedFrontVelocity = darcyVelocity_*storageW/storageTotal/porosity;
140
2/4
✓ Branch 2 taken 86 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 86 times.
✗ Branch 6 not taken.
 172 std::cout << "retarded velocity: " << retardedFrontVelocity << '\n';
141
2/4
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
 172 auto fvGeometry = localView(this->gridGeometry());
142
5/10
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 6966 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 6880 times.
✗ Branch 14 not taken.
 14104 for (const auto& element : elements(this->gridGeometry().gridView()))
143 {
144
1/2
✓ Branch 1 taken 6880 times.
✗ Branch 2 not taken.
 6880 fvGeometry.bindElement(element);
145
4/4
✓ Branch 0 taken 17200 times.
✓ Branch 1 taken 6880 times.
✓ Branch 2 taken 13760 times.
✓ Branch 3 taken 3440 times.
 41280 for (auto&& scv : scvs(fvGeometry))
146 {
147
2/2
✓ Branch 0 taken 271 times.
✓ Branch 1 taken 3169 times.
 17200 auto dofIdxGlobal = scv.dofIndex();
148 17200 auto dofPosition = scv.dofPosition();
149
4/4
✓ Branch 0 taken 1361 times.
✓ Branch 1 taken 15839 times.
✓ Branch 2 taken 1361 times.
✓ Branch 3 taken 15839 times.
 34400 temperatureExact_[dofIdxGlobal] = (dofPosition[0] < retardedFrontVelocity*time) ? temperatureHigh_ : temperatureLow_;
150 }
151 }
152 86 }
153
154 /*!
155 * \name Problem parameters
156 */
157 // \{
158
159 /*!
160 * \brief The problem name.
161 *
162 * This is used as a prefix for files generated by the simulation.
163 */
164 const std::string& name() const
165 {
166
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return name_;
167 }
168
169 // \}
170
171 /*!
172 * \name Boundary conditions
173 */
174 // \{
175
176 /*!
177 * \brief Specifies which kind of boundary condition should be
178 * used for which equation on a given boundary segment.
179 *
180 * \param globalPos The position for which the bc type should be evaluated
181 */
182 146988 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
183 {
184 146988 BoundaryTypes values;
185
10/10
✓ Branch 0 taken 146078 times.
✓ Branch 1 taken 910 times.
✓ Branch 2 taken 146078 times.
✓ Branch 3 taken 910 times.
✓ Branch 4 taken 146078 times.
✓ Branch 5 taken 910 times.
✓ Branch 6 taken 146078 times.
✓ Branch 7 taken 910 times.
✓ Branch 8 taken 146078 times.
✓ Branch 9 taken 910 times.
 734940 if(globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
186 {
187 values.setAllDirichlet();
188 }
189 else
190 {
191 values.setAllNeumann();
192 }
193 146988 return values;
194 }
195
196 /*!
197 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
198 *
199 * \param globalPos The position for which the bc type should be evaluated
200 *
201 */
202 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
203 {
204 330 return initialAtPos(globalPos);
205 }
206
207 /*!
208 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
209 *
210 * \param element The finite element
211 * \param fvGeometry The finite-volume geometry in the box scheme
212 * \param elemVolVars The element volume variables
213 * \param elemFluxVarsCache Flux variables caches for all faces in stencil
214 * \param scvf The subcontrolvolume face
215 * Negative values mean influx.
216 */
217 676070 NumEqVector neumann(const Element &element,
218 const FVElementGeometry& fvGeometry,
219 const ElementVolumeVariables& elemVolVars,
220 const ElementFluxVariablesCache& elemFluxVarsCache,
221 const SubControlVolumeFace& scvf) const
222 {
223 676070 NumEqVector values(0.0);
224 676070 const auto globalPos = scvf.ipGlobal();
225
4/4
✓ Branch 0 taken 3672 times.
✓ Branch 1 taken 583848 times.
✓ Branch 2 taken 3672 times.
✓ Branch 3 taken 583848 times.
 1352140 const auto& volVars = elemVolVars[scvf.insideScvIdx()];
226
227
4/4
✓ Branch 0 taken 4222 times.
✓ Branch 1 taken 671848 times.
✓ Branch 2 taken 4222 times.
✓ Branch 3 taken 671848 times.
 1352140 if(globalPos[0] < eps_)
228 {
229 8444 values[conti0EqIdx] = -darcyVelocity_*volVars.density(wPhaseIdx);
230 4222 values[energyEqIdx] = -darcyVelocity_*volVars.density(wPhaseIdx)
231 8444 *IapwsH2O::liquidEnthalpy(temperatureHigh_, volVars.pressure(wPhaseIdx));
232 }
233 676070 return values;
234 }
235
236 // \}
237
238 /*!
239 * \name Volume terms
240 */
241 // \{
242
243 /*!
244 * \brief Evaluates the initial value for a control volume.
245 *
246 * \param globalPos The position for which the initial condition should be evaluated
247 *
248 */
249 ✗ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
250 {
251 PrimaryVariables values;
252
2/3
✓ Branch 2 taken 43 times.
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
 658 values[pressureIdx] = pressureLow_; // initial condition for the pressure
253
2/3
✓ Branch 2 taken 43 times.
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
 658 values[swIdx] = 1.0; // initial condition for the wetting phase saturation
254
2/3
✓ Branch 2 taken 43 times.
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
 658 values[snIdx] = 1e-10; // initial condition for the nonwetting phase saturation
255
4/7
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 43 times.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 43 times.
✗ Branch 9 not taken.
 1316 values[temperatureIdx] = temperatureLow_;
256 ✗ return values;
257 }
258
259 // \}
260
261 private:
262 Scalar temperatureHigh_;
263 Scalar temperatureLow_;
264 Scalar pressureHigh_;
265 Scalar pressureLow_;
266 Scalar darcyVelocity_;
267 static constexpr Scalar eps_ = 1e-6;
268 std::string name_;
269 int outputInterval_;
270 std::vector<Scalar> temperatureExact_;
271 };
272
273 } // end namespace Dumux
274
275 #endif
276