GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 61 64 95.3%
Functions: 9 9 100.0%
Branches: 59 79 74.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-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 NavierStokesTests
10 * \brief Test for the staggered grid Navier-Stokes model with analytical solution (Kovasznay 1948, \cite Kovasznay1948)
11 */
12
13 #ifndef DUMUX_KOVASZNAY_TEST_PROBLEM_NEW_HH
14 #define DUMUX_KOVASZNAY_TEST_PROBLEM_NEW_HH
15
16
17 #include <dumux/common/properties.hh>
18 #include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
19 #include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
20 #include <dumux/freeflow/navierstokes/mass/1p/advectiveflux.hh>
21
22 namespace Dumux {
23
24 /*!
25 * \ingroup NavierStokesTests
26 * \brief Test problem for the staggered grid (Kovasznay 1948, \cite Kovasznay1948)
27 *
28 * A two-dimensional Navier-Stokes flow with a periodicity in one direction
29 * is considered. The set-up represents a wake behind a two-dimensional grid
30 * and is chosen in a way such that an exact solution is available.
31 */
32 template <class TypeTag, class BaseProblem>
33 class FreeFlowOnePTestProblem : public BaseProblem
34 {
35 using ParentType = BaseProblem;
36
37 using BoundaryTypes = typename ParentType::BoundaryTypes;
38 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
39 using FVElementGeometry = typename GridGeometry::LocalView;
40 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
41 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
42 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
43 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
44 using InitialValues = typename ParentType::InitialValues;
45 using Sources = typename ParentType::Sources;
46 using DirichletValues = typename ParentType::DirichletValues;
47 using BoundaryFluxes = typename ParentType::BoundaryFluxes;
48 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
49 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
50
51 static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
52 using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
53 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
54 using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
55
56 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
57
58 // static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
59
60 public:
61 12 FreeFlowOnePTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
62 : ParentType(gridGeometry, couplingManager, "FreeFlow")
63
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)
64 {
65
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");
66
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 initialPressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InitialPressure", 0.0);
67
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference", 0.0);
68
69 // determine whether to simulate a vertical or horizontal flow configuration
70 12 verticalFlow_ = problemName_.find("vertical") != std::string::npos;
71 12 }
72
73 /*!
74 * \name Problem parameters
75 */
76 // \{
77
78 3 const std::string& name() const
79
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 { return problemName_; }
80
81 // \}
82
83 /*!
84 * \name Boundary conditions
85 */
86 // \{
87
88 /*!
89 * \brief Specifies which kind of boundary condition should be
90 * used for which equation on a given boundary segment.
91 *
92 * \param element The finite element
93 * \param scvf The sub control volume face
94 */
95
2/2
✓ Branch 0 taken 2168 times.
✓ Branch 1 taken 42000 times.
 123412 BoundaryTypes boundaryTypes(const Element& element,
96 const SubControlVolumeFace& scvf) const
97 {
98
3/3
✓ Branch 0 taken 2168 times.
✓ Branch 1 taken 42966 times.
✓ Branch 2 taken 16572 times.
 123412 BoundaryTypes values;
99 123412 const auto& globalPos = scvf.center(); //avoid ambiguities at corners
100
101
2/2
✓ Branch 0 taken 3134 times.
✓ Branch 1 taken 58572 times.
 123412 if (verticalFlow_)
102 {
103 if constexpr (ParentType::isMomentumProblem())
104 {
105
2/2
✓ Branch 1 taken 726 times.
✓ Branch 2 taken 240 times.
 1932 if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
106 {
107 480 values.setCouplingNeumann(Indices::momentumYBalanceIdx);
108 480 values.setCouplingNeumann(Indices::momentumXBalanceIdx);
109 }
110 else
111 35076 values.setAllDirichlet();
112 }
113 else
114 {
115
2/2
✓ Branch 1 taken 1568 times.
✓ Branch 2 taken 600 times.
 4336 if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
116 values.setAllCouplingNeumann();
117 else
118 values.setAllNeumann();
119 }
120 }
121 else // horizontal flow
122 {
123 if constexpr (ParentType::isMomentumProblem())
124 {
125
4/4
✓ Branch 0 taken 12450 times.
✓ Branch 1 taken 4122 times.
✓ Branch 2 taken 4122 times.
✓ Branch 3 taken 8328 times.
 33144 if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
126 35076 values.setAllNeumann();
127
2/2
✓ Branch 1 taken 4164 times.
✓ Branch 2 taken 4164 times.
 16656 else if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
128 {
129 8328 values.setCouplingNeumann(Indices::momentumYBalanceIdx);
130 8328 values.setCouplingNeumann(Indices::momentumXBalanceIdx);
131 }
132 else
133 35076 values.setAllDirichlet();
134 }
135 else
136 {
137
2/2
✓ Branch 1 taken 30340 times.
✓ Branch 2 taken 11660 times.
 84000 if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
138 values.setAllCouplingNeumann();
139 else
140 values.setAllNeumann();
141 }
142 }
143
144 123412 return values;
145 }
146
147 /*!
148 * \brief Returns Dirichlet boundary values at a given position.
149 *
150 * \param globalPos The global position
151 */
152 3050 DirichletValues dirichletAtPos(const GlobalPosition& globalPos) const
153 {
154 3050 DirichletValues values(0.0);
155 if constexpr (ParentType::isMomentumProblem())
156 {
157
2/2
✓ Branch 0 taken 486 times.
✓ Branch 1 taken 2564 times.
 3050 if (verticalFlow_)
158 {
159
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 485 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 486 static const Scalar inletVelocity = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity");
160 486 values[Indices::velocityYIdx] = inletVelocity * globalPos[0] * (this->gridGeometry().bBoxMax()[0] - globalPos[0]);
161 }
162 }
163 3050 return values;
164 }
165
166 /*!
167 * \brief Evaluates the boundary conditions for a Neumann control volume.
168 *
169 * \param element The element for which the Neumann boundary condition is set
170 * \param fvGeometry The fvGeometry
171 * \param elemVolVars The element volume variables
172 * \param elemFaceVars The element face variables
173 * \param scvf The boundary sub control volume face
174 */
175 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
176 311656 BoundaryFluxes neumann(const Element& element,
177 const FVElementGeometry& fvGeometry,
178 const ElementVolumeVariables& elemVolVars,
179 const ElementFluxVariablesCache& elemFluxVarsCache,
180 const SubControlVolumeFace& scvf) const
181 {
182 311656 BoundaryFluxes values(0.0);
183 311656 const auto& globalPos = scvf.ipGlobal();
184 using SlipVelocityPolicy = NavierStokesSlipVelocity<typename GridGeometry::DiscretizationMethod, NavierStokes::SlipConditions::BJS>;
185 using FluxHelper = NavierStokesMomentumBoundaryFlux<typename GridGeometry::DiscretizationMethod, SlipVelocityPolicy>;
186
187 if constexpr (ParentType::isMomentumProblem())
188 {
189
190
2/2
✓ Branch 1 taken 46580 times.
✓ Branch 2 taken 73800 times.
 240760 if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
191 {
192 93160 values[scvf.normalAxis()] += couplingManager_->momentumCouplingCondition(
193 CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex,
194 fvGeometry, scvf, elemVolVars
195 );
196
197 186320 values += FluxHelper::slipVelocityMomentumFlux(
198 *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache
199 );
200 }
201 else
202 {
203
2/2
✓ Branch 0 taken 36900 times.
✓ Branch 1 taken 36900 times.
 147600 const Scalar pressure = onLeftBoundary_(globalPos) ? deltaP_ : 0.0;
204 147600 values = FluxHelper::fixedPressureMomentumFlux(
205 *this, fvGeometry, scvf, elemVolVars,
206 elemFluxVarsCache, pressure, true /*zeroNormalVelocityGradient*/
207 );
208 }
209 }
210 else
211 {
212
2/2
✓ Branch 1 taken 10080 times.
✓ Branch 2 taken 25368 times.
 70896 if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
213 {
214 20160 values = couplingManager_->massCouplingCondition(
215 CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex,
216 fvGeometry, scvf, elemVolVars
217 );
218 }
219 else
220 {
221 using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
222 50736 values = FluxHelper::scalarOutflowFlux(
223 *this, element, fvGeometry, scvf, elemVolVars
224 );
225 }
226 }
227
228 311656 return values;
229 }
230
231 // The following function defines the initial conditions
232 400 InitialValues initialAtPos(const GlobalPosition &globalPos) const
233 {
234 400 InitialValues values(0.0); //velocity is 0.0
235
236 if constexpr (!ParentType::isMomentumProblem())
237 {
238 400 values[Indices::pressureIdx] = initialPressure_;
239 }
240 return values;
241 }
242
243 /*!
244 * \brief Returns true if the scvf lies on a porous slip boundary
245 */
246
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 67944 times.
 67944 bool onSlipBoundary(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
247 {
248
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 67944 times.
 67944 assert(scvf.isFrontal());
249
3/4
✓ Branch 0 taken 67944 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 28840 times.
✓ Branch 4 taken 39104 times.
 67944 return scvf.boundary() && couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf);
250 }
251
252 /*!
253 * \brief Returns the beta value which is the alpha value divided by the square root of the (scalar-valued) interface permeability.
254 */
255 38264 Scalar betaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
256 {
257 38264 const auto& K = couplingManager_->darcyPermeability(fvGeometry, scvf);
258
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 38264 times.
 38264 const auto alphaBJ = couplingManager_->problem(CouplingManager::porousMediumIndex).spatialParams().beaversJosephCoeffAtPos(scvf.ipGlobal());
259 38264 const auto interfacePermeability = vtmv(tangentialVector, K, tangentialVector);
260 using std::sqrt;
261 38264 return alphaBJ / sqrt(interfacePermeability);
262 }
263
264 /*!
265 * \brief Returns the pressure difference between inlet and outlet for the horizontal flow scenario
266 */
267 ✗ Scalar pressureDifference() const
268 {
269 ✗ assert(!verticalFlow_);
270 ✗ return deltaP_;
271 }
272
273 /*!
274 * \brief Returns true if a vertical flow scenario is considered
275 */
276 2 bool verticalFlow() const
277
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 { return verticalFlow_; }
278
279 // \}
280
281 private:
282
283
4/4
✓ Branch 0 taken 36900 times.
✓ Branch 1 taken 36900 times.
✓ Branch 2 taken 12450 times.
✓ Branch 3 taken 4122 times.
 90372 bool onLeftBoundary_(const GlobalPosition &globalPos) const
284
4/4
✓ Branch 0 taken 36900 times.
✓ Branch 1 taken 36900 times.
✓ Branch 2 taken 12450 times.
✓ Branch 3 taken 4122 times.
 90372 { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
285
286
2/2
✓ Branch 0 taken 4122 times.
✓ Branch 1 taken 8328 times.
 12450 bool onRightBoundary_(const GlobalPosition &globalPos) const
287
2/2
✓ Branch 0 taken 4122 times.
✓ Branch 1 taken 8328 times.
 12450 { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
288
289 bool onLowerBoundary_(const GlobalPosition &globalPos) const
290 { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; }
291
292 bool onUpperBoundary_(const GlobalPosition &globalPos) const
293 { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
294
295 std::string problemName_;
296 static constexpr Scalar eps_ = 1e-6;
297 Scalar initialPressure_;
298 Scalar deltaP_;
299 bool verticalFlow_;
300
301 std::shared_ptr<CouplingManager> couplingManager_;
302 };
303 } // end namespace Dumux
304
305 #endif
306