GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/freeflow/navierstokes/unstructured/problem.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 86 92 93.5%
Functions: 21 33 63.6%
Branches: 94 172 54.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-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 #ifndef DUMUX_TEST_FREEFLOW_NAVIERSTOKES_3D_CHANNEL_PROBLEM_HH
8 #define DUMUX_TEST_FREEFLOW_NAVIERSTOKES_3D_CHANNEL_PROBLEM_HH
9
10 #include <cmath>
11 #include <numeric>
12
13 #include <dune/common/float_cmp.hh>
14 #include <dune/common/fmatrix.hh>
15 #include <dune/geometry/quadraturerules.hh>
16
17 #include <dumux/common/math.hh>
18 #include <dumux/common/properties.hh>
19 #include <dumux/common/parameters.hh>
20 #include <dumux/io/vtk/function.hh>
21 #include <dumux/io/grid/griddata.hh>
22 #include <dumux/discretization/evalsolution.hh>
23 #include <dumux/discretization/evalgradients.hh>
24 #include <dumux/discretization/elementsolution.hh>
25
26 namespace Dumux {
27
28 /*!
29 * \brief Test problem for the (Navier-) Stokes model in a 3D channel
30 *
31 * Benchmark case from
32 * Turek, Schaefer et a (1996) Benchmark computations of laminar flow around cylinder.
33 * Flow Simulation with High-Performance Computers II,
34 * Notes on Numerical Fluid Mechanics 52, 547-566, Vieweg 1996
35 * https://doi.org/10.1007/978-3-322-89849-4_39
36 */
37 template <class TypeTag, class BaseProblem>
38 12 class DFGChannelTestProblem : public BaseProblem
39 {
40 using ParentType = BaseProblem;
41
42 using BoundaryTypes = typename ParentType::BoundaryTypes;
43 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
44 using FVElementGeometry = typename GridGeometry::LocalView;
45 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
46 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
47 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
48 using BoundaryFluxes = typename ParentType::BoundaryFluxes;
49 using DirichletValues = typename ParentType::DirichletValues;
50 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
51 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
52 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
53
54 static constexpr int dim = GridGeometry::GridView::dimension;
55 static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
56
57 using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
58 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
59 using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
60
61 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
62
63 public:
64 24 DFGChannelTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
65
7/20
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 12 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 12 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 12 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 12 times.
✓ Branch 18 taken 12 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.
 96 : ParentType(gridGeometry, couplingManager, ParentType::isMomentumProblem() ? "Momentum" : "Mass")
66 {
67
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 24 density_ = getParam<Scalar>("Component.LiquidDensity");
68
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 24 viscosity_ = getParam<Scalar>("Component.LiquidDynamicViscosity");
69
70
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 12 times.
 72 for (int i = 0; i < dimWorld; ++i)
71 384 domainSize_[i] = gridGeometry->bBoxMax()[i] - gridGeometry->bBoxMin()[i];
72 24 }
73
74 /*!
75 * \brief Specifies which kind of boundary condition should be
76 * used for which equation on a given boundary control volume.
77 *
78 * \param globalPos The position of the center of the finite volume
79 */
80 104146 BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
81 {
82
3/4
✓ Branch 0 taken 32772 times.
✓ Branch 1 taken 49104 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 16368 times.
 202390 BoundaryTypes values;
83
84 if constexpr (ParentType::isMomentumProblem())
85 {
86 208292 if (isOutlet_(globalPos))
87 values.setAllNeumann();
88 else
89 values.setAllDirichlet();
90 }
91 else
92
6/8
✓ Branch 0 taken 32772 times.
✓ Branch 1 taken 49104 times.
✓ Branch 2 taken 32772 times.
✓ Branch 3 taken 49104 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 16368 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 16368 times.
 196488 values.setNeumann(Indices::conti0EqIdx);
93
94 104146 return values;
95 }
96
97 /*!
98 * \brief Evaluates the boundary conditions for a Dirichlet boundary face
99 */
100 46284 DirichletValues dirichletAtPos(const GlobalPosition& globalPos) const
101 {
102
10/40
✓ Branch 0 taken 3000 times.
✓ Branch 1 taken 43284 times.
✓ Branch 2 taken 3000 times.
✓ Branch 3 taken 43284 times.
✓ Branch 4 taken 3000 times.
✓ Branch 5 taken 43284 times.
✓ Branch 6 taken 3000 times.
✓ Branch 7 taken 43284 times.
✓ Branch 8 taken 3000 times.
✓ Branch 9 taken 43284 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ 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.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
 231420 if (globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_)
103 3000 return inflowVelocity_(globalPos[1]);
104 else
105 86568 return DirichletValues(0.0);
106 }
107
108 /*!
109 * \brief Evaluates the boundary conditions for a Neumann boundary face
110 *
111 * \param element The element for which the Neumann boundary condition is set
112 * \param fvGeometry The fvGeometry
113 * \param elemVolVars The element volume variables
114 * \param elemFaceVars The element face variables
115 * \param scvf The boundary sub control volume face
116 */
117 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
118 214056 BoundaryFluxes neumann(const Element& element,
119 const FVElementGeometry& fvGeometry,
120 const ElementVolumeVariables& elemVolVars,
121 const ElementFluxVariablesCache& elemFluxVarsCache,
122 const SubControlVolumeFace& scvf) const
123 {
124
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 214056 BoundaryFluxes values(0.0);
125 if constexpr (ParentType::isMomentumProblem())
126 {
127
2/2
✓ Branch 0 taken 17370 times.
✓ Branch 1 taken 3474 times.
 41688 if (this->enableInertiaTerms())
128 {
129
2/2
✓ Branch 0 taken 17010 times.
✓ Branch 1 taken 360 times.
 34740 if (isOutlet_(scvf.ipGlobal()))
130 {
131 // advective term: vv*n
132 34020 const auto elemSol = elementSolution(element, elemVolVars, fvGeometry);
133 136080 const auto v = evalSolution(element, element.geometry(), fvGeometry.gridGeometry(), elemSol, scvf.ipGlobal());
134 102060 values.axpy(density_*(v*scvf.unitOuterNormal()), v);
135 }
136 }
137 }
138 else
139 {
140 172368 values[Indices::conti0EqIdx]
141 = this->faceVelocity(element, fvGeometry, scvf)
142
0/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 344736 * density_ * scvf.unitOuterNormal();
143 }
144
145 214056 return values;
146 }
147
148 template<class ScvOrScvf>
149 Scalar density(const Element&,
150 const FVElementGeometry&,
151 const ScvOrScvf&,
152 const bool isPreviousTimeStep = false) const
153 {
154 53920312 return density_;
155 }
156
157 template<class ScvOrScvf>
158 ✗ Scalar effectiveViscosity(const Element&,
159 const FVElementGeometry&,
160 const ScvOrScvf&) const
161 {
162 ✗ return viscosity_;
163 }
164
165 //! Computes pressure difference benchmark indicator
166 //! (only works for correct domain size 2.2 x 0.41 and boundary conditions, Re=20)
167 template<class SolutionVector, class GridVariables, class P = ParentType, typename std::enable_if_t<!P::isMomentumProblem(), int> = 0>
168 3 Scalar evalPressureDifference(const GridVariables& gridVariables, const SolutionVector& p) const
169 {
170 3 const auto& gg = gridVariables.gridGeometry();
171 3 const auto& tree = gg.boundingBoxTree();
172 3 GlobalPosition evalPoint1({0.15, 0.2});
173 3 GlobalPosition evalPoint2({0.25, 0.2});
174 3 const Scalar stepSize = 0.001;
175
176 9 const auto evalPressure = [&](auto pos, bool backwards)
177 {
178 6 auto entities = intersectingEntities(pos, tree);
179
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
 12 while (entities.empty())
180 {
181 ✗ pos[0] += backwards ? -stepSize : stepSize;
182 ✗ entities = intersectingEntities(pos, tree);
183 }
184
185
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 const auto element = gg.element(entities[0]);
186
2/6
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 12 const auto fvGeometry = localView(gg).bindElement(element);
187
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 22 const auto elemVolVars = localView(gridVariables.curGridVolVars()).bindElement(element, fvGeometry, p);
188 6 const auto elemSol = elementSolution(element, elemVolVars, fvGeometry);
189
4/10
✓ 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 9 taken 4 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 24 return evalSolution(element, element.geometry(), fvGeometry.gridGeometry(), elemSol, pos);
190 };
191
192 3 return evalPressure(evalPoint1, true) - evalPressure(evalPoint2, false);
193 }
194
195 //! Computes drag and lift coefficient benchmark indicator
196 //! (only works for correct domain size 2.2 x 0.41 and boundary conditions, Re=20)
197 template<class SolutionVector, class GridVariables, class P = ParentType, typename std::enable_if_t<P::isMomentumProblem(), int> = 0>
198 3 auto evalDragAndLiftCoefficient(const GridVariables& gridVariables, const SolutionVector& v) const
199 {
200 6 auto fvGeometry = localView(this->gridGeometry());
201 6 auto elemVolVars = localView(gridVariables.curGridVolVars());
202
203 3 BoundaryFluxes forces(0.0);
204 3 Scalar cylinderSurface = 0.0;
205
4/4
✓ Branch 3 taken 71722 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 71722 times.
✓ Branch 6 taken 3 times.
 71731 for (const auto& element : elements(this->gridGeometry().gridView()))
206 {
207 71722 fvGeometry.bind(element);
208
4/4
✓ Branch 0 taken 2044 times.
✓ Branch 1 taken 69678 times.
✓ Branch 2 taken 2044 times.
✓ Branch 3 taken 69678 times.
 143444 if (fvGeometry.hasBoundaryScvf())
209 {
210
4/4
✓ Branch 0 taken 14992 times.
✓ Branch 1 taken 2044 times.
✓ Branch 2 taken 14992 times.
✓ Branch 3 taken 2044 times.
 19080 for (const auto& scvf : scvfs(fvGeometry))
211 {
212
4/4
✓ Branch 0 taken 3414 times.
✓ Branch 1 taken 11578 times.
✓ Branch 2 taken 800 times.
✓ Branch 3 taken 2614 times.
 14992 if (scvf.boundary() && onCylinderBoundary_(scvf.ipGlobal()))
213 {
214 800 elemVolVars.bind(element, fvGeometry, v);
215 1600 Dune::FieldVector<GlobalPosition, 2> stress(GlobalPosition(0.0));
216
217 800 const auto geometry = element.geometry();
218 800 const auto isgeometry = fvGeometry.geometry(scvf);
219 800 const auto& quad = Dune::QuadratureRules<Scalar, std::decay_t<decltype(isgeometry)>::mydimension>::rule(isgeometry.type(), 5);
220
4/4
✓ Branch 0 taken 2400 times.
✓ Branch 1 taken 800 times.
✓ Branch 2 taken 2400 times.
✓ Branch 3 taken 800 times.
 7200 for (auto&& qp : quad)
221 {
222 4800 const auto ipGlobal = isgeometry.global(qp.position());
223 2400 const auto elemSol = elementSolution(element, elemVolVars, fvGeometry);
224 4800 auto stressLocal = evalGradients(element, geometry, fvGeometry.gridGeometry(), elemSol, ipGlobal);
225 4800 stressLocal *= qp.weight()*isgeometry.integrationElement(qp.position());
226 4800 stress = stress + stressLocal;
227 }
228 800 stress *= viscosity_;
229 800 stress /= scvf.area();
230
231 800 BoundaryFluxes normalStress(0.0);
232
2/2
✓ Branch 0 taken 1600 times.
✓ Branch 1 taken 800 times.
 2400 for (int dir = 0; dir < dim; ++dir)
233 {
234 // Add pressure contribuition
235 1600 stress[dir][dir] -= this->pressure(element, fvGeometry, scvf);
236 6400 normalStress[dir] = stress[dir]*scvf.unitOuterNormal();
237 }
238 800 forces.axpy(scvf.area(), normalStress);
239
240 800 cylinderSurface += scvf.area();
241 }
242 }
243 }
244 }
245
246 3 const Scalar dragForce = -forces[0];
247 3 const Scalar liftForce = -forces[1];
248
249 3 const Scalar uMean = 0.2;
250 3 const Scalar lChar = 0.1;
251 3 const Scalar coefficientFactor = 2.0/(uMean*uMean*lChar);
252
253 // sanity check
254 6 std::cout << "Reynolds number: " << uMean*lChar/viscosity_ << std::endl;
255 6 std::cout << "Cylinder surface: " << cylinderSurface
256 3 << " (reference: 0.31415926535)"
257 3 << std::endl;
258
259 6 return std::make_pair( coefficientFactor*dragForce, coefficientFactor*liftForce );
260 }
261
262 private:
263 bool isInlet_(const GlobalPosition& globalPos) const
264 { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
265
266 bool isOutlet_(const GlobalPosition& globalPos) const
267
20/20
✓ Branch 0 taken 17010 times.
✓ Branch 1 taken 360 times.
✓ Branch 2 taken 17010 times.
✓ Branch 3 taken 360 times.
✓ Branch 4 taken 17010 times.
✓ Branch 5 taken 360 times.
✓ Branch 6 taken 17010 times.
✓ Branch 7 taken 360 times.
✓ Branch 8 taken 17010 times.
✓ Branch 9 taken 360 times.
✓ Branch 10 taken 48847 times.
✓ Branch 11 taken 3226 times.
✓ Branch 12 taken 48847 times.
✓ Branch 13 taken 3226 times.
✓ Branch 14 taken 48847 times.
✓ Branch 15 taken 3226 times.
✓ Branch 16 taken 48847 times.
✓ Branch 17 taken 3226 times.
✓ Branch 18 taken 48847 times.
✓ Branch 19 taken 3226 times.
 347215 { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
268
269 ✗ bool onCylinderBoundary_(const GlobalPosition& globalPos) const
270
6/6
✓ Branch 0 taken 800 times.
✓ Branch 1 taken 2614 times.
✓ Branch 2 taken 800 times.
✓ Branch 3 taken 2614 times.
✓ Branch 4 taken 800 times.
✓ Branch 5 taken 2614 times.
 10242 { return std::hypot(globalPos[0] - 0.2, globalPos[1] - 0.2) < 0.06; }
271
272 ✗ DirichletValues inflowVelocity_(Scalar y) const
273 {
274 3000 constexpr Scalar maxVelocity = 0.3;
275 15000 return { 4*maxVelocity*y*(domainSize_[1]-y)/(domainSize_[1]*domainSize_[1]), 0.0 };
276 }
277
278 static constexpr Scalar eps_ = 1e-10;
279 Scalar density_, viscosity_;
280 std::array<Scalar, dimWorld> domainSize_;
281 };
282
283 } // end namespace Dumux
284
285 #endif
286