GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/1p/nonisothermal/problem_convection.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 57 61 93.4%
Functions: 8 14 57.1%
Branches: 59 108 54.6%
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 OnePTests
10 * \brief Test for the OnePModel in combination with the NI model for a convection problem.
11 *
12 * The simulation domain is a tube where water with an elevated temperature is injected
13 * at a constant rate on the left hand side.
14 */
15
16 #ifndef DUMUX_1PNI_CONVECTION_PROBLEM_HH
17 #define DUMUX_1PNI_CONVECTION_PROBLEM_HH
18
19 #include <cmath>
20
21 #include <dumux/common/properties.hh>
22 #include <dumux/common/parameters.hh>
23
24 #include <dumux/common/boundarytypes.hh>
25 #include <dumux/common/numeqvector.hh>
26 #include <dumux/porousmediumflow/problem.hh>
27
28 #include <dumux/material/components/h2o.hh>
29 namespace Dumux {
30
31 /*!
32 * \ingroup OnePTests
33 * \brief Test for the OnePModel in combination with the NI model for a convection problem.
34 *
35 * The simulation domain is a tube where water with an elevated temperature is injected
36 * at a constant rate on the left hand side.
37 *
38 * Initially the domain is fully saturated with water at a constant temperature.
39 * On the left hand side water is injected at a constant rate and on the right hand side
40 * a Dirichlet boundary with constant pressure, saturation and temperature is applied.
41 *
42 * The result of the analytical solution is written into the vtu files.
43 * This problem uses the \ref OnePModel and \ref NIModel model.
44 */
45 template <class TypeTag>
46 class OnePNIConvectionProblem : public PorousMediumFlowProblem<TypeTag>
47 {
48 using ParentType = PorousMediumFlowProblem<TypeTag>;
49 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
50 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
51 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
52 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
53 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
54 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
55 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
56
57 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
58 using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
59 using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
60 using VolumeVariables = typename GridVariables::GridVolumeVariables::VolumeVariables;
61 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
62 using IapwsH2O = Components::H2O<Scalar>;
63
64 enum { dimWorld = GridView::dimensionworld };
65
66 // copy some indices for convenience
67 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
68 enum {
69 // indices of the primary variables
70 pressureIdx = Indices::pressureIdx,
71 temperatureIdx = Indices::temperatureIdx
72 };
73 enum {
74 // index of the transport equation
75 conti0EqIdx = Indices::conti0EqIdx,
76 energyEqIdx = Indices::energyEqIdx
77 };
78
79 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
80 using Element = typename GridView::template Codim<0>::Entity;
81 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
82 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
83
84 public:
85 2 OnePNIConvectionProblem(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup)
86
4/14
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 6 : ParentType(gridGeometry, paramGroup)
87 {
88 //initialize fluid system
89 2 FluidSystem::init();
90
91
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");
92
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 darcyVelocity_ = getParam<Scalar>("Problem.DarcyVelocity");
93
94 2 temperatureHigh_ = 291.0;
95 2 temperatureLow_ = 290.0;
96 2 pressureHigh_ = 2e5;
97 2 pressureLow_ = 1e5;
98
99
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());
100 2 }
101
102 //! Get exact temperature vector for output
103 const std::vector<Scalar>& getExactTemperature()
104 {
105
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return temperatureExact_;
106 }
107
108 /*!
109 * \brief Update the analytical temperature
110 * The results are compared to an analytical solution where a retarded front velocity is calculated as follows:
111 \f[
112 v_{Front}=\frac{q S_{water}}{\phi S_{total}}
113 \f]
114 */
115 88 void updateExactTemperature(const SolutionVector& curSol, Scalar time)
116 {
117 264 const auto someElement = *(elements(this->gridGeometry().gridView()).begin());
118
119 176 auto someElemSol = elementSolution(someElement, curSol, this->gridGeometry());
120 88 const auto someInitSol = initialAtPos(someElement.geometry().center());
121
122
2/4
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
 352 const auto someFvGeometry = localView(this->gridGeometry()).bindElement(someElement);
123
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 const auto someScv = *(scvs(someFvGeometry).begin());
124
125
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 VolumeVariables volVars;
126
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 volVars.update(someElemSol, *this, someElement, someScv);
127
128
2/4
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
 176 const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
129
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 const auto densityW = volVars.density();
130
3/6
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 88 times.
✗ Branch 8 not taken.
 264 const auto heatCapacityW = IapwsH2O::liquidHeatCapacity(someInitSol[temperatureIdx], someInitSol[pressureIdx]);
131 88 const auto storageW = densityW*heatCapacityW*porosity;
132 88 const auto densityS = volVars.solidDensity();
133 88 const auto heatCapacityS = volVars.solidHeatCapacity();
134 88 const auto storageTotal = storageW + densityS*heatCapacityS*(1 - porosity);
135
2/4
✓ Branch 2 taken 88 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 88 times.
✗ Branch 6 not taken.
 176 std::cout << "storage: " << storageTotal << '\n';
136
137 using std::max;
138
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 88 times.
 88 time = max(time, 1e-10);
139 88 const Scalar retardedFrontVelocity = darcyVelocity_*storageW/storageTotal/porosity;
140
2/4
✓ Branch 2 taken 88 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 88 times.
✗ Branch 6 not taken.
 176 std::cout << "retarded velocity: " << retardedFrontVelocity << '\n';
141
2/4
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
 176 auto fvGeometry = localView(this->gridGeometry());
142
5/10
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 88 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 7128 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 7040 times.
✗ Branch 14 not taken.
 14432 for (const auto& element : elements(this->gridGeometry().gridView()))
143 {
144
1/2
✓ Branch 1 taken 7040 times.
✗ Branch 2 not taken.
 7040 fvGeometry.bindElement(element);
145
4/4
✓ Branch 0 taken 10560 times.
✓ Branch 1 taken 7040 times.
✓ Branch 2 taken 7040 times.
✓ Branch 3 taken 3520 times.
 28160 for (auto&& scv : scvs(fvGeometry))
146 {
147
2/2
✓ Branch 0 taken 280 times.
✓ Branch 1 taken 3240 times.
 10560 auto dofIdxGlobal = scv.dofIndex();
148 10560 auto dofPosition = scv.dofPosition();
149
2/2
✓ Branch 0 taken 842 times.
✓ Branch 1 taken 9718 times.
 10560 temperatureExact_[dofIdxGlobal] = (dofPosition[0] < retardedFrontVelocity*time) ? temperatureHigh_ : temperatureLow_;
150 }
151 }
152 88 }
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 1346 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
183 {
184 1346 BoundaryTypes bcTypes;
185
186
6/6
✓ Branch 0 taken 673 times.
✓ Branch 1 taken 673 times.
✓ Branch 2 taken 673 times.
✓ Branch 3 taken 673 times.
✓ Branch 4 taken 673 times.
✓ Branch 5 taken 673 times.
 4038 if(globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
187 bcTypes.setAllDirichlet();
188 else
189 bcTypes.setAllNeumann();
190
191 1346 return bcTypes;
192 }
193
194 /*!
195 * \brief Evaluates the boundary conditions for a Dirichlet control volume.
196 *
197 * \param globalPos The center of the finite volume which ought to be set.
198 *
199 * For this method, the \a values parameter stores primary variables.
200 */
201 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
202 {
203 538 return initial_(globalPos);
204 }
205
206 /*!
207 * \brief Evaluates the boundary conditions for a Neumann
208 * boundary segment in dependency on the current solution.
209 *
210 * \param element The finite element
211 * \param fvGeometry The finite volume geometry of the element
212 * \param elemVolVars All volume variables for the element
213 * \param elemFluxVarsCache Flux variables caches for all faces in stencil
214 * \param scvf The sub-control volume face
215 *
216 * This method is used for cases, when the Neumann condition depends on the
217 * solution and requires some quantities that are specific to the fully-implicit method.
218 * The \a values store the mass flux of each phase normal to the boundary.
219 * Negative values indicate an inflow.
220 */
221 1049 NumEqVector neumann(const Element& element,
222 const FVElementGeometry& fvGeometry,
223 const ElementVolumeVariables& elemVolVars,
224 const ElementFluxVariablesCache& elemFluxVarsCache,
225 const SubControlVolumeFace& scvf) const
226 {
227 1049 NumEqVector values(0.0);
228 1049 const auto globalPos = scvf.ipGlobal();
229
2/4
✓ Branch 0 taken 655 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 655 times.
✗ Branch 3 not taken.
 2098 const auto& volVars = elemVolVars[scvf.insideScvIdx()];
230
231
1/2
✓ Branch 0 taken 1049 times.
✗ Branch 1 not taken.
 1049 if(globalPos[0] < eps_)
232 {
233 2098 values[conti0EqIdx] = -darcyVelocity_*volVars.density();
234 3147 values[energyEqIdx] = values[conti0EqIdx]*IapwsH2O::liquidEnthalpy(temperatureHigh_, volVars.pressure());
235 }
236 1049 return values;
237 }
238
239 // \}
240
241 /*!
242 * \name Volume terms
243 */
244 // \{
245
246 /*!
247 * \brief Evaluates the initial value for a control volume.
248 *
249 * \param globalPos The position for which the initial condition should be evaluated
250 *
251 * For this method, the \a values parameter stores primary
252 * variables.
253 */
254 ✗ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
255 {
256 498 return initial_(globalPos);
257 }
258
259 // \}
260
261 private:
262 // the internal method for the initial condition
263 ✗ PrimaryVariables initial_(const GlobalPosition &globalPos) const
264 {
265 518 PrimaryVariables priVars(0.0);
266
2/3
✓ Branch 2 taken 44 times.
✓ Branch 3 taken 44 times.
✗ Branch 4 not taken.
 518 priVars[pressureIdx] = pressureLow_; // initial condition for the pressure
267
4/7
✓ Branch 3 taken 44 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 44 times.
✓ Branch 6 taken 44 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 44 times.
✗ Branch 9 not taken.
 1036 priVars[temperatureIdx] = temperatureLow_;
268 ✗ return priVars;
269 }
270
271 Scalar temperatureHigh_;
272 Scalar temperatureLow_;
273 Scalar pressureHigh_;
274 Scalar pressureLow_;
275 Scalar darcyVelocity_;
276 static constexpr Scalar eps_ = 1e-6;
277 std::string name_;
278 std::vector<Scalar> temperatureExact_;
279 };
280
281 } // end namespace Dumux
282
283 #endif
284