GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/richards/nonisothermal/conduction/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 49 56 87.5%
Functions: 6 16 37.5%
Branches: 74 152 48.7%
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 conduction problem:
11 * The simulation domain is a tube with an elevated temperature on the left hand side.
12 */
13
14 #ifndef DUMUX_RICHARDS_CONDUCTION_PROBLEM_HH
15 #define DUMUX_RICHARDS_CONDUCTION_PROBLEM_HH
16
17 #include <cmath>
18
19 #include <dumux/common/properties.hh>
20 #include <dumux/common/parameters.hh>
21 #include <dumux/common/boundarytypes.hh>
22 #include <dumux/common/numeqvector.hh>
23 #include <dumux/discretization/elementsolution.hh>
24
25 #include <dumux/porousmediumflow/problem.hh>
26 #include <dumux/material/components/h2o.hh>
27
28 namespace Dumux {
29
30 /*!
31 * \ingroup RichardsTests
32 *
33 * \brief Test for the RichardsModel in combination with the NI model for a conduction problem:
34 * The simulation domain is a tube with an elevated temperature on the left hand side.
35 *
36 * Initially the domain is fully saturated with water at a constant temperature.
37 * On the left hand side there is a Dirichlet boundary condition with an increased temperature
38 * and on the right hand side a Dirichlet boundary with constant pressure, saturation
39 * and temperature is applied.
40 *
41 * The results are compared to an analytical solution for a diffusion process:
42 \f[
43 T =T_{high} + (T_{init} - T_{high})erf \left(0.5\sqrt{\frac{x^2 S_{total}}{t \lambda_{eff}}}\right)
44 \f]
45 *
46 * This problem uses the \ref RichardsModel and \ref NIModel model.
47 */
48 template <class TypeTag>
49 class RichardsNIConductionProblem :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 GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
56 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
57 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
58 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
59 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
60 using ThermalConductivityModel = GetPropType<TypeTag, Properties::ThermalConductivityModel>;
61 using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
62 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
63 using IapwsH2O = Components::H2O<Scalar>;
64
65 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
66 enum { dimWorld = GridView::dimensionworld };
67
68 enum {
69 pressureIdx = Indices::pressureIdx,
70 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
71 temperatureIdx = Indices::temperatureIdx
72 };
73
74 using Element = typename GridView::template Codim<0>::Entity;
75
76 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
77
78 public:
79 2 RichardsNIConductionProblem(std::shared_ptr<const GridGeometry> gridGeometry)
80
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)
81 {
82 //initialize fluid system
83
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 FluidSystem::init();
84
85
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");
86 2 temperatureHigh_ = 300.;
87
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());
88 2 }
89
90 //! Get the analytical temperature
91 const std::vector<Scalar>& getExactTemperature()
92 {
93
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return temperatureExact_;
94 }
95
96 //! Update the analytical temperature
97 48 void updateExactTemperature(const SolutionVector& curSol, Scalar time)
98 {
99 144 const auto someElement = *(elements(this->gridGeometry().gridView()).begin());
100
101 96 const auto someElemSol = elementSolution(someElement, curSol, this->gridGeometry());
102
1/2
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
 48 const auto someInitSol = initialAtPos(someElement.geometry().center());
103
104
2/4
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
 192 const auto someFvGeometry = localView(this->gridGeometry()).bindElement(someElement);
105
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 const auto someScv = *(scvs(someFvGeometry).begin());
106
107
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 VolumeVariables volVars;
108
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 volVars.update(someElemSol, *this, someElement, someScv);
109
110
2/4
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
 96 const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
111
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 const auto densityW = volVars.density(liquidPhaseIdx);
112
3/6
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 48 times.
✗ Branch 8 not taken.
 144 const auto heatCapacityW = IapwsH2O::liquidHeatCapacity(someInitSol[temperatureIdx], someInitSol[pressureIdx]);
113
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 const auto densityS =volVars.solidDensity();
114
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 const auto heatCapacityS = volVars.solidHeatCapacity();
115 48 const auto storage = densityW*heatCapacityW*porosity + densityS*heatCapacityS*(1 - porosity);
116
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 const auto effectiveThermalConductivity = ThermalConductivityModel::effectiveThermalConductivity(volVars);
117 using std::max;
118
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
 48 time = max(time, 1e-10);
119
2/4
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
 96 auto fvGeometry = localView(this->gridGeometry());
120
5/10
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 48 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 9648 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 9600 times.
✗ Branch 14 not taken.
 19392 for (const auto& element : elements(this->gridGeometry().gridView()))
121 {
122
1/2
✓ Branch 1 taken 9600 times.
✗ Branch 2 not taken.
 9600 fvGeometry.bindElement(element);
123
4/4
✓ Branch 0 taken 24000 times.
✓ Branch 1 taken 9600 times.
✓ Branch 2 taken 19200 times.
✓ Branch 3 taken 4800 times.
 57600 for (auto&& scv : scvs(fvGeometry))
124 {
125 24000 auto globalIdx = scv.dofIndex();
126 24000 const auto& globalPos = scv.dofPosition();
127 using std::erf;
128 using std::sqrt;
129 48000 temperatureExact_[globalIdx] = temperatureHigh_ + (someInitSol[temperatureIdx] - temperatureHigh_)
130 72000 *erf(0.5*sqrt(globalPos[0]*globalPos[0]*storage/time/effectiveThermalConductivity));
131
132 }
133 }
134 48 }
135 /*!
136 * \name Problem parameters
137 */
138 // \{
139
140 /*!
141 * \brief The problem name.
142 *
143 * This is used as a prefix for files generated by the simulation.
144 */
145 const std::string& name() const
146 {
147
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return name_;
148 }
149 // \}
150
151 /*!
152 * \name Boundary conditions
153 */
154 // \{
155
156 /*!
157 * \brief Specifies which kind of boundary condition should be
158 * used for which equation on a given boundary segment.
159 *
160 * \param globalPos The position for which the boundary type is set
161 */
162 189078 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
163 {
164 189078 BoundaryTypes values;
165
14/14
✓ Branch 0 taken 188607 times.
✓ Branch 1 taken 471 times.
✓ Branch 2 taken 188607 times.
✓ Branch 3 taken 471 times.
✓ Branch 4 taken 188136 times.
✓ Branch 5 taken 471 times.
✓ Branch 6 taken 188136 times.
✓ Branch 7 taken 471 times.
✓ Branch 8 taken 188136 times.
✓ Branch 9 taken 471 times.
✓ Branch 10 taken 188136 times.
✓ Branch 11 taken 471 times.
✓ Branch 12 taken 188136 times.
✓ Branch 13 taken 471 times.
 378156 if(globalPos[0] < eps_ || globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
166 {
167 values.setAllDirichlet();
168 }
169 else
170 {
171 values.setAllNeumann();
172 }
173 189078 return values;
174 }
175
176 /*!
177 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
178 *
179 * \param globalPos The position for which the bc type should be evaluated
180 *
181 * For this method, the \a values parameter stores primary variables.
182 */
183 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
184 {
185 446 PrimaryVariables values(0.0);
186 892 values = initial_(globalPos);
187
188
2/8
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 223 times.
✓ Branch 5 taken 223 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 446 if (globalPos[0] < eps_)
189 {
190 446 values[temperatureIdx] = temperatureHigh_;
191 }
192 446 return values;
193 }
194
195 /*!
196 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
197 *
198 * \param globalPos The global position where we evaluate
199 *
200 * For this method, the \a priVars parameter stores the mass flux
201 * in normal direction of each component. Negative values mean
202 * influx.
203 */
204 ✗ NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
205 {
206 562024 NumEqVector values(0.0);
207
1/2
✓ Branch 0 taken 79200 times.
✗ Branch 1 not taken.
 562024 return values;
208 }
209
210 // \}
211
212 /*!
213 * \name Volume terms
214 */
215 // \{
216
217 /*!
218 * \brief Returns the reference pressure [Pa] of the nonwetting
219 * fluid phase within a finite volume.
220 *
221 * This problem assumes a constant reference pressure of 1 bar.
222 */
223 ✗ Scalar nonwettingReferencePressure() const
224 ✗ { return 1e5; };
225
226 /*!
227 * \brief Evaluates the initial value for a control volume.
228 *
229 * \param globalPos The position for which the initial condition should be evaluated
230 *
231 * For this method, the \a values parameter stores primary
232 * variables.
233 */
234 ✗ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
235 {
236 1300 return initial_(globalPos);
237 }
238
239 // \}
240
241 private:
242 ✗ PrimaryVariables initial_(const GlobalPosition &globalPos) const
243 {
244
1/2
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
 1096 PrimaryVariables priVars(0.0);
245
4/12
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 223 times.
✓ Branch 5 taken 223 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 24 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 24 times.
✗ Branch 11 not taken.
 1096 priVars[pressureIdx] = 1e5; // initial condition for the pressure
246
247
8/22
✗ 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 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 223 times.
✓ Branch 9 taken 223 times.
✓ Branch 10 taken 223 times.
✓ Branch 11 taken 223 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 24 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 24 times.
✓ Branch 19 taken 24 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 24 times.
✗ Branch 23 not taken.
 2192 priVars[temperatureIdx] = 290.;
248 ✗ return priVars;
249 }
250
251 Scalar temperatureHigh_;
252 static constexpr Scalar eps_ = 1e-6;
253 std::string name_;
254 int outputInterval_;
255 std::vector<Scalar> temperatureExact_;
256 };
257
258 } // end namespace Dumux
259 #endif
260