GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/3p3c/kuevette/problem.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 36 38 94.7%
Functions: 8 10 80.0%
Branches: 65 94 69.1%
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 gas injection problem where a gas (e.g. steam/air)
11 * is injected into a unsaturated porous medium with a residually
12 * trapped NAPL contamination.
13 */
14
15 #ifndef DUMUX_KUEVETTE3P3CNIPROBLEM_HH
16 #define DUMUX_KUEVETTE3P3CNIPROBLEM_HH
17
18 #include <dune/common/float_cmp.hh>
19
20 #include <dumux/common/parameters.hh>
21 #include <dumux/common/properties.hh>
22 #include <dumux/common/boundarytypes.hh>
23 #include <dumux/common/numeqvector.hh>
24
25 #include <dumux/porousmediumflow/problem.hh>
26
27 namespace Dumux {
28
29 /*!
30 * \ingroup ThreePThreeCTests
31 * \brief Non-isothermal gas injection problem where a gas (e.g. steam/air)
32 * is injected into a unsaturated porous medium with a residually
33 * trapped NAPL contamination.
34 *
35 * The domain is a quasi-two-dimensional container (kuevette). Its dimensions
36 * are 1.5m x 0.74m. The top and bottom boundaries are closed, the right
37 * boundary is a Dirichlet boundary allowing fluids to escape. From the left,
38 * an injection of a hot water-air mixture is applied (Neumann boundary condition
39 * for the mass components and the enthalpy), aimed at remediating an initial
40 * NAPL (Non-Aquoeus Phase Liquid) contamination in the heterogeneous domain.
41 * The contamination is initially placed partly into the coarse sand
42 * and partly into a fine sand lens.
43 *
44 * This simulation can be varied through assigning different boundary conditions
45 * at the left boundary as described in \cite Class2001.
46 *
47 * This problem uses the \ref ThreePThreeCModel and \ref NIModel model.
48 *
49 * To see the basic effect and the differences to scenarios with pure steam or
50 * pure air injection, it is sufficient to simulate for about 2-3 hours (10000 s).
51 * Complete remediation of the domain requires much longer (about 10 days simulated time).
52 */
53 template <class TypeTag >
54 class KuevetteProblem : public PorousMediumFlowProblem<TypeTag>
55 {
56 using ParentType = PorousMediumFlowProblem<TypeTag>;
57 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
58 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
59
60 // copy some indices for convenience
61 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
62 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
63 enum {
64
65 pressureIdx = Indices::pressureIdx,
66 switch1Idx = Indices::switch1Idx,
67 switch2Idx = Indices::switch2Idx,
68 temperatureIdx = Indices::temperatureIdx,
69 contiWEqIdx = Indices::conti0EqIdx + FluidSystem::wCompIdx,
70 contiGEqIdx = Indices::conti0EqIdx + FluidSystem::gCompIdx,
71 contiNEqIdx = Indices::conti0EqIdx + FluidSystem::nCompIdx,
72 energyEqIdx = Indices::energyEqIdx,
73
74 // phase states
75 threePhases = Indices::threePhases,
76 wgPhaseOnly = Indices::wgPhaseOnly,
77
78 // world dimension
79 dimWorld = GridView::dimensionworld
80 };
81
82 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
83 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
84 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
85 using Element = typename GridView::template Codim<0>::Entity;
86 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
87 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
88 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
89 using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
90 using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
91 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
92
93 using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
94
95 public:
96 2 KuevetteProblem(std::shared_ptr<const GridGeometry> gridGeometry)
97
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)
98 {
99
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 FluidSystem::init();
100
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");
101 2 }
102
103 /*!
104 * \name Problem parameters
105 */
106 // \{
107
108 /*!
109 * \brief The problem name.
110 *
111 * This is used as a prefix for files generated by the simulation.
112 */
113 const std::string name() const
114
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_; }
115
116 // \}
117
118 /*!
119 * \name Boundary conditions
120 */
121 // \{
122
123 /*!
124 * \brief Specifies which kind of boundary condition should be
125 * used for which equation on a given boundary segment.
126 *
127 * \param globalPos The position for which the bc type should be evaluated
128 */
129 99585 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
130 {
131 99585 BoundaryTypes bcTypes;
132
10/10
✓ Branch 0 taken 81025 times.
✓ Branch 1 taken 18560 times.
✓ Branch 2 taken 81025 times.
✓ Branch 3 taken 18560 times.
✓ Branch 4 taken 81025 times.
✓ Branch 5 taken 18560 times.
✓ Branch 6 taken 81025 times.
✓ Branch 7 taken 18560 times.
✓ Branch 8 taken 81025 times.
✓ Branch 9 taken 18560 times.
 497925 if(globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
133 bcTypes.setAllDirichlet();
134 else
135 bcTypes.setAllNeumann();
136 99585 return bcTypes;
137 }
138
139 /*!
140 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
141 *
142 * \param globalPos The position for which the bc type should be evaluated
143 *
144 * For this method, the \a values parameter stores primary variables.
145 */
146 PrimaryVariables dirichletAtPos( const GlobalPosition &globalPos) const
147 {
148 4680 return initial_(globalPos);
149 }
150
151 /*!
152 * \brief Evaluates the boundary conditions for a N eumann boundary segment.
153 *
154 * \param element The finite element
155 * \param fvGeometry The finite-volume geometry in the box scheme
156 * \param elemVolVars The element volume variables
157 * \param elemFluxVarsCache Flux variables caches for all faces in stencil
158 * \param scvf The sub-control volume face
159 *
160 * For this method, the \a values parameter stores the mass flux
161 * in normal direction of each phase. Negative values mean influx.
162 */
163 ✗ NumEqVector neumann(const Element& element,
164 const FVElementGeometry& fvGeometry,
165 const ElementVolumeVariables& elemVolVars,
166 const ElementFluxVariablesCache& elemFluxVarsCache,
167 const SubControlVolumeFace& scvf) const
168 {
169 446856 NumEqVector values(0.0);
170
4/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 85008 times.
✓ Branch 3 taken 308154 times.
✓ Branch 4 taken 11304 times.
✓ Branch 5 taken 42390 times.
 446856 const auto& globalPos = scvf.ipGlobal();
171
172 // negative values for injection
173
8/12
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 85008 times.
✓ Branch 5 taken 308154 times.
✓ Branch 6 taken 85008 times.
✓ Branch 7 taken 308154 times.
✓ Branch 8 taken 11304 times.
✓ Branch 9 taken 42390 times.
✓ Branch 10 taken 11304 times.
✓ Branch 11 taken 42390 times.
 893712 if (globalPos[0] < eps_)
174 {
175 96312 values[contiWEqIdx] = -0.1435; // 0.3435 [mol/(s m)] in total
176 96312 values[contiGEqIdx] = -0.2;
177 96312 values[contiNEqIdx] = 0.0;
178 192624 values[energyEqIdx] = -6929.;
179 }
180 ✗ return values;
181 }
182
183 // \}
184
185 /*!
186 * \name Volume terms
187 */
188 // \{
189
190 /*!
191 * \brief Evaluates the initial value for a control volume.
192 *
193 * \param globalPos The position for which the initial condition should be evaluated
194 *
195 * For this method, the \a values parameter stores primary
196 * variables.
197 */
198 PrimaryVariables initialAtPos( const GlobalPosition &globalPos) const
199 {
200 264 return initial_(globalPos);
201 }
202
203 /*!
204 * \brief Appends all quantities of interest which can be derived
205 * from the solution of the current time step to the VTK
206 * writer.
207 *
208 * Adjust this in case of anisotropic permeabilities.
209 */
210 template<class VTKWriter>
211 void addVtkFields(VTKWriter& vtk)
212 {
213 const auto& gg = this->gridGeometry();
214 Kxx_.resize(gg.numDofs());
215 vtk.addField(Kxx_, "permeability");
216 auto fvGeometry = localView(gg);
217 for (const auto& element : elements(this->gridView()))
218 {
219 fvGeometry.bindElement(element);
220 for (const auto& scv : scvs(fvGeometry))
221 Kxx_[scv.dofIndex()] = this->spatialParams().intrinsicPermeabilityAtPos(scv.dofPosition());
222 }
223 }
224
225 private:
226 // checks, whether a point is located inside the contamination zone
227 9888 bool isInContaminationZone(const GlobalPosition &globalPos) const
228 {
229
4/4
✓ Branch 0 taken 86 times.
✓ Branch 1 taken 9802 times.
✓ Branch 2 taken 86 times.
✓ Branch 3 taken 9802 times.
 19776 return (Dune::FloatCmp::ge<Scalar>(globalPos[0], 0.2)
230
4/4
✓ Branch 0 taken 9616 times.
✓ Branch 1 taken 204 times.
✓ Branch 2 taken 9616 times.
✓ Branch 3 taken 204 times.
 19640 && Dune::FloatCmp::le<Scalar>(globalPos[0], 0.8)
231
4/4
✓ Branch 0 taken 118 times.
✓ Branch 1 taken 104 times.
✓ Branch 2 taken 118 times.
✓ Branch 3 taken 104 times.
 444 && Dune::FloatCmp::ge<Scalar>(globalPos[1], 0.4)
232
4/4
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 78 times.
✓ Branch 2 taken 26 times.
✓ Branch 3 taken 78 times.
 208 && Dune::FloatCmp::le<Scalar>(globalPos[1], 0.65));
233 }
234
235 // internal method for the initial condition (reused for the
236 // dirichlet conditions!)
237 4944 PrimaryVariables initial_(const GlobalPosition &globalPos) const
238 {
239
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 4905 times.
 4944 PrimaryVariables values;
240
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 4905 times.
 4944 if (isInContaminationZone(globalPos))
241 39 values.setState(threePhases);
242 else
243 4905 values.setState(wgPhaseOnly);
244
245
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 4905 times.
 4944 values[pressureIdx] = 1e5 ;
246
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 4905 times.
 4944 values[switch1Idx] = 0.12;
247
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 4905 times.
 4944 values[switch2Idx] = 1.e-6;
248
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 4905 times.
 4944 values[temperatureIdx] = 293.0;
249
250
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 4905 times.
 4944 if (isInContaminationZone(globalPos))
251 {
252 78 values[switch2Idx] = 0.07;
253 }
254 4944 return values;
255 }
256
257 static constexpr Scalar eps_ = 1e-6;
258 std::string name_;
259 std::vector<Scalar> Kxx_;
260 };
261
262 } // end namespace Dumux
263
264 #endif
265