GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/1p/rootbenchmark/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 67 68 98.5%
Functions: 4 5 80.0%
Branches: 88 176 50.0%
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 OnePTests
10 * \brief Root benchmark case Schnepf et al 2020 M3.1 doi: 10.3389/fpls.2020.00316
11 */
12
13 #ifndef DUMUX_ROOTBENCHMARK_TEST_PROBLEM_HH
14 #define DUMUX_ROOTBENCHMARK_TEST_PROBLEM_HH
15
16 #include <dune/common/fvector.hh>
17 #include <dune/common/fmatrix.hh>
18
19 #include <dumux/common/properties.hh>
20 #include <dumux/common/parameters.hh>
21 #include <dumux/common/boundarytypes.hh>
22 #include <dumux/common/numeqvector.hh>
23
24 #include <dumux/io/format.hh>
25
26 #include <dumux/porousmediumflow/problem.hh>
27
28 #include <dumux/material/components/constant.hh>
29
30 namespace Dumux {
31
32 /*!
33 * \ingroup OnePTests
34 * \brief Root benchmark case Schnepf et al 2020 M3.1 doi: 10.3389/fpls.2020.00316
35 */
36 template <class TypeTag>
37 class RootBenchmarkProblem : public PorousMediumFlowProblem<TypeTag>
38 {
39 using ParentType = PorousMediumFlowProblem<TypeTag>;
40 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
41 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
42 using GridView = typename GridGeometry::GridView;
43 using FVElementGeometry = typename GridGeometry::LocalView;
44 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
45 using Element = typename GridView::template Codim<0>::Entity;
46 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
47
48 static constexpr int dimWorld = GridView::dimensionworld;
49
50 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
51 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
52 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
53
54 public:
55 12 RootBenchmarkProblem(std::shared_ptr<const GridGeometry> gridGeometry)
56
10/28
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 12 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 12 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 12 times.
✓ Branch 15 taken 12 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 12 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 12 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 12 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 12 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
 36 : ParentType(gridGeometry)
57 {
58
2/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 12 times.
 12 name_ = getParam<std::string>("Problem.Name");
59
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 rootCollarPressure_ = getParam<double>("Problem.RootCollarPressure");
60
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 soilPressure_ = getParam<double>("Problem.SoilPressure");
61
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 radius_ = getParam<double>("Problem.Radius");
62
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 density_ = getParam<double>("Component.LiquidDensity");
63
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 Kr_ = this->spatialParams().radialConductivity();
64
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 Kx_ = this->spatialParams().axialConductivity();
65
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 const auto enableGravity = getParam<bool>("Problem.EnableGravity");
66
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 6 times.
 12 const Scalar gravity = enableGravity ? 9.81 : 0.0;
67
68 // cf. Schnepf et al 2020 doi: 10.3389/fpls.2020.00316
69 // Equation is 2*pi*R*Kr(ps - p0) = d/ds(-Kx (dp/ds + rho*g))
70 // Analytical solution is p = ps + A*exp(Cs) + B*exp(-Cs)
71 // where C = sqrt(2*pi*R*Kr/Kx)
72 // Two boundary conditions are needed to determine A and B
73 // Dirichlet at root collar (s=0, p=p0)
74 // (1) A + B = p0 - ps
75 // Neumann no-flow at tip (s=-L, dp/ds + rho*g = 0)
76 // (2) A*C*exp(-CL) - B*C*exp(CL) = -rho*g
77 // Assemble (1) and (2) and solve for A and B
78 Dune::FieldMatrix<double, 2, 2> M;
79 12 Dune::FieldVector<double, 2> b;
80 using std::exp; using std::sqrt;
81 12 const auto C = sqrt(2.0*M_PI*radius_*Kr_/Kx_); C_ = C;
82
4/8
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12 times.
✗ Branch 11 not taken.
 48 const auto L = this->gridGeometry().bBoxMax()[dimWorld-1] - this->gridGeometry().bBoxMin()[dimWorld-1];
83
5/10
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 12 times.
✗ Branch 14 not taken.
 60 M[0][0] = 1.0; M[0][1] = 1.0; b[0] = rootCollarPressure_ - soilPressure_;
84
5/10
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 12 times.
✗ Branch 14 not taken.
 60 M[1][0] = C*exp(-C*L); M[1][1] = -C*exp(C*L); b[1] = -density_*gravity;
85
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 M.solve(d_, b);
86
87
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 std::cout << "Computed analytical solution with:\n";
88
3/8
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 12 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 36 std::cout << Fmt::format("A or d1 = {}, B or d2 = {}, c = {}, √c = {}\n", d_[0], d_[1], C*C, C);
89
3/8
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 12 std::cout << Fmt::format("Kx = {}, Kr = {}\n", Kx_, Kr_);
90
91 12 const auto psiCollar = 100*(rootCollarPressure_ - 1e5)/(1000*9.81);
92 12 const auto psiSoil = 100*(soilPressure_ - 1e5)/(1000*9.81);
93 36 const auto psi0 = 100*(analyticalSolution(gridGeometry->bBoxMax()) - 1e5)/(1000*9.81);
94 36 const auto psiTip = 100*(analyticalSolution(gridGeometry->bBoxMin()) - 1e5)/(1000*9.81);
95
3/10
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 12 std::cout << Fmt::format("ψ(0) = {}, ψ(-0.5) = {}, ψ0 = {}, ψs = {} cm\n", psi0, psiTip, psiCollar, psiSoil);
96
97
4/8
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12 times.
✗ Branch 11 not taken.
 36 analyticalPressures_.resize(gridGeometry->gridView().size(0));
98
4/8
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12 times.
✗ Branch 11 not taken.
 36 analyticalHead_.resize(gridGeometry->gridView().size(0));
99
4/8
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1272 times.
✗ Branch 10 not taken.
 2556 for (const auto& element : elements(gridGeometry->gridView()))
100 {
101 3780 const auto eIdx = gridGeometry->elementMapper().index(element);
102 1260 const auto& globalPos = element.geometry().center();
103 1260 analyticalPressures_[eIdx] = analyticalSolution(globalPos);
104 2520 analyticalHead_[eIdx] = 100*(analyticalPressures_[eIdx] - 1e5)/(1000*9.81);
105 }
106 12 }
107
108 /*!
109 * \brief The problem name.
110 */
111 const std::string& name() const
112
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 { return name_; }
113
114 /*!
115 * \brief Specifies which kind of boundary condition should be
116 * used for which equation on a given boundary control volume.
117 */
118 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
119 {
120
4/6
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 24 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 72 BoundaryTypes bcTypes;
121
12/18
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 12 times.
✓ Branch 5 taken 12 times.
✓ Branch 6 taken 24 times.
✓ Branch 7 taken 24 times.
✓ Branch 8 taken 24 times.
✓ Branch 9 taken 24 times.
✓ Branch 10 taken 24 times.
✓ Branch 11 taken 24 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
 216 if (globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps_)
122 bcTypes.setAllDirichlet(); // Dirichlet at root collar
123 else
124 bcTypes.setAllNeumann(); // No-flow at root tips
125 return bcTypes;
126 }
127
128 /*!
129 * \brief Evaluates the boundary conditions for a Dirichlet control volume.
130 */
131 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
132 24 { return PrimaryVariables(rootCollarPressure_); }
133
134 /*!
135 * \brief Evaluates the source term for all phases within a given sub-control volume.
136 * Positive values mean that the conserved quantity is created, negative ones mean that it vanishes.
137 * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
138 */
139 template<class ElementVolumeVariables>
140 2520 NumEqVector source(const Element& element,
141 const FVElementGeometry& fvGeometry,
142 const ElementVolumeVariables& elemVolVars,
143 const SubControlVolume& scv) const
144 {
145 2520 NumEqVector values(0.0);
146 2520 const auto& volVars = elemVolVars[scv];
147 // source term with static soil pressure
148 5040 values[0] = density_*2.0*M_PI*radius_*Kr_*(soilPressure_ - volVars.pressure());
149 // make this volume-specific source
150 2520 values[0] /= volVars.extrusionFactor();
151 2520 return values;
152 }
153
154 template<class SolutionVector>
155 12 void outputL2Norm(const SolutionVector& solution) const
156 {
157 // calculate the discrete L2-Norm and maximum element size
158 12 Scalar l2Norm = 0.0;
159 12 Scalar l2NormNormalization = 0.0;
160 12 Scalar hMax = 0.0;
161
162 12 const auto& gg = this->gridGeometry();
163
1/2
✓ Branch 2 taken 1272 times.
✗ Branch 3 not taken.
 2544 for (const auto& element : elements(gg.gridView()))
164 {
165 2520 const auto eIdx = gg.elementMapper().index(element);
166 1260 const auto geo = element.geometry();
167 1260 const auto globalPos = geo.center();
168 1260 const auto pe = analyticalSolution(globalPos);
169
2/2
✓ Branch 0 taken 670 times.
✓ Branch 1 taken 590 times.
 1260 const auto p = solution[eIdx][0];
170
2/2
✓ Branch 0 taken 670 times.
✓ Branch 1 taken 590 times.
 1260 const auto dV = geo.volume();
171
172 1260 l2Norm += (p - pe)*(p - pe)*dV;
173 1260 l2NormNormalization += pe*pe*dV;
174
2/2
✓ Branch 0 taken 670 times.
✓ Branch 1 taken 590 times.
 1930 hMax = std::max(dV, hMax);
175 }
176
177 12 l2Norm = std::sqrt(l2Norm/l2NormNormalization);
178
179 // write the norm into a log file
180
1/2
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
 12 std::ofstream logFile(this->name() + ".log", std::ios::app);
181
3/6
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 12 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 12 times.
✗ Branch 10 not taken.
 36 logFile << "[ConvergenceTest] L2-norm(pressure) = " << l2Norm << " hMax = " << hMax << std::endl;
182 12 }
183
184 2544 Scalar analyticalSolution(const GlobalPosition& globalPos) const
185 {
186 using std::exp;
187 2544 const auto s = globalPos[dimWorld-1];
188 7632 return soilPressure_ + d_[0]*exp(s*C_) + d_[1]*exp(-s*C_);
189 }
190
191 const std::vector<Scalar>& analyticalPressure() const
192
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 { return analyticalPressures_; }
193
194 const std::vector<Scalar>& analyticalHead() const
195
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 { return analyticalHead_; }
196
197 private:
198 Scalar rootCollarPressure_, soilPressure_;
199 Scalar radius_;
200 Scalar density_;
201 Scalar Kr_, Kx_;
202
203 static constexpr Scalar eps_ = 1e-8;
204 std::string name_;
205
206 // analytical solution
207 Dune::FieldVector<Scalar, 2> d_;
208 Scalar C_;
209
210 std::vector<Scalar> analyticalPressures_;
211 std::vector<Scalar> analyticalHead_;
212 };
213
214 } // end namespace Dumux
215
216 #endif
217