GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/3p3c/infiltration/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 50 52 96.2%
Functions: 8 10 80.0%
Branches: 75 118 63.6%
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 Isothermal NAPL infiltration problem: LNAPL contaminates
11 * the unsaturated and the saturated groundwater zone.
12 */
13
14 #ifndef DUMUX_INFILTRATION_THREEPTHREEC_PROBLEM_HH
15 #define DUMUX_INFILTRATION_THREEPTHREEC_PROBLEM_HH
16
17 #include <dumux/common/properties.hh>
18 #include <dumux/common/parameters.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 Isothermal NAPL infiltration problem: LNAPL contaminates
29 * the unsaturated and the saturated groundwater zone.
30 *
31 * The 2D domain of this test problem is 500 m long and 10 m deep, where
32 * the lower part represents a slightly inclined groundwater table, and the
33 * upper part is the vadose zone.
34 * A LNAPL (Non-Aqueous Phase Liquid which is lighter than water) infiltrates
35 * (modelled with a Neumann boundary condition) into the vadose zone. Upon
36 * reaching the water table, it spreads (since lighter than water) and migrates
37 * on top of the water table in the direction of the slope.
38 * On its way through the vadose zone, it leaves a trace of residually trapped
39 * immobile NAPL, which can in the following dissolve and evaporate slowly,
40 * and eventually be transported by advection and diffusion.
41 *
42 * Left and right boundaries are constant head boundaries (Dirichlet),
43 * Top and bottom are Neumann boundaries, all no-flow except for the small
44 * infiltration zone in the upper left part.
45 *
46 * This problem uses the \ref ThreePThreeCModel.
47 * */
48 template <class TypeTag >
49 class InfiltrationThreePThreeCProblem : public PorousMediumFlowProblem<TypeTag>
50 {
51 using ParentType = PorousMediumFlowProblem<TypeTag>;
52
53 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
54 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
55 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
56 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
57
58 // copy some indices for convenience
59 enum {
60 pressureIdx = Indices::pressureIdx,
61 switch1Idx = Indices::switch1Idx,
62 switch2Idx = Indices::switch2Idx,
63
64 // phase state
65 wgPhaseOnly = Indices::wgPhaseOnly,
66
67 contiWEqIdx = Indices::conti0EqIdx + FluidSystem::wCompIdx, //!< Index of the mass conservation equation for the water component
68 contiNEqIdx = Indices::conti0EqIdx + FluidSystem::nCompIdx,//!< Index of the mass conservation equation for the contaminant component
69 contiAEqIdx = Indices::conti0EqIdx + FluidSystem::gCompIdx,//!< Index of the mass conservation equation for the gas component
70
71 // world dimension
72 dimWorld = GridView::dimensionworld
73 };
74
75 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
76 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
77 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
78 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
79
80 using Element = typename GridView::template Codim<0>::Entity;
81 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
82
83 public:
84 2 InfiltrationThreePThreeCProblem(std::shared_ptr<const GridGeometry> gridGeometry)
85
7/20
✓ 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 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
 6 : ParentType(gridGeometry)
86 {
87
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 FluidSystem::init(/*tempMin=*/282.15,
88 /*tempMax=*/284.15,
89 /*nTemp=*/3,
90 /*pressMin=*/0.8*1e5,
91 /*pressMax=*/3*1e5,
92 /*nPress=*/200);
93
94
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");
95 2 }
96
97 /*!
98 * \brief The problem name.
99 *
100 * This is used as a prefix for files generated by the simulation.
101 */
102 const std::string& name() const
103
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 { return name_; }
104
105 // \}
106
107 /*!
108 * \name Boundary conditions
109 */
110 // \{
111
112 /*!
113 * \brief Specifies which kind of boundary condition should be
114 * used for which equation on a given boundary segment.
115 *
116 * \param globalPos The position for which the bc type should be evaluated
117 */
118 22998 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
119 {
120 22998 BoundaryTypes values;
121
10/10
✓ Branch 0 taken 1144 times.
✓ Branch 1 taken 21854 times.
✓ Branch 2 taken 1144 times.
✓ Branch 3 taken 21854 times.
✓ Branch 4 taken 1144 times.
✓ Branch 5 taken 21854 times.
✓ Branch 6 taken 1144 times.
✓ Branch 7 taken 21854 times.
✓ Branch 8 taken 1144 times.
✓ Branch 9 taken 21854 times.
 114990 if(globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
122 values.setAllDirichlet();
123
10/10
✓ Branch 0 taken 20710 times.
✓ Branch 1 taken 1144 times.
✓ Branch 2 taken 20710 times.
✓ Branch 3 taken 1144 times.
✓ Branch 4 taken 20710 times.
✓ Branch 5 taken 1144 times.
✓ Branch 6 taken 20710 times.
✓ Branch 7 taken 1144 times.
✓ Branch 8 taken 20710 times.
✓ Branch 9 taken 1144 times.
 109270 else if(globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_)
124 values.setAllDirichlet();
125 else
126 values.setAllNeumann();
127 22998 return values;
128 }
129
130 /*!
131 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
132 *
133 * \param globalPos The position for which the bc type should be evaluated
134 *
135 * For this method, the \a values parameter stores primary variables.
136 */
137 PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
138 656 { return initial_(globalPos); }
139
140 /*!
141 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
142 *
143 * \param globalPos The position for which the bc type should be evaluated
144 *
145 * For this method, the \a values parameter stores the mass flux
146 * in normal direction of each phase. Negative values mean influx.
147 */
148 ✗ NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
149 {
150 109820 NumEqVector values(0.0);
151
152 // negative values for injection
153
24/48
✗ 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 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 15000 times.
✓ Branch 25 taken 82500 times.
✓ Branch 26 taken 15000 times.
✓ Branch 27 taken 82500 times.
✓ Branch 28 taken 5000 times.
✓ Branch 29 taken 10000 times.
✓ Branch 30 taken 5000 times.
✓ Branch 31 taken 10000 times.
✓ Branch 32 taken 2500 times.
✓ Branch 33 taken 2500 times.
✓ Branch 34 taken 2500 times.
✓ Branch 35 taken 2500 times.
✓ Branch 36 taken 1848 times.
✓ Branch 37 taken 10472 times.
✓ Branch 38 taken 1848 times.
✓ Branch 39 taken 10472 times.
✓ Branch 40 taken 616 times.
✓ Branch 41 taken 1232 times.
✓ Branch 42 taken 616 times.
✓ Branch 43 taken 1232 times.
✓ Branch 44 taken 308 times.
✓ Branch 45 taken 308 times.
✓ Branch 46 taken 308 times.
✓ Branch 47 taken 308 times.
 219640 if ((globalPos[0] < 80.0 + eps_) && (globalPos[0] > 55.0 - eps_) && (globalPos[1] > 10.0 - eps_))
154 {
155 2808 values[contiWEqIdx] = -0.0;
156 //mole flow conversion to mass flow with molar mass M(Mesit.)=0,120 kg/mol --> 1.2e-4 kg/(sm)
157 //the 3p3c model uses mole fractions
158 2808 values[contiNEqIdx] = -0.001;
159 5616 values[contiAEqIdx] = -0.0;
160 }
161
162 ✗ return values;
163 }
164
165 // \}
166
167 /*!
168 * \name Volume terms
169 */
170 // \{
171
172 /*!
173 * \brief Evaluates the initial value for a control volume.
174 *
175 * \param globalPos The position for which the initial condition should be evaluated
176 *
177 * For this method, the \a values parameter stores primary
178 * variables.
179 */
180 PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
181 365 { return initial_(globalPos); }
182
183 private:
184 // internal method for the initial condition
185 1021 PrimaryVariables initial_(const GlobalPosition &globalPos) const
186 {
187 1021 PrimaryVariables values(0.0);
188
2/2
✓ Branch 0 taken 591 times.
✓ Branch 1 taken 430 times.
 1021 values.setState(wgPhaseOnly);
189
190
2/2
✓ Branch 0 taken 591 times.
✓ Branch 1 taken 430 times.
 1021 Scalar y = globalPos[1];
191
2/2
✓ Branch 0 taken 591 times.
✓ Branch 1 taken 430 times.
 1021 Scalar x = globalPos[0];
192 1021 Scalar sw, swr=0.12, sgr=0.03;
193
194
2/2
✓ Branch 0 taken 591 times.
✓ Branch 1 taken 430 times.
 1021 if(y >(-1.E-3*x+5) - eps_)
195 {
196 591 Scalar pc = 9.81 * 1000.0 * (y - (-5E-4*x+5));
197
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 591 times.
 591 if (pc < 0.0) pc = 0.0;
198
199 1773 sw = invertPcgw_(pc, this->spatialParams().fluidMatrixInteractionAtPos(globalPos));
200
2/2
✓ Branch 0 taken 406 times.
✓ Branch 1 taken 185 times.
 591 if (sw < swr) sw = swr;
201
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 515 times.
 591 if (sw > 1.-sgr) sw = 1.-sgr;
202
203 591 values[pressureIdx] = 1e5 ;
204 591 values[switch1Idx] = sw;
205 1182 values[switch2Idx] = 1.e-6;
206 }else {
207 430 values[pressureIdx] = 1e5 + 9.81 * 1000.0 * ((-5E-4*x+5) - y);
208 430 values[switch1Idx] = 1.-sgr;
209 860 values[switch2Idx] = 1.e-6;
210 }
211 1021 return values;
212 }
213
214 template<class FluidMatrixInteraction>
215 591 static Scalar invertPcgw_(Scalar pcIn, const FluidMatrixInteraction& fluidmatrixinteraction)
216 {
217 Scalar lower,upper;
218 int k;
219 591 int maxIt = 50;
220 591 Scalar bisLimit = 1.;
221 Scalar sw, pcgw;
222 591 lower=0.0; upper=1.0;
223
1/2
✓ Branch 0 taken 9175 times.
✗ Branch 1 not taken.
 9175 for (k=1; k<=25; k++)
224 {
225 9175 sw = 0.5*(upper+lower);
226 9175 pcgw = fluidmatrixinteraction.pcgw(sw, 0.0/*sn*/);
227 9175 Scalar delta = pcgw-pcIn;
228
2/2
✓ Branch 0 taken 4574 times.
✓ Branch 1 taken 4601 times.
 9175 if (delta<0.) delta*=-1.;
229
2/2
✓ Branch 0 taken 591 times.
✓ Branch 1 taken 8584 times.
 9175 if (delta<bisLimit)
230 {
231 591 return(sw);
232 }
233 if (k==maxIt) {
234 return(sw);
235 }
236
2/2
✓ Branch 0 taken 4378 times.
✓ Branch 1 taken 4206 times.
 8584 if (pcgw>pcIn) lower=sw;
237 4378 else upper=sw;
238 }
239 return(sw);
240 }
241
242 static constexpr Scalar eps_ = 1e-6;
243 std::string name_;
244 };
245
246 } //end namespace Dumux
247
248 #endif
249