GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/3p3c/columnxylol/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 66 70 94.3%
Functions: 4 6 66.7%
Branches: 93 130 71.5%
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 ThreePThreeCTests
10 * \brief Non-isothermal injection problem where water is injected into a
11 * sand column with a NAPL contamination.
12 */
13
14 #ifndef DUMUX_COLUMNXYLOLPROBLEM_HH
15 #define DUMUX_COLUMNXYLOLPROBLEM_HH
16
17 #include <dumux/common/parameters.hh>
18 #include <dumux/common/properties.hh>
19 #include <dumux/common/boundarytypes.hh>
20 #include <dumux/common/numeqvector.hh>
21
22 #include <dumux/porousmediumflow/problem.hh>
23
24 namespace Dumux {
25
26 /*!
27 * \ingroup ThreePThreeCTests
28 * \brief Non-isothermal injection problem where a water is injected into a
29 * sand column with a NAPL contamination.
30 *
31 * The 2D domain of this test problem is 0.1m long and 1.2m high.
32 * Initially the column is filled with NAPL, gas and water, the NAPL saturation
33 * increases to the bottom of the columns, the water saturation is constant.
34 * Then water is injected from the top with a rate of 0.395710 mol/(s m) and
35 * an enthalpy of 17452.97 [J/(s m)]. *
36 *
37 * Left, right and top boundaries are Neumann boundaries. Left and right are
38 * no-flow boundaries, on top the injection takes places.
39 * The bottom is a Dirichlet boundary.
40 *
41 * This problem uses the \ref ThreePThreeCModel and \ref NIModel model.
42 */
43 template <class TypeTag >
44 class ColumnProblem : public PorousMediumFlowProblem<TypeTag>
45 {
46 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
47 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
48 using ParentType = PorousMediumFlowProblem<TypeTag>;
49
50 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
51 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
52 using Indices = typename ModelTraits::Indices;
53
54 enum
55 {
56 // primary variable indices
57 pressureIdx = Indices::pressureIdx,
58 switch1Idx = Indices::switch1Idx,
59 switch2Idx = Indices::switch2Idx,
60 temperatureIdx = Indices::temperatureIdx,
61
62 // equation indices
63 energyEqIdx = Indices::energyEqIdx,
64 contiWEqIdx = Indices::conti0EqIdx + FluidSystem::wCompIdx, //!< Index of the mass conservation equation for the water component,
65 contiGEqIdx = Indices::conti0EqIdx + FluidSystem::gCompIdx, //!< Index of the mass conservation equation for the gas component,
66 contiNEqIdx = Indices::conti0EqIdx + FluidSystem::nCompIdx, //!< Index of the mass conservation equation for the contaminant component
67
68 // Phase State
69 threePhases = Indices::threePhases
70 };
71
72 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
73 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
74 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
75 using Element = typename GridView::template Codim<0>::Entity;
76 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
77 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
78 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
79 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
80 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
81
82 public:
83 2 ColumnProblem(std::shared_ptr<const GridGeometry> gridGeometry)
84
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)
85 {
86
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 FluidSystem::init();
87
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");
88 2 }
89
90 /*!
91 * \name Problem parameters
92 */
93 // \{
94
95 /*!
96 * \brief The problem name.
97 *
98 * This is used as a prefix for files generated by the simulation.
99 */
100 const std::string name() const
101
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2 { return name_; }
102
103 // \}
104
105 /*!
106 * \name Boundary conditions
107 */
108 // \{
109
110 /*!
111 * \brief Specifies which kind of boundary condition should be
112 * used for which equation on a given boundary segment.
113 *
114 * \param globalPos The position for which the bc type should be evaluated
115 */
116 ✗ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
117 {
118 ✗ BoundaryTypes bcTypes;
119 ✗ if (globalPos[1] < eps_)
120 bcTypes.setAllDirichlet();
121 else
122 bcTypes.setAllNeumann();
123 ✗ return bcTypes;
124 }
125
126 /*!
127 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
128 *
129 * \param globalPos The position for which the bc type should be evaluated
130 *
131 * For this method, the \a values parameter stores primary variables.
132 */
133 PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
134 {
135 1417 return initial_(globalPos);
136 }
137
138 /*!
139 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
140 *
141 * \param globalPos The position for which the bc type should be evaluated
142 *
143 * For this method, the \a values parameter stores the mass flux
144 * in normal direction of each phase. Negative values mean influx.
145 */
146 NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
147 {
148 8641287 NumEqVector values(0.0);
149
150 // negative values for injection
151
20/30
✗ 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 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 33990 times.
✓ Branch 11 taken 8123610 times.
✓ Branch 12 taken 33990 times.
✓ Branch 13 taken 8123610 times.
✓ Branch 14 taken 33990 times.
✓ Branch 15 taken 8123610 times.
✓ Branch 16 taken 33990 times.
✓ Branch 17 taken 8123610 times.
✓ Branch 18 taken 33990 times.
✓ Branch 19 taken 8123610 times.
✓ Branch 20 taken 2007 times.
✓ Branch 21 taken 481680 times.
✓ Branch 22 taken 2007 times.
✓ Branch 23 taken 481680 times.
✓ Branch 24 taken 2007 times.
✓ Branch 25 taken 481680 times.
✓ Branch 26 taken 2007 times.
✓ Branch 27 taken 481680 times.
✓ Branch 28 taken 2007 times.
✓ Branch 29 taken 481680 times.
 43206435 if (globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_)
152 {
153 35997 values[contiWEqIdx] = -0.395710;
154 35997 values[contiGEqIdx] = -0.000001;
155 35997 values[contiNEqIdx] = -0.00;
156 71994 values[energyEqIdx] = -17452.97;
157 }
158 return values;
159 }
160
161 // \}
162
163
164 /*!
165 * \brief Evaluates the initial value for a control volume.
166 *
167 * \param globalPos The position for which the initial condition should be evaluated
168 *
169 * For this method, the \a values parameter stores primary
170 * variables.
171 */
172 PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
173 {
174 362 return initial_(globalPos);
175 }
176
177
178 /*!
179 * \brief Appends all quantities of interest which can be derived
180 * from the solution of the current time step to the VTK
181 * writer.
182 *
183 * Adjust this in case of anisotropic permeabilities.
184 */
185 template<class VTKWriter>
186 void addVtkFields(VTKWriter& vtk)
187 {
188 const auto& gg = this->gridGeometry();
189 Kxx_.resize(gg.numDofs());
190 vtk.addField(Kxx_, "permeability");
191 auto fvGeometry = localView(gg);
192 for (const auto& element : elements(this->gridView()))
193 {
194 fvGeometry.bindElement(element);
195 for (const auto& scv : scvs(fvGeometry))
196 Kxx_[scv.dofIndex()] = this->spatialParams().intrinsicPermeabilityAtPos(scv.dofPosition());
197 }
198 }
199
200 private:
201 // internal method for the initial condition (reused for the
202 // dirichlet conditions!)
203 1779 PrimaryVariables initial_(const GlobalPosition &globalPos) const
204 {
205
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1777 times.
 1779 PrimaryVariables values;
206
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1777 times.
 1779 values.setState(threePhases);
207
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1777 times.
 1779 const auto y = globalPos[1];
208
6/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1777 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 1777 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 1777 times.
 5337 const auto yMax = this->gridGeometry().bBoxMax()[1];
209
210
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1777 times.
 1779 values[temperatureIdx] = 296.15;
211
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1777 times.
 1779 values[pressureIdx] = 1.e5;
212
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1777 times.
 1779 values[switch1Idx] = 0.005;
213
214
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1777 times.
 1779 if (y > yMax - eps_)
215 4 values[switch2Idx] = 0.112;
216
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1774 times.
 1777 else if (y > yMax - 0.0148 - eps_)
217 6 values[switch2Idx] = 0 + ((yMax - y)/0.0148)*0.112;
218
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1770 times.
 1774 else if (y > yMax - 0.0296 - eps_)
219 8 values[switch2Idx] = 0.112 + (((yMax - y) - 0.0148)/0.0148)*(0.120 - 0.112);
220
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1765 times.
 1770 else if (y > yMax - 0.0444 - eps_)
221 10 values[switch2Idx] = 0.120 + (((yMax - y) - 0.0296)/0.0148)*(0.125 - 0.120);
222
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1761 times.
 1765 else if (y > yMax - 0.0592 - eps_)
223 8 values[switch2Idx] = 0.125 + (((yMax - y) - 0.0444)/0.0148)*(0.137 - 0.125);
224
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1756 times.
 1761 else if (y > yMax - 0.0740 - eps_)
225 10 values[switch2Idx] = 0.137 + (((yMax - y) - 0.0592)/0.0148)*(0.150 - 0.137);
226
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1752 times.
 1756 else if (y > yMax - 0.0888 - eps_)
227 8 values[switch2Idx] = 0.150 + (((yMax - y) - 0.0740)/0.0148)*(0.165 - 0.150);
228
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1747 times.
 1752 else if (y > yMax - 0.1036 - eps_)
229 10 values[switch2Idx] = 0.165 + (((yMax - y) - 0.0888)/0.0148)*(0.182 - 0.165);
230
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1743 times.
 1747 else if (y > yMax - 0.1184 - eps_)
231 8 values[switch2Idx] = 0.182 + (((yMax - y) - 0.1036)/0.0148)*(0.202 - 0.182);
232
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1738 times.
 1743 else if (y > yMax - 0.1332 - eps_)
233 10 values[switch2Idx] = 0.202 + (((yMax - y) - 0.1184)/0.0148)*(0.226 - 0.202);
234
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1734 times.
 1738 else if (y > yMax - 0.1480 - eps_)
235 8 values[switch2Idx] = 0.226 + (((yMax - y) - 0.1332)/0.0148)*(0.257 - 0.226);
236
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1729 times.
 1734 else if (y > yMax - 0.1628 - eps_)
237 10 values[switch2Idx] = 0.257 + (((yMax - y) - 0.1480)/0.0148)*(0.297 - 0.257);
238
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1725 times.
 1729 else if (y > yMax - 0.1776 - eps_)
239 8 values[switch2Idx] = 0.297 + (((yMax - y) - 0.1628)/0.0148)*(0.352 - 0.297);
240
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1720 times.
 1725 else if (y > yMax - 0.1924 - eps_)
241 10 values[switch2Idx] = 0.352 + (((yMax - y) - 0.1776)/0.0148)*(0.426 - 0.352);
242
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1716 times.
 1720 else if (y > yMax - 0.2072 - eps_)
243 8 values[switch2Idx] = 0.426 + (((yMax - y) - 0.1924)/0.0148)*(0.522 - 0.426);
244
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1711 times.
 1716 else if (y > yMax - 0.2220 - eps_)
245 10 values[switch2Idx] = 0.522 + (((yMax - y) - 0.2072)/0.0148)*(0.640 - 0.522);
246
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1707 times.
 1711 else if (y > yMax - 0.2368 - eps_)
247 8 values[switch2Idx] = 0.640 + (((yMax - y) - 0.2220)/0.0148)*(0.767 - 0.640);
248
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1702 times.
 1707 else if (y > yMax - 0.2516 - eps_)
249 10 values[switch2Idx] = 0.767 + (((yMax - y) - 0.2368)/0.0148)*(0.878 - 0.767);
250
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1698 times.
 1702 else if (y > yMax - 0.2664 - eps_)
251 8 values[switch2Idx] = 0.878 + (((yMax - y) - 0.2516)/0.0148)*(0.953 - 0.878);
252
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1693 times.
 1698 else if (y > yMax - 0.2812 - eps_)
253 10 values[switch2Idx] = 0.953 + (((yMax - y) - 0.2664)/0.0148)*(0.988 - 0.953);
254
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 1687 times.
 1693 else if (y > yMax - 0.3000 - eps_)
255 12 values[switch2Idx] = 0.988;
256 else
257 3374 values[switch2Idx] = 1.e-4;
258 1779 return values;
259 }
260
261 static constexpr Scalar eps_ = 1e-6;
262 std::string name_;
263 std::vector<Scalar> Kxx_;
264 };
265
266 } // end namespace Dumux
267
268 #endif
269