GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/1pnc/1p2c/nonisothermal/convection/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 60 66 90.9%
Functions: 12 24 50.0%
Branches: 85 143 59.4%
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 OnePNCTests
10 * \brief Test for the OnePNCModel in combination with the NI model for a conduction problem.
11 *
12 * The simulation domain is a tube with an elevated temperature on the left hand side.
13 */
14
15 #ifndef DUMUX_1P2CNI_CONVECTION_TEST_PROBLEM_HH
16 #define DUMUX_1P2CNI_CONVECTION_TEST_PROBLEM_HH
17 #include <dumux/common/properties.hh>
18 #include <dumux/common/parameters.hh>
19
20 #include <dumux/common/boundarytypes.hh>
21 #include <dumux/common/numeqvector.hh>
22 #include <dumux/porousmediumflow/problem.hh>
23
24 #include <dumux/material/components/h2o.hh>
25 namespace Dumux {
26
27 /*!
28 * \ingroup OnePNCTests
29 * \brief Test for the OnePTwoCModel in combination with the NI model for a convection problem.
30 *
31 * The simulation domain is a tube where water with an elevated temperature is injected
32 * at a constant rate on the left hand side.
33 *
34 * Initially, the domain is fully saturated with water at a constant temperature.
35 * On the left hand side water is injected at a constant rate and on the right hand side
36 * a Dirichlet boundary with constant pressure, saturation and temperature is applied.
37 *
38 * The result of the analytical solution is written into the vtu files.
39 *
40 * This problem uses the \ref OnePModel and \ref NIModel model.
41 */
42
43 template <class TypeTag>
44 class OnePTwoCNIConvectionProblem : public PorousMediumFlowProblem<TypeTag>
45 {
46 using ParentType = PorousMediumFlowProblem<TypeTag>;
47
48 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
49 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
50 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
51 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
52 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
53 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
54 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
55 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
56 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
57
58 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
59 using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
60 using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
61
62 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
63 using Element = typename GridView::template Codim<0>::Entity;
64 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
65 using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
66 using IapwsH2O = Components::H2O<Scalar>;
67
68 // copy some indices for convenience
69 enum
70 {
71 // indices of the primary variables
72 pressureIdx = Indices::pressureIdx,
73 temperatureIdx = Indices::temperatureIdx,
74
75 // component indices
76 H2OIdx = FluidSystem::compIdx(FluidSystem::MultiPhaseFluidSystem::H2OIdx),
77 N2Idx = FluidSystem::compIdx(FluidSystem::MultiPhaseFluidSystem::N2Idx),
78
79 // indices of the equations
80 contiH2OEqIdx = Indices::conti0EqIdx + H2OIdx,
81 contiN2EqIdx = Indices::conti0EqIdx + N2Idx,
82 energyEqIdx = Indices::energyEqIdx
83 };
84
85 //! Property that defines whether mole or mass fractions are used
86 static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
87 static const int dimWorld = GridView::dimensionworld;
88 using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
89
90 public:
91 3 OnePTwoCNIConvectionProblem(std::shared_ptr<const GridGeometry> gridGeometry)
92
7/18
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 3 times.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 3 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
 9 : ParentType(gridGeometry)
93 {
94 //initialize fluid system
95
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 FluidSystem::init();
96
97 // stating in the console whether mole or mass fractions are used
98 if(useMoles)
99
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 6 std::cout<<"problem uses mole fractions"<<std::endl;
100 else
101 std::cout<<"problem uses mass fractions"<<std::endl;
102
103
4/10
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 7 temperatureExact_.resize(gridGeometry->numDofs(), 290.0);
104
105
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 darcyVelocity_ = getParam<Scalar>("Problem.DarcyVelocity");
106
107 3 temperatureHigh_ = 291.;
108 3 temperatureLow_ = 290.;
109 3 pressureHigh_ = 2e5;
110 3 pressureLow_ = 1e5;
111 3 }
112
113 //! Get the analytical temperature
114 const std::vector<Scalar>& getExactTemperature()
115 {
116
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 return temperatureExact_;
117 }
118
119 /*!
120 * \brief Update the analytical temperature
121 * The results are compared to an analytical solution where a retarded front velocity is calculated as follows:
122 \f[
123 v_{Front}=\frac{q S_{water}}{\phi S_{total}}
124 \f]
125 */
126 130 void updateExactTemperature(const SolutionVector& curSol, Scalar time)
127 {
128 390 const auto someElement = *(elements(this->gridGeometry().gridView()).begin());
129
130 260 auto someElemSol = elementSolution(someElement, curSol, this->gridGeometry());
131 260 const auto someInitSol = initialAtPos(someElement.geometry().center());
132
133
2/4
✓ Branch 1 taken 130 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 130 times.
✗ Branch 5 not taken.
 520 const auto someFvGeometry = localView(this->gridGeometry()).bindElement(someElement);
134
1/2
✓ Branch 1 taken 130 times.
✗ Branch 2 not taken.
 130 const auto someScv = *(scvs(someFvGeometry).begin());
135
136
1/2
✓ Branch 1 taken 130 times.
✗ Branch 2 not taken.
 130 VolumeVariables volVars;
137
1/2
✓ Branch 1 taken 130 times.
✗ Branch 2 not taken.
 130 volVars.update(someElemSol, *this, someElement, someScv);
138
139
2/4
✓ Branch 1 taken 130 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 130 times.
✗ Branch 5 not taken.
 260 const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
140
1/2
✓ Branch 1 taken 130 times.
✗ Branch 2 not taken.
 130 const auto densityW = volVars.density();
141
3/6
✓ Branch 1 taken 130 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 130 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 130 times.
✗ Branch 8 not taken.
 390 const auto heatCapacityW = IapwsH2O::liquidHeatCapacity(someInitSol[temperatureIdx], someInitSol[pressureIdx]);
142 130 const auto storageW = densityW*heatCapacityW*porosity;
143 130 const auto densityS = volVars.solidDensity();
144 130 const auto heatCapacityS = volVars.solidHeatCapacity();
145 130 const auto storageTotal = storageW + densityS*heatCapacityS*(1 - porosity);
146
2/4
✓ Branch 2 taken 130 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 130 times.
✗ Branch 6 not taken.
 260 std::cout << "storage: " << storageTotal << '\n';
147
148 using std::max;
149
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 130 times.
 130 time = max(time, 1e-10);
150 130 const Scalar retardedFrontVelocity = darcyVelocity_*storageW/storageTotal/porosity;
151
2/4
✓ Branch 2 taken 130 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 130 times.
✗ Branch 6 not taken.
 260 std::cout << "retarded velocity: " << retardedFrontVelocity << '\n';
152
153
2/4
✓ Branch 1 taken 130 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 130 times.
✗ Branch 5 not taken.
 390 auto fvGeometry = localView(this->gridGeometry());
154
9/14
✓ Branch 1 taken 130 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 130 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 130 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 10400 times.
✓ Branch 10 taken 130 times.
✓ Branch 11 taken 10400 times.
✓ Branch 12 taken 130 times.
✓ Branch 14 taken 10400 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 10400 times.
✗ Branch 18 not taken.
 10790 for (const auto& element : elements(this->gridGeometry().gridView()))
155 {
156
1/2
✓ Branch 1 taken 10400 times.
✗ Branch 2 not taken.
 10400 fvGeometry.bindElement(element);
157
4/4
✓ Branch 0 taken 20960 times.
✓ Branch 1 taken 10400 times.
✓ Branch 2 taken 14080 times.
✓ Branch 3 taken 3520 times.
 48960 for (auto&& scv : scvs(fvGeometry))
158 {
159
2/2
✓ Branch 0 taken 542 times.
✓ Branch 1 taken 6338 times.
 20960 auto dofIdxGlobal = scv.dofIndex();
160 20960 auto dofPosition = scv.dofPosition();
161
4/4
✓ Branch 0 taken 1682 times.
✓ Branch 1 taken 19278 times.
✓ Branch 2 taken 1682 times.
✓ Branch 3 taken 19278 times.
 41920 temperatureExact_[dofIdxGlobal] = (dofPosition[0] < retardedFrontVelocity*time) ? temperatureHigh_ : temperatureLow_;
162 }
163 }
164 130 }
165
166 /*!
167 * \name Problem parameters
168 */
169 // \{
170
171 /*!
172 * \brief Returns the temperature within the domain [K].
173 *
174 * This problem assumes a temperature of 20 degrees Celsius.
175 */
176 Scalar temperature() const
177 { return 273.15 + 20; } // in [K]
178
179 // \}
180
181 /*!
182 * \name Boundary conditions
183 */
184 // \{
185
186 /*!
187 * \brief Specifies which kind of boundary condition should be
188 * used for which equation on a given boundary segment.
189 *
190 * \param globalPos The position for which the bc type should be evaluated
191 */
192 580244 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
193 {
194 580244 BoundaryTypes values;
195
196
10/10
✓ Branch 0 taken 578138 times.
✓ Branch 1 taken 2106 times.
✓ Branch 2 taken 578138 times.
✓ Branch 3 taken 2106 times.
✓ Branch 4 taken 578138 times.
✓ Branch 5 taken 2106 times.
✓ Branch 6 taken 578138 times.
✓ Branch 7 taken 2106 times.
✓ Branch 8 taken 578138 times.
✓ Branch 9 taken 2106 times.
 2901220 if(globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
197 {
198 values.setAllDirichlet();
199 }
200 else
201 {
202 values.setAllNeumann();
203 }
204
205 580244 return values;
206 }
207
208 /*!
209 * \brief Evaluates the boundary conditions for a Dirichlet Sboundary segment.
210 *
211 * \param globalPos The position for which the bc type should be evaluated
212 */
213 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
214 {
215 570 PrimaryVariables values = initial_(globalPos);
216 ✗ return values;
217 }
218
219 /*!
220 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
221 *
222 * This is the method for the case where the Neumann condition is
223 * potentially solution dependent and requires some quantities that
224 * are specific to the fully-implicit method.
225 *
226 * \param element The finite element
227 * \param fvGeometry The finite-volume geometry
228 * \param elemVolVars All volume variables for the element
229 * \param elemFluxVarsCache The cache related to flux computation
230 * \param scvf The sub-control volume face
231 *
232 * For this method, the \a values parameter stores the flux
233 * in normal direction of each phase. Negative values mean influx.
234 * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
235 */
236 811454 NumEqVector neumann(const Element& element,
237 const FVElementGeometry& fvGeometry,
238 const ElementVolumeVariables& elemVolVars,
239 const ElementFluxVariablesCache& elemFluxVarsCache,
240 const SubControlVolumeFace& scvf) const
241 {
242 811454 NumEqVector flux(0.0);
243
2/2
✓ Branch 0 taken 208722 times.
✓ Branch 1 taken 602732 times.
 811454 const auto& globalPos = scvf.ipGlobal();
244
4/4
✓ Branch 0 taken 208722 times.
✓ Branch 1 taken 602732 times.
✓ Branch 2 taken 71550 times.
✓ Branch 3 taken 471276 times.
 1354280 const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
245
246
4/4
✓ Branch 0 taken 4242 times.
✓ Branch 1 taken 807212 times.
✓ Branch 2 taken 4242 times.
✓ Branch 3 taken 807212 times.
 1622908 if(globalPos[0] < eps_)
247 {
248 10170 flux[contiH2OEqIdx] = -darcyVelocity_*elemVolVars[scv].molarDensity();
249 16524 flux[contiN2EqIdx] = -darcyVelocity_*elemVolVars[scv].molarDensity()*elemVolVars[scv].moleFraction(0, N2Idx);
250 4242 flux[energyEqIdx] = -darcyVelocity_
251 7206 *elemVolVars[scv].density()
252 10170 *IapwsH2O::liquidEnthalpy(temperatureHigh_, elemVolVars[scv].pressure());
253 }
254
255 811454 return flux;
256 }
257
258 // \}
259
260 /*!
261 * \name Volume terms
262 */
263 // \{
264
265 /*!
266 * \brief Evaluates the source term for all phases within a given
267 * sub-control volume.
268 *
269 * For this method, the \a priVars parameter stores the rate mass
270 * of a component is generated or annihilated per volume
271 * unit. Positive values mean that mass is created, negative ones
272 * mean that it vanishes.
273 *
274 * The units must be according to either using mole or mass fractions (mole/(m^3*s) or kg/(m^3*s)).
275 */
276 ✗ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
277 540800 { return NumEqVector(0.0); }
278
279 /*!
280 * \brief Evaluates the initial value for a control volume.
281 *
282 * \param globalPos The position for which the initial condition should be evaluated
283 *
284 * For this method, the \a values parameter stores primary
285 * variables.
286 */
287 ✗ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
288 452 { return initial_(globalPos); }
289
290 // \}
291 private:
292
293 // the internal method for the initial condition
294 ✗ PrimaryVariables initial_(const GlobalPosition &globalPos) const
295 {
296 PrimaryVariables priVars;
297
3/6
✗ Branch 1 not taken.
✓ Branch 2 taken 44 times.
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
 1022 priVars[pressureIdx] = pressureLow_; // initial condition for the pressure
298
3/6
✗ Branch 1 not taken.
✓ Branch 2 taken 44 times.
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
 1022 priVars[N2Idx] = 1e-10; // initial condition for the N2 molefraction
299
6/13
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 44 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 43 times.
✓ Branch 6 taken 44 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 43 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 43 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 43 times.
✗ Branch 15 not taken.
 2044 priVars[temperatureIdx] = temperatureLow_;
300 ✗ return priVars;
301 }
302 static constexpr Scalar eps_ = 1e-6;
303 Scalar temperatureHigh_;
304 Scalar temperatureLow_;
305 Scalar pressureHigh_;
306 Scalar pressureLow_;
307 Scalar darcyVelocity_;
308
309 std::vector<Scalar> temperatureExact_;
310 };
311
312 } // end namespace Dumux
313
314 #endif
315