GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/multidomain/embedded/1d3d/1p_1p/problem_bloodflow.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 58 65 89.2%
Functions: 13 21 61.9%
Branches: 89 179 49.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 A test problem for the one-phase blood flow model:
11 * Blood is flowing through a 1d network grid.
12 */
13
14 #ifndef DUMUX_BLOOD_FLOW_PROBLEM_HH
15 #define DUMUX_BLOOD_FLOW_PROBLEM_HH
16
17 #include <dumux/common/boundarytypes.hh>
18 #include <dumux/common/parameters.hh>
19 #include <dumux/common/properties.hh>
20
21 #include <dumux/porousmediumflow/problem.hh>
22
23 #include <dumux/multidomain/embedded/couplingmanager1d3d.hh> // for coupling mode
24
25 namespace Dumux {
26
27 /*!
28 * \ingroup EmbeddedTests
29 * \brief Exact solution 1D-3D.
30 */
31 template <class TypeTag>
32 class BloodFlowProblem : public PorousMediumFlowProblem<TypeTag>
33 {
34 using ParentType = PorousMediumFlowProblem<TypeTag>;
35 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
36 using PointSource = GetPropType<TypeTag, Properties::PointSource>;
37 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
38 using PrimaryVariables = GetPropType<TypeTag, Properties::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 SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
45 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
46 using Element = typename GridView::template Codim<0>::Entity;
47 using GlobalPosition = typename GridGeometry::GlobalCoordinate;
48 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
49
50 public:
51 19 BloodFlowProblem(std::shared_ptr<const GridGeometry> gridGeometry,
52 std::shared_ptr<CouplingManager> couplingManager)
53 : ParentType(gridGeometry, "Vessel")
54
5/18
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 19 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 19 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 19 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
 57 , couplingManager_(couplingManager)
55 {
56 //read parameters from input file
57
9/26
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 19 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 19 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 19 times.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 19 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 19 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 19 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 19 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.
 38 name_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
58
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 19 p_in_ = getParam<Scalar>("BoundaryConditions1D.PressureInput");
59
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 19 delta_p_ = getParam<Scalar>("BoundaryConditions1D.DeltaPressure");
60
4/7
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 17 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 40 exactPressure_.resize(this->gridGeometry().numDofs());
61 38 auto fvGeometry = localView(this->gridGeometry());
62
4/4
✓ Branch 3 taken 365 times.
✓ Branch 4 taken 19 times.
✓ Branch 5 taken 365 times.
✓ Branch 6 taken 19 times.
 422 for (const auto& element : elements(this->gridGeometry().gridView()))
63 {
64 365 fvGeometry.bindElement(element);
65
6/6
✓ Branch 0 taken 80 times.
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 405 times.
✓ Branch 3 taken 365 times.
✓ Branch 4 taken 325 times.
✓ Branch 5 taken 325 times.
 1215 for (auto&& scv : scvs(fvGeometry))
66 1540 exactPressure_[scv.dofIndex()] = exactSolution(scv.dofPosition());
67 }
68 19 }
69
70
71
72 /*!
73 * \name Problem parameters
74 */
75 // \{
76
77 /*!
78 * \brief The problem name.
79 *
80 * This is used as a prefix for files generated by the simulation.
81 */
82 const std::string& name() const
83
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 19 { return name_; }
84
85 // \}
86 /*!
87 * \name Boundary conditions
88 */
89 // \{
90
91 /*!
92 * \brief Specifies which kind of boundary condition should be
93 * used for which equation on a given boundary segment.
94 *
95 * \param globalPos The global position
96 */
97 ✗ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
98 {
99
2/4
✓ Branch 0 taken 72 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 68 times.
✗ Branch 3 not taken.
 140 BoundaryTypes values;
100
2/4
✓ Branch 0 taken 72 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 68 times.
✗ Branch 3 not taken.
 140 values.setAllDirichlet();
101 ✗ return values;
102 }
103
104 /*!
105 * \brief Evaluates the boundary conditions for a Dirichlet control volume.
106 *
107 * \param globalPos The global position
108 *
109 * For this method, the \a values parameter stores primary variables.
110 */
111 ✗ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
112 {
113 72 PrimaryVariables values;
114
4/12
✗ 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 taken 36 times.
✓ Branch 9 taken 36 times.
✓ Branch 10 taken 36 times.
✓ Branch 11 taken 36 times.
 144 if(globalPos[2] > 0.5)
115 ✗ values[Indices::pressureIdx] = p_in_;
116 else
117 36 values[Indices::pressureIdx] = p_in_-delta_p_;
118 ✗ return values;
119 }
120
121
122 // \}
123
124 /*!
125 * \name Volume terms
126 */
127 // \{
128
129 /*!
130 * \brief Applies a vector of point sources which are possibly solution dependent.
131 *
132 * \param pointSources A vector of PointSource s that contain
133 source values for all phases and space positions.
134 *
135 * For this method, the \a values method of the point source
136 * has to return the absolute mass rate in kg/s. Positive values mean
137 * that mass is created, negative ones mean that it vanishes.
138 */
139 void addPointSources(std::vector<PointSource>& pointSources) const
140
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 19 { pointSources = this->couplingManager().lowDimPointSources(); }
141
142 /*!
143 * \brief Evaluates the point sources (added by addPointSources)
144 * for all phases within a given sub-control-volume.
145 *
146 * This is the method for the case where the point source is
147 * solution dependent and requires some quantities that
148 * are specific to the fully-implicit method.
149 *
150 * \param source A single point source
151 * \param element The finite element
152 * \param fvGeometry The finite-volume geometry
153 * \param elemVolVars All volume variables for the element
154 * \param scv The sub-control volume within the element
155 *
156 * For this method, the \a values() method of the point sources returns
157 * the absolute rate mass generated or annihilated in kg/s. Positive values mean
158 * that mass is created, negative ones mean that it vanishes.
159 */
160 template<class ElementVolumeVariables>
161 2097520 void pointSource(PointSource& source,
162 const Element &element,
163 const FVElementGeometry& fvGeometry,
164 const ElementVolumeVariables& elemVolVars,
165 const SubControlVolume &scv) const
166 {
167 // compute source at every integration point
168 6292560 const Scalar pressure1D = this->couplingManager().lowDimPriVars(source.id())[Indices::pressureIdx];
169 6292560 const Scalar pressure3D = this->couplingManager().bulkPriVars(source.id())[Indices::pressureIdx];
170
171 // calculate the source
172 6292560 const Scalar radius = this->couplingManager().radius(source.id());
173 2097520 const Scalar beta = 2*M_PI/(2*M_PI + std::log(radius));
174 2097520 Scalar sourceValue = beta*(pressure3D - pressure1D);
175 if constexpr(CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel)
176 12540 sourceValue *= this->couplingManager().fluxScalingFactor(source.id());
177
178 4195040 source = sourceValue*source.quadratureWeight()*source.integrationElement();
179 2097520 }
180
181 /*!
182 * \brief Evaluates the initial value for a control volume.
183 *
184 * For this method, the \a priVars parameter stores primary
185 * variables.
186 */
187 ✗ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
188 367 { return PrimaryVariables(0.0); }
189
190 // \}
191
192 //! The exact pressure solution
193 ✗ Scalar exactSolution(const GlobalPosition &globalPos) const
194 1540 { return 1 + globalPos[2]; }
195
196 //! Called after every time step
197 //! Output the total global exchange term
198 19 Scalar computeSourceIntegral(const SolutionVector& sol, const GridVariables& gridVars)
199 {
200
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 19 SolutionVector source(sol.size());
201
2/6
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 19 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 38 SolutionVector sourceExact(sol.size());
202 19 Scalar volume = 0.0;
203 38 auto fvGeometry = localView(this->gridGeometry());
204
1/2
✓ Branch 2 taken 17 times.
✗ Branch 3 not taken.
 55 auto elemVolVars = localView(gridVars.curGridVolVars());
205
206
4/4
✓ Branch 3 taken 365 times.
✓ Branch 4 taken 19 times.
✓ Branch 5 taken 365 times.
✓ Branch 6 taken 19 times.
 422 for (const auto& element : elements(this->gridGeometry().gridView()))
207 {
208 365 fvGeometry.bindElement(element);
209 365 elemVolVars.bindElement(element, fvGeometry, sol);
210
211
7/8
✓ Branch 0 taken 80 times.
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 405 times.
✓ Branch 3 taken 365 times.
✓ Branch 4 taken 325 times.
✓ Branch 5 taken 325 times.
✓ Branch 7 taken 325 times.
✗ Branch 8 not taken.
 1540 for (auto&& scv : scvs(fvGeometry))
212 {
213
1/2
✓ Branch 1 taken 405 times.
✗ Branch 2 not taken.
 405 auto pointSources = this->scvPointSources(element, fvGeometry, elemVolVars, scv);
214
3/5
✓ Branch 1 taken 80 times.
✓ Branch 2 taken 325 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 80 times.
✗ Branch 5 not taken.
 485 pointSources *= scv.volume()*elemVolVars[scv].extrusionFactor();
215
3/6
✓ Branch 1 taken 405 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 405 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 325 times.
✗ Branch 8 not taken.
 1135 source[scv.dofIndex()] += pointSources;
216
3/5
✓ Branch 1 taken 405 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 325 times.
✓ Branch 5 taken 80 times.
✗ Branch 6 not taken.
 730 auto a = fvGeometry.geometry(scv).corner(0)[2];
217
3/5
✓ Branch 1 taken 405 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 365 times.
✓ Branch 5 taken 40 times.
 730 auto b = fvGeometry.geometry(scv).corner(1)[2];
218
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 365 times.
 405 if (a > b) std::swap(a, b);
219 1135 sourceExact[scv.dofIndex()] += a*(1.0+0.5*a) - b*(1.0+0.5*b);
220 730 volume += scv.volume();
221 }
222 }
223
2/4
✓ Branch 5 taken 19 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 19 times.
✗ Branch 9 not taken.
 76 std::cout << "Global integrated source (1D): " << std::accumulate(source.begin(), source.end(), 0.0) << '\n';
224
2/6
✓ Branch 5 taken 19 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 19 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 76 std::cout << "Global integrated source exact (1D): " << std::accumulate(sourceExact.begin(), sourceExact.end(), 0.0) << '\n';
225
226 19 source -= sourceExact;
227 19 const auto norm = source.two_norm()/sourceExact.two_norm();
228
2/4
✓ Branch 2 taken 19 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 19 times.
✗ Branch 6 not taken.
 38 std::cout << "relative L2 Norm source (1D): " << norm << '\n';
229
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 38 return source.two_norm()/volume;
230 }
231
232 /*!
233 * \brief Adds additional VTK output data to the VTKWriter.
234 *
235 * Function is called by the output module on every write.
236 */
237 template<class VtkOutputModule>
238 19 void addVtkOutputFields(VtkOutputModule& vtk) const
239 {
240
4/10
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 19 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 19 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 38 vtk.addField(exactPressure_, "exact pressure");
241 19 }
242
243 //! Get the coupling manager
244 const CouplingManager& couplingManager() const
245
9/17
✗ Branch 0 not taken.
✓ Branch 1 taken 5940 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5940 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2097520 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2097520 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 5953 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 5940 times.
✓ Branch 12 taken 13 times.
✓ Branch 13 taken 6 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 6 times.
✗ Branch 17 not taken.
 8344302 { return *couplingManager_; }
246
247 private:
248 std::vector<Scalar> exactPressure_;
249
250 static constexpr Scalar eps_ = 1.5e-7;
251 std::string name_;
252
253 Scalar p_in_;
254 Scalar delta_p_;
255
256 std::shared_ptr<CouplingManager> couplingManager_;
257 };
258
259 } // end namespace Dumux
260
261 #endif
262