GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/richards/analytical/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 60 82 73.2%
Functions: 4 10 40.0%
Branches: 64 154 41.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 RichardsTests
10 * \brief A one-dimensional infiltration problem with a smooth, given solution.
11 *
12 * The source term is calculated analytically. Thus, this example can be used
13 * to calculate the L2 error and to show convergence for grid and time-step
14 * refinement.
15 */
16
17 #ifndef DUMUX_RICHARDS_ANALYTICALPROBLEM_HH
18 #define DUMUX_RICHARDS_ANALYTICALPROBLEM_HH
19
20 #include <cmath>
21 #include <dune/common/math.hh>
22 #include <dune/geometry/quadraturerules.hh>
23
24 #include <dumux/common/properties.hh>
25 #include <dumux/common/parameters.hh>
26 #include <dumux/common/boundarytypes.hh>
27 #include <dumux/common/numeqvector.hh>
28
29 #include <dumux/porousmediumflow/problem.hh>
30
31 namespace Dumux {
32
33 /*!
34 * \ingroup RichardsTests
35 *
36 * \brief A water infiltration problem using Richards model and comparing
37 * to an analytical solution. Implemented by using the source term
38 * defined as the analytical solution.
39 *
40 * The domain is box shaped. Top and bottom boundaries are Dirichlet
41 * boundaries with fixed water pressure (fixed Saturation \f$S_w = 0\f$),
42 * left and right boundary are closed (Neumann 0 boundary).
43 * This problem uses the \ref RichardsModel
44 *
45 * The L2 error is decreasing with decreasing time and space discretization.
46 */
47 template <class TypeTag>
48 class RichardsAnalyticalProblem : public PorousMediumFlowProblem<TypeTag>
49 {
50 using ParentType = PorousMediumFlowProblem<TypeTag>;
51 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
52 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
53 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
54 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
55 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
56 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
57 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
58 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
59 // copy pressure index for convenience
60 enum { pwIdx = Indices::pressureIdx };
61 // Grid and world dimension
62 static const int dimWorld = GridView::dimensionworld;
63 static const int dim = GridView::dimension;
64 using Element = typename GridView::template Codim<0>::Entity;
65 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
66 using Geometry = typename GridView::template Codim<0>::Entity::Geometry;
67
68 public:
69 1 RichardsAnalyticalProblem(std::shared_ptr<const GridGeometry> gridGeometry)
70
7/20
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 1 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.
 3 : ParentType(gridGeometry)
71 {
72 1 pnRef_ = 1e5; // air pressure
73
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 1 name_ = getParam<std::string>("Problem.Name");
74 1 time_ = 0.0;
75 1 }
76
77 /*!
78 * \name Problem parameters
79 */
80 // \{
81
82 /*!
83 * \brief The problem name
84 *
85 * This is used as a prefix for files generated by the simulation.
86 */
87 const std::string name() const
88
2/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 1 { return name_; }
89
90 ✗ void setTime(Scalar time)
91
1/2
✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
 1000 { time_ = time; }
92
93 /*!
94 * \brief Returns the reference pressure [Pa] of the nonwetting
95 * fluid phase within a finite volume
96 *
97 * This problem assumes a constant reference pressure of 1 bar.
98 */
99 ✗ Scalar nonwettingReferencePressure() const
100 ✗ { return pnRef_; }
101
102 /*!
103 * \brief Evaluates the source values for a control volume.
104 *
105 * For this method, the \a values parameter stores primary
106 * variables. For this test case, the analytical solution is
107 * used to calculate the source term. See the Matlab script
108 * Richards.m which uses Matlab's Symbolic Toolbox to calculate
109 * the source term.
110 *
111 * \param globalPos The position for which the source term is set
112 */
113 2869248 NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
114 {
115 2869248 NumEqVector values(0.0);
116 2869248 const Scalar time = time_;
117 2869248 const Scalar pwTop = 98942.8;
118 2869248 const Scalar pwBottom = 95641.1;
119
120 // linear model with complex solution
121 // calculated with Matlab script "Richards.m"
122 using Dune::power;
123 using std::tanh;
124
125 5738496 values = (power(tanh(globalPos[1]*5.0+time*(1.0/1.0E1)-1.5E1),2)*(1.0/1.0E1)
126 5738496 -1.0/1.0E1)*(pwBottom*(1.0/2.0)-pwTop*(1.0/2.0))*4.0E-8-((power(tanh(globalPos[1]
127 2869248 *5.0+time*(1.0/1.0E1)-1.5E1),2)*5.0-5.0)*(pwBottom*(1.0/2.0)-pwTop*(1.0/2.0))-1.0E3)
128 5738496 *(power(tanh(globalPos[1]*5.0+time*(1.0/1.0E1)-1.5E1),2)*5.0-5.0)*(pwBottom
129 2869248 *(1.0/2.0)-pwTop*(1.0/2.0))*5.0E-16+tanh(globalPos[1]*5.0+time*(1.0/1.0E1)-1.5E1)
130 5738496 *(power(tanh(globalPos[1]*5.0+time*(1.0/1.0E1)-1.5E1),2)*5.0-5.0)*(pwBottom
131 11476992 *(1.0/2.0)-pwTop*(1.0/2.0))*(pwBottom*5.0E-16-(tanh(globalPos[1]*5.0+time*(1.0/1.0E1)
132 8607744 -1.5E1)+1.0)*(pwBottom*(1.0/2.0)-pwTop*(1.0/2.0))*5.0E-16+4.99995E-6)*1.0E1;
133 2869248 return values;
134 }
135
136 // \}
137
138 /*!
139 * \name Boundary conditions
140 */
141 // \{
142
143 /*!
144 * \brief Specifies which kind of boundary condition should be
145 * used for which equation on a given boundary segment.
146 *
147 * \param globalPos The position for which the boundary type is set
148 */
149 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
150 {
151
4/10
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2872050 times.
✓ Branch 5 taken 2802 times.
✓ Branch 6 taken 5744100 times.
✓ Branch 7 taken 5604 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 8624556 BoundaryTypes bcTypes;
152 34481412 if (onLowerBoundary_(globalPos) ||
153 onUpperBoundary_(globalPos))
154 {
155 bcTypes.setAllDirichlet();
156 }
157 else
158 bcTypes.setAllNeumann();
159 return bcTypes;
160 }
161
162 /*!
163 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
164 *
165 * \param globalPos The position for which the Dirichlet value is set
166 *
167 * For this method, the \a values parameter stores primary variables.
168 */
169 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
170 {
171 ✗ PrimaryVariables values(0.0);
172 ✗ const Scalar time = time_;
173 ✗ const Scalar pwTop = 98942.8;
174 ✗ const Scalar pwBottom = 95641.1;
175 using std::tanh;
176 ✗ Scalar pw = pwBottom
177 ✗ + 0.5 * (tanh( (5.0 * globalPos[1]) - 15.0 + time/10.0) + 1.0) * (pwTop - pwBottom);
178
179 ✗ values[pwIdx] = pw;
180 ✗ return values;
181 }
182
183 /*!
184 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
185 *
186 * For this method, the \a values parameter stores the mass flux
187 * in normal direction of each phase. Negative values mean influx.
188 *
189 * \param globalPos The position for which the Neumann value is set
190 */
191 ✗ NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
192 {
193
2/4
✓ Branch 0 taken 5738496 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5738496 times.
✗ Branch 3 not taken.
 11476992 NumEqVector values(0.0);
194 ✗ return values;
195 }
196
197 /*!
198 * \name Volume terms
199 */
200 // \{
201
202 /*!
203 * \brief Evaluates the initial values for a control volume.
204 *
205 * For this method, the \a values parameter stores primary
206 * variables.
207 *
208 * \param globalPos The position for which the boundary type is set
209 */
210 ✗ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
211 {
212 512 PrimaryVariables values(0.0);
213 512 analyticalSolution(values, time_, globalPos);
214 512 return values;
215 }
216
217 // \}
218
219 /*!
220 * \brief Evaluates the analytical solution.
221 *
222 * \param values The Dirichlet values for the primary variables
223 * \param time The time at which the solution should be evaluated
224 * \param globalPos The position for which the Dirichlet value is set
225 *
226 * For this method, the \a values parameter stores primary variables.
227 */
228 ✗ void analyticalSolution(PrimaryVariables &values,
229 const Scalar time,
230 const GlobalPosition &globalPos) const
231 {
232
233 ✗ const Scalar pwTop = 98942.8;
234 ✗ const Scalar pwBottom = 95641.1;
235 using std::tanh;
236 ✗ Scalar pw = pwBottom
237 ✗ + 0.5 * (tanh( (5.0 * globalPos[1]) - 15.0 + time/10.0) + 1.0) * (pwTop - pwBottom);
238
239 ✗ values[pwIdx] = pw;
240 ✗ }
241
242 /*!
243 * \brief Calculate the L2 error between the solution given by
244 * dirichletAtPos and the numerical approximation.
245 *
246 * \param curSol The current solution vector
247 * \note Works for cell-centered FV only because the numerical
248 * approximation is only evaluated in the cell center (once).
249 * To extend this function to the box method the evaluation
250 * has to be extended to box' sub-volumes.
251 */
252 1000 Scalar calculateL2Error(const SolutionVector& curSol)
253 {
254 1000 const unsigned int qOrder = 4;
255 1000 Scalar l2error = 0.0;
256 1000 Scalar l2analytic = 0.0;
257 1000 const Scalar time = time_;
258
259
1/2
✓ Branch 3 taken 513000 times.
✗ Branch 4 not taken.
 1027000 for (const auto& element :elements(this->gridGeometry().gridView()))
260 {
261 1536000 int eIdx = this->gridGeometry().elementMapper().index(element);
262 // value from numerical approximation
263 512000 Scalar numericalSolution = curSol[eIdx];
264
265 // integrate over element using a quadrature rule
266 512000 Geometry geometry = element.geometry();
267 512000 Dune::GeometryType gt = geometry.type();
268
1/2
✓ Branch 0 taken 512000 times.
✗ Branch 1 not taken.
 512000 Dune::QuadratureRule<Scalar, dim> rule =
269 1024000 Dune::QuadratureRules<Scalar, dim>::rule(gt, qOrder);
270
271
4/4
✓ Branch 0 taken 4608000 times.
✓ Branch 1 taken 512000 times.
✓ Branch 2 taken 4608000 times.
✓ Branch 3 taken 512000 times.
 10752000 for (const auto& qp : rule)
272 {
273 // evaluate analytical solution
274 9216000 Dune::FieldVector<Scalar, dim> globalPos = geometry.global(qp.position());
275 4608000 PrimaryVariables values(0.0);
276 4608000 analyticalSolution(values, time, globalPos);
277 // add contributino of current quadrature point
278 9216000 l2error += (numericalSolution - values[0]) * (numericalSolution - values[0]) *
279 9216000 qp.weight() * geometry.integrationElement(qp.position());
280 9216000 l2analytic += values[0] * values[0] *
281 9216000 qp.weight() * geometry.integrationElement(qp.position());
282 }
283 }
284 using std::sqrt;
285 1000 return sqrt(l2error/l2analytic);
286 }
287
288 /*!
289 * \brief Writes the relevant secondary variables of the current
290 * solution into an VTK output file.
291 */
292 1000 void writeOutput(const SolutionVector& curSol)
293 {
294
295 1000 Scalar l2error = calculateL2Error(curSol);
296
297 // compute L2 error if analytical solution is available
298 1000 std::cout.precision(8);
299 std::cout << "L2 error for "
300 5000 << std::setw(6) << this->gridGeometry().gridView().size(0)
301 1000 << " elements: "
302 1000 << std::scientific
303 1000 << l2error
304 1000 << std::endl;
305 1000 }
306
307 private:
308
309 // evaluates if global position is at lower boundary
310 bool onLowerBoundary_(const GlobalPosition &globalPos) const
311 {
312
20/50
✗ 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 taken 2872050 times.
✓ Branch 21 taken 2802 times.
✓ Branch 22 taken 2872050 times.
✓ Branch 23 taken 2802 times.
✓ Branch 24 taken 2872050 times.
✓ Branch 25 taken 2802 times.
✓ Branch 26 taken 2872050 times.
✓ Branch 27 taken 2802 times.
✓ Branch 28 taken 2872050 times.
✓ Branch 29 taken 2802 times.
✓ Branch 30 taken 5744100 times.
✓ Branch 31 taken 5604 times.
✓ Branch 32 taken 5744100 times.
✓ Branch 33 taken 5604 times.
✓ Branch 34 taken 5744100 times.
✓ Branch 35 taken 5604 times.
✓ Branch 36 taken 5744100 times.
✓ Branch 37 taken 5604 times.
✓ Branch 38 taken 5744100 times.
✓ Branch 39 taken 5604 times.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
 43122780 return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_;
313 }
314
315 // evaluates if global position is at upper boundary
316 bool onUpperBoundary_(const GlobalPosition &globalPos) const
317 {
318
20/50
✗ 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 taken 2802 times.
✓ Branch 21 taken 2869248 times.
✓ Branch 22 taken 2802 times.
✓ Branch 23 taken 2869248 times.
✓ Branch 24 taken 2802 times.
✓ Branch 25 taken 2869248 times.
✓ Branch 26 taken 2802 times.
✓ Branch 27 taken 2869248 times.
✓ Branch 28 taken 2802 times.
✓ Branch 29 taken 2869248 times.
✓ Branch 30 taken 5604 times.
✓ Branch 31 taken 5738496 times.
✓ Branch 32 taken 5604 times.
✓ Branch 33 taken 5738496 times.
✓ Branch 34 taken 5604 times.
✓ Branch 35 taken 5738496 times.
✓ Branch 36 taken 5604 times.
✓ Branch 37 taken 5738496 times.
✓ Branch 38 taken 5604 times.
✓ Branch 39 taken 5738496 times.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
 43080750 return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_;
319 }
320
321 static constexpr Scalar eps_ = 3e-6;
322 Scalar pnRef_;
323 std::string name_;
324 Scalar time_;
325 };
326 } // end namespace Dumux
327
328 #endif
329