GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/test/multidomain/boundary/freeflowporenetwork/1p2c_1p2c/problem_freeflow.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 68 68 100.0%
Functions: 8 8 100.0%
Branches: 66 88 75.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 BoundaryTests
10 * \brief Free-flow sub-problem for the coupled 1p2c_1p2c free-flow/pore-network-model test
11 */
12
13 #ifndef DUMUX_TEST_MULTIDOMAIN_BOUNDARY_FREEFLOW_PORE_NETWORK_ONEPTWOC_PROBLEM_FREEFLOW_HH
14 #define DUMUX_TEST_MULTIDOMAIN_BOUNDARY_FREEFLOW_PORE_NETWORK_ONEPTWOC_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/1pnc/advectiveflux.hh>
20
21 namespace Dumux {
22
23 /*!
24 * \ingroup BoundaryTests
25 * \brief Free-flow sub-problem for the coupled 1p2c_1p2c 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 FreeFlowOnePNCTestProblem : 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 static constexpr bool useMoles = ModelTraits::useMoles();
55
56 public:
57 12 FreeFlowOnePNCTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
58 : ParentType(gridGeometry, couplingManager, "FreeFlow")
59
2/6
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 48 , couplingManager_(couplingManager)
60 {
61
3/6
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
 24 problemName_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
62
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 initialPressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InitialPressure", 1e5);
63
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 initialMoleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InitialMoleFraction", 0.001);
64
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 onlyDiffusion_ = getParam<bool>("Problem.OnlyDiffusion", false);
65
3/4
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 12 pressureDifference_ = onlyDiffusion_ ? 0.0 : getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference", 10);
66
67 #if !ISOTHERMAL
68
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 initialTemperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InitialTemperature", 273.15 + 20.0);
69 #endif
70 12 }
71
72 /*!
73 * \name Problem parameters
74 */
75 // \{
76
77 3 const std::string& name() const
78
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 { return problemName_; }
79
80 // \}
81
82 /*!
83 * \name Boundary conditions
84 */
85 // \{
86
87 /*!
88 * \brief Specifies which kind of boundary condition should be
89 * used for which equation on a given boundary segment.
90 *
91 * \param element The finite element
92 * \param scvf The sub control volume face
93 */
94 459244 BoundaryTypes boundaryTypes(const Element& element,
95 const SubControlVolumeFace& scvf) const
96 {
97 459244 BoundaryTypes values;
98 459244 const auto& globalPos = scvf.center(); //avoid ambiguities at corners
99
100 if constexpr (ParentType::isMomentumProblem())
101 {
102
2/2
✓ Branch 1 taken 2031 times.
✓ Branch 2 taken 56747 times.
 117556 if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::poreNetworkIndex, scvf))
103 {
104 4062 values.setCouplingNeumann(Indices::momentumYBalanceIdx);
105 4062 values.setCouplingNeumann(Indices::momentumXBalanceIdx);
106 }
107
6/6
✓ Branch 0 taken 43207 times.
✓ Branch 1 taken 13540 times.
✓ Branch 2 taken 13540 times.
✓ Branch 3 taken 29667 times.
✓ Branch 4 taken 29667 times.
✓ Branch 5 taken 27080 times.
 167654 else if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
108 117556 values.setAllNeumann();
109 else
110 117556 values.setAllDirichlet();
111 }
112 else
113 {
114 356190 if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::poreNetworkIndex, scvf))
115 341688 values.setAllCouplingNeumann();
116
6/6
✓ Branch 0 taken 115453 times.
✓ Branch 1 taken 48140 times.
✓ Branch 2 taken 87633 times.
✓ Branch 3 taken 27820 times.
✓ Branch 4 taken 135773 times.
✓ Branch 5 taken 27820 times.
 598732 else if (!onlyDiffusion_ && onLeftBoundary_(globalPos))
