GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/test/freeflow/navierstokes/channel/2d/problem.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 112 113 99.1%
Functions: 10 10 100.0%
Branches: 102 158 64.6%
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 Channel flow test for the staggered grid (Navier-)Stokes model.
11 */
12
13 #ifndef DUMUX_CHANNEL_TEST_PROBLEM_HH
14 #define DUMUX_CHANNEL_TEST_PROBLEM_HH
15
16 #include <dumux/common/parameters.hh>
17 #include <dumux/common/properties.hh>
18 #include <dumux/common/timeloop.hh>
19
20 #include <dumux/freeflow/navierstokes/boundarytypes.hh>
21
22 #include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
23 #include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
24 #include <dumux/freeflow/navierstokes/mass/1p/advectiveflux.hh>
25
26 namespace Dumux {
27
28 /*!
29 * \ingroup NavierStokesTests
30 * \brief Test problem for the one-phase (Navier-) Stokes problem in a channel.
31 *
32 * Flow from left to right in a two-dimensional channel is considered. At the inlet (left),
33 * fixed values for velocity are set, while at the outlet (right), a fixed pressure
34 * boundary condition is used. The channel is confined by solid walls at the top and bottom
35 * of the domain which corresponds to no-slip/no-flow conditions.
36 * For the non-isothermal test, water of increased temperature is injected at the inlet
37 * while the walls are fully isolating.
38 */
39 template <class TypeTag, class BaseProblem>
40 18 class ChannelTestProblem : public BaseProblem
41 {
42 using ParentType = BaseProblem;
43
44 using BoundaryTypes = typename ParentType::BoundaryTypes;
45 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
46 using FVElementGeometry = typename GridGeometry::LocalView;
47 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
48 using InitialValues = typename ParentType::InitialValues;
49 using Sources = typename ParentType::Sources;
50 using DirichletValues = typename ParentType::DirichletValues;
51 using BoundaryFluxes = typename ParentType::BoundaryFluxes;
52 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
53 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
54 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
55 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
56 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
57
58 static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
59 using Element = typename FVElementGeometry::Element;
60 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
61
62 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
63
64 // the types of outlet boundary conditions
65 enum class OutletCondition
66 {
67 outflow, unconstrainedOutflow, neumannXdirichletY, neumannXneumannY
68 };
69
70 public:
71 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
72
73 36 ChannelTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
74
2/6
✓ Branch 3 taken 18 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 18 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 144 : ParentType(gridGeometry, couplingManager)
75 {
76
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
77
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 dynamicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0)
78
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 * getParam<Scalar>("Component.LiquidDensity", 1.0);
79
80
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 const auto outletBC = getParam<std::string>("Problem.OutletCondition", "Outflow");
81
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 6 times.
 36 if (outletBC == "Outflow")
82 24 outletCondition_ = OutletCondition::outflow;
83
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 4 times.
 12 else if (outletBC == "UnconstrainedOutflow")
84 4 outletCondition_ = OutletCondition::unconstrainedOutflow;
85
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 8 else if (outletBC == "NeumannX_DirichletY")
86 4 outletCondition_ = OutletCondition::neumannXdirichletY;
87
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 else if (outletBC == "NeumannX_NeumannY")
88 4 outletCondition_ = OutletCondition::neumannXneumannY;
89 else
90 ✗ DUNE_THROW(Dune::InvalidStateException, outletBC + " is not a valid outlet boundary condition");
91
92
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 useVelocityProfile_ = getParam<bool>("Problem.UseVelocityProfile", false);
93
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 enableGravity_ = getParam<bool>("Problem.EnableGravity", false);
94
95
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 outletPressure_ = getParam<Scalar>("Problem.OutletPressure", 1.1e5);
96
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 referencePressure_ = getParam<Scalar>("Problem.RefPressure", 1.1e5);
97
98 36 FluidState fluidState;
99
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 fluidState.setPressure(0, outletPressure_);
100
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 const auto upperRight = getParam<GlobalPosition>("Grid.UpperRight");
101
2/4
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
 36 fluidState.setTemperature(this->spatialParams().temperatureAtPos(upperRight));
102
103 36 outletDensity_ = FluidSystem::density(fluidState, 0);
104 36 }
105
106 /*!
107 * \name Boundary conditions
108 */
109 // \{
110
111 /*!
112 * \brief Specifies which kind of boundary condition should be
113 * used for which equation on a given boundary control volume.
114 *
115 * \param globalPos The position of the center of the finite volume
116 */
117
4/4
✓ Branch 0 taken 13500 times.
✓ Branch 1 taken 67500 times.
✓ Branch 2 taken 3875 times.
✓ Branch 3 taken 19375 times.
 388726 BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
118 {
119
6/6
✓ Branch 0 taken 197700 times.
✓ Branch 1 taken 182298 times.
✓ Branch 2 taken 31315 times.
✓ Branch 3 taken 19375 times.
✓ Branch 5 taken 23330 times.
✓ Branch 6 taken 38442 times.
 818898 BoundaryTypes values;
120
121 if constexpr (ParentType::isMomentumProblem())
122 {
123 162048 values.setDirichlet(Indices::velocityXIdx);
124 162048 values.setDirichlet(Indices::velocityYIdx);
125
126
2/2
✓ Branch 0 taken 46028 times.
✓ Branch 1 taken 65882 times.
 162048 if (isOutlet_(globalPos))
127 {
128
2/2
✓ Branch 0 taken 44878 times.
✓ Branch 1 taken 1150 times.
 68726 if (outletCondition_ == OutletCondition::neumannXdirichletY)
129 {
130 1150 values.setNeumann(Indices::momentumXBalanceIdx);
131 1150 values.setDirichlet(Indices::velocityYIdx);
132 }
133 else
134 162048 values.setAllNeumann();
135 }
136 }
137 else
138 {
139 288450 values.setAllNeumann();
140
141 288450 if (isInlet_(globalPos))
142 {
143 48075 values.setDirichlet(Indices::pressureIdx);
144 #if NONISOTHERMAL
145 30700 values.setDirichlet(Indices::temperatureIdx);
146 #endif
147 }
148 }
149
150 346248 return values;
151 }
152
153
154 /*!
155 * \brief Evaluates the boundary conditions for a Dirichlet control volume.
156 *
157 * \param globalPos The center of the finite volume which ought to be set.
158 */
159
4/5
✓ Branch 0 taken 121308 times.
✓ Branch 1 taken 2766 times.
✓ Branch 2 taken 2050 times.
✓ Branch 3 taken 5625 times.
✗ Branch 4 not taken.
 131749 DirichletValues dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
160 {
161 131749 const auto& globalPos = scvf.ipGlobal();
162 127874 DirichletValues values = initialAtPos(globalPos);
163
164 if constexpr (ParentType::isMomentumProblem())
165 {
166 124074 values[Indices::velocityXIdx] = parabolicProfile(globalPos[1], inletVelocity_);
167 245382 if (!useVelocityProfile_ && isInlet_(globalPos))
168 values[Indices::velocityXIdx] = inletVelocity_;
169 }
170 else
171 {
172 7675 values[Indices::pressureIdx] = this->couplingManager().cellPressure(element, scvf);
173
174 #if NONISOTHERMAL
175 // give the system some time so that the pressure can equilibrate, then start the injection of the hot liquid
176 3800 if (time() >= 200.0)
177 2050 values[Indices::temperatureIdx] = 293.15;
178 #endif
179 }
180
181 124074 return values;
182 }
183
184 /*!
185 * \brief Evaluates the boundary conditions for a Neumann control volume.
186 *
187 * \param element The element for which the Neumann boundary condition is set
188 * \param fvGeometry The fvGeometry
189 * \param elemVolVars The element volume variables
190 * \param elemFaceVars The element face variables
191 * \param scvf The boundary sub control volume face
192 */
193 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
194 468856 BoundaryFluxes neumann(const Element& element,
195 const FVElementGeometry& fvGeometry,
196 const ElementVolumeVariables& elemVolVars,
197 const ElementFluxVariablesCache& elemFluxVarsCache,
198 const SubControlVolumeFace& scvf) const
199 {
200 468856 BoundaryFluxes values(0.0);
201
202
2/2
✓ Branch 0 taken 25050 times.
✓ Branch 1 taken 100200 times.
 468856 Scalar outletPressure = outletPressure_;
203
2/2
✓ Branch 0 taken 9296 times.
✓ Branch 1 taken 99882 times.
 218356 if (enableGravity_)
204 18592 outletPressure += hydrostaticPressure_(scvf.center());
205
206 if constexpr (ParentType::isMomentumProblem())
207 {
208 using FluxHelper = NavierStokesMomentumBoundaryFlux<typename GridGeometry::DiscretizationMethod>;
209
210
2/2
✓ Branch 0 taken 4596 times.
✓ Branch 1 taken 104582 times.
 218356 if (outletCondition_ == OutletCondition::unconstrainedOutflow)
211 9192 values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, outletPressure, false /*zeroNormalVelocityGradient*/);
212
2/2
✓ Branch 0 taken 100492 times.
✓ Branch 1 taken 4090 times.
 209164 else if (outletCondition_ == OutletCondition::outflow)
213 200984 values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, outletPressure, true /*zeroNormalVelocityGradient*/);
214 else
215 {
216
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4090 times.
 8180 assert(outletCondition_ == OutletCondition::neumannXneumannY || outletCondition_ == OutletCondition::neumannXdirichletY);
217 8180 values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, outletPressure, false /*zeroNormalVelocityGradient*/);
218
219 8180 Dune::FieldMatrix<Scalar, dimWorld, dimWorld> shearRate(0.0); // gradV + gradV^T
220 8180 shearRate[0][1] = dudy(scvf.ipGlobal()[1], inletVelocity_); // we assume du/dx = dv/dy = dv/dx = 0
221 8180 shearRate[1][0] = shearRate[0][1];
222
223 8180 const auto normal = scvf.unitOuterNormal();
224 8180 BoundaryFluxes normalGradient(0.0);
225 8180 shearRate.mv(normal, normalGradient);
226 24540 values -= this->effectiveViscosity(element, fvGeometry, scvf)*normalGradient;
227 }
228 }
229 else
230 {
231 using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
232
233 250500 if (isOutlet_(scvf.ipGlobal()))
234 50100 values = FluxHelper::scalarOutflowFlux(*this, element, fvGeometry, scvf, elemVolVars);
235 }
236
237 468856 return values;
238 }
239
240 /*!
241 * \brief A parabolic velocity profile.
242 *
243 * \param y The position where the velocity is evaluated.
244 * \param vMax The profile's maximum velocity.
245 */
246 136724 Scalar parabolicProfile(const Scalar y, const Scalar vMax) const
247 {
248
2/4
✓ Branch 0 taken 121308 times.
✓ Branch 1 taken 2766 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 136724 const Scalar yMin = this->gridGeometry().bBoxMin()[1];
249
2/4
✓ Branch 0 taken 121308 times.
✓ Branch 1 taken 2766 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 136724 const Scalar yMax = this->gridGeometry().bBoxMax()[1];
250
2/4
✓ Branch 0 taken 121308 times.
✓ Branch 1 taken 2766 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 136724 return vMax * (y - yMin)*(yMax - y) / (0.25*(yMax - yMin)*(yMax - yMin));
251 }
252
253 /*!
254 * \brief The partial derivative of the horizontal velocity (following a parabolic profile for
255 * Stokes flow) w.r.t. to the y-coordinate (du/dy).
256 *
257 * \param y The position where the derivative is evaluated.
258 * \param vMax The profile's maximum velocity.
259 */
260 4090 Scalar dudy(const Scalar y, const Scalar vMax) const
261 {
262 4090 const Scalar yMin = this->gridGeometry().bBoxMin()[1];
263 4090 const Scalar yMax = this->gridGeometry().bBoxMax()[1];
264 4090 return vMax * (4.0*yMin + 4*yMax - 8.0*y) / ((yMin-yMax)*(yMin-yMax));
265 }
266
267 /*!
268 * \brief Return the analytical solution of the problem at a given position
269 *
270 * \param globalPos The global position
271 */
272 8750 DirichletValues analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
273 {
274 8750 DirichletValues values;
275
276 if constexpr (ParentType::isMomentumProblem())
277 {
278
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 7500 values[Indices::velocityXIdx] = parabolicProfile(globalPos[1], inletVelocity_);
279 values[Indices::velocityYIdx] = 0.0;
280 }
281 else
282 {
283 1250 const Scalar yMin = this->gridGeometry().bBoxMin()[1];
284 1250 const Scalar yMax = this->gridGeometry().bBoxMax()[1];
285 1250 const Scalar velocityQuadraticCoefficient = - inletVelocity_ / (0.25*(yMax - yMin)*(yMax - yMin));
286 1250 const Scalar pressureLinearCoefficient = 2.0 * velocityQuadraticCoefficient * dynamicViscosity_;
287 1250 const Scalar pressureConstant = -pressureLinearCoefficient * this->gridGeometry().bBoxMax()[0] + outletPressure_;
288 1250 values[Indices::pressureIdx] = pressureLinearCoefficient * globalPos[0] + pressureConstant;
289 }
290
291 1250 return values;
292 }
293
294 // \}
295
296 /*!
297 * \name Volume terms
298 */
299 // \{
300
301 /*!
302 * \brief Evaluates the initial value for a control volume.
303 *
304 * \param globalPos The global position
305 */
306 237774 InitialValues initialAtPos(const GlobalPosition& globalPos) const
307 {
308
4/4
✓ Branch 0 taken 121308 times.
✓ Branch 1 taken 2766 times.
✓ Branch 2 taken 5150 times.
✓ Branch 3 taken 71050 times.
 200274 InitialValues values;
309
310 if constexpr (ParentType::isMomentumProblem())
311 {
312 200274 values[Indices::velocityYIdx] = 0.0;
313
2/2
✓ Branch 0 taken 5150 times.
✓ Branch 1 taken 71050 times.
 76200 if (useVelocityProfile_)
314 5150 values[Indices::velocityXIdx] = parabolicProfile(globalPos[1], inletVelocity_);
315 else
316
2/2
✓ Branch 0 taken 121308 times.
✓ Branch 1 taken 2766 times.
 195124 values[Indices::velocityXIdx] = inletVelocity_;
317 }
318 else
319 {
320
3/4
✓ Branch 0 taken 2050 times.
✓ Branch 1 taken 5625 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10000 times.
 17675 values[Indices::pressureIdx] = outletPressure_;
321
2/2
✓ Branch 0 taken 5000 times.
✓ Branch 1 taken 32500 times.
 37500 if (enableGravity_)
322 5000 values[Indices::pressureIdx] += hydrostaticPressure_(globalPos);
323
324 #if NONISOTHERMAL
325 values[Indices::temperatureIdx] = 283.15;
326 #endif
327 }
328
329
330
331 return values;
332 }
333
334 // \}
335
336 /*!
337 * \brief Returns a reference pressure at a given sub control volume face.
338 * This pressure is subtracted from the actual pressure for the momentum balance
339 * which potentially helps to improve numerical accuracy by avoiding issues related do floating point arithmetic.
340 */
341 10502004 Scalar referencePressure(const Element& element,
342 const FVElementGeometry& fvGeometry,
343 const SubControlVolumeFace& scvf) const
344
2/2
✓ Branch 0 taken 7568 times.
✓ Branch 1 taken 14980 times.
 10502004 { return referencePressure_; }
