GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/multidomain/boundary/freeflowporenetwork/1p_1p/problem_freeflow.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 49 54 90.7%
Functions: 7 12 58.3%
Branches: 83 130 63.8%

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 BoundaryTests
10 * \brief Free-flow sub-problem for the coupled 1p_1p free-flow/pore-network-model test
11 */
12
13 #ifndef DUMUX_TEST_MULTIDOMAIN_BOUNDARY_FREEFLOW_PORE_NETWORK_PROBLEM_FREEFLOW_HH
14 #define DUMUX_TEST_MULTIDOMAIN_BOUNDARY_FREEFLOW_PORE_NETWORK_PROBLEM_FREEFLOW_HH
15
16 #include <dumux/common/properties.hh>
17 #include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
18 #include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
19 #include <dumux/freeflow/navierstokes/mass/1p/advectiveflux.hh>
20
21 namespace Dumux {
22
23 /*!
24 * \ingroup BoundaryTests
25 * \brief Free-flow sub-problem for the coupled 1p_1p free-flow/pore-network-model test
26 * A two-dimensional Stokes flow region coupled to a pore-network model.
27 */
28 template <class TypeTag, class BaseProblem>
29 class FreeFlowOnePTestProblem : public BaseProblem
30 {
31 using ParentType = BaseProblem;
32
33 using BoundaryTypes = typename ParentType::BoundaryTypes;
34 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
35 using FVElementGeometry = typename GridGeometry::LocalView;
36 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
37 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
38 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
39 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
40 using InitialValues = typename ParentType::InitialValues;
41 using Sources = typename ParentType::Sources;
42 using DirichletValues = typename ParentType::DirichletValues;
43 using BoundaryFluxes = typename ParentType::BoundaryFluxes;
44 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
45 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
46
47 static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
48 using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
49 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
50 using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
51
52 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
53
54 public:
55 8 FreeFlowOnePTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
56 : ParentType(gridGeometry, couplingManager, "FreeFlow")
57
6/20
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 4 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 4 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 4 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.
32 , couplingManager_(couplingManager)
58 {
59
9/24
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4 times.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 4 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 4 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 4 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 4 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.
16 problemName_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
60
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
8 deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference", 0.0);
61
1/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
8 singleThroatTest_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.SingleThroatTest", true);
62 8 enablePseudoThreeDWallFriction_ = !singleThroatTest_;
63 8 }
64
65 /*!
66 * \name Problem parameters
67 */
68 // \{
69
70 const std::string& name() const
71
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
2 { return problemName_; }
72
73 // \}
74
75 /*!
76 * \name Boundary conditions
77 */
78 // \{
79
80 /*!
81 * \brief Specifies which kind of boundary condition should be
82 * used for which equation on a given boundary segment.
83 *
84 * \param element The finite element
85 * \param scvf The sub control volume face
86 */
87 14080 BoundaryTypes boundaryTypes(const Element& element,
88 const SubControlVolumeFace& scvf) const
89 {
90
2/2
✓ Branch 0 taken 1995 times.
✓ Branch 1 taken 2697 times.
14080 BoundaryTypes values;
91
2/2
✓ Branch 0 taken 2955 times.
✓ Branch 1 taken 4085 times.
14080 const auto& globalPos = scvf.center(); //avoid ambiguities at corners
92
93
2/2
✓ Branch 0 taken 2955 times.
✓ Branch 1 taken 4085 times.
14080 if (singleThroatTest_) // vertical flow
94 {
95 if constexpr (ParentType::isMomentumProblem())
96 {
97
2/2
✓ Branch 2 taken 95 times.
✓ Branch 3 taken 865 times.
3840 if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::poreNetworkIndex, scvf))
98 {
99 190 values.setCouplingNeumann(Indices::momentumYBalanceIdx);
100 190 values.setCouplingNeumann(Indices::momentumXBalanceIdx);
101 }
102 3460 else if (onUpperBoundary_(globalPos))
103 values.setAllNeumann();
104 else
105 values.setAllDirichlet();
106 }
107 else
108 {
109
2/2
✓ Branch 2 taken 225 times.
✓ Branch 3 taken 1770 times.
7980 if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::poreNetworkIndex, scvf))
110 values.setAllCouplingNeumann();
111 else
112 values.setAllNeumann();
113 }
114 }
115 else // horizontal flow
116 {
117 if constexpr (ParentType::isMomentumProblem())
118 {
119 10768 if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
120 values.setAllNeumann();
121
2/2
✓ Branch 2 taken 126 times.
✓ Branch 3 taken 1094 times.
4880 else if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::poreNetworkIndex, scvf))
122 {
123 252 values.setCouplingNeumann(Indices::momentumXBalanceIdx);
124 252 values.setCouplingNeumann(Indices::momentumYBalanceIdx);
125 }
126 else
127 values.setAllDirichlet();
128 }
129 else
130 {
131
2/2
✓ Branch 2 taken 288 times.
✓ Branch 3 taken 2409 times.
10788 if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::poreNetworkIndex, scvf))
132 values.setAllCouplingNeumann();
133 else
134 values.setAllNeumann();
135 }
136 }
137
138 14080 return values;
139 }
140
141 /*!
142 * \brief Returns Dirichlet boundary values at a given position.
143 *
144 * \param globalPos The global position
145 */
146 DirichletValues dirichletAtPos(const GlobalPosition& globalPos) const
147 {
148 2372 return DirichletValues(0.0);
149 }
150
151 /*!
152 * \brief Evaluates the boundary conditions for a Neumann control volume.
153 *
154 * \param element The element for which the Neumann boundary condition is set
155 * \param fvGeometry The fvGeometry
156 * \param elemVolVars The element volume variables
157 * \param elemFaceVars The element face variables
158 * \param scvf The boundary sub control volume face
159 */
160 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
161 22216 BoundaryFluxes neumann(const Element& element,
162 const FVElementGeometry& fvGeometry,
163 const ElementVolumeVariables& elemVolVars,
164 const ElementFluxVariablesCache& elemFluxVarsCache,
165 const SubControlVolumeFace& scvf) const
166 {
167 22216 BoundaryFluxes values(0.0);
168 22216 const auto& globalPos = scvf.ipGlobal();
169 using SlipVelocityPolicy = NavierStokesSlipVelocity<typename GridGeometry::DiscretizationMethod, NavierStokes::SlipConditions::BJ>;
170 using FluxHelper = NavierStokesMomentumBoundaryFlux<typename GridGeometry::DiscretizationMethod, SlipVelocityPolicy>;
171
172 if constexpr (ParentType::isMomentumProblem())
173 {
174
175
2/2
✓ Branch 2 taken 4491 times.
✓ Branch 3 taken 3061 times.
30208 if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::poreNetworkIndex, scvf))
176 {
177 17964 values += couplingManager_->momentumCouplingCondition(
178 CouplingManager::freeFlowMomentumIndex, CouplingManager::poreNetworkIndex,
179 fvGeometry, scvf, elemVolVars
180 );
181
182 17964 values += FluxHelper::slipVelocityMomentumFlux(
183 *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache
184 );
185 }
186 else
187 {
188 1072 const Scalar pressure = onLeftBoundary_(globalPos) ? deltaP_ : 0.0;
189 6122 values = FluxHelper::fixedPressureMomentumFlux(
190 *this, fvGeometry, scvf, elemVolVars,
191 elemFluxVarsCache, pressure, true /*zeroNormalVelocityGradient*/
192 );
193 }
194 }
195 else
196 {
197
2/2
✓ Branch 2 taken 405 times.
✓ Branch 3 taken 3151 times.
14224 if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::poreNetworkIndex, scvf))
198 {
199 1620 values = couplingManager_->massCouplingCondition(
200 CouplingManager::freeFlowMassIndex, CouplingManager::poreNetworkIndex,
201 fvGeometry, scvf, elemVolVars
202 );
203 }
204 else
205 {
206 using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
207 6302 values = FluxHelper::scalarOutflowFlux(
208 *this, element, fvGeometry, scvf, elemVolVars
209 );
210 }
211 }
212
213 22216 return values;
214 }
215
216 /*!
217 * \brief Evaluates the source term for all phases within a given
218 * sub-control volume
219 */
220 template<class ElementVolumeVariables>
221 87520 Sources source(const Element& element,
222 const FVElementGeometry& fvGeometry,
223 const ElementVolumeVariables& elemVolVars,
224 const SubControlVolume& scv) const
225 {
226
0/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
106972 auto source = Sources(0.0);
227
228 if constexpr (ParentType::isMomentumProblem())
229 {
230
2/2
✓ Branch 0 taken 23360 times.
✓ Branch 1 taken 20400 times.
87520 if (enablePseudoThreeDWallFriction_)
231 {
232 46720 const Scalar height = elemVolVars[scv].extrusionFactor();
233
5/8
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 23359 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
46720 static const Scalar factor = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PseudoWallFractionFactor", 8.0);
234 46720 source[scv.dofAxis()] = this->pseudo3DWallFriction(element, fvGeometry, elemVolVars, scv, height, factor);
235 }
236 }
237
238 87520 return source;
239 }
240
241 /*!
242 * \brief Returns true if the scvf lies on a porous slip boundary
243 */
244 bool onSlipBoundary(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
245 {
246 assert(scvf.isFrontal());
247 return scvf.boundary() && couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::poreNetworkIndex, scvf);
248 }
249
250 /*!
251 * \brief Returns the beta value
252 */
253 Scalar betaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
254 {
255 7212 const Scalar radius = couplingManager_->coupledPoreInscribedRadius(fvGeometry, scvf);
256 3606 return 5.73 / radius; // this values is only an approximation of wall friction is considered
257 }
258
259 /*!
260 * \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann).
261 */
262 VelocityVector porousMediumVelocity(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
263 {
264 7212 return couplingManager_->interfaceThroatVelocity(fvGeometry, scvf);
265 }
266
267 // \}
268
269 private:
270
271 bool onLeftBoundary_(const GlobalPosition &globalPos) const
272
20/20
✓ Branch 0 taken 536 times.
✓ Branch 1 taken 2525 times.
✓ Branch 2 taken 536 times.
✓ Branch 3 taken 2525 times.
✓ Branch 4 taken 536 times.
✓ Branch 5 taken 2525 times.
✓ Branch 6 taken 536 times.
✓ Branch 7 taken 2525 times.
✓ Branch 8 taken 536 times.
✓ Branch 9 taken 2525 times.
✓ Branch 10 taken 1304 times.
✓ Branch 11 taken 84 times.
✓ Branch 12 taken 1304 times.
✓ Branch 13 taken 84 times.
✓ Branch 14 taken 1304 times.
✓ Branch 15 taken 84 times.
✓ Branch 16 taken 1304 times.
✓ Branch 17 taken 84 times.
✓ Branch 18 taken 1304 times.
✓ Branch 19 taken 84 times.
22245 { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
273
274 bool onRightBoundary_(const GlobalPosition &globalPos) const
275
10/10
✓ Branch 0 taken 84 times.
✓ Branch 1 taken 1220 times.
✓ Branch 2 taken 84 times.
✓ Branch 3 taken 1220 times.
✓ Branch 4 taken 84 times.
✓ Branch 5 taken 1220 times.
✓ Branch 6 taken 84 times.
✓ Branch 7 taken 1220 times.
✓ Branch 8 taken 84 times.
✓ Branch 9 taken 1220 times.
6520 { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
276
277 bool onUpperBoundary_(const GlobalPosition &globalPos) const
278
10/10
✓ Branch 0 taken 285 times.
✓ Branch 1 taken 580 times.
✓ Branch 2 taken 285 times.
✓ Branch 3 taken 580 times.
✓ Branch 4 taken 285 times.
✓ Branch 5 taken 580 times.
✓ Branch 6 taken 285 times.
✓ Branch 7 taken 580 times.
✓ Branch 8 taken 285 times.
✓ Branch 9 taken 580 times.
4325 { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
279
280 std::string problemName_;
281 static constexpr Scalar eps_ = 1e-6;
282 Scalar deltaP_;
283 bool singleThroatTest_;
284 bool enablePseudoThreeDWallFriction_;
285
286 std::shared_ptr<CouplingManager> couplingManager_;
287 };
288 } // end namespace Dumux
289
290 #endif
291