117 341688 values.setAllDirichlet();
118 else
119 341688 values.setAllNeumann();
120 }
121 459244 return values;
122 }
123
124 /*!
125 * \brief Returns Dirichlet boundary values at a given position.
126 *
127 * \param globalPos The global position
128 */
129 6420 DirichletValues dirichletAtPos(const GlobalPosition& globalPos) const
130 {
131 18229 return initialAtPos(globalPos);
132 }
133
134 /*!
135 * \brief Evaluates the boundary conditions for a Neumann control volume.
136 *
137 * \param element The element for which the Neumann boundary condition is set
138 * \param fvGeometry The fvGeometry
139 * \param elemVolVars The element volume variables
140 * \param elemFaceVars The element face variables
141 * \param scvf The boundary sub control volume face
142 */
143 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
144 1227846 BoundaryFluxes neumann(const Element& element,
145 const FVElementGeometry& fvGeometry,
146 const ElementVolumeVariables& elemVolVars,
147 const ElementFluxVariablesCache& elemFluxVarsCache,
148 const SubControlVolumeFace& scvf) const
149 {
150 1227846 BoundaryFluxes values(0.0);
151 1227846 const auto& globalPos = scvf.ipGlobal();
152 using SlipVelocityPolicy = NavierStokesSlipVelocity<typename GridGeometry::DiscretizationMethod, NavierStokes::SlipConditions::BJ>;
153 using FluxHelper = NavierStokesMomentumBoundaryFlux<typename GridGeometry::DiscretizationMethod, SlipVelocityPolicy>;
154
155 if constexpr (ParentType::isMomentumProblem())
156 {
157
158
2/2
✓ Branch 1 taken 90137 times.
✓ Branch 2 taken 413988 times.
 1008250 if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::poreNetworkIndex, scvf))
159 {
160 180274 values[scvf.normalAxis()] += couplingManager_->momentumCouplingCondition(
161 CouplingManager::freeFlowMomentumIndex, CouplingManager::poreNetworkIndex,
162 fvGeometry, scvf, elemVolVars
163 );
164
165 360548 values += FluxHelper::slipVelocityMomentumFlux(
166 *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache
167 );
168 }
169
2/2
✓ Branch 0 taken 206994 times.
✓ Branch 1 taken 206994 times.
 827976 else if (onLeftBoundary_(globalPos))
170 {
171 413988 values = FluxHelper::fixedPressureMomentumFlux(
172 *this, fvGeometry, scvf, elemVolVars,
173 413988 elemFluxVarsCache, initialPressure_ + pressureDifference_, true /*zeroNormalVelocityGradient*/
174 );
175 }
176
1/2
✓ Branch 0 taken 206994 times.
✗ Branch 1 not taken.
 413988 else if (onRightBoundary_(globalPos))
177 {
178 413988 values = FluxHelper::fixedPressureMomentumFlux(
179 *this, fvGeometry, scvf, elemVolVars,
180 413988 elemFluxVarsCache, initialPressure_, true /*zeroNormalVelocityGradient*/
181 );
182 }
183 }
184 else
185 {
186 219596 if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::poreNetworkIndex, scvf))
187 {
188 11736 const auto& massCouplingCondition = couplingManager_->massCouplingCondition(
189 CouplingManager::freeFlowMassIndex, CouplingManager::poreNetworkIndex,
190 fvGeometry, scvf, elemVolVars
191 );
192 11736 values[Indices::conti0EqIdx] = massCouplingCondition[Indices::conti0EqIdx];
193 11736 values[Indices::conti0EqIdx + 1] = massCouplingCondition[Indices::conti0EqIdx + 1];
194
195 #if !ISOTHERMAL
196 4512 values[Indices::energyEqIdx] = couplingManager_->energyCouplingCondition(CouplingManager::freeFlowMassIndex, CouplingManager::poreNetworkIndex,
197 fvGeometry, scvf, elemVolVars);
198 #endif
199 }
200
4/4
✓ Branch 0 taken 67410 times.
✓ Branch 1 taken 36520 times.
✓ Branch 2 taken 21400 times.
✓ Branch 3 taken 46010 times.
 207860 else if (!onlyDiffusion_ && onRightBoundary_(globalPos))
