GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/richardsnc/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 59 62 95.2%
Functions: 13 17 76.5%
Branches: 64 110 58.2%
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 RichardsNCTests
10 * \brief A water infiltration problem with a low-permeability lens
11 * embedded into a high-permeability domain which uses the
12 * Richards box model.
13 */
14
15 #ifndef DUMUX_RICHARDS_NC_WELL_TRACER_PROBLEM_HH
16 #define DUMUX_RICHARDS_NC_WELL_TRACER_PROBLEM_HH
17
18 #include <dumux/common/properties.hh>
19 #include <dumux/common/parameters.hh>
20 #include <dumux/common/boundarytypes.hh>
21 #include <dumux/common/numeqvector.hh>
22
23 #include <dumux/porousmediumflow/problem.hh>
24
25 namespace Dumux {
26
27 /*!
28 * \ingroup RichardsNCTests
29 *
30 * \brief A water infiltration problem with a low-permeability lens
31 * embedded into a high-permeability domain which uses the
32 * Richards model.
33 *
34 * The domain is box shaped. Left and right boundaries are Dirichlet
35 * boundaries with fixed water pressure (hydrostatic, gradient from right to left),
36 * bottom boundary is closed (Neumann 0 boundary), the top boundary
37 * (Neumann 0 boundary) is also closed. Water is extracted at a point in
38 * the middle of the domain.
39 * This problem is very similar to the LensProblem
40 * which uses the TwoPBoxModel, with the main difference being that
41 * the domain is initially fully saturated by gas instead of water and
42 * water instead of a %DNAPL infiltrates from the top.
43 *
44 * This problem uses the \ref RichardsNCModel
45 */
46 template <class TypeTag>
47 class RichardsWellTracerProblem : public PorousMediumFlowProblem<TypeTag>
48 {
49 using ParentType = PorousMediumFlowProblem<TypeTag>;
50 using Problem = GetPropType<TypeTag, Properties::Problem>;
51 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
52 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
53 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
54 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
55 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
56 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
57 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
58 using PointSource = GetPropType<TypeTag, Properties::PointSource>;
59 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
60 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
61 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
62 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
63 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
64 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
65 enum {
66 pressureIdx = Indices::pressureIdx,
67 compIdx = Indices::compMainIdx + 1,
68 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
69 dimWorld = GridView::dimensionworld
70 };
71 using Element = typename GridView::template Codim<0>::Entity;
72 using GlobalPosition = typename SubControlVolume::GlobalPosition;
73
74 public:
75 2 RichardsWellTracerProblem(std::shared_ptr<const GridGeometry> gridGeometry)
76
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)
77 {
78
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");
79
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 contaminantMoleFraction_ = getParam<Scalar>("Problem.ContaminantMoleFraction");
80
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 pumpRate_ = getParam<Scalar>("Problem.PumpRate"); // in kg/s
81
82 // for initial conditions
83 2 const Scalar sw = 0.4; // start with 40% saturation on top
84 10 pcTop_ = this->spatialParams().fluidMatrixInteractionAtPos(this->gridGeometry().bBoxMax()).pc(sw);
85
86 // for post time step mass balance
87 2 accumulatedSource_ = 0.0;
88 2 }
89
90 68 void printTracerMass(const SolutionVector& curSol,
91 const GridVariables& gridVariables,
92 const Scalar timeStepSize)
93
94 {
95 // compute the mass in the entire domain to make sure the tracer is conserved
96 68 Scalar tracerMass = 0.0;
97
98
2/4
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 68 times.
✗ Branch 5 not taken.
 169 auto fvGeometry = localView(this->gridGeometry());
99
2/4
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 68 times.
✗ Branch 5 not taken.
 171 auto elemVolVars = localView(gridVariables.curGridVolVars());
100
101 // bulk elements
102
5/10
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 68 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 68 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 68068 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 68000 times.
✗ Branch 13 not taken.
 136204 for (const auto& element : elements(this->gridGeometry().gridView()))
103 {
104
1/2
✓ Branch 1 taken 68000 times.
✗ Branch 2 not taken.
 68000 fvGeometry.bindElement(element);
105
1/2
✓ Branch 1 taken 68000 times.
✗ Branch 2 not taken.
 68000 elemVolVars.bindElement(element, fvGeometry, curSol);
106
107
4/4
✓ Branch 0 taken 173000 times.
✓ Branch 1 taken 68000 times.
✓ Branch 2 taken 140000 times.
✓ Branch 3 taken 35000 times.
 416000 for (auto&& scv : scvs(fvGeometry))
108 {
109
1/2
✓ Branch 1 taken 140000 times.
✗ Branch 2 not taken.
 173000 const auto& volVars = elemVolVars[scv];
110
2/4
✓ Branch 1 taken 173000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 173000 times.
✗ Branch 5 not taken.
 173000 tracerMass += volVars.massFraction(liquidPhaseIdx, compIdx)*volVars.density(liquidPhaseIdx)
111
2/4
✓ Branch 1 taken 173000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 173000 times.
✗ Branch 5 not taken.
 346000 * scv.volume() * volVars.saturation(liquidPhaseIdx) * volVars.porosity() * volVars.extrusionFactor();
112
113
2/6
✓ Branch 1 taken 173000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 173000 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 173000 accumulatedSource_ += this->scvPointSources(element, fvGeometry, elemVolVars, scv)[compIdx]
114 173000 * scv.volume() * volVars.extrusionFactor()
115
1/2
✓ Branch 1 taken 173000 times.
✗ Branch 2 not taken.
 173000 * FluidSystem::molarMass(compIdx)
116 173000 * timeStepSize;
117 }
118 }
119
120
1/2
✓ Branch 3 taken 68 times.
✗ Branch 4 not taken.
 204 std::cout << "\033[1;33m" << "The domain contains " << tracerMass*1e9 << " µg tracer, "
121
1/2
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
 68 << accumulatedSource_*1e9 << " µg ("<< int(std::round(-accumulatedSource_/(tracerMass - accumulatedSource_)*100))
122
1/2
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
 68 <<"%) was already extracted (balanced: "
123
3/6
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 68 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 35 times.
✗ Branch 8 not taken.
 136 << (tracerMass - accumulatedSource_)*1e9 << " µg)\033[0m" << '\n';
124
125 68 }
126
127 /*!
128 * \name Problem parameters
129 */
130 // \{
131
132 /*!
133 * \brief The problem name
134 *
135 * This is used as a prefix for files generated by the simulation.
136 */
137 const std::string& name() const
138
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 { return name_; }
139
140 /*!
141 * \brief Returns the reference pressure [Pa] of the nonwetting
142 * fluid phase within a finite volume
143 *
144 * This problem assumes a constant reference pressure of 1 bar.
145 */
146 ✗ Scalar nonwettingReferencePressure() const
147 ✗ { return 1.0e5; };
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 boundary type is set
161 */
162 163576 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
163 {
164 163576 BoundaryTypes bcTypes;
165 607344 if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
166 bcTypes.setAllDirichlet();
167 else
168 bcTypes.setAllNeumann();
169 163576 return bcTypes;
170 }
171
172 /*!
173 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
174 *
175 * \param globalPos The position for which the Dirichlet value is set
176 *
177 * For this method, the \a values parameter stores primary variables.
178 */
179 PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
180 23840 { return initial_(globalPos); }
181
182 /*!
183 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
184 *
185 * For this method, the \a values parameter stores the mass flux
186 * in normal direction of each phase. Negative values mean influx.
187 *
188 * \param globalPos The position for which the Neumann value is set
189 */
190 ✗ NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
191
1/2
✓ Branch 0 taken 51000 times.
✗ Branch 1 not taken.
 800288 { return NumEqVector(0.0); }
192
193 /*!
194 * \name Volume terms
195 */
196 // \{
197
198 /*!
199 * \brief Applies a vector of point sources which are possibly solution dependent.
200 *
201 * \param pointSources A vector of PointSource s that contain
202 source values for all phases and space positions.
203 *
204 * For this method, the \a values method of the point source
205 * has to return the absolute rate values in units
206 * \f$ [ \textnormal{unit of conserved quantity} / s ] \f$.
207 * Positive values mean that the conserved quantity is created, negative ones mean that it vanishes.
208 * E.g. for the mass balance that would be a mass rate in \f$ [ kg / s ] \f$.
209 */
210 2 void addPointSources(std::vector<PointSource>& pointSources) const
211 {
212 10 auto globalPos = this->gridGeometry().bBoxMax()-this->gridGeometry().bBoxMin();
213 2 globalPos *= 0.5;
214 //! Add point source in middle of domain
215 2 pointSources.emplace_back(globalPos,
216 [this](const Problem &problem,
217 const Element &element,
218 const FVElementGeometry &fvGeometry,
219 const ElementVolumeVariables &elemVolVars,
220 2172 const SubControlVolume &scv)
221 {
222 16564 const auto& volVars = elemVolVars[scv];
223 //! convert pump rate from kg/s to mol/s
224 //! We assume we can't keep up the pump rate if the saturation sinks
225 37472 const Scalar value = pumpRate_*volVars.molarDensity(liquidPhaseIdx)/volVars.density(liquidPhaseIdx)*volVars.saturation(liquidPhaseIdx);
226 28104 return PrimaryVariables({-value, -value*volVars.moleFraction(liquidPhaseIdx, compIdx)});
227 });
228 2 }
229
230 /*!
231 * \brief Evaluates the initial values for a control volume.
232 *
233 * For this method, the \a values parameter stores primary
234 * variables.
235 *
236 * \param globalPos The position for which the boundary type is set
237 */
238 PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
239 2071 { return initial_(globalPos); };
240
241 // \}
242
243 private:
244 25911 PrimaryVariables initial_(const GlobalPosition &globalPos) const
245 {
246 25911 const auto xTracer = [&,this]()
247 {
248 155466 const GlobalPosition contaminationPos({0.2*this->gridGeometry().bBoxMax()[0], 0.5*this->gridGeometry().bBoxMax()[1]});
249
2/2
✓ Branch 0 taken 215 times.
✓ Branch 1 taken 25696 times.
 207288 if ((globalPos - contaminationPos).two_norm() < 0.1*(this->gridGeometry().bBoxMax()-this->gridGeometry().bBoxMin()).two_norm() + eps_)
250 215 return contaminantMoleFraction_;
251 else
252 return 0.0;
253 77948 }();
254
255 25911 PrimaryVariables values(0.0);
256 //! Hydrostatic pressure profile
257 77733 values[pressureIdx] = (nonwettingReferencePressure() - pcTop_)
258 129555 - 9.81*1000*(globalPos[dimWorld-1] - this->gridGeometry().bBoxMax()[dimWorld-1]);
259 51822 values[compIdx] = xTracer;
260 25911 return values;
261 }
262
263 bool onLeftBoundary_(const GlobalPosition &globalPos) const
264 {
265
10/10
✓ Branch 0 taken 140096 times.
✓ Branch 1 taken 23480 times.
✓ Branch 2 taken 140096 times.
✓ Branch 3 taken 23480 times.
✓ Branch 4 taken 140096 times.
✓ Branch 5 taken 23480 times.
✓ Branch 6 taken 140096 times.
✓ Branch 7 taken 23480 times.
✓ Branch 8 taken 140096 times.
✓ Branch 9 taken 23480 times.
 817880 return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_;
266 }
267
268 bool onRightBoundary_(const GlobalPosition &globalPos) const
269 {
270
10/10
✓ Branch 0 taken 116616 times.
✓ Branch 1 taken 23480 times.
✓ Branch 2 taken 116616 times.
✓ Branch 3 taken 23480 times.
✓ Branch 4 taken 116616 times.
✓ Branch 5 taken 23480 times.
✓ Branch 6 taken 116616 times.
✓ Branch 7 taken 23480 times.
✓ Branch 8 taken 116616 times.
✓ Branch 9 taken 23480 times.
 700480 return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_;
271 }
272
273 bool onLowerBoundary_(const GlobalPosition &globalPos) const
274 {
275 return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_;
276 }
277
278 bool onUpperBoundary_(const GlobalPosition &globalPos) const
279 {
280 return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_;
281 }
282
283 static constexpr Scalar eps_ = 1.5e-7;
284 std::string name_;
285 Scalar contaminantMoleFraction_;
286 Scalar pumpRate_;
287 Scalar pcTop_;
288 Scalar accumulatedSource_;
289 };
290
291 } // end namespace Dumux
292
293 #endif
294