GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/richards/nonisothermal/convection/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 55 61 90.2%
Functions: 8 16 50.0%
Branches: 64 116 55.2%
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 RichardsTests
10 * \brief Test for the RichardsModel in combination with the NI model for a convection problem:
11 * The simulation domain is a tube where water with an elevated temperature is injected
12 * at a constant rate on the left hand side.
13 */
14
15 #ifndef DUMUX_RICHARDS_CONVECTION_PROBLEM_HH
16 #define DUMUX_RICHARDS_CONVECTION_PROBLEM_HH
17
18 #include <cmath>
19
20 #include <dumux/common/properties.hh>
21 #include <dumux/common/parameters.hh>
22 #include <dumux/common/boundarytypes.hh>
23 #include <dumux/common/numeqvector.hh>
24
25 #include <dumux/porousmediumflow/problem.hh>
26
27 #include <dumux/material/components/h2o.hh>
28
29 namespace Dumux {
30
31 /*!
32 * \ingroup RichardsTests
33 *
34 * \brief Test for the RichardsModel in combination with the NI model for a convection problem:
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 results are compared to an analytical solution where a retarded front velocity is calculated as follows:
43 \f[
44 v_{Front}=\frac{q S_{water}}{\phi S_{total}}
45 \f]
46 *
47 * The result of the analytical solution is written into the vtu files.
48 * This problem uses the \ref RichardsModel and \ref NIModel model.
49 *
50 * To run the simulation execute the following line in shell: <br>
51 * <tt>./test_boxrichardsniconvection -ParameterFile ./test_boxrichardsniconvection.input</tt> or <br>
52 * <tt>./test_ccrichardsniconvection -ParameterFile ./test_ccrichardsniconvection.input</tt>
53 */
54 template <class TypeTag>
55 class RichardsNIConvectionProblem : public PorousMediumFlowProblem<TypeTag>
56 {
57 using ParentType = PorousMediumFlowProblem<TypeTag>;
58
59 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
60 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
61 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
62 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
63 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
64 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
65 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
66 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
67 using ThermalConductivityModel = GetPropType<TypeTag, Properties::ThermalConductivityModel>;
68 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
69 using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
70 using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
71 using VolumeVariables = typename GridVariables::GridVolumeVariables::VolumeVariables;
72 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
73 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
74 using IapwsH2O = Components::H2O<Scalar>;
75
76 // copy some indices for convenience
77 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
78 enum { dimWorld = GridView::dimensionworld };
79
80 enum {
81 pressureIdx = Indices::pressureIdx,
82 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
83 temperatureIdx = Indices::temperatureIdx
84 };
85
86 using Element = typename GridView::template Codim<0>::Entity;
87
88 using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
89
90 public:
91 2 RichardsNIConvectionProblem(std::shared_ptr<const GridGeometry> gridGeometry)
92
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)
93 {
94 // initialize fluid system
95
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 FluidSystem::init();
96
97
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");
98
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 darcyVelocity_ = getParam<Scalar>("Problem.DarcyVelocity");
99 2 temperatureHigh_ = 291.;
100 2 temperatureLow_ = 290.;
101 2 pressureHigh_ = 2e5;
102 2 pressureLow_ = 1e5;
103
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(gridGeometry->numDofs());
104 2 }
105
106 //! Get the analytical temperature
107 const std::vector<Scalar>& getExactTemperature()
108 {
109
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return temperatureExact_;
110 }
111
112 //! Update the analytical temperature
113 87 void updateExactTemperature(const SolutionVector& curSol, Scalar time)
114 {
115 261 const auto someElement = *(elements(this->gridGeometry().gridView()).begin());
116
117 174 const auto someElemSol = elementSolution(someElement, curSol, this->gridGeometry());
118 87 const auto someInitSol = initialAtPos(someElement.geometry().center());
119
120
2/4
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87 times.
✗ Branch 5 not taken.
 348 const auto someFvGeometry = localView(this->gridGeometry()).bindElement(someElement);
121
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 87 const auto someScv = *(scvs(someFvGeometry).begin());
122
123
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 87 VolumeVariables volVars;
124
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 87 volVars.update(someElemSol, *this, someElement, someScv);
125
126
2/4
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87 times.
✗ Branch 5 not taken.
 174 const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
127
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 87 const auto densityW = volVars.density(liquidPhaseIdx);
128
3/6
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 87 times.
✗ Branch 8 not taken.
 261 const auto heatCapacityW = IapwsH2O::liquidHeatCapacity(someInitSol[temperatureIdx], someInitSol[pressureIdx]);
129
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 87 const auto densityS = volVars.solidDensity();
130
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 87 const auto heatCapacityS = volVars.solidHeatCapacity();
131 87 const auto storage = densityW*heatCapacityW*porosity + densityS*heatCapacityS*(1 - porosity);
132
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 87 const auto effectiveThermalConductivity = ThermalConductivityModel::effectiveThermalConductivity(volVars);
133 using std::max;
134
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 87 times.
 87 time = max(time, 1e-10);
135
2/4
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87 times.
✗ Branch 5 not taken.
 174 auto fvGeometry = localView(this->gridGeometry());
136
5/10
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 87 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 7047 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 6960 times.
✗ Branch 14 not taken.
 14268 for (const auto& element : elements(this->gridGeometry().gridView()))
137 {
138
1/2
✓ Branch 1 taken 6960 times.
✗ Branch 2 not taken.
 6960 fvGeometry.bindElement(element);
139
4/4
✓ Branch 0 taken 17520 times.
✓ Branch 1 taken 6960 times.
✓ Branch 2 taken 14080 times.
✓ Branch 3 taken 3520 times.
 42080 for (auto&& scv : scvs(fvGeometry))
140 {
141 17520 auto globalIdx = scv.dofIndex();
142 17520 const auto& globalPos = scv.dofPosition();
143 using std::erf;
144 using std::sqrt;
145 35040 temperatureExact_[globalIdx] = temperatureHigh_ + (someInitSol[temperatureIdx] - temperatureHigh_)
146 52560 *erf(0.5*sqrt(globalPos[0]*globalPos[0]*storage/time/effectiveThermalConductivity));
147
148 }
149 }
150 87 }
151 /*!
152 * \name Problem parameters
153 */
154 // \{
155
156 /*!
157 * \brief The problem name.
158 *
159 * This is used as a prefix for files generated by the simulation.
160 */
161 const std::string& name() const
162 {
163
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return name_;
164 }
165
166 // \}
167
168 /*!
169 * \name Boundary conditions
170 */
171 // \{
172
173 /*!
174 * \brief Specifies which kind of boundary condition should be
175 * used for which equation on a given boundary segment.
176 *
177 * \param globalPos The position for which the boundary type is set
178 */
179 104848 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
180 {
181 104848 BoundaryTypes values;
182
10/10
✓ Branch 0 taken 104198 times.
✓ Branch 1 taken 650 times.
✓ Branch 2 taken 104198 times.
✓ Branch 3 taken 650 times.
✓ Branch 4 taken 104198 times.
✓ Branch 5 taken 650 times.
✓ Branch 6 taken 104198 times.
✓ Branch 7 taken 650 times.
✓ Branch 8 taken 104198 times.
✓ Branch 9 taken 650 times.
 524240 if(globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
183 {
184 values.setAllDirichlet();
185 }
186 else
187 {
188 values.setAllNeumann();
189 }
190 104848 return values;
191 }
192
193 /*!
194 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
195 *
196 * \param globalPos The position for which the bc type should be evaluated
197 *
198 * For this method, the \a values parameter stores primary variables.
199 */
200 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
201 {
202 664 return initial_(globalPos);
203 }
204
205 /*!
206 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
207 *
208 * \param element The finite element
209 * \param fvGeometry The finite-volume geometry in the box scheme
210 * \param elemVolVars The element volume variables
211 * \param elemFluxVarsCache Flux variables caches for all faces in stencil
212 * \param scvf The sub-control volume face
213 * Negative values mean influx.
214 */
215 376638 NumEqVector neumann(const Element &element,
216 const FVElementGeometry& fvGeometry,
217 const ElementVolumeVariables& elemVolVars,
218 const ElementFluxVariablesCache& elemFluxVarsCache,
219 const SubControlVolumeFace& scvf) const
220 {
221 376638 NumEqVector values(0.0);
222 376638 const auto globalPos = scvf.ipGlobal();
223
224
4/4
✓ Branch 0 taken 2352 times.
✓ Branch 1 taken 374286 times.
✓ Branch 2 taken 2352 times.
✓ Branch 3 taken 374286 times.
 753276 if(globalPos[0] < eps_)
225 {
226 8772 values[pressureIdx] = -darcyVelocity_*elemVolVars[scvf.insideScvIdx()].density(liquidPhaseIdx);
227 7056 values[temperatureIdx] = -darcyVelocity_*elemVolVars[scvf.insideScvIdx()].density(liquidPhaseIdx)
228 8772 *IapwsH2O::liquidEnthalpy(temperatureHigh_, elemVolVars[scvf.insideScvIdx()].pressure(liquidPhaseIdx));
229 }
230 376638 return values;
231 }
232
233 // \}
234
235 /*!
236 * \name Volume terms
237 */
238 // \{
239
240
241 /*!
242 * \brief Returns the reference pressure [Pa] of the nonwetting
243 * fluid phase within a finite volume.
244 *
245 * This problem assumes a constant reference pressure of 1 bar.
246 */
247 ✗ Scalar nonwettingReferencePressure() const
248 ✗ { return 1e5; };
249
250 /*!
251 * \brief Evaluates the initial value for a control volume.
252 *
253 * \param globalPos The position for which the initial condition should be evaluated
254 *
255 * For this method, the \a values parameter stores primary
256 * variables.
257 */
258 ✗ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
259 {
260 658 return initial_(globalPos);
261 }
262
263 // \}
264
265 private:
266 ✗ PrimaryVariables initial_(const GlobalPosition &globalPos) const
267 {
268 661 PrimaryVariables priVars(0.0);
269
2/3
✓ Branch 1 taken 44 times.
✓ Branch 2 taken 43 times.
✗ Branch 3 not taken.
 419 priVars[pressureIdx] = pressureLow_; // initial condition for the pressure
270
4/7
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 43 times.
✓ Branch 4 taken 44 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
 838 priVars[temperatureIdx] = temperatureLow_;
271 ✗ return priVars;
272 }
273
274 Scalar temperatureHigh_;
275 Scalar temperatureLow_;
276 Scalar pressureHigh_;
277 Scalar pressureLow_;
278 Scalar darcyVelocity_;
279 static constexpr Scalar eps_ = 1e-6;
280 std::string name_;
281 std::vector<Scalar> temperatureExact_;
282 };
283
284 } // end namespace
285 #endif
286