345
346 30 void setTime(Scalar time)
347 {
348
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
 30 time_ = time;
349 }
350
351 3800 Scalar time() const
352 {
353
2/2
✓ Branch 0 taken 2050 times.
✓ Branch 1 taken 1750 times.
 3800 return time_;
354 }
355
356 10 bool hasAnalyticalSolution() const
357
3/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 10 { return outletCondition_ == OutletCondition::neumannXneumannY; }
358
359 private:
360
6/6
✓ Branch 0 taken 61892 times.
✓ Branch 1 taken 79520 times.
✓ Branch 2 taken 64248 times.
✓ Branch 3 taken 88748 times.
✓ Branch 4 taken 3875 times.
✓ Branch 5 taken 19375 times.
 317658 bool isInlet_(const GlobalPosition& globalPos) const
361 {
362
6/6
✓ Branch 0 taken 61892 times.
✓ Branch 1 taken 79520 times.
✓ Branch 2 taken 64248 times.
✓ Branch 3 taken 88748 times.
✓ Branch 4 taken 3875 times.
✓ Branch 5 taken 19375 times.
 317658 return globalPos[0] < eps_;
363 }
364
365
4/4
✓ Branch 0 taken 25050 times.
✓ Branch 1 taken 100200 times.
✓ Branch 2 taken 46028 times.
✓ Branch 3 taken 65882 times.
 237160 bool isOutlet_(const GlobalPosition& globalPos) const
