GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/freeflow/navierstokes/channel/2d/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 82 97 84.5%
Functions: 7 20 35.0%
Branches: 103 206 50.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 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 14 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
56 static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
57 using Element = typename FVElementGeometry::Element;
58 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
59
60 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
61
62 // the types of outlet boundary conditions
63 enum class OutletCondition
64 {
65 outflow, unconstrainedOutflow, neumannXdirichletY, neumannXneumannY
66 };
67
68 public:
69 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
70
71 28 ChannelTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
72
7/20
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 14 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 14 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 14 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 14 times.
✓ Branch 18 taken 14 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
 112 : ParentType(gridGeometry, couplingManager)
73 {
74
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 28 inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
75
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 28 dynamicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0)
76
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 28 * getParam<Scalar>("Component.LiquidDensity", 1.0);
77
78
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 28 const auto outletBC = getParam<std::string>("Problem.OutletCondition", "Outflow");
79
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 6 times.
 28 if (outletBC == "Outflow")
80 16 outletCondition_ = OutletCondition::outflow;
81
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 4 times.
 12 else if (outletBC == "UnconstrainedOutflow")
82 4 outletCondition_ = OutletCondition::unconstrainedOutflow;
83
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 8 else if (outletBC == "NeumannX_DirichletY")
84 4 outletCondition_ = OutletCondition::neumannXdirichletY;
85
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 else if (outletBC == "NeumannX_NeumannY")
86 4 outletCondition_ = OutletCondition::neumannXneumannY;
87 else
88 ✗ DUNE_THROW(Dune::InvalidStateException, outletBC + " is not a valid outlet boundary condition");
89
90
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 28 useVelocityProfile_ = getParam<bool>("Problem.UseVelocityProfile", false);
91
3/6
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 28 outletPressure_ = getParam<Scalar>("Problem.OutletPressure", 1.1e5);
92 28 }
93
94 /*!
95 * \name Boundary conditions
96 */
97 // \{
98
99 /*!
100 * \brief Specifies which kind of boundary condition should be
101 * used for which equation on a given boundary control volume.
102 *
103 * \param globalPos The position of the center of the finite volume
104 */
105 196180 BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
106 {
107
4/6
✓ Branch 0 taken 9500 times.
✓ Branch 1 taken 47500 times.
✓ Branch 2 taken 2675 times.
✓ Branch 3 taken 13375 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 269230 BoundaryTypes values;
108
109 if constexpr (ParentType::isMomentumProblem())
110 {
111
2/2
✓ Branch 0 taken 43728 times.
✓ Branch 1 taken 54362 times.
 196180 values.setDirichlet(Indices::velocityXIdx);
112
2/2
✓ Branch 0 taken 43728 times.
✓ Branch 1 taken 54362 times.
 196180 values.setDirichlet(Indices::velocityYIdx);
113
114 392360 if (isOutlet_(globalPos))
115 {
116
2/2
✓ Branch 0 taken 42578 times.
✓ Branch 1 taken 1150 times.
 87456 if (outletCondition_ == OutletCondition::neumannXdirichletY)
117 {
118 2300 values.setNeumann(Indices::momentumXBalanceIdx);
119 2300 values.setDirichlet(Indices::velocityYIdx);
120 }
121 else
122 values.setAllNeumann();
123 }
124 }
125 else
126 {
127 73050 values.setAllNeumann();
128
129 146100 if (isInlet_(globalPos))
130 {
131 12175 values.setDirichlet(Indices::pressureIdx);
132 #if NONISOTHERMAL
133 ✗ values.setDirichlet(Indices::temperatureIdx);
134 #endif
135 }
136 }
137
138 196180 return values;
139 }
140
141
142 /*!
143 * \brief Evaluates the boundary conditions for a Dirichlet control volume.
144 *
145 * \param globalPos The center of the finite volume which ought to be set.
146 */
147 234108 DirichletValues dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
148 {
149
2/2
✓ Branch 0 taken 116338 times.
✓ Branch 1 taken 4516 times.
 240583 const auto& globalPos = scvf.ipGlobal();
150 481166 DirichletValues values = initialAtPos(globalPos);
151
152 if constexpr (ParentType::isMomentumProblem())
153 {
154 936432 values[Indices::velocityXIdx] = parabolicProfile(globalPos[1], inletVelocity_);
155
2/2
✓ Branch 0 taken 114288 times.
✓ Branch 1 taken 2766 times.
 234108 if (!useVelocityProfile_ && isInlet_(globalPos))
156 values[Indices::velocityXIdx] = inletVelocity_;
157 }
158 else
159 {
160
4/4
✓ Branch 0 taken 2050 times.
✓ Branch 1 taken 1750 times.
✓ Branch 2 taken 2050 times.
✓ Branch 3 taken 1750 times.
 12950 values[Indices::pressureIdx] = this->couplingManager().cellPressure(element, scvf);
161
162 #if NONISOTHERMAL
163 // give the system some time so that the pressure can equilibrate, then start the injection of the hot liquid
164
2/2
✓ Branch 0 taken 2050 times.
✓ Branch 1 taken 1750 times.
 3800 if (time() >= 200.0)
165 4100 values[Indices::temperatureIdx] = 293.15;
166 #endif
167 }
168
169 234108 return values;
170 }
171
172 /*!
173 * \brief Evaluates the boundary conditions for a Neumann control volume.
174 *
175 * \param element The element for which the Neumann boundary condition is set
176 * \param fvGeometry The fvGeometry
177 * \param elemVolVars The element volume variables
178 * \param elemFaceVars The element face variables
179 * \param scvf The boundary sub control volume face
180 */
181 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
182 181172 BoundaryFluxes neumann(const Element& element,
183 const FVElementGeometry& fvGeometry,
184 const ElementVolumeVariables& elemVolVars,
185 const ElementFluxVariablesCache& elemFluxVarsCache,
186 const SubControlVolumeFace& scvf) const
187 {
188
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 181172 BoundaryFluxes values(0.0);
189
190 if constexpr (ParentType::isMomentumProblem())
191 {
192 using FluxHelper = NavierStokesMomentumBoundaryFlux<typename GridGeometry::DiscretizationMethod>;
193
194
2/2
✓ Branch 0 taken 4596 times.
✓ Branch 1 taken 85990 times.
 181172 if (outletCondition_ == OutletCondition::unconstrainedOutflow)
195 9192 values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, outletPressure_, false /*zeroNormalVelocityGradient*/);
196
2/2
✓ Branch 0 taken 81900 times.
✓ Branch 1 taken 4090 times.
 171980 else if (outletCondition_ == OutletCondition::outflow)
197 163800 values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, outletPressure_, true /*zeroNormalVelocityGradient*/);
198 else
199 {
200
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4090 times.
 8180 assert(outletCondition_ == OutletCondition::neumannXneumannY || outletCondition_ == OutletCondition::neumannXdirichletY);
201 8180 values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, outletPressure_, false /*zeroNormalVelocityGradient*/);
202
203 8180 Dune::FieldMatrix<Scalar, dimWorld, dimWorld> shearRate(0.0); // gradV + gradV^T
204 24540 shearRate[0][1] = dudy(scvf.ipGlobal()[1], inletVelocity_); // we assume du/dx = dv/dy = dv/dx = 0
205 32720 shearRate[1][0] = shearRate[0][1];
206
207 8180 const auto normal = scvf.unitOuterNormal();
208 8180 BoundaryFluxes normalGradient(0.0);
209 8180 shearRate.mv(normal, normalGradient);
210 24540 values -= this->effectiveViscosity(element, fvGeometry, scvf)*normalGradient;
211 }
212 }
213 else
214 {
215 using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
216
217 ✗ if (isOutlet_(scvf.ipGlobal()))
218 ✗ values = FluxHelper::scalarOutflowFlux(*this, element, fvGeometry, scvf, elemVolVars);
219 }
220
221 181172 return values;
222 }
223
224 /*!
225 * \brief A parabolic velocity profile.
226 *
227 * \param y The position where the velocity is evaluated.
228 * \param vMax The profile's maximum velocity.
229 */
230 Scalar parabolicProfile(const Scalar y, const Scalar vMax) const
231 {
232
6/10
✓ Branch 0 taken 114288 times.
✓ Branch 1 taken 2766 times.
✓ Branch 2 taken 114288 times.
✓ Branch 3 taken 2766 times.
✓ Branch 4 taken 114288 times.
✓ Branch 5 taken 2766 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 376462 const Scalar yMin = this->gridGeometry().bBoxMin()[1];
233
6/12
✓ Branch 0 taken 114288 times.
✓ Branch 1 taken 2766 times.
✓ Branch 2 taken 114288 times.
✓ Branch 3 taken 2766 times.
✓ Branch 4 taken 114288 times.
✓ Branch 5 taken 2766 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 389112 const Scalar yMax = this->gridGeometry().bBoxMax()[1];
234
2/2
✓ Branch 0 taken 114288 times.
✓ Branch 1 taken 2766 times.
 117054 return vMax * (y - yMin)*(yMax - y) / (0.25*(yMax - yMin)*(yMax - yMin));
235 }
236
237 /*!
238 * \brief The partial derivative of the horizontal velocity (following a parabolic profile for
239 * Stokes flow) w.r.t. to the y-coordinate (du/dy).
240 *
241 * \param y The position where the derivative is evaluated.
242 * \param vMax The profile's maximum velocity.
243 */
244 Scalar dudy(const Scalar y, const Scalar vMax) const
245 {
246 12270 const Scalar yMin = this->gridGeometry().bBoxMin()[1];
247 12270 const Scalar yMax = this->gridGeometry().bBoxMax()[1];
248 4090 return vMax * (4.0*yMin + 4*yMax - 8.0*y) / ((yMin-yMax)*(yMin-yMax));
249 }
250
251 /*!
252 * \brief Return the analytical solution of the problem at a given position
253 *
254 * \param globalPos The global position
255 */
256 2500 DirichletValues analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
257 {
258 10000 DirichletValues values;
259
260 if constexpr (ParentType::isMomentumProblem())
261 {
262
0/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 22500 values[Indices::velocityXIdx] = parabolicProfile(globalPos[1], inletVelocity_);
263
0/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 15000 values[Indices::velocityYIdx] = 0.0;
264 }
265 else
266 {
267 10000 const Scalar yMin = this->gridGeometry().bBoxMin()[1];
268 10000 const Scalar yMax = this->gridGeometry().bBoxMax()[1];
269 2500 const Scalar velocityQuadraticCoefficient = - inletVelocity_ / (0.25*(yMax - yMin)*(yMax - yMin));
270 2500 const Scalar pressureLinearCoefficient = 2.0 * velocityQuadraticCoefficient * dynamicViscosity_;
271 10000 const Scalar pressureConstant = -pressureLinearCoefficient * this->gridGeometry().bBoxMax()[0] + outletPressure_;
272 5000 values[Indices::pressureIdx] = pressureLinearCoefficient * globalPos[0] + pressureConstant;
273 }
274
275 2500 return values;
276 }
277
278 // \}
279
280 /*!
281 * \name Volume terms
282 */
283 // \{
284
285 /*!
286 * \brief Evaluates the initial value for a control volume.
287 *
288 * \param globalPos The global position
289 */
290 ✗ InitialValues initialAtPos(const GlobalPosition& globalPos) const
291 {
292 206929 InitialValues values;
293
294 if constexpr (ParentType::isMomentumProblem())
295 {
296
4/4
✓ Branch 0 taken 114288 times.
✓ Branch 1 taken 2766 times.
✓ Branch 2 taken 5150 times.
✓ Branch 3 taken 50750 times.
 172954 values[Indices::velocityYIdx] = 0.0;
297
2/2
✓ Branch 0 taken 5150 times.
✓ Branch 1 taken 50750 times.
 55900 if (useVelocityProfile_)
298 20600 values[Indices::velocityXIdx] = parabolicProfile(globalPos[1], inletVelocity_);
299 else
300
4/4
✓ Branch 0 taken 114288 times.
✓ Branch 1 taken 2766 times.
✓ Branch 2 taken 114288 times.
✓ Branch 3 taken 2766 times.
 335608 values[Indices::velocityXIdx] = inletVelocity_;
301 }
302 else
303 {
304
2/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 2050 times.
✓ Branch 3 taken 1750 times.
 33975 values[Indices::pressureIdx] = outletPressure_;
305 #if NONISOTHERMAL
306
4/8
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2050 times.
✓ Branch 5 taken 1750 times.
✓ Branch 6 taken 2050 times.
✓ Branch 7 taken 1750 times.
 27600 values[Indices::temperatureIdx] = 283.15;
307 #endif
308 }
309
310
311
312 ✗ return values;
313 }
314
315 // \}
316
317 /*!
318 * \brief Returns a reference pressure at a given sub control volume face.
319 * This pressure is subtracted from the actual pressure for the momentum balance
320 * which potentially helps to improve numerical accuracy by avoiding issues related do floating point arithmetic.
321 */
322 ✗ Scalar referencePressure(const Element& element,
323 const FVElementGeometry& fvGeometry,
324 const SubControlVolumeFace& scvf) const
325 ✗ { return outletPressure_; }
326
327 ✗ void setTime(Scalar time)
328 {
329
4/8
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 17 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 17 times.
✗ Branch 11 not taken.
 48 time_ = time;
330 ✗ }
331
332 ✗ Scalar time() const
333 {
334 ✗ return time_;
335 }
336
337 ✗ bool hasAnalyticalSolution() const
338 ✗ { return outletCondition_ == OutletCondition::neumannXneumannY; }
339
340 private:
341 ✗ bool isInlet_(const GlobalPosition& globalPos) const
342 {
343
12/20
✓ Branch 0 taken 21520 times.
✓ Branch 1 taken 95892 times.
✓ Branch 2 taken 21520 times.
✓ Branch 3 taken 95892 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 10578 times.
✓ Branch 9 taken 43298 times.
✓ Branch 10 taken 10578 times.
✓ Branch 11 taken 43298 times.
✓ Branch 12 taken 2675 times.
✓ Branch 13 taken 13375 times.
✓ Branch 14 taken 2675 times.
✓ Branch 15 taken 13375 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 374676 return globalPos[0] < eps_;
344 }
345
346 bool isOutlet_(const GlobalPosition& globalPos) const
347 {
348
10/20
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 43728 times.
✓ Branch 11 taken 54362 times.
✓ Branch 12 taken 43728 times.
✓ Branch 13 taken 54362 times.
✓ Branch 14 taken 43728 times.
✓ Branch 15 taken 54362 times.
✓ Branch 16 taken 43728 times.
✓ Branch 17 taken 54362 times.
✓ Branch 18 taken 43728 times.
✓ Branch 19 taken 54362 times.
 490450 return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_;
349 }
350
351 static constexpr Scalar eps_=1e-6;
352 Scalar inletVelocity_;
353 Scalar dynamicViscosity_;
354 Scalar outletPressure_;
355 OutletCondition outletCondition_;
356 bool useVelocityProfile_;
357 Scalar time_;
358 };
359 } // end namespace Dumux
360
361 #endif
362