GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 46 54 85.2%
Functions: 8 12 66.7%
Branches: 75 122 61.5%
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 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 8 FreeFlowOnePTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
62 : ParentType(gridGeometry, couplingManager, "FreeFlow")
63
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)
64 {
65
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");
66
1/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 8 deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference", 0.0);
67
68 // determine whether to simulate a vertical or horizontal flow configuration
69 8 verticalFlow_ = problemName_.find("vertical") != std::string::npos;
70 8 }
71
72 /*!
73 * \name Problem parameters
74 */
75 // \{
76
77 const std::string& name() const
78
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 { 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 12560 BoundaryTypes boundaryTypes(const Element& element,
95 const SubControlVolumeFace& scvf) const
96 {
97
2/2
✓ Branch 0 taken 2168 times.
✓ Branch 1 taken 2168 times.
 12560 BoundaryTypes values;
98
2/2
✓ Branch 0 taken 3134 times.
✓ Branch 1 taken 3146 times.
 12560 const auto& globalPos = scvf.center(); //avoid ambiguities at corners
99
100
2/2
✓ Branch 0 taken 3134 times.
✓ Branch 1 taken 3146 times.
 12560 if (verticalFlow_)
101 {
102 if constexpr (ParentType::isMomentumProblem())
103 {
104
2/2
✓ Branch 2 taken 726 times.
✓ Branch 3 taken 240 times.
 3864 if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
105 {
106 480 values.setCouplingNeumann(Indices::momentumYBalanceIdx);
107 480 values.setCouplingNeumann(Indices::momentumXBalanceIdx);
108 }
109 else
110 values.setAllDirichlet();
111 }
112 else
113 {
114
2/2
✓ Branch 2 taken 600 times.
✓ Branch 3 taken 1568 times.
 8672 if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
115 values.setAllCouplingNeumann();
116 else
117 values.setAllNeumann();
118 }
119 }
120 else // horizontal flow
121 {
122 if constexpr (ParentType::isMomentumProblem())
123 {
124 6852 if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
125 values.setAllNeumann();
126
2/2
✓ Branch 2 taken 246 times.
✓ Branch 3 taken 246 times.
 1968 else if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
127 {
128 492 values.setCouplingNeumann(Indices::momentumYBalanceIdx);
129 492 values.setCouplingNeumann(Indices::momentumXBalanceIdx);
130 }
131 else
132 values.setAllDirichlet();
133 }
134 else
135 {
136
2/2
✓ Branch 2 taken 600 times.
✓ Branch 3 taken 1568 times.
 8672 if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
137 values.setAllCouplingNeumann();
138 else
139 values.setAllNeumann();
140 }
141 }
142
143 12560 return values;
144 }
145
146 /*!
147 * \brief Returns Dirichlet boundary values at a given position.
148 *
149 * \param globalPos The global position
150 */
151 1304 DirichletValues dirichletAtPos(const GlobalPosition& globalPos) const
152 {
153 1304 DirichletValues values(0.0);
154 if constexpr (ParentType::isMomentumProblem())
155 {
156
2/2
✓ Branch 0 taken 486 times.
✓ Branch 1 taken 166 times.
 1304 if (verticalFlow_)
157 {
158
6/8
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 484 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
 972 static const Scalar inletVelocity = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity");
159 6804 values[Indices::velocityYIdx] = inletVelocity * globalPos[0] * (this->gridGeometry().bBoxMax()[0] - globalPos[0]);
160 }
161 }
162 1304 return values;
163 }
164
165 /*!
166 * \brief Evaluates the boundary conditions for a Neumann control volume.
167 *
168 * \param element The element for which the Neumann boundary condition is set
169 * \param fvGeometry The fvGeometry
170 * \param elemVolVars The element volume variables
171 * \param elemFaceVars The element face variables
172 * \param scvf The boundary sub control volume face
173 */
174 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
175 23008 BoundaryFluxes neumann(const Element& element,
176 const FVElementGeometry& fvGeometry,
177 const ElementVolumeVariables& elemVolVars,
178 const ElementFluxVariablesCache& elemFluxVarsCache,
179 const SubControlVolumeFace& scvf) const
180 {
181 23008 BoundaryFluxes values(0.0);
182 23008 const auto& globalPos = scvf.ipGlobal();
183 using SlipVelocityPolicy = NavierStokesSlipVelocity<typename GridGeometry::DiscretizationMethod, NavierStokes::SlipConditions::BJS>;
184 using FluxHelper = NavierStokesMomentumBoundaryFlux<typename GridGeometry::DiscretizationMethod, SlipVelocityPolicy>;
185
186 if constexpr (ParentType::isMomentumProblem())
187 {
188
189
2/2
✓ Branch 2 taken 4400 times.
✓ Branch 3 taken 3728 times.
 32512 if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
190 {
191 17600 values += couplingManager_->momentumCouplingCondition(
192 CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex,
193 fvGeometry, scvf, elemVolVars
194 );
195
196 17600 values += FluxHelper::slipVelocityMomentumFlux(
197 *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache
198 );
199 }
200 else
201 {
202 3728 const Scalar pressure = onLeftBoundary_(globalPos) ? deltaP_ : 0.0;
203 7456 values = FluxHelper::fixedPressureMomentumFlux(
204 *this, fvGeometry, scvf, elemVolVars,
205 elemFluxVarsCache, pressure, true /*zeroNormalVelocityGradient*/
206 );
207 }
208 }
209 else
210 {
211
2/2
✓ Branch 2 taken 960 times.
✓ Branch 3 taken 2416 times.
 13504 if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
212 {
213 3840 values = couplingManager_->massCouplingCondition(
214 CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex,
215 fvGeometry, scvf, elemVolVars
216 );
217 }
218 else
219 {
220 using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
221 4832 values = FluxHelper::scalarOutflowFlux(
222 *this, element, fvGeometry, scvf, elemVolVars
223 );
224 }
225 }
226
227 23008 return values;
228 }
229
230 /*!
231 * \brief Returns true if the scvf lies on a porous slip boundary
232 */
233 ✗ bool onSlipBoundary(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
234 {
235 ✗ assert(scvf.isFrontal());
236 ✗ return scvf.boundary() && couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf);
237 }
238
239 /*!
240 * \brief Returns the beta value which is the alpha value divided by the square root of the (scalar-valued) interface permeability.
241 */
242 3608 Scalar betaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
243 {
244 7216 const auto& K = couplingManager_->darcyPermeability(fvGeometry, scvf);
245
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3608 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3608 times.
 7216 const auto alphaBJ = couplingManager_->problem(CouplingManager::porousMediumIndex).spatialParams().beaversJosephCoeffAtPos(scvf.ipGlobal());
246 3608 const auto interfacePermeability = vtmv(tangentialVector, K, tangentialVector);
247 using std::sqrt;
248 3608 return alphaBJ / sqrt(interfacePermeability);
249 }
250
251 /*!
252 * \brief Returns the pressure difference between inlet and outlet for the horizontal flow scenario
253 */
254 ✗ Scalar pressureDifference() const
255 {
256 ✗ assert(!verticalFlow_);
257 ✗ return deltaP_;
258 }
259
260 /*!
261 * \brief Returns true if a vertical flow scenario is considered
262 */
263 ✗ bool verticalFlow() const
264 ✗ { return verticalFlow_; }
265
266 // \}
267
268 private:
269
270 bool onLeftBoundary_(const GlobalPosition &globalPos) const
271
20/20
✓ Branch 0 taken 1864 times.
✓ Branch 1 taken 1864 times.
✓ Branch 2 taken 1864 times.
✓ Branch 3 taken 1864 times.
✓ Branch 4 taken 1864 times.
✓ Branch 5 taken 1864 times.
✓ Branch 6 taken 1864 times.
✓ Branch 7 taken 1864 times.
✓ Branch 8 taken 1864 times.
✓ Branch 9 taken 1864 times.
✓ Branch 10 taken 735 times.
✓ Branch 11 taken 243 times.
✓ Branch 12 taken 735 times.
✓ Branch 13 taken 243 times.
✓ Branch 14 taken 735 times.
✓ Branch 15 taken 243 times.
✓ Branch 16 taken 735 times.
✓ Branch 17 taken 243 times.
✓ Branch 18 taken 735 times.
✓ Branch 19 taken 243 times.
 23530 { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
272
273 bool onRightBoundary_(const GlobalPosition &globalPos) const
274
10/10
✓ Branch 0 taken 243 times.
✓ Branch 1 taken 492 times.
✓ Branch 2 taken 243 times.
✓ Branch 3 taken 492 times.
✓ Branch 4 taken 243 times.
✓ Branch 5 taken 492 times.
✓ Branch 6 taken 243 times.
✓ Branch 7 taken 492 times.
✓ Branch 8 taken 243 times.
✓ Branch 9 taken 492 times.
 3675 { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
275
276 bool onLowerBoundary_(const GlobalPosition &globalPos) const
277 { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; }
278
279 bool onUpperBoundary_(const GlobalPosition &globalPos) const
280 { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
281
282 std::string problemName_;
283 static constexpr Scalar eps_ = 1e-6;
284 Scalar deltaP_;
285 bool verticalFlow_;
286
287 std::shared_ptr<CouplingManager> couplingManager_;
288 };
289 } // end namespace Dumux
290
291 #endif
292