GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/multidomain/embedded/1d3d/root_soil_benchmark/spatialparams_root.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 50 68 73.5%
Functions: 1 5 20.0%
Branches: 136 334 40.7%
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 EmbeddedTests
10 * \brief The spatial parameters class for the root problem
11 */
12 #ifndef DUMUX_TEST_ROOT_SOIL_BENCHMARK_ROOT_SPATIALPARAMS_HH
13 #define DUMUX_TEST_ROOT_SOIL_BENCHMARK_ROOT_SPATIALPARAMS_HH
14
15 #include <vector>
16 #include <utility>
17 #include <memory>
18
19 #include <dumux/common/math.hh>
20 #include <dumux/common/parameters.hh>
21 #include <dumux/io/grid/griddata.hh>
22 #include <dumux/material/components/simpleh2o.hh>
23 #include <dumux/porousmediumflow/fvspatialparams1p.hh>
24
25 namespace Dumux {
26
27 /*!
28 * \ingroup EmbeddedTests
29 * \brief The spatial parameters class for the root problem
30 */
31 template<class GridGeometry, class Scalar>
32 class RootSpatialParams
33 : public FVPorousMediumFlowSpatialParamsOneP<GridGeometry, Scalar, RootSpatialParams<GridGeometry, Scalar>>
34 {
35 using ThisType = RootSpatialParams<GridGeometry, Scalar>;
36 using ParentType = FVPorousMediumFlowSpatialParamsOneP<GridGeometry, Scalar, ThisType>;
37 using Grid = typename GridGeometry::Grid;
38 using GridView = typename GridGeometry::GridView;
39 using Element = typename GridView::template Codim<0>::Entity;
40 using SubControlVolume = typename GridGeometry::SubControlVolume;
41 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
42
43 public:
44 // export permeability type
45 using PermeabilityType = Scalar;
46
47 1 RootSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry,
48 std::shared_ptr<const GridData<Grid>> gridData)
49
9/32
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 1 times.
✗ Branch 24 not taken.
✓ Branch 26 taken 1 times.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
 2 : ParentType(gridGeometry), gridData_(gridData)
50 {
51
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 porosity_ = getParam<Scalar>("Root.SpatialParams.Porosity", 0.4);
52
4/12
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 2 const auto scenarioParamPrefix = "Root.SpatialParams." + getParam<std::string>("Root.SpatialParams.Scenario", "C12a") + ".";
53
54 // read the tabularized root conductivities
55
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 axialRootConductivity_.resize(2);
56
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 radialRootConductivity_.resize(2);
57
8/22
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
 6 axialRootConductivity_[0] = { getParam<std::vector<Scalar>>(scenarioParamPrefix + "AxialConductivites.Order0.Age"),
58 getParam<std::vector<Scalar>>(scenarioParamPrefix + "AxialConductivites.Order0.Kx") };
59
8/22
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
 6 axialRootConductivity_[1] = { getParam<std::vector<Scalar>>(scenarioParamPrefix + "AxialConductivites.Order1.Age"),
60 getParam<std::vector<Scalar>>(scenarioParamPrefix + "AxialConductivites.Order1.Kx") };
61
8/22
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
 6 radialRootConductivity_[0] = { getParam<std::vector<Scalar>>(scenarioParamPrefix + "RadialConductivites.Order0.Age"),
62 getParam<std::vector<Scalar>>(scenarioParamPrefix + "RadialConductivites.Order0.Kr") };
63
8/24
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
 6 radialRootConductivity_[1] = { getParam<std::vector<Scalar>>(scenarioParamPrefix + "RadialConductivites.Order1.Age"),
64 getParam<std::vector<Scalar>>(scenarioParamPrefix + "RadialConductivites.Order1.Kr") };
65
66 // sanity checks
67
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 1 times.
 5 for (const auto& k : axialRootConductivity_)
68
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 6 if (k.first.size() != k.second.size())
69 ✗ DUNE_THROW(Dune::IOError, "AxialConductivites.Age and AxialConductivites.Kx have to have the same length!");
70
71
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 1 times.
 5 for (const auto& k : radialRootConductivity_)
72
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 6 if (k.first.size() != k.second.size())
73 ✗ DUNE_THROW(Dune::IOError, "RadialConductivites.Age and RadialConductivites.Kr have to have the same length!");
74
75 2 const auto& gv = gridGeometry->gridView();
76
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 radii_.resize(gv.size(0));
77
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 age_.resize(gv.size(0));
78
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 order_.resize(gv.size(0));
79
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Kx_.resize(gv.size(0));
80
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Kr_.resize(gv.size(0));
81
3/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2 const auto paramType = getParam<std::string>("Root.Grid.DGFParam", "Radius");
82
4/4
✓ Branch 1 taken 1170 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 1170 times.
✓ Branch 4 taken 1 times.
 2341 for (const auto& element : elements(gv))
83 {
84 3510 const auto eIdx = gridGeometry->elementMapper().index(element);
85 1170 auto level0element = element;
86
4/4
✓ Branch 0 taken 1170 times.
✓ Branch 1 taken 1170 times.
✓ Branch 2 taken 1170 times.
✓ Branch 3 taken 1170 times.
 4680 while (level0element.hasFather())
87 3510 level0element = level0element.father();
88
89
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1170 times.
 1170 if (paramType == "Surface")
90 {
91 // see e.g. lupine.dgf
92 ✗ const auto rootLength = element.geometry().volume();
93 ✗ const auto rootSurface = gridData->parameters(level0element)[2]/(1 << element.level());
94 ✗ radii_[eIdx] = rootSurface / ( rootLength * 2.0 * M_PI );
95 }
96
1/2
✓ Branch 1 taken 1170 times.
✗ Branch 2 not taken.
 1170 else if (paramType == "Radius")
97 {
98
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1169 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1170 static const auto radiusIdx = getParam<int>("Root.Grid.DGFRadiusIndex", 4);
99
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1169 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1170 static const auto radiusScaling = getParam<double>("Root.Grid.DGFRadiusScaling", 1.0);
100
2/4
✓ Branch 1 taken 1170 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1170 times.
✗ Branch 5 not taken.
 2340 radii_[eIdx] = gridData->parameters(level0element)[radiusIdx]*radiusScaling; // read radius from grid file
101 }
102 else
103 ✗ DUNE_THROW(Dune::NotImplemented, "Unknown DGF Parameter " << paramType);
104
105 // read age
106
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1169 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1170 static const auto ageIdx = getParam<int>("Root.Grid.DGFAgeIndex", 7);
107
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1169 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1170 static const auto ageScaling = getParam<int>("Root.Grid.DGFAgeScaling", 1.0);
108
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1169 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1170 static const auto currentPlantAge = getParam<int>("Root.Grid.DGFPlantAge", 8.0);
109
6/10
✓ Branch 1 taken 1170 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1170 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1170 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1170 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 1169 times.
 3510 age_[eIdx] = std::max(0.0, currentPlantAge - gridData->parameters(level0element)[ageIdx]*ageScaling);
110
111 // read order
112
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1169 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1170 static const auto orderIdx = getParam<int>("Root.Grid.DGFOrderIndex", 0);
113
6/8
✓ Branch 1 taken 1170 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1170 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 906 times.
✓ Branch 7 taken 264 times.
✓ Branch 8 taken 906 times.
✓ Branch 9 taken 264 times.
 2340 order_[eIdx] = gridData->parameters(level0element)[orderIdx];
114
115 // permeability
116
4/4
✓ Branch 0 taken 906 times.
✓ Branch 1 taken 264 times.
✓ Branch 2 taken 82 times.
✓ Branch 3 taken 1088 times.
 2076 const auto order = static_cast<std::size_t>(std::min(std::max(0.0, order_[eIdx]), 1.0));
117 1170 Kx_[eIdx] = computeAxialConductivity(order, age_[eIdx]);
118 1170 Kr_[eIdx] = computeRadialConductivity(order, age_[eIdx]);
119 }
120 1 }
121
122 /*!
123 * \brief Returns how much the domain is extruded at a given sub-control volume.
124 *
125 * The extrusion factor here extrudes the 1d line to a circular tube with
126 * cross-section area pi*r^2.
127 */
128 template<class ElementSolution>
129 ✗ Scalar extrusionFactor(const Element &element,
130 const SubControlVolume &scv,
131 const ElementSolution& elemSol) const
132 {
133 ✗ const auto eIdx = this->gridGeometry().elementMapper().index(element);
134 ✗ const auto radius = this->radius(eIdx);
135 ✗ return M_PI*radius*radius;
136 }
137
138 /*!
139 * \brief Returns the intrinsic permeability for the current sub-control volume in [m^2].
140 *
141 * \note Kx has units [m^4/(Pa*s)] so we have to divide by the cross-section area
142 * and multiply with a characteristic viscosity
143 */
144 template<class ElementSolution>
145 ✗ PermeabilityType permeability(const Element& element,
146 const SubControlVolume& scv,
147 const ElementSolution& elemSol) const
148 {
149 ✗ const auto eIdx = this->gridGeometry().elementMapper().index(element);
150 ✗ const Scalar r = radius(eIdx);
151 ✗ return Kx(eIdx) / (M_PI*r*r) * Components::SimpleH2O<Scalar>::liquidViscosity(285.15, 1e5);
152 }
153
154 /*!
155 * \brief Returns the radius of the circular pipe for the current sub-control volume in [m].
156 *
157 * \param eIdxGlobal the index of the element
158 */
159 Scalar radius(std::size_t eIdxGlobal) const
160 {
161
2/4
✓ Branch 1 taken 1329 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1329 times.
✗ Branch 5 not taken.
 5316 return radii_[eIdxGlobal];
162 }
163
164 /*!
165 * \brief Returns the radial permeability.
166 *
167 * \param eIdxGlobal the index of the element
168 */
169 Scalar Kr(std::size_t eIdxGlobal) const
170 {
171 9458934 return Kr_[eIdxGlobal];
172 }
173
174 /*!
175 * \brief Returns the radial permeability.
176 *
177 * \param eIdxGlobal the index of the element
178 */
179 Scalar Kx(std::size_t eIdxGlobal) const
180 {
181
8/8
✓ Branch 0 taken 1681 times.
✓ Branch 1 taken 150 times.
✓ Branch 2 taken 1681 times.
✓ Branch 3 taken 150 times.
✓ Branch 4 taken 199 times.
✓ Branch 5 taken 18 times.
✓ Branch 6 taken 199 times.
✓ Branch 7 taken 18 times.
 4096 return Kx_[eIdxGlobal];
182 }
183
184 // get radii for output
185
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const std::vector<Scalar>& getRadii() const { return radii_; }
186 // get ages for output
187 const std::vector<Scalar>& getAges() const { return age_; }
188 // get orders for output
189
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const std::vector<Scalar>& getOrders() const { return order_; }
190 // get kx for output
191
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const std::vector<Scalar>& getKx() const { return Kx_; }
192 // get kr for output
193
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const std::vector<Scalar>& getKr() const { return Kr_; }
194
195 /*!
196 * \brief Returns the porosity \f$[-]\f$.
197 *
198 * \param globalPos the scv center
199 */
200 ✗ Scalar porosityAtPos(const GlobalPosition& globalPos) const
201 {
202 ✗ return porosity_;
203 }
204
205 //! compute the radial conductivity (m/Pa/s) given the segment age in days
206 double computeRadialConductivity(int order, double age) const
207 {
208 2340 return interpolate<InterpolationPolicy::LinearTable>(age, radialRootConductivity_[order]);
209 }
210
211 //! compute the axial conductivity in (m^4/Pa/s) given the segment age in days
212 double computeAxialConductivity(int order, double age) const
213 {
214 2340 return interpolate<InterpolationPolicy::LinearTable>(age, axialRootConductivity_[order]);
215 }
216
217 /*!
218 * \brief Returns the temperature within the domain in [K].
219 */
220 ✗ Scalar temperatureAtPos(const GlobalPosition& globalPos) const
221 ✗ { return 273.15 + 10.0; }
222
223 private:
224 std::shared_ptr<const GridData<Grid>> gridData_;
225 Scalar porosity_;
226
227 //! Tabularized root conductivity functions (pairs of x, y) for each order
228 std::vector<std::pair<std::vector<double>, std::vector<double>>> axialRootConductivity_;
229 std::vector<std::pair<std::vector<double>, std::vector<double>>> radialRootConductivity_;
230
231 std::vector<Scalar> radii_;
232 std::vector<Scalar> age_;
233 std::vector<Scalar> order_;
234 std::vector<Scalar> Kx_;
235 std::vector<Scalar> Kr_;
236 };
237
238 } // end namespace Dumux
239
240 #endif
241