GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/1p/network1d3d/problem.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 49 54 90.7%
Functions: 6 14 42.9%
Branches: 66 112 58.9%
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 A test problem for the 1p model. A pipe system with circular cross-section
11 * and a branching point embedded in a three-dimensional world
12 */
13
14 #ifndef DUMUX_ONEP_TUBES_TEST_PROBLEM_HH
15 #define DUMUX_ONEP_TUBES_TEST_PROBLEM_HH
16
17 #include <dune/geometry/quadraturerules.hh>
18 #include <dune/localfunctions/lagrange/lagrangelfecache.hh>
19
20 #include <dumux/common/properties.hh>
21 #include <dumux/common/parameters.hh>
22 #include <dumux/common/boundarytypes.hh>
23 #include <dumux/common/numeqvector.hh>
24
25 #include <dumux/discretization/elementsolution.hh>
26 #include <dumux/porousmediumflow/problem.hh>
27
28 namespace Dumux {
29
30 /*!
31 * \ingroup OnePTests
32 * \brief A test problem for the 1p model: A pipe system with circular cross-section
33 * and a branching point embedded in a three-dimensional world
34 */
35 template <class TypeTag>
36 class TubesTestProblem : public PorousMediumFlowProblem<TypeTag>
37 {
38 using ParentType = PorousMediumFlowProblem<TypeTag>;
39 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
40 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
41 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
42
43 // Grid and world dimension
44 static const int dim = GridView::dimension;
45 static const int dimWorld = GridView::dimensionworld;
46
47 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
48 enum {
49 // indices of the primary variables
50 conti0EqIdx = Indices::conti0EqIdx,
51 pressureIdx = Indices::pressureIdx
52 };
53
54 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
55 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
56 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
57 using Element = typename GridView::template Codim<0>::Entity;
58 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
59 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
60 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
61 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
62 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
63
64 enum { isBox = GetPropType<TypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethods::box };
65
66 public:
67 12 TubesTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
68
8/22
✓ 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 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
 36 : ParentType(gridGeometry)
69 {
70
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");
71
72 //get hMax_ of the grid
73 12 hMax_ = 0.0;
74
5/6
✓ Branch 3 taken 4032 times.
✓ Branch 4 taken 12 times.
✓ Branch 5 taken 4032 times.
✓ Branch 6 taken 12 times.
✓ Branch 8 taken 4032 times.
✗ Branch 9 not taken.
 8100 for (const auto& element : elements(gridGeometry->gridView()))
75
5/6
✓ Branch 1 taken 4032 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3528 times.
✓ Branch 4 taken 504 times.
✓ Branch 5 taken 3528 times.
✓ Branch 6 taken 504 times.
 7560 hMax_ = std::max(element.geometry().volume(), hMax_);
76 12 }
77
78 /*!
79 * \name Problem parameters
80 */
81 // \{
82
83 /*!
84 * \brief The problem name.
85 *
86 * This is used as a prefix for files generated by the simulation.
87 */
88 const std::string& name() const
89 {
90
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 return name_;
91 }
92
93 // \}
94 /*!
95 * \name Boundary conditions
96 */
97 // \{
98
99 /*!
100 * \brief Specifies which kind of boundary condition should be
101 * used for which equation on a given boundary control volume.
102 *
103 * \param globalPos The position of the center of the finite volume
104 */
105 ✗ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
106 {
107
3/8
✓ Branch 0 taken 85 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 90 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 180 times.
✗ Branch 7 not taken.
 355 BoundaryTypes bcTypes;
108
3/8
✓ Branch 0 taken 85 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 90 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 180 times.
✗ Branch 7 not taken.
 355 bcTypes.setAllDirichlet();
109 ✗ return bcTypes;
110 }
111
112 /*!
113 * \brief Evaluates the boundary conditions for a Dirichlet control volume.
114 *
115 * \param globalPos The center of the finite volume which ought to be set.
116 *
117 * For this method, the \a values parameter stores primary variables.
118 */
119 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
120 {
121 175 return PrimaryVariables(0.0);
122 }
123
124 // \}
125
126 /*!
127 * \name Volume terms
128 */
129 // \{
130
131 /*!
132 * \brief Evaluates the source term for all phases within a given
133 * sub-control volume.
134 *
135 * This is the method for the case where the source term is
136 * potentially solution dependent and requires some quantities that
137 * are specific to the fully-implicit method.
138 *
139 * \param element The finite element
140 * \param fvGeometry The finite-volume geometry
141 * \param elemVolVars All volume variables for the element
142 * \param scv The sub-control volume
143 *
144 * For this method, the \a values parameter stores the conserved quantity rate
145 * generated or annihilated per volume unit. Positive values mean
146 * that the conserved quantity is created, negative ones mean that it vanishes.
147 * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
148 */
149 45312 NumEqVector source(const Element &element,
150 const FVElementGeometry& fvGeometry,
151 const ElementVolumeVariables& elemVolVars,
152 const SubControlVolume &scv) const
153 {
154
2/2
✓ Branch 0 taken 4152 times.
✓ Branch 1 taken 29064 times.
 45312 NumEqVector source(0.0);
155
156
2/2
✓ Branch 0 taken 4152 times.
✓ Branch 1 taken 29064 times.
 45312 const auto& globalPos = scv.center();
157
2/2
✓ Branch 0 taken 4152 times.
✓ Branch 1 taken 29064 times.
 45312 const auto& volVars = elemVolVars[scv];
158
4/4
✓ Branch 0 taken 5664 times.
✓ Branch 1 taken 39648 times.
✓ Branch 2 taken 5664 times.
✓ Branch 3 taken 39648 times.
 90624 const auto K = this->spatialParams().permeability(element, scv, EmptyElementSolution{});
159
160 using std::sin;
161
4/4
✓ Branch 0 taken 5664 times.
✓ Branch 1 taken 39648 times.
✓ Branch 2 taken 5664 times.
✓ Branch 3 taken 39648 times.
 90624 if (globalPos[2] > 0.5 - eps_)
162 11328 source[conti0EqIdx] = K/volVars.viscosity()*volVars.density()
163 11328 *4.0*4.0*M_PI*M_PI*sin(4.0*M_PI*globalPos[2]);
164 else
165 // the "/3.0" stems from the coordinate transformations on the lower branches
166 79296 source[conti0EqIdx] = K/volVars.viscosity()*volVars.density()
167 79296 *4.0*4.0*M_PI*M_PI*sin(4.0*M_PI*globalPos[2])/3.0;
168
169 45312 return source;
170 }
171
172 /*!
173 * \brief Evaluates the initial value for a control volume.
174 *
175 * For this method, the \a priVars parameter stores primary
176 * variables.
177 */
178 ✗ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
179 {
180 4038 return PrimaryVariables(0.0);
181 }
182
183 // \}
184
185 12 void outputL2Norm(const SolutionVector solution) const
186 {
187 // calculate the discrete L2-Norm
188 12 Scalar lTwoNorm = 0.0;
189
190 // get the Gaussian quadrature rule for intervals
191 12 const auto& quad = Dune::QuadratureRules<Scalar, dim>::rule(Dune::GeometryTypes::line, 1);
192
193 12 const auto& gg = this->gridGeometry();
194
4/4
✓ Branch 2 taken 4032 times.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 4032 times.
✓ Branch 5 taken 12 times.
 4056 for (const auto& element : elements(gg.gridView()))
195 {
196 8064 const auto eIdx = gg.elementMapper().index(element);
197 4032 const auto geometry = element.geometry();
198
4/4
✓ Branch 0 taken 4032 times.
✓ Branch 1 taken 4032 times.
✓ Branch 2 taken 4032 times.
✓ Branch 3 taken 4032 times.
 20160 for(auto&& qp : quad)
199 {
200 8064 const auto globalPos = geometry.global(qp.position());
201 8064 const Scalar pe = exactPressure_(globalPos);
202 4032 const Scalar integrationElement = geometry.integrationElement(qp.position());
203 4032 Scalar p = 0.0;
204 if (!isBox)
205 2016 p = solution[eIdx][pressureIdx];
206 else
207 {
208 // do interpolation with ansatz functions
209
1/2
✓ Branch 1 taken 2016 times.
✗ Branch 2 not taken.
 4032 std::vector<Dune::FieldVector<Scalar, 1> > shapeValues;
210
1/4
✓ Branch 2 taken 2016 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2016 const auto& localFiniteElement = feCache_.get(geometry.type());
211
2/4
✓ Branch 1 taken 2016 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2016 times.
✗ Branch 5 not taken.
 4032 localFiniteElement.localBasis().evaluateFunction(qp.position(), shapeValues);
212
4/4
✓ Branch 0 taken 4032 times.
✓ Branch 1 taken 2016 times.
✓ Branch 2 taken 4032 times.
✓ Branch 3 taken 2016 times.
 10080 for (unsigned int i = 0; i < shapeValues.size(); ++i)
213
2/4
✓ Branch 1 taken 4032 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4032 times.
✗ Branch 5 not taken.
 8064 p += shapeValues[i]*solution[gg.dofMapper().subIndex(element, i, dim)][pressureIdx];
214 }
215 4032 lTwoNorm += (p - pe)*(p - pe)*qp.weight()*integrationElement;
216 }
217 }
218 12 lTwoNorm = std::sqrt(lTwoNorm);
219
220 // write the norm into a log file
221 24 std::ofstream logFile;
222
3/8
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 12 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 12 logFile.open(this->name() + ".log", std::ios::app);
223
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) = " << lTwoNorm << " hMax = " << hMax_ << std::endl;
224
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 logFile.close();
225 12 }
226
227 private:
228
229 ✗ Scalar exactPressure_(const GlobalPosition& globalPos) const
230 {
231 using std::sin;
232 8064 return sin(4.0*M_PI*globalPos[2]);
233 }
234
235 static constexpr Scalar eps_ = 1e-8;
236 std::string name_;
237
238 Scalar hMax_;
239
240 typename Dune::LagrangeLocalFiniteElementCache<Scalar, Scalar, dim, 1> feCache_;
241 };
242
243 } // end namespace Dumux
244
245 #endif
246