GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: dumux/test/multidomain/embedded/1d3d/1p_1p/problem_tissue.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 97 97 100.0%
Functions: 21 21 100.0%
Branches: 80 114 70.2%

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-FileCopyrightText: 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 Definition of a problem, for the 1p2c problem:
11 * Component transport of oxygen in interstitial fluid.
12 */
13
14 #ifndef DUMUX_TISSUE_PROBLEM_HH
15 #define DUMUX_TISSUE_PROBLEM_HH
16
17 #include <dune/geometry/quadraturerules.hh>
18 #include <dune/localfunctions/lagrange/pqkfactory.hh>
19
20 #include <dumux/common/boundarytypes.hh>
21 #include <dumux/common/math.hh>
22 #include <dumux/common/parameters.hh>
23 #include <dumux/common/properties.hh>
24 #include <dumux/common/numeqvector.hh>
25
26 #include <dumux/porousmediumflow/problem.hh>
27
28 #include <dumux/multidomain/embedded/couplingmanager1d3d.hh> // for coupling mode
29
30 namespace Dumux {
31
32 /*!
33 * \ingroup EmbeddedTests
34 * \brief Definition of a problem, for the 1p2c problem:
35 * Component transport of oxygen in interstitial fluid.
36 */
37 template <class TypeTag>
38 class TissueProblem : public PorousMediumFlowProblem<TypeTag>
39 {
40 using ParentType = PorousMediumFlowProblem<TypeTag>;
41 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
42 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
43 using FVElementGeometry = typename GridGeometry::LocalView;
44 using SubControlVolume = typename GridGeometry::SubControlVolume;
45 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
46 using GridView = typename GridGeometry::GridView;
47 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
48 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
49 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
50 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
51 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
52 using PointSource = GetPropType<TypeTag, Properties::PointSource>;
53 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
54 using Element = typename GridView::template Codim<0>::Entity;
55 using GlobalPosition = typename GridGeometry::GlobalCoordinate;
56
57 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
58
59 public:
60 19 TissueProblem(std::shared_ptr<const GridGeometry> gridGeometry,
61 std::shared_ptr<CouplingManager> couplingManager)
62 : ParentType(gridGeometry, "Tissue")
63
2/4
✓ Branch 2 taken 19 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
57 , couplingManager_(couplingManager)
64 {
65 // read parameters from input file
66
3/6
✓ 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.
38 name_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
67
68
2/4
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 17 times.
✗ Branch 5 not taken.
19 exactPressure_.resize(this->gridGeometry().numDofs());
69
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
19 auto fvGeometry = localView(this->gridGeometry());
70
2/4
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 931237 times.
✗ Branch 5 not taken.
2793673 for (const auto& element : elements(this->gridGeometry().gridView()))
71 {
72 931218 fvGeometry.bindElement(element);
73
74
3/4
✓ Branch 1 taken 1027244 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1027244 times.
✓ Branch 4 taken 931218 times.
1958462 for (auto&& scv : scvs(fvGeometry))
75
1/2
✓ Branch 1 taken 1027244 times.
✗ Branch 2 not taken.
1027244 exactPressure_[scv.dofIndex()] = exactSolution(scv.dofPosition());
76 }
77 19 }
78
79 /*!
80 * \brief The problem name.
81 *
82 * This is used as a prefix for files generated by the simulation.
83 */
84 19 const std::string& name() const
85
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
19 { return name_; }
86
87 // \}
88
89 /*!
90 * \name Boundary conditions
91 */
92 // \{
93
94 /*!
95 * \brief Specifies which kind of boundary condition should be
96 * used for which equation on a given boundary segment.
97 *
98 * \param globalPos The position for which the bc type should be evaluated
99 */
100
2/2
✓ Branch 0 taken 359544 times.
✓ Branch 1 taken 112888 times.
472432 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
101 {
102
2/2
✓ Branch 0 taken 359544 times.
✓ Branch 1 taken 112888 times.
472432 BoundaryTypes values;
103
6/6
✓ Branch 0 taken 359544 times.
✓ Branch 1 taken 112888 times.
✓ Branch 2 taken 112888 times.
✓ Branch 3 taken 246656 times.
✓ Branch 4 taken 225776 times.
✓ Branch 5 taken 246656 times.
698208 if (globalPos[2] > this->gridGeometry().bBoxMax()[2] - eps_ || globalPos[2] < this->gridGeometry().bBoxMin()[2] + eps_)
104 values.setAllNeumann();
105 else
106 values.setAllDirichlet();
107 472432 return values;
108 }
109
110 /*!
111 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
112 *
113 * \param globalPos The position for which the bc type should be evaluated
114 *
115 * For this method, the \a values parameter stores primary variables.
116 */
117 128656 PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
118 {
119 128656 PrimaryVariables values;
120
2/4
✓ Branch 1 taken 128656 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 118000 times.
✗ Branch 5 not taken.
128656 values = exactSolution(globalPos);
121 return values;
122 }
123
124 /*!
125 * \brief Evaluate the boundary conditions for a neumann
126 * boundary segment.
127 *
128 * This is the method for the case where the Neumann condition is
129 * potentially solution dependent
130 *
131 * \param element The finite element
132 * \param fvGeometry The finite-volume geometry
133 * \param elemVolVars All volume variables for the element
134 * \param scvf The sub control volume face
135 *
136 * Negative values mean influx.
137 * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
138 */
139 template<class ElementVolumeVariables, class ElemFluxVarsCache>
140 167456 NumEqVector neumann(const Element& element,
141 const FVElementGeometry& fvGeometry,
142 const ElementVolumeVariables& elemVolVars,
143 const ElemFluxVarsCache&,
144 const SubControlVolumeFace& scvf) const
145 {
146 167456 NumEqVector flux(0.0);
147 // integrate over the scvf to compute the flux
148 167456 const auto geometry = fvGeometry.geometry(scvf);
149 167456 Scalar derivative = 0.0;
150 167456 const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension-1>::rule(geometry.type(), 4);
151
2/2
✓ Branch 1 taken 1507104 times.
✓ Branch 2 taken 167456 times.
1674560 for(auto&& qp : quad)
152 {
153 1507104 auto globalPos = geometry.global(qp.position());
154 1507104 globalPos[2] = 0; // the derivative in z-direction is the exact solution evaluated with z=0
155 2497104 derivative += exactSolution(globalPos)*qp.weight()*geometry.integrationElement(qp.position());
156 }
157
158
2/2
✓ Branch 0 taken 80992 times.
✓ Branch 1 taken 86464 times.
167456 const auto globalPos = scvf.ipGlobal();
159
2/2
✓ Branch 0 taken 80992 times.
✓ Branch 1 taken 86464 times.
167456 if (globalPos[2] > this->gridGeometry().bBoxMax()[2] - eps_ )
160 80992 flux[0] = -derivative/scvf.area();
161 else
162 86464 flux[0] = derivative/scvf.area();
163 167456 return flux;
164 }
165
166
167 // \}
168
169 /*!
170 * \name Volume terms
171 */
172 // \{
173
174 /*!
175 * \brief Applies a vector of point sources which are possibly solution dependent.
176 *
177 * \param pointSources A vector of Dumux::PointSource s that contain
178 source values for all phases and space positions.
179 *
180 * For this method, the \a values method of the point source
181 * has to return the absolute mass rate in kg/s. Positive values mean
182 * that mass is created, negative ones mean that it vanishes.
183 */
184 19 void addPointSources(std::vector<PointSource>& pointSources) const
185
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
19 { pointSources = this->couplingManager().bulkPointSources(); }
186
187 /*!
188 * \brief Evaluates the point sources (added by addPointSources)
189 * for all phases within a given sub control volume.
190 *
191 * This is the method for the case where the point source is
192 * solution dependent and requires some quantities that
193 * are specific to the fully-implicit method.
194 *
195 * \param source A single point source
196 * \param element The finite element
197 * \param fvGeometry The finite-volume geometry
198 * \param elemVolVars All volume variables for the element
199 * \param scv The sub-control volume within the element
200 *
201 * For this method, the \a values() method of the point sources returns
202 * the absolute rate mass generated or annihilated in kg/s. Positive values mean
203 * that mass is created, negative ones mean that it vanishes.
204 */
205 template<class ElementVolumeVariables>
206 1267416 void pointSource(PointSource& source,
207 const Element &element,
208 const FVElementGeometry& fvGeometry,
209 const ElementVolumeVariables& elemVolVars,
210 const SubControlVolume &scv) const
211 {
212 if constexpr (CouplingManager::couplingMode != Embedded1d3dCouplingMode::kernel)
213 {
214 // compute source at every integration point
215
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1199320 times.
1267416 const Scalar pressure3D = this->couplingManager().bulkPriVars(source.id())[Indices::pressureIdx];
216
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1267416 times.
1267416 const Scalar pressure1D = this->couplingManager().lowDimPriVars(source.id())[Indices::pressureIdx];
217
218 // calculate the source
219
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1267416 times.
1267416 const Scalar radius = this->couplingManager().radius(source.id());
220 1267416 const Scalar beta = 2*M_PI/(2*M_PI + std::log(radius));
221 1267416 const Scalar sourceValue = beta*(pressure1D - pressure3D);
222 1267416 source = sourceValue*source.quadratureWeight()*source.integrationElement();
223 }
224 1267416 }
225
226 /*!
227 * \brief Evaluates the source term for all phases within a given
228 * sub control volume.
229 *
230 * This is the method for the case where the source term is
231 * potentially solution dependent and requires some quantities that
232 * are specific to the fully-implicit method.
233 *
234 * \param element The finite element
235 * \param fvGeometry The finite-volume geometry
236 * \param elemVolVars All volume variables for the element
237 * \param scv The sub control volume
238 *
239 * For this method, the return parameter stores the conserved quantity rate
240 * generated or annihilated per volume unit. Positive values mean
241 * that the conserved quantity is created, negative ones mean that it vanishes.
242 * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
243 */
244 template<class ElementVolumeVariables>
245 925620 NumEqVector source(const Element &element,
246 const FVElementGeometry& fvGeometry,
247 const ElementVolumeVariables& elemVolVars,
248 const SubControlVolume &scv) const
249 {
250 925620 NumEqVector source(0.0);
251
252 if constexpr(CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel)
253 {
254 925620 const auto eIdx = this->gridGeometry().elementMapper().index(element);
255 925620 const auto& sourceIds = this->couplingManager().bulkSourceIds(eIdx, scv.indexInElement());
256 925620 const auto& sourceWeights = this->couplingManager().bulkSourceWeights(eIdx, scv.indexInElement());
257
258
2/2
✓ Branch 0 taken 138480 times.
✓ Branch 1 taken 925620 times.
1064100 for (int i = 0; i < sourceIds.size(); ++i)
259 {
260
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 138480 times.
138480 const auto id = sourceIds[i];
261 138480 const auto weight = sourceWeights[i];
262
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 138480 times.
138480 const auto xi = this->couplingManager().fluxScalingFactor(id);
263
264
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 138480 times.
138480 const Scalar radius = this->couplingManager().radius(id);
265 138480 const Scalar pressure3D = this->couplingManager().bulkPriVars(id)[Indices::pressureIdx];
266
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 138475 times.
138480 const Scalar pressure1D = this->couplingManager().lowDimPriVars(id)[Indices::pressureIdx];
267
268 // calculate the source
269
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 138475 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
138480 static const Scalar beta = 2*M_PI/(2*M_PI + std::log(radius));
270 138480 const Scalar sourceValue = beta*(pressure1D - pressure3D)*xi;
271 138480 source[Indices::conti0EqIdx] += sourceValue*weight;
272 }
273
274 925620 const auto volume = scv.volume()*elemVolVars[scv].extrusionFactor();
275 925620 source[Indices::conti0EqIdx] /= volume;
276 }
277
278 925620 return source;
279 }
280
281 //! compute the flux scaling factor (xi) for a distance r for a vessel with radius R and kernel width rho
282 380 Scalar fluxScalingFactor(const Scalar r, const Scalar R, const Scalar rho) const
283 {
284 using std::log;
285
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 375 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
380 static const Scalar beta = 2*M_PI/(2*M_PI + log(R));
286
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 375 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
380 static const Scalar kernelWidthFactorLog = log(rho/R);
287
2/2
✓ Branch 0 taken 360 times.
✓ Branch 1 taken 20 times.
380 return r < rho ? 1.0/(1.0 + R*beta/(2*M_PI*R)*(r*r/(2*rho*rho) + kernelWidthFactorLog - 0.5))
288 20 : 1.0/(1.0 + R*beta/(2*M_PI*R)*log(r/R));
289 }
290
291 /*!
292 * \brief Evaluates the initial value for a control volume.
293 *
294 * \param globalPos The position for which the initial condition should be evaluated
295 *
296 * For this method, the \a values parameter stores primary
297 * variables.
298 */
299 933500 PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
300 { return PrimaryVariables(0.0); }
301
302 3594222 Scalar p1DExact(const GlobalPosition &globalPos) const
303 3594222 { return 1.0 + globalPos[2]; }
304
305 //! The exact solution
306
2/2
✓ Branch 0 taken 19 times.
✓ Branch 1 taken 3594203 times.
3594222 Scalar exactSolution(const GlobalPosition &globalPos) const
307 {
308 3594222 const auto r = std::hypot(globalPos[0], globalPos[1]);
309
4/6
✓ Branch 0 taken 19 times.
✓ Branch 1 taken 3594203 times.
✓ Branch 3 taken 19 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 19 times.
✗ Branch 7 not taken.
3594222 static const auto R = getParam<Scalar>("SpatialParams.Radius");
310
311 if (CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel)
312 {
313
4/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 958595 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
958600 static const auto rho = getParam<Scalar>("MixedDimension.KernelWidthFactor")*R;
314
2/2
✓ Branch 0 taken 949208 times.
✓ Branch 1 taken 9392 times.
958600 if (r > rho)
315 949208 return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(r);
316 else
317 9392 return -1.0*p1DExact(globalPos)/(2*M_PI)*(r*r/(2.0*rho*rho) + std::log(rho) - 0.5);
318 }
319 else if (CouplingManager::couplingMode == Embedded1d3dCouplingMode::surface)
320 {
321
2/2
✓ Branch 0 taken 949208 times.
✓ Branch 1 taken 9392 times.
958600 if (r > R)
322 949208 return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(r);
323 else
324 9392 return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(R);
325 }
326
327 1677022 return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(r);
328 }
329
330 //! Called after every time step
331 //! Output the total global exchange term
332 19 void computeSourceIntegral(const SolutionVector& sol, const GridVariables& gridVars)
333 {
334 19 PrimaryVariables source(0.0);
335
2/3
✓ Branch 2 taken 13720 times.
✓ Branch 3 taken 917517 times.
✗ Branch 4 not taken.
2793692 for (const auto& element : elements(this->gridGeometry().gridView()))
336 {
337 1862436 const auto fvGeometry = localView(this->gridGeometry()).bindElement(element);
338 931218 const auto elemVolVars = localView(gridVars.curGridVolVars()).bindElement(element, fvGeometry, sol);
339
340
2/2
✓ Branch 0 taken 1027244 times.
✓ Branch 1 taken 931218 times.
1958462 for (auto&& scv : scvs(fvGeometry))
341 {
342
1/2
✓ Branch 1 taken 917500 times.
✗ Branch 2 not taken.
1027244 auto pointSources = this->scvPointSources(element, fvGeometry, elemVolVars, scv);
343 1027244 pointSources += this->source(element, fvGeometry, elemVolVars, scv);
344 2054488 pointSources *= scv.volume()*elemVolVars[scv].extrusionFactor();
345 1027244 source += pointSources;
346 }
347 }
348
349 19 std::cout << "Global integrated source (3D): " << source << '\n';
350 19 }
351
352 /*!
353 * \brief Adds additional VTK output data to the VTKWriter.
354 *
355 * Function is called by the output module on every write.
356 */
357 template<class VtkOutputModule>
358 19 void addVtkOutputFields(VtkOutputModule& vtk) const
359 {
360
1/2
✓ Branch 2 taken 19 times.
✗ Branch 3 not taken.
19 vtk.addField(exactPressure_, "exact pressure");
361 19 }
362
363 //! Get the coupling manager
364 2331535 const CouplingManager& couplingManager() const
365
5/7
✗ Branch 0 not taken.
✓ Branch 1 taken 1337800 times.
✓ Branch 2 taken 68101 times.
✓ Branch 3 taken 138487 times.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
✓ Branch 4 taken 2 times.
2331535 { return *couplingManager_; }
366
367 private:
368 std::vector<Scalar> exactPressure_;
369
370 static constexpr Scalar eps_ = 1.5e-7;
371 std::string name_;
372
373 std::shared_ptr<CouplingManager> couplingManager_;
374 };
375
376 } // end namespace Dumux
377
378 #endif
379