GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: dumux/test/freeflow/rans/problem.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 103 103 100.0%
Functions: 62 62 100.0%
Branches: 132 154 85.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 RANSTests
10 * \brief Pipe flow test for the staggered grid RANS model
11 *
12 * This test simulates pipe flow experiments performed by John Laufer in 1954
13 * \cite Laufer1954a.
14 */
15 #ifndef DUMUX_PIPE_LAUFER_PROBLEM_HH
16 #define DUMUX_PIPE_LAUFER_PROBLEM_HH
17
18 #include <dumux/common/properties.hh>
19 #include <dumux/common/parameters.hh>
20 #include <dumux/common/timeloop.hh>
21 #include <dumux/common/numeqvector.hh>
22
23 #include <dumux/freeflow/rans/boundarytypes.hh>
24 #include <dumux/freeflow/rans/problem.hh>
25 #include <dumux/freeflow/turbulencemodel.hh>
26 #include <dumux/freeflow/turbulenceproperties.hh>
27
28 namespace Dumux {
29
30 /*!
31 * \ingroup NavierStokesTests
32 * \brief Test problem for the one-phase (Navier-) Stokes problem in a channel.
33 *
34 * This test simulates is based on pipe flow experiments by
35 * John Laufers experiments in 1954 \cite Laufer1954a.
36 */
37 template <class TypeTag>
38 class PipeLauferProblem : public RANSProblem<TypeTag>
39 {
40 using ParentType = RANSProblem<TypeTag>;
41
42 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
43 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
44 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
45 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
46 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
47 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
48
49 static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
50 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
51 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
52 using BoundaryTypes = Dumux::RANSBoundaryTypes<ModelTraits, ModelTraits::numEq()>;
53 using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
54 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
55 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
56 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
57 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
58
59 using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
60
61
62 public:
63 15 PipeLauferProblem(std::shared_ptr<const GridGeometry> gridGeometry)
64
2/4
✓ Branch 2 taken 15 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 15 times.
✗ Branch 7 not taken.
45 : ParentType(gridGeometry), eps_(1e-6)
65 {
66
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
15 inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
67
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
15 inletTemperature_ = getParam<Scalar>("Problem.InletTemperature", 283.15);
68
2/4
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
15 wallTemperature_ = getParam<Scalar>("Problem.WallTemperature", 303.15);
69
70 FluidSystem::init();
71 Dumux::TurbulenceProperties<Scalar, dimWorld, true> turbulenceProperties;
72
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
15 FluidState fluidState;
73 15 fluidState.setPressure(0, 1e5);
74
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
15 fluidState.setTemperature(this->spatialParams().temperatureAtPos({}));
75 15 Scalar density = FluidSystem::density(fluidState, 0);
76 15 Scalar kinematicViscosity = FluidSystem::viscosity(fluidState, 0) / density;
77
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
15 Scalar diameter = this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1];
78
79 // ideally the viscosityTilde parameter as inflow for the Spalart-Allmaras model should be zero
80
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
15 viscosityTilde_ = 1e-3 * turbulenceProperties.viscosityTilde(inletVelocity_, diameter, kinematicViscosity);
81
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
15 turbulentKineticEnergy_ = turbulenceProperties.turbulentKineticEnergy(inletVelocity_, diameter, kinematicViscosity);
82
83 if (ModelTraits::turbulenceModel() == TurbulenceModel::komega || ModelTraits::turbulenceModel() == TurbulenceModel::sst)
84
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
5 dissipation_ = turbulenceProperties.dissipationRate(inletVelocity_, diameter, kinematicViscosity);
85 else
86
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
10 dissipation_ = turbulenceProperties.dissipation(inletVelocity_, diameter, kinematicViscosity);
87
88 if (ModelTraits::turbulenceModel() == TurbulenceModel::oneeq)
89
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
2 initializationTime_ = getParam<Scalar>("TimeLoop.Initialization", 1.0);
90 else
91
2/4
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
13 initializationTime_ = getParam<Scalar>("TimeLoop.Initialization", -1.0);
92
93
1/2
✓ Branch 3 taken 15 times.
✗ Branch 4 not taken.
15 turbulenceModelName_ = turbulenceModelToString(ModelTraits::turbulenceModel());
94
2/4
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
15 std::cout << "Using the "<< turbulenceModelName_ << " Turbulence Model. \n";
95 15 std::cout << std::endl;
96 15 }
97
98 /*!
99 * \name Boundary conditions
100 */
101 // \{
102
103 /*!
104 * \brief Specifies which kind of boundary condition should be
105 * used for which equation on a given boundary control volume.
106 *
107 * \param globalPos The position of the center of the finite volume
108 */
109 242886163 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
110 {
111 242886163 BoundaryTypes values;
112
113 // common boundary types for all turbulence models
114
2/2
✓ Branch 0 taken 19594712 times.
✓ Branch 1 taken 223291451 times.
242886163 if(isOutlet_(globalPos))
115 19594712 values.setDirichlet(Indices::pressureIdx);
116 else
117 {
118 223291451 values.setDirichlet(Indices::velocityXIdx);
119 223291451 values.setDirichlet(Indices::velocityYIdx);
120 }
121
122
4/4
✓ Branch 0 taken 233808525 times.
✓ Branch 1 taken 9077638 times.
✓ Branch 2 taken 9077623 times.
✓ Branch 3 taken 224730902 times.
242886163 if (isLowerWall_(globalPos) || isUpperWall_(globalPos)) // All walls
123
2/2
✓ Branch 0 taken 11456862 times.
✓ Branch 1 taken 94010613 times.
242886163 values.setWall();
124
125 #if NONISOTHERMAL
126
2/2
✓ Branch 0 taken 11456862 times.
✓ Branch 1 taken 94010613 times.
105467475 if(isOutlet_(globalPos))
127 11456862 values.setOutflow(Indices::energyEqIdx);
128 else
129 94010613 values.setDirichlet(Indices::temperatureIdx);
130 #endif
131 // turbulence model-specific boundary types
132 240721791 setBcTypes_(values, globalPos);
133
134 242886163 return values;
135 }
136
137 /*!
138 * \brief Returns whether a fixed Dirichlet value shall be used inside a given cell.
139 *
140 * \param element The finite element
141 * \param fvGeometry The finite-volume geometry
142 * \param scv The sub control volume
143 * \param pvIdx The primary variable index in the solution vector
144 */
145 template<class Element, class FVElementGeometry, class SubControlVolume>
146 57385236 bool isDirichletCell(const Element& element,
147 const FVElementGeometry& fvGeometry,
148 const SubControlVolume& scv,
149 int pvIdx) const
150
4/4
✓ Branch 0 taken 720384 times.
✓ Branch 1 taken 26905216 times.
✓ Branch 2 taken 726676 times.
✓ Branch 3 taken 649040 times.
57385236 { return isDirichletCell_(element, fvGeometry, scv, pvIdx); }
151
152 /*!
153 * \brief Evaluate the boundary conditions for dirichlet values at the given boundary.
154 *
155 * \param element The finite element
156 * \param scvf the sub control volume face
157 * \note used for cell-centered discretization schemes
158 */
159 16790462 PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
160 {
161 16790462 const auto globalPos = scvf.ipGlobal();
162
1/2
✓ Branch 6 taken 39076 times.
✗ Branch 7 not taken.
16110130 PrimaryVariables values(initialAtPos(globalPos));
163
164 #if NONISOTHERMAL
165
4/4
✓ Branch 0 taken 6069906 times.
✓ Branch 1 taken 1166786 times.
✓ Branch 2 taken 1166786 times.
✓ Branch 3 taken 4903120 times.
7236692 values[Indices::temperatureIdx] = (isLowerWall_(globalPos) || isUpperWall_(globalPos)) ? wallTemperature_ : inletTemperature_;
166 #endif
167
168 7236692 return values;
169 }
170
171 /*!
172 * \brief Evaluate the boundary conditions for fixed values at cell centers
173 *
174 * \param element The finite element
175 * \param scv the sub control volume
176 * \note used for cell-centered discretization schemes
177 */
178 3831222 PrimaryVariables dirichlet([[maybe_unused]] const Element& element, const SubControlVolume& scv) const
179 {
180 if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::kepsilon
181 || ModelTraits::turbulenceModel() == TurbulenceModel::komega
182 || ModelTraits::turbulenceModel() == TurbulenceModel::sst)
183 3831222 return dirichletTurbulentTwoEq_(element, scv);
184 else
185 {
186 const auto globalPos = scv.center();
187 PrimaryVariables values(initialAtPos(globalPos));
188 return values;
189 }
190 }
191
192 /*!
193 * \brief Evaluate the initial value at the given global position.
194 *
195 * \param globalPos The global position
196 */
197 103181226 PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
198 {
199 // common initial conditions for all turbulence models
200
2/2
✓ Branch 0 taken 98546698 times.
✓ Branch 1 taken 4634528 times.
103181226 PrimaryVariables values(0.0);
201 103181226 values[Indices::pressureIdx] = 1.0e+5;
202
2/2
✓ Branch 0 taken 99160050 times.
✓ Branch 1 taken 4021176 times.
103181226 values[Indices::velocityXIdx] = time() > initializationTime_
203
2/2
✓ Branch 0 taken 98546698 times.
✓ Branch 1 taken 4634528 times.
103181226 ? inletVelocity_
204 4634528 : time() / initializationTime_ * inletVelocity_;
205
4/4
✓ Branch 0 taken 99160050 times.
✓ Branch 1 taken 4021176 times.
✓ Branch 2 taken 4021176 times.
✓ Branch 3 taken 95138874 times.
103181226 if (isLowerWall_(globalPos) || isUpperWall_(globalPos))
206 8042352 values[Indices::velocityXIdx] = 0.0;
207 #if NONISOTHERMAL
208
6/6
✓ Branch 0 taken 38625684 times.
✓ Branch 1 taken 1166844 times.
✓ Branch 2 taken 1166844 times.
✓ Branch 3 taken 37458840 times.
✓ Branch 4 taken 35131518 times.
✓ Branch 5 taken 1060578 times.
39792528 values[Indices::temperatureIdx] = (isLowerWall_(globalPos) || isUpperWall_(globalPos)) ? wallTemperature_ : inletTemperature_;
209 #endif
210 // turbulence model-specific initial conditions
211 96940442 setInitialAtPos_(values, globalPos);
212 103181226 return values;
213 }
214
215 // \}
216
217 15 void setTimeLoop(TimeLoopPtr timeLoop)
218
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
15 { timeLoop_ = timeLoop; }
219
220 103181226 Scalar time() const
221
2/2
✓ Branch 0 taken 98546698 times.
✓ Branch 1 taken 4634528 times.
103181226 { return timeLoop_->time(); }
222
223 private:
224 bool isInlet_(const GlobalPosition& globalPos) const
225 { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
226
227
2/2
✓ Branch 0 taken 19594712 times.
✓ Branch 1 taken 223291451 times.
242886163 bool isOutlet_(const GlobalPosition& globalPos) const
228
6/6
✓ Branch 0 taken 19594712 times.
✓ Branch 1 taken 223291451 times.
✓ Branch 2 taken 19361500 times.
✓ Branch 3 taken 222849944 times.
✓ Branch 4 taken 10723528 times.
✓ Branch 5 taken 93254294 times.
589075429 { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
229
230
6/6
✓ Branch 0 taken 141227660 times.
✓ Branch 1 taken 7529561 times.
✓ Branch 2 taken 157345928 times.
✓ Branch 3 taken 5479251 times.
✓ Branch 4 taken 40464893 times.
✓ Branch 5 taken 1256788 times.
353304081 bool isLowerWall_(const GlobalPosition& globalPos) const
231
10/10
✓ Branch 0 taken 141227660 times.
✓ Branch 1 taken 7529561 times.
✓ Branch 2 taken 161514948 times.
✓ Branch 3 taken 5412072 times.
✓ Branch 4 taken 97121218 times.
✓ Branch 5 taken 3943527 times.
✓ Branch 6 taken 33999980 times.
✓ Branch 7 taken 1166836 times.
✓ Branch 8 taken 36970727 times.
✓ Branch 9 taken 1150522 times.
490037051 { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; }
232
233
10/10
✓ Branch 0 taken 7529554 times.
✓ Branch 1 taken 133698106 times.
✓ Branch 2 taken 5412066 times.
✓ Branch 3 taken 156102882 times.
✓ Branch 4 taken 3943526 times.
✓ Branch 5 taken 93177692 times.
✓ Branch 6 taken 1166836 times.
✓ Branch 7 taken 32833144 times.
✓ Branch 8 taken 1150521 times.
✓ Branch 9 taken 35820206 times.
470834533 bool isUpperWall_(const GlobalPosition& globalPos) const
234
10/10
✓ Branch 0 taken 7529554 times.
✓ Branch 1 taken 133698106 times.
✓ Branch 2 taken 5412066 times.
✓ Branch 3 taken 156102882 times.
✓ Branch 4 taken 3943526 times.
✓ Branch 5 taken 93177692 times.
✓ Branch 6 taken 1166836 times.
✓ Branch 7 taken 32833144 times.
✓ Branch 8 taken 1150521 times.
✓ Branch 9 taken 35820206 times.
470834533 { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
235
236 //! Initial conditions for the komega, kepsilon and lowrekepsilon turbulence models
237 96940442 void setInitialAtPos_([[maybe_unused]] PrimaryVariables& values,
238 [[maybe_unused]] const GlobalPosition &globalPos) const
239 {
240 if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 0) // zero equation models
241 return;
242 else if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 1) // one equation models
243 {
244
2/2
✓ Branch 0 taken 10098046 times.
✓ Branch 1 taken 432814 times.
10530860 values[Indices::viscosityTildeIdx] = viscosityTilde_;
245
4/4
✓ Branch 0 taken 10098046 times.
✓ Branch 1 taken 432814 times.
✓ Branch 2 taken 432814 times.
✓ Branch 3 taken 9665232 times.
10530860 if (isLowerWall_(globalPos) || isUpperWall_(globalPos))
246 865628 values[Indices::viscosityTildeIdx] = 0.0;
247 }
248 else // two equation models
249 {
250 static_assert(numTurbulenceEq(ModelTraits::turbulenceModel()) == 2, "Only reached by 2eq models");
251
2/2
✓ Branch 0 taken 53493590 times.
✓ Branch 1 taken 2445588 times.
86409582 values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_;
252
2/2
✓ Branch 0 taken 83072322 times.
✓ Branch 1 taken 3337260 times.
86409582 values[Indices::dissipationIdx] = dissipation_;
253
4/4
✓ Branch 0 taken 83072322 times.
✓ Branch 1 taken 3337260 times.
✓ Branch 2 taken 3337260 times.
✓ Branch 3 taken 79735062 times.
86409582 if (isLowerWall_(globalPos) || isUpperWall_(globalPos))
254 {
255 6674520 values[Indices::turbulentKineticEnergyIdx] = 0.0;
256 6674520 values[Indices::dissipationIdx] = 0.0;
257 }
258 }
259 }
260
261 //! Boundary condition types for the one-eq turbulence model
262
2/2
✓ Branch 0 taken 18628166 times.
✓ Branch 1 taken 222093625 times.
240721791 void setBcTypes_([[maybe_unused]] BoundaryTypes& values,
263 [[maybe_unused]] const GlobalPosition& pos) const
264 {
265 if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 0) // zero equation models
266 return;
267 else if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 1) // one equation models
268 {
269
2/2
✓ Branch 0 taken 1836194 times.
✓ Branch 1 taken 11147690 times.
12983884 if(isOutlet_(pos))
270 1836194 values.setOutflow(Indices::viscosityTildeIdx);
271 else // walls and inflow
272 11147690 values.setDirichlet(Indices::viscosityTildeIdx);
273 }
274 else // two equation models
275 {
276 static_assert(numTurbulenceEq(ModelTraits::turbulenceModel()) == 2, "Only reached by 2eq models");
277
2/2
✓ Branch 0 taken 16791972 times.
✓ Branch 1 taken 210945935 times.
227737907 if(isOutlet_(pos))
278 {
279 16791972 values.setOutflow(Indices::turbulentKineticEnergyEqIdx);
280 16791972 values.setOutflow(Indices::dissipationEqIdx);
281 }
282 else
283 {
284 // walls and inflow
285 210945935 values.setDirichlet(Indices::turbulentKineticEnergyIdx);
286 210945935 values.setDirichlet(Indices::dissipationIdx);
287 }
288 }
289 227737907 }
290
291 template<class Element, class FVElementGeometry, class SubControlVolume>
292 57385236 bool isDirichletCell_([[maybe_unused]] const Element& element,
293 const FVElementGeometry& fvGeometry,
294 [[maybe_unused]] const SubControlVolume& scv,
295 const int& pvIdx) const
296 {
297 if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::kepsilon)
298 {
299 29089300 const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(element);
300 // For the kepsilon model we set fixed values within the matching point and at the wall
301
2/2
✓ Branch 1 taken 4241312 times.
✓ Branch 2 taken 24847988 times.
29089300 if (this->inNearWallRegion(eIdx))
302 4241312 return pvIdx == Indices::turbulentKineticEnergyEqIdx || pvIdx == Indices::dissipationEqIdx;
303
2/2
✓ Branch 1 taken 2885116 times.
✓ Branch 2 taken 21962872 times.
24847988 if (this->isMatchingPoint(eIdx))
304 2885116 return pvIdx == Indices::dissipationEqIdx;
305 return false;
306 }
307 else if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::komega ||
308 ModelTraits::turbulenceModel() == TurbulenceModel::sst)
309 {
310 // For the komega model we set a fixed dissipation (omega) for all cells at the wall
311
2/2
✓ Branch 0 taken 112812904 times.
✓ Branch 1 taken 27554256 times.
140367160 for (const auto& scvf : scvfs(fvGeometry))
312
4/4
✓ Branch 0 taken 2424080 times.
✓ Branch 1 taken 110388824 times.
✓ Branch 2 taken 741680 times.
✓ Branch 3 taken 1682400 times.
112812904 if (this->boundaryTypes(element, scvf).hasWall() && pvIdx == Indices::dissipationIdx)
313 741680 return true;
314 return false;
315 }
316 else
317 return ParentType::isDirichletCell(element, fvGeometry, scv, pvIdx);
318 }
319
320 //! Specialization for the kepsilon and komega
321 template<class Element, class SubControlVolume>
322 3831222 PrimaryVariables dirichletTurbulentTwoEq_(const Element& element,
323 const SubControlVolume& scv) const
324 {
325 3831222 const auto globalPos = scv.center();
326 3831222 PrimaryVariables values(initialAtPos(globalPos));
327 3831222 unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
328
329 if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::kepsilon)
330 {
331 // For the kepsilon model we set a fixed value for the turbulent kinetic energy and the dissipation
332 3110838 values[Indices::turbulentKineticEnergyEqIdx] = this->turbulentKineticEnergyWallFunction(elementIdx);
333 3110838 values[Indices::dissipationEqIdx] = this->dissipationWallFunction(elementIdx);
334 3110838 return values;
335 }
336 else
337 {
338 static_assert(ModelTraits::turbulenceModel() == TurbulenceModel::komega
339 || ModelTraits::turbulenceModel() == TurbulenceModel::sst, "Only valid for Komega");
340 // For the komega model we set a fixed value for the dissipation
341 720384 const auto wallDistance = ParentType::wallDistance(elementIdx);
342 using std::pow;
343 720384 values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity(elementIdx)
344 720384 / (ParentType::betaOmega() * wallDistance * wallDistance);
345 720384 return values;
346 }
347 }
348
349 Scalar eps_;
350 Scalar inletVelocity_;
351 Scalar inletTemperature_;
352 Scalar wallTemperature_;
353 Scalar initializationTime_;
354 Scalar viscosityTilde_;
355 Scalar turbulentKineticEnergy_;
356 Scalar dissipation_;
357 std::string turbulenceModelName_;
358 TimeLoopPtr timeLoop_;
359 };
360 } // end namespace Dumux
361
362 #endif
363