366 {
367
4/4
✓ Branch 0 taken 25050 times.
✓ Branch 1 taken 100200 times.
✓ Branch 2 taken 46028 times.
✓ Branch 3 taken 65882 times.
 237160 return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_;
368 }
369
370 14296 Scalar hydrostaticPressure_(const GlobalPosition& globalPos) const
371 {
372
2/3
✓ Branch 0 taken 2050 times.
✓ Branch 1 taken 5625 times.
✗ Branch 2 not taken.
 21971 Scalar height = globalPos[dimWorld - 1] - this->gridGeometry().bBoxMax()[dimWorld - 1];
373
2/3
✓ Branch 0 taken 2050 times.
✓ Branch 1 taken 5625 times.
✗ Branch 2 not taken.
 21971 Scalar g = this->spatialParams().gravity(globalPos)[dimWorld - 1];
374
375
2/3
✓ Branch 0 taken 2050 times.
✓ Branch 1 taken 5625 times.
✗ Branch 2 not taken.
 21971 return outletDensity_ * g * height;
376 }
377
378 static constexpr Scalar eps_=1e-6;
379 Scalar inletVelocity_;
380 Scalar dynamicViscosity_;
381 Scalar outletPressure_;
382 OutletCondition outletCondition_;
383 bool useVelocityProfile_;
384 bool enableGravity_;
385 Scalar time_;
386 Scalar outletDensity_;
387 Scalar referencePressure_;
388 };
389 } // end namespace Dumux
390
391 #endif
392