201 {
202 using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
203 42800 DirichletValues outsideBoundaryPriVars = initialAtPos(globalPos);
204 42800 values = FluxHelper::scalarOutflowFlux(
205 *this, element, fvGeometry, scvf, elemVolVars, std::move(outsideBoundaryPriVars)
206 );
207 }
208 }
209
210 1227846 return values;
211 }
212
213 /*!
214 * \brief Evaluates the source term for all phases within a given
215 * sub-control volume
216 */
217 template<class ElementVolumeVariables>
218 260004 Sources source(const Element& element,
219 const FVElementGeometry& fvGeometry,
220 const ElementVolumeVariables& elemVolVars,
221 const SubControlVolume& scv) const
222 {
223 auto source = Sources(0.0);
224
225 260004 return source;
226 }
227
228 // The following function defines the initial conditions
229 29200 InitialValues initialAtPos(const GlobalPosition &globalPos) const
230 {
231 14220 InitialValues values(0.0); //velocity is 0.0
232
233 if constexpr (!ParentType::isMomentumProblem())
234 {
235
1/2
✓ Branch 2 taken 6420 times.
✗ Branch 3 not taken.
 29200 values[Indices::pressureIdx] = initialPressure_;
236
1/2
✓ Branch 2 taken 3180 times.
✗ Branch 3 not taken.
 29200 values[Indices::conti0EqIdx +1] = initialMoleFraction_;
237 #if !ISOTHERMAL
238
1/2
✓ Branch 2 taken 3240 times.
✗ Branch 3 not taken.
 14980 values[Indices::temperatureIdx] = initialTemperature_;
239 #endif
240 }
241
242 return values;
243 }
244
245 /*!
246 * \brief Returns true if the scvf lies on a porous slip boundary
247 */
248
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 248622 times.
 248622 bool onSlipBoundary(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
249 {
250
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 248622 times.
 248622 assert(scvf.isFrontal());
251
3/4
✓ Branch 0 taken 248622 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 176888 times.
✓ Branch 4 taken 71734 times.
 248622 return scvf.boundary() && couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::poreNetworkIndex, scvf);
252 }
253
254 /*!
255 * \brief Returns the beta value
256 */
257 71734 Scalar betaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
258 {
259 71734 const Scalar radius = couplingManager_->coupledPoreInscribedRadius(fvGeometry, scvf);
260 71734 return 5.73 / radius; // this value only an approximation of wall friction is considered
261 }
262
263 /*!
264 * \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann).
265 */
266 71734 VelocityVector porousMediumVelocity(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
267 {
268 71734 return couplingManager_->interfaceThroatVelocity(fvGeometry, scvf);
269 }
270
271 // \}
272
273 private:
274
275
6/6
✓ Branch 0 taken 206994 times.
✓ Branch 1 taken 206994 times.
✓ Branch 2 taken 87633 times.
✓ Branch 3 taken 27820 times.
✓ Branch 4 taken 43207 times.
✓ Branch 5 taken 13540 times.
 586188 bool onLeftBoundary_(const GlobalPosition &globalPos) const
276
6/6
✓ Branch 0 taken 206994 times.
✓ Branch 1 taken 206994 times.
✓ Branch 2 taken 87633 times.
✓ Branch 3 taken 27820 times.
✓ Branch 4 taken 43207 times.
✓ Branch 5 taken 13540 times.
 586188 { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
277
278
5/6
✓ Branch 0 taken 206994 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 21400 times.
✓ Branch 3 taken 46010 times.
✓ Branch 4 taken 13540 times.
✓ Branch 5 taken 29667 times.
 317611 bool onRightBoundary_(const GlobalPosition &globalPos) const
279
5/6
✓ Branch 0 taken 206994 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 21400 times.
✓ Branch 3 taken 46010 times.
✓ Branch 4 taken 13540 times.
✓ Branch 5 taken 29667 times.
 317611 { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
280
281 bool onUpperBoundary_(const GlobalPosition &globalPos) const
282 { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
283
284 std::string problemName_;
285 static constexpr Scalar eps_ = 1e-6;
286 Scalar deltaP_;
287 Scalar initialPressure_;
288 Scalar pressureDifference_;
289 Scalar initialMoleFraction_;
290 bool onlyDiffusion_;
291 #if !ISOTHERMAL
292 Scalar initialTemperature_;
293 #endif
294 std::shared_ptr<CouplingManager> couplingManager_;
295 };
296 } // end namespace Dumux
297
298 #endif
299