GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/test/porousmediumflow/3p/conduction/problem.hh
Date: 2026-05-16 22:02:50
Exec Total Coverage
Lines: 54 54 100.0%
Functions: 6 6 100.0%
Branches: 47 80 58.8%
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-FileCopyrightText: 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 ThreePTests
10 * \brief Definition of a 3pni problem:
11 * Component transport of nitrogen dissolved in the water phase.
12 */
13
14 #ifndef DUMUX_3PNI_CONDUCTION_PROBLEM_HH
15 #define DUMUX_3PNI_CONDUCTION_PROBLEM_HH
16
17 #include <cmath>
18 #include <algorithm>
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 #include <dumux/material/components/h2o.hh>
27
28 namespace Dumux {
29
30 /*!
31 * \ingroup ThreePTests
32 *
33 * \brief Test for the ThreePModel in combination with the NI model for a conduction problem.
34 *
35 * The simulation domain is a tube where with an elevated temperature on the left hand side.
36 *
37 * Initially the domain is fully saturated with water at a constant temperature.
38 * On the left hand side there is a Dirichlet boundary condition with an increased
39 * temperature and on the right hand side a Dirichlet boundary with constant
40 * pressure, saturation and temperature is applied.
41 *
42 * The results are compared to an analytical solution for a diffusion process:
43 \f[
44 T =T_{high} + (T_{init} - T_{high})erf \left(0.5\sqrt{\frac{x^2 S_{total}}{t \lambda_{eff}}}\right)
45 \f]
46 * This problem uses the \ref ThreePModel and \ref NIModel model.
47 */
48 template <class TypeTag>
49 class ThreePNIConductionProblem : 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 ThermalConductivityModel = GetPropType<TypeTag, Properties::ThermalConductivityModel>;
60 using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
61 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
62 using IapwsH2O = Components::H2O<Scalar>;
63 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
64
65 // copy some indices for convenience
66 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
67 enum {
68 // index of the primary variables
69 pressureIdx = Indices::pressureIdx,
70 swIdx = Indices::swIdx,
71 snIdx = Indices::snIdx,
72 temperatureIdx = Indices::temperatureIdx,
73 wPhaseIdx = FluidSystem::wPhaseIdx
74 };
75
76 enum { dimWorld = GridView::dimensionworld };
77
78 using Element = typename GridView::template Codim<0>::Entity;
79 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
80
81 public:
82 2 ThreePNIConductionProblem(std::shared_ptr<const GridGeometry> gridGeometry)
83
3/6
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
 6 : ParentType(gridGeometry)
84 {
85 //initialize fluid system
86 FluidSystem::init();
87
88
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 name_ = getParam<std::string>("Problem.Name");
89 2 temperatureHigh_ = 300.0;
90
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 temperatureExact_.resize(gridGeometry->numDofs());
91 2 }
92
93 //! Get the analytical temperature
94 2 const std::vector<Scalar>& getExactTemperature()
95 {
96
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return temperatureExact_;
97 }
98
99 //! Update the analytical temperature
100 48 void updateExactTemperature(const SolutionVector& curSol, Scalar time)
101 {
102 48 const auto someElement = *(elements(this->gridGeometry().gridView()).begin());
103
104
1/2
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
 48 const auto someElemSol = elementSolution(someElement, curSol, this->gridGeometry());
105
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 const auto someInitSol = initialAtPos(someElement.geometry().center());
106
107
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 96 const auto someFvGeometry = localView(this->gridGeometry()).bindElement(someElement);
108 48 const auto someScv = *(scvs(someFvGeometry).begin());
109
110 48 VolumeVariables volVars;
111
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 volVars.update(someElemSol, *this, someElement, someScv);
112
113
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 const auto porosity = this->spatialParams().porosity(someElement, someScv, someElemSol);
114
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 const auto densityW = volVars.density(wPhaseIdx);
115
2/4
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
 48 const auto heatCapacityW = IapwsH2O::liquidHeatCapacity(someInitSol[temperatureIdx], someInitSol[pressureIdx]);
116
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 const auto densityS = volVars.solidDensity();
117 48 const auto heatCapacityS = volVars.solidHeatCapacity();
118 48 const auto storage = densityW*heatCapacityW*porosity + densityS*heatCapacityS*(1 - porosity);
119
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 const auto effectiveThermalConductivity = ThermalConductivityModel::effectiveThermalConductivity(volVars);
120 using std::max;
121
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
 48 time = max(time, 1e-10);
122
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 auto fvGeometry = localView(this->gridGeometry());
123
3/6
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9648 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9600 times.
✗ Branch 8 not taken.
 28848 for (const auto& element : elements(this->gridGeometry().gridView()))
124 {
125
1/2
✓ Branch 1 taken 4800 times.
✗ Branch 2 not taken.
 14400 fvGeometry.bindElement(element);
126
2/2
✓ Branch 0 taken 24000 times.
✓ Branch 1 taken 9600 times.
 33600 for (auto&& scv : scvs(fvGeometry))
127 {
128 43200 auto globalIdx = scv.dofIndex();
129 24000 const auto& globalPos = scv.dofPosition();
130 using std::erf;
131 using std::sqrt;
132 24000 temperatureExact_[globalIdx] = temperatureHigh_ + (someInitSol[temperatureIdx] - temperatureHigh_)
133 24000 *erf(0.5*sqrt(globalPos[0]*globalPos[0]*storage/time/effectiveThermalConductivity));
134
135 }
136 }
137 48 }
138
139 /*!
140 * \name Problem parameters
141 */
142 // \{
143
144 /*!
145 * \brief The problem name.
146 *
147 * This is used as a prefix for files generated by the simulation.
148 */
149 2 const std::string& name() const
150 {
151
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return name_;
152 }
153 // \}
154
155 /*!
156 * \name Boundary conditions
157 */
158 // \{
159
160 /*!
161 * \brief Specifies which kind of boundary condition should be
162 * used for which equation on a given boundary segment.
163 *
164 * \param globalPos The position for which the bc type should be evaluated
165 */
166 228876 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
167 {
168
2/2
✓ Branch 0 taken 228306 times.
✓ Branch 1 taken 570 times.
 228876 BoundaryTypes values;
169
6/6
✓ Branch 0 taken 228306 times.
✓ Branch 1 taken 570 times.
✓ Branch 2 taken 570 times.
✓ Branch 3 taken 227736 times.
✓ Branch 4 taken 227736 times.
✓ Branch 5 taken 1140 times.
 230016 if(globalPos[0] < eps_ || globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
170 {
171 228876 values.setAllDirichlet();
172 }
173 else
174 {
175 228876 values.setAllNeumann();
176 }
177 228876 return values;
178 }
179
180 /*!
181 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
182 *
183 * \param globalPos The position for which the bc type should be evaluated
184 */
185 408 PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
186 {
187
2/2
✓ Branch 0 taken 72 times.
✓ Branch 1 taken 72 times.
 408 PrimaryVariables values = initialAtPos(globalPos);
188
189
2/4
✓ Branch 0 taken 204 times.
✓ Branch 1 taken 204 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 408 if (globalPos[0] < eps_)
190 204 values[temperatureIdx] = temperatureHigh_;
191 return values;
192 }
193
194 /*!
195 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
196 *
197 * \param globalPos The position of the integration point of the boundary segment.
198 *
199 * For this method, the \a values parameter stores the mass flux
200 * in normal direction of each phase. Negative values mean influx.
201 */
202 1032312 NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
203 {
204
1/2
✓ Branch 0 taken 132000 times.
✗ Branch 1 not taken.
 2064624 return NumEqVector(0.0);
205 }
206
207 // \}
208
209 /*!
210 * \name Volume terms
211 */
212 // \{
213
214 /*!
215 * \brief Evaluates the source term for all phases within a given
216 * sub-controlvolume.
217 *
218 * \param globalPos The position for which the source should be evaluated
219 *
220 * Returns the rate mass of a component is generated or annihilated
221 * per volume unit. Positive values mean that mass is created,
222 * negative ones mean that it vanishes.
223 *
224 * The units must be according to either using mole or mass fractions (mole/(m^3*s) or kg/(m^3*s)).
225 */
226 963600 NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
227 {
228
2/2
✓ Branch 0 taken 3854400 times.
✓ Branch 1 taken 963600 times.
 5781600 return NumEqVector(0.0);
229 }
230
231 /*!
232 * \brief Evaluates the initial value for a control volume.
233 *
234 * \param globalPos The position for which the initial condition should be evaluated
235 *
236 */
237
4/8
✓ Branch 0 taken 204 times.
✓ Branch 1 taken 204 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 6 taken 24 times.
✗ Branch 7 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
 1058 PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
238 {
239 PrimaryVariables values;
240 264 values[pressureIdx] = 1e5; // initial condition for the pressure
241 264 values[swIdx] = 1.0; // initial condition for the wetting phase saturation
242 264 values[snIdx] = 1e-5; // initial condition for the nonwetting phase saturation
243
2/4
✓ Branch 0 taken 132 times.
✓ Branch 1 taken 132 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 264 values[temperatureIdx] = 290;
244 return values;
245 }
246
247 // \}
248
249 private:
250 Scalar temperatureHigh_;
251 static constexpr Scalar eps_ = 1e-6;
252 std::string name_;
253 std::vector<Scalar> temperatureExact_;
254 };
255
256 } // end namespace Dumux
257 #endif
258