GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/1p/nonisothermal/problem_conduction.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 45 50 90.0%
Functions: 9 18 50.0%
Branches: 75 138 54.3%
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 conduction problem.
11 *
12 * The simulation domain is a tube with an elevated temperature on the left hand side.
13 */
14
15 #ifndef DUMUX_1PNI_CONDUCTION_PROBLEM_HH
16 #define DUMUX_1PNI_CONDUCTION_PROBLEM_HH
17
18 #include <cmath>
19
20 #include <dumux/common/properties.hh>
21 #include <dumux/common/parameters.hh>
22
23 #include <dumux/common/boundarytypes.hh>
24 #include <dumux/porousmediumflow/problem.hh>
25
26 #include <dumux/material/components/h2o.hh>
27 namespace Dumux {
28
29 /*!
30 * \ingroup OnePTests
31 * \brief Test for the OnePModel in combination with the NI model for a conduction problem.
32 *
33 * The simulation domain is a tube with an elevated temperature on the left hand side.
34 *
35 * Initially the domain is fully saturated with water at a constant temperature.
36 * On the left hand side there is a Dirichlet boundary condition with an increased
37 * temperature and on the right hand side a Dirichlet boundary with constant pressure,
38 * saturation and temperature is applied.
39 *
40 * The results are compared to an analytical solution for a diffusion process.
41 * This problem uses the \ref OnePModel and \ref NIModel model.
42 */
43 template <class TypeTag>
44 class OnePNIConductionProblem : public PorousMediumFlowProblem<TypeTag>
45 {
46 using ParentType = PorousMediumFlowProblem<TypeTag>;
47 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
48 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
49 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
50 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
51 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
52 using ThermalConductivityModel = GetPropType<TypeTag, Properties::ThermalConductivityModel>;
53 using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
54 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
55 using IapwsH2O = Components::H2O<Scalar>;
56
57 enum { dimWorld = GridView::dimensionworld };
58
59 // copy some indices for convenience
60 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
61 enum {
62 // indices of the primary variables
63 pressureIdx = Indices::pressureIdx,
64 temperatureIdx = Indices::temperatureIdx
65 };
66
67 using Element = typename GridView::template Codim<0>::Entity;
68 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
69 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
70
71 public:
72 3 OnePNIConductionProblem(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup)
73
4/14
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 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.
 9 : ParentType(gridGeometry, paramGroup)
74 {
75 //initialize fluid system
76 3 FluidSystem::init();
77
78
2/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
 3 name_ = getParam<std::string>("Problem.Name");
79 3 temperatureHigh_ = 300.0;
80
3/6
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
 7 temperatureExact_.resize(gridGeometry->numDofs());
81 3 }
82
83 //! Get the analytical temperature
84 const std::vector<Scalar>& getExactTemperature()
85 {
86
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 return temperatureExact_;
87 }
88
89 /*!
90 * \brief Update the analytical temperature
91 * The results are compared to an analytical solution for a diffusion process:
92 \f[
93 T =T_{high} + (T_{init} - T_{high})erf \left(0.5\sqrt{\frac{x^2 S_{total}}{t \lambda_{eff}}}\right)
94 \f]
95 */
96 72 void updateExactTemperature(const SolutionVector& curSol, Scalar time)
97 {
98 216 const auto someElement = *(elements(this->gridGeometry().gridView()).begin());
99
100 144 auto someElemSol = elementSolution(someElement, curSol, this->gridGeometry());
101
1/2
✓ Branch 2 taken 72 times.
✗ Branch 3 not taken.
 72 const auto someInitSol = initialAtPos(someElement.geometry().center());
102
103
2/4
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
 288 const auto someFvGeometry = localView(this->gridGeometry()).bindElement(someElement);
104
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 const auto someScv = *(scvs(someFvGeometry).begin());
105
106
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 VolumeVariables volVars;
107
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 volVars.update(someElemSol, *this, someElement, someScv);
108
109
2/4
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
 144 const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
110
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 const auto densityW = volVars.density();
111
3/6
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
 216 const auto heatCapacityW = IapwsH2O::liquidHeatCapacity(someInitSol[temperatureIdx], someInitSol[pressureIdx]);
112
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 const auto densityS = volVars.solidDensity();
113
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 const auto heatCapacityS = volVars.solidHeatCapacity();
114 72 const auto storage = densityW*heatCapacityW*porosity + densityS*heatCapacityS*(1 - porosity);
115
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 const auto effectiveThermalConductivity = ThermalConductivityModel::effectiveThermalConductivity(volVars);
116 using std::max;
117
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
 72 time = max(time, 1e-10);
118
2/4
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
 144 auto fvGeometry = localView(this->gridGeometry());
119
5/10
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 72 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 14472 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 14400 times.
✗ Branch 14 not taken.
 29088 for (const auto& element : elements(this->gridGeometry().gridView()))
120 {
121
1/2
✓ Branch 1 taken 14400 times.
✗ Branch 2 not taken.
 14400 fvGeometry.bindElement(element);
122
4/4
✓ Branch 0 taken 28800 times.
✓ Branch 1 taken 14400 times.
✓ Branch 2 taken 19200 times.
✓ Branch 3 taken 4800 times.
 67200 for (auto&& scv : scvs(fvGeometry))
123 {
124 28800 auto globalIdx = scv.dofIndex();
125 28800 const auto& globalPos = scv.dofPosition();
126 using std::erf;
127 using std::sqrt;
128 57600 temperatureExact_[globalIdx] = temperatureHigh_ + (someInitSol[temperatureIdx] - temperatureHigh_)
129 86400 *erf(0.5*sqrt(globalPos[0]*globalPos[0]*storage/time/effectiveThermalConductivity));
130
131 }
132 }
133 72 }
134
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 3 times.
✗ Branch 2 not taken.
 3 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 bc type should be evaluated
161 */
162 768308 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
163 {
164 768308 BoundaryTypes bcTypes;
165
166
14/14
✓ Branch 0 taken 767204 times.
✓ Branch 1 taken 1104 times.
✓ Branch 2 taken 767204 times.
✓ Branch 3 taken 1104 times.
✓ Branch 4 taken 766100 times.
✓ Branch 5 taken 1104 times.
✓ Branch 6 taken 766100 times.
✓ Branch 7 taken 1104 times.
✓ Branch 8 taken 766100 times.
✓ Branch 9 taken 1104 times.
✓ Branch 10 taken 766100 times.
✓ Branch 11 taken 1104 times.
✓ Branch 12 taken 766100 times.
✓ Branch 13 taken 1104 times.
 1536616 if(globalPos[0] < eps_ || globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
167 bcTypes.setAllDirichlet();
168 else
169 bcTypes.setAllNeumann();
170
171 768308 return bcTypes;
172 }
173
174 /*!
175 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
176 *
177 * \param globalPos The position for which the bc type should be evaluated
178 *
179 * For this method, the \a values parameter stores primary variables.
180 */
181 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
182 {
183 1376 PrimaryVariables priVars(initial_());
184
4/8
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 208 times.
✓ Branch 5 taken 208 times.
✓ Branch 6 taken 136 times.
✓ Branch 7 taken 136 times.
 688 if (globalPos[0] < eps_)
185 688 priVars[temperatureIdx] = temperatureHigh_;
186 ✗ return priVars;
187 }
188
189 // \}
190
191 /*!
192 * \name Volume terms
193 */
194 // \{
195
196 /*!
197 * \brief Evaluates the initial value for a control volume.
198 *
199 * \param globalPos The position for which the initial condition should be evaluated
200 *
201 * For this method, the \a values parameter stores primary
202 * variables.
203 */
204 ✗ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
205 {
206 1748 return initial_();
207 }
208
209 // \}
210
211 private:
212 // the internal method for the initial condition
213 ✗ PrimaryVariables initial_() const
214 {
215
1/2
✓ Branch 2 taken 72 times.
✗ Branch 3 not taken.
 1562 PrimaryVariables priVars(0.0);
216
6/12
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 208 times.
✓ Branch 5 taken 208 times.
✓ Branch 6 taken 136 times.
✓ Branch 7 taken 136 times.
✓ Branch 8 taken 24 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 48 times.
✗ Branch 11 not taken.
 1562 priVars[pressureIdx] = 1.0e5;
217
11/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 208 times.
✓ Branch 9 taken 208 times.
✓ Branch 10 taken 208 times.
✓ Branch 11 taken 208 times.
✓ Branch 12 taken 136 times.
✓ Branch 13 taken 136 times.
✓ Branch 14 taken 136 times.
✓ Branch 15 taken 160 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 24 times.
✓ Branch 19 taken 48 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 48 times.
✗ Branch 23 not taken.
 3124 priVars[temperatureIdx] = 290.0;
218 ✗ return priVars;
219 }
220
221 Scalar temperatureHigh_;
222 static constexpr Scalar eps_ = 1e-6;
223 std::string name_;
224 std::vector<Scalar> temperatureExact_;
225 };
226
227 } // end namespace Dumux
228
229 #endif
230