GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/multidomain/embedded/1d3d/1p_richards/problem_root.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 43 50 86.0%
Functions: 6 10 60.0%
Branches: 57 124 46.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 EmbeddedTests
10 * \brief A test problem for the one-phase root model:
11 * Sap is flowing through a 1d network root xylem.
12 */
13
14 #ifndef DUMUX_ROOT_PROBLEM_HH
15 #define DUMUX_ROOT_PROBLEM_HH
16
17 #include <dumux/common/boundarytypes.hh>
18 #include <dumux/common/parameters.hh>
19 #include <dumux/common/properties.hh>
20 #include <dumux/common/numeqvector.hh>
21
22 #include <dumux/porousmediumflow/problem.hh>
23 #include <dumux/multidomain/embedded/couplingmanager1d3d_projection.hh>
24 namespace Dumux {
25
26 /*!
27 * \ingroup EmbeddedTests
28 * \brief Exact solution 1D-3D.
29 */
30 template <class TypeTag>
31 class RootProblem : public PorousMediumFlowProblem<TypeTag>
32 {
33 using ParentType = PorousMediumFlowProblem<TypeTag>;
34 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
35 using PointSource = GetPropType<TypeTag, Properties::PointSource>;
36 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
37 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
38 using NeumannFluxes = Dumux::NumEqVector<PrimaryVariables>;
39 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
40 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
41 using GridView = typename GridGeometry::GridView;
42 using FVElementGeometry = typename GridGeometry::LocalView;
43 using SubControlVolume = typename GridGeometry::SubControlVolume;
44 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
45 using GlobalPosition = typename GridGeometry::GlobalCoordinate;
46 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
47 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
48 using Element = typename GridView::template Codim<0>::Entity;
49 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
50
51 public:
52
53 template<class SpatialParams>
54 3 RootProblem(std::shared_ptr<const GridGeometry> gridGeometry,
55 std::shared_ptr<SpatialParams> spatialParams,
56 std::shared_ptr<CouplingManager> couplingManager)
57 : ParentType(gridGeometry, spatialParams, "Root")
58
6/20
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 3 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 3 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
 12 , couplingManager_(couplingManager)
59 {
60 // read parameters from input file
61
9/26
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 3 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 3 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 3 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 3 times.
✗ 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.
✗ Branch 31 not taken.
 6 name_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
62
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 transpirationRate_ = getParam<Scalar>("BoundaryConditions.TranspirationRate");
63 3 }
64
65 /*!
66 * \name Problem parameters
67 */
68 // \{
69
70 /*!
71 * \brief The problem name.
72 *
73 * This is used as a prefix for files generated by the simulation.
74 */
75 const std::string& name() const
76
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 { return name_; }
77
78 // \}
79 /*!
80 * \name Boundary conditions
81 */
82 // \{
83
84 /*!
85 * \brief Specifies which kind of boundary condition should be
86 * used for which equation on a given boundary segment.
87 *
88 * \param globalPos The global position
89 */
90 ✗ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
91 {
92
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 7208 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3827 times.
 11035 BoundaryTypes values;
93
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 7208 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3827 times.
 11035 values.setAllNeumann();
94 ✗ return values;
95 }
96
97 /*!
98 * \brief Evaluates the boundary conditions for a Dirichlet control volume.
99 *
100 * \param globalPos The global position
101 *
102 * For this method, the \a values parameter stores primary variables.
103 */
104 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
105 ✗ { return initialAtPos(globalPos); }
106
107
108 /*!
109 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
110 *
111 * For this method, the \a priVars parameter stores the mass flux
112 * in normal direction of each component. Negative values mean
113 * influx.
114 */
115 template<class ElementVolumeVariables, class ElementFluxVarsCache>
116 ✗ NeumannFluxes neumann(const Element& element,
117 const FVElementGeometry& fvGeometry,
118 const ElementVolumeVariables& elemVolvars,
119 const ElementFluxVarsCache& elemFluxVarsCache,
120 const SubControlVolumeFace& scvf) const
121 {
122
2/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 232 times.
✓ Branch 3 taken 6976 times.
 7208 NeumannFluxes values(0.0);
123
12/24
✗ 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 taken 232 times.
✓ Branch 13 taken 6976 times.
✓ Branch 14 taken 232 times.
✓ Branch 15 taken 6976 times.
✓ Branch 16 taken 232 times.
✓ Branch 17 taken 6976 times.
✓ Branch 18 taken 232 times.
✓ Branch 19 taken 6976 times.
✓ Branch 20 taken 232 times.
✓ Branch 21 taken 6976 times.
✓ Branch 22 taken 232 times.
✓ Branch 23 taken 6976 times.
 43248 if (scvf.center()[2] + eps_ > this->gridGeometry().bBoxMax()[2])
124 {
125 696 const auto r = this->spatialParams().radius(scvf.insideScvIdx());
126 232 values[Indices::conti0EqIdx] = transpirationRate_/(M_PI*r*r)/scvf.area();
127 }
128 ✗ return values;
129
130 }
131
132 // \}
133
134 /*!
135 * \name Volume terms
136 */
137 // \{
138
139 /*!
140 * \brief Applies a vector of point sources which are possibly solution dependent.
141 *
142 * \param pointSources A vector of PointSource s that contain
143 source values for all phases and space positions.
144 *
145 * For this method, the \a values method of the point source
146 * has to return the absolute mass rate in kg/s. Positive values mean
147 * that mass is created, negative ones mean that it vanishes.
148 */
149 void addPointSources(std::vector<PointSource>& pointSources) const
150
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 { pointSources = this->couplingManager().lowDimPointSources(); }
151
152 /*!
153 * \brief Evaluates the point sources (added by addPointSources)
154 * for all phases within a given sub control volume.
155 *
156 * This is the method for the case where the point source is
157 * solution dependent and requires some quantities that
158 * are specific to the fully-implicit method.
159 *
160 * \param source A single point source
161 * \param element The finite element
162 * \param fvGeometry The finite-volume geometry
163 * \param elemVolVars All volume variables for the element
164 * \param scv The sub control volume within the element
165 *
166 * For this method, the \a values() method of the point sources returns
167 * the absolute rate mass generated or annihilated in kg/s. Positive values mean
168 * that mass is created, negative ones mean that it vanishes.
169 */
170 template<class ElementVolumeVariables>
171 14167512 void pointSource(PointSource& source,
172 const Element &element,
173 const FVElementGeometry& fvGeometry,
174 const ElementVolumeVariables& elemVolVars,
175 const SubControlVolume &scv) const
176 {
177 // compute source at every integration point
178 42502536 const Scalar pressure3D = this->couplingManager().bulkPriVars(source.id())[Indices::pressureIdx];
179 42502536 const Scalar pressure1D = this->couplingManager().lowDimPriVars(source.id())[Indices::pressureIdx];
180
181 42502536 const auto lowDimElementIdx = this->couplingManager().pointSourceData(source.id()).lowDimElementIdx();
182 28335024 const Scalar Kr = this->spatialParams().Kr(lowDimElementIdx);
183
184 // sink defined as radial flow Jr * density [m^2 s-1]* [kg m-3]
185 14167512 const auto density = 1000;
186 14167512 Scalar sourceValue = Kr *(pressure3D - pressure1D)*density;
187
188 // For the projection method, we are integrating over the two-dimensional root surface
189 // so surface is included in the weight/integration element.
190 // All other currently implemented schemes are implicit interface schemes
191 // that assume cylindrical segments, so we multiply with the cylinder surface here
192 if constexpr(CouplingManager::couplingMode != Embedded1d3dCouplingMode::projection)
193 {
194 5156844 const Scalar rootRadius = this->spatialParams().radius(lowDimElementIdx);
195 1718948 sourceValue *= 2*M_PI*rootRadius;
196 }
197
198 28335024 source = sourceValue*source.quadratureWeight()*source.integrationElement();
199 14167512 }
200
201 /*!
202 * \brief Evaluates the initial value for a control volume.
203 *
204 * For this method, the \a priVars parameter stores primary
205 * variables.
206 */
207 ✗ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
208 3520 { return PrimaryVariables(0.0); }
209
210 // \}
211
212 //! Called after every time step
213 //! Output the total global exchange term
214 21 void computeSourceIntegral(const SolutionVector& sol, const GridVariables& gridVars)
215 {
216 21 PrimaryVariables source(0.0);
217 42 auto fvGeometry = localView(this->gridGeometry());
218 42 auto elemVolVars = localView(gridVars.curGridVolVars());
219
4/4
✓ Branch 3 taken 24640 times.
✓ Branch 4 taken 21 times.
✓ Branch 5 taken 24640 times.
✓ Branch 6 taken 21 times.
 24703 for (const auto& element : elements(this->gridGeometry().gridView()))
220 {
221 24640 fvGeometry.bindElement(element);
222 24640 elemVolVars.bindElement(element, fvGeometry, sol);
223
5/6
✓ Branch 2 taken 24640 times.
✓ Branch 3 taken 24640 times.
✓ Branch 4 taken 24640 times.
✓ Branch 5 taken 24640 times.
✓ Branch 7 taken 24640 times.
✗ Branch 8 not taken.
 98560 for (auto&& scv : scvs(fvGeometry))
224 {
225
1/2
✓ Branch 1 taken 24640 times.
✗ Branch 2 not taken.
 24640 auto pointSources = this->scvPointSources(element, fvGeometry, elemVolVars, scv);
226 24640 pointSources *= scv.volume()*elemVolVars[scv].extrusionFactor();
227 24640 source += pointSources;
228 }
229 }
230
231
1/2
✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
 42 std::cout << "Global integrated source (root): " << source << " (kg/s) / "
232
2/4
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 21 times.
✗ Branch 6 not taken.
 42 << source*3600*24*1000 << " (g/day)" << '\n';
233 21 }
234
235 /*!
236 * \brief Adds additional VTK output data to the VTKWriter.
237 *
238 * Function is called by the output module on every write.
239 */
240 template<class VtkOutputModule>
241 3 void addVtkOutputFields(VtkOutputModule& vtk) const
242 {
243
6/14
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 3 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 6 vtk.addField(this->spatialParams().getRadii(), "radius");
244 3 }
245
246 //! Get the coupling manager
247 const CouplingManager& couplingManager() const
248
2/4
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 28335030 { return *couplingManager_; }
249
250 private:
251 Scalar transpirationRate_;
252
253 static constexpr Scalar eps_ = 1.5e-7;
254 std::string name_;
255
256 std::shared_ptr<CouplingManager> couplingManager_;
257 };
258
259 } // end namespace Dumux
260
261 #endif
262