GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/freeflow/navierstokes/sincos/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 89 137 65.0%
Functions: 26 47 55.3%
Branches: 83 190 43.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 /*!
8 * \file
9 * \ingroup NavierStokesTests
10 * \brief Test for the staggered grid Navier-Stokes model with analytical solution.
11 */
12 #ifndef DUMUX_SINCOS_TEST_PROBLEM_HH
13 #define DUMUX_SINCOS_TEST_PROBLEM_HH
14
15 #include <dune/common/fmatrix.hh>
16
17 #include <dumux/common/parameters.hh>
18 #include <dumux/common/properties.hh>
19
20 #include <dumux/freeflow/navierstokes/boundarytypes.hh>
21
22 namespace Dumux {
23
24 /*!
25 * \ingroup NavierStokesTests
26 * \brief Test problem for the staggered grid.
27 *
28 * The 2D, incompressible Navier-Stokes equations for zero gravity and a Newtonian
29 * flow is solved and compared to an analytical solution (sums/products of trigonometric functions).
30 * For the instationary case, the velocities and pressures are periodical in time. The Dirichlet boundary conditions are
31 * consistent with the analytical solution and in the instationary case time-dependent.
32 */
33 template <class TypeTag, class BaseProblem>
34 8 class SincosTestProblem : public BaseProblem
35 {
36 using ParentType = BaseProblem;
37
38 using BoundaryTypes = typename ParentType::BoundaryTypes;
39 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
40 using InitialValues = typename ParentType::InitialValues;
41 using Sources = typename ParentType::Sources;
42 using DirichletValues = typename ParentType::DirichletValues;
43 using BoundaryFluxes = typename ParentType::BoundaryFluxes;
44 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
45 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
46 using SubControlVolume = typename GridGeometry::SubControlVolume;
47 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
48 using FVElementGeometry = typename GridGeometry::LocalView;
49
50 static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
51 using Element = typename FVElementGeometry::Element;
52 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
53 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
54
55 public:
56 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
57 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
58
59 16 SincosTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
60
7/20
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 8 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 8 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 8 times.
✓ Branch 18 taken 8 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.
 64 : ParentType(gridGeometry, couplingManager), time_(0.0)
61 {
62
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 isStationary_ = getParam<bool>("Problem.IsStationary");
63
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 rho_ = getParam<Scalar>("Component.LiquidDensity");
64
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 Scalar nu = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
65 16 mu_ = rho_*nu; // dynamic viscosity
66
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 useNeumann_ = getParam<bool>("Problem.UseNeumann", false);
67 16 }
68
69 /*!
70 * \brief Return the sources within the domain.
71 *
72 * \param globalPos The global position
73 */
74 67694200 Sources sourceAtPos(const GlobalPosition &globalPos) const
75 {
76
2/8
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 7755000 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 7755000 times.
 83204200 Sources source(0.0);
77 if constexpr (ParentType::isMomentumProblem())
78 {
79
1/2
✓ Branch 0 taken 33847100 times.
✗ Branch 1 not taken.
 67694200 const Scalar x = globalPos[0];
80
1/2
✓ Branch 0 taken 33847100 times.
✗ Branch 1 not taken.
 67694200 const Scalar y = globalPos[1];
81 67694200 const Scalar t = time_;
82
83
1/2
✓ Branch 0 taken 33847100 times.
✗ Branch 1 not taken.
 67694200 source[Indices::momentumXBalanceIdx] = rho_*dtU_(x,y,t) - 2.0*mu_*dxxU_(x,y,t) - mu_*dyyU_(x,y,t) - mu_*dxyV_(x,y,t) + dxP_(x,y,t);
84
1/2
✓ Branch 0 taken 33847100 times.
✗ Branch 1 not taken.
 67694200 source[Indices::momentumYBalanceIdx] = rho_*dtV_(x,y,t) - 2.0*mu_*dyyV_(x,y,t) - mu_*dxyU_(x,y,t) - mu_*dxxV_(x,y,t) + dyP_(x,y,t);
85
86
1/2
✓ Branch 0 taken 33847100 times.
✗ Branch 1 not taken.
 67694200 if (this->enableInertiaTerms())
87 {
88 67694200 source[Indices::momentumXBalanceIdx] += rho_*dxUU_(x,y,t) + rho_*dyUV_(x,y,t);
89 135388400 source[Indices::momentumYBalanceIdx] += rho_*dxUV_(x,y,t) + rho_*dyVV_(x,y,t);
90 }
91 }
92
93 67694200 return source;
94 }
95
96 // \}
97 /*!
98 * \name Boundary conditions
99 */
100 // \{
101
102 /*!
103 * \brief Specifies which kind of boundary condition should be
104 * used for which equation on a given boundary control volume.
105 *
106 * \param globalPos The position of the center of the finite volume
107 */
108 560800 BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
109 {
110
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 303600 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 123400 times.
 987800 BoundaryTypes values;
111
112 if constexpr (ParentType::isMomentumProblem())
113 {
114
1/2
✓ Branch 0 taken 280400 times.
✗ Branch 1 not taken.
 560800 if (useNeumann_)
115 {
116 ✗ if (globalPos[1] > this->gridGeometry().bBoxMax()[1] - 1e-6
117 ✗ || globalPos[0] > this->gridGeometry().bBoxMax()[0] - 1e-6
118 ✗ || globalPos[1] < this->gridGeometry().bBoxMin()[1] + 1e-6)
119 {
120 ✗ values.setNeumann(Indices::velocityXIdx);
121 ✗ values.setNeumann(Indices::velocityYIdx);
122 }
123 else
124 {
125 ✗ values.setDirichlet(Indices::velocityXIdx);
126 ✗ values.setDirichlet(Indices::velocityYIdx);
127 }
128 }
129 else
130 {
131 // set Dirichlet values for the velocity everywhere
132 values.setAllDirichlet();
133 }
134 }
135 else
136
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 303600 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 123400 times.
 427000 values.setAllNeumann();
137
138 560800 return values;
139 }
140
141 //! Enable internal Dirichlet constraints
142 static constexpr bool enableInternalDirichletConstraints()
143 { return !ParentType::isMomentumProblem(); }
144
145 /*!
146 * \brief Tag a degree of freedom to carry internal Dirichlet constraints.
147 * If true is returned for a dof, the equation for this dof is replaced
148 * by the constraint that its primary variable values must match the
149 * user-defined values obtained from the function internalDirichlet(),
150 * which must be defined in the problem.
151 *
152 * \param element The finite element
153 * \param scv The sub-control volume
154 */
155 ✗ std::bitset<DirichletValues::dimension> hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
156 {
157 // set fixed pressure in one cell
158 7022400 std::bitset<DirichletValues::dimension> values;
159
160
20/30
✓ Branch 0 taken 2820000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 440 times.
✓ Branch 3 taken 2819560 times.
✓ Branch 4 taken 440 times.
✓ Branch 5 taken 2819560 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.
✓ Branch 12 taken 705000 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 110 times.
✓ Branch 15 taken 704890 times.
✓ Branch 16 taken 110 times.
✓ Branch 17 taken 704890 times.
✓ Branch 18 taken 2792400 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 220 times.
✓ Branch 21 taken 2792180 times.
✓ Branch 22 taken 220 times.
✓ Branch 23 taken 2792180 times.
✓ Branch 24 taken 705000 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 110 times.
✓ Branch 27 taken 704890 times.
✓ Branch 28 taken 110 times.
✓ Branch 29 taken 704890 times.
 7022400 if (!useNeumann_ && scv.dofIndex() == 0)
161 values.set(Indices::pressureIdx);
162
163 ✗ return values;
164 }
165
166 /*!
167 * \brief Define the values of internal Dirichlet constraints for a degree of freedom.
168 * \param element The finite element
169 * \param scv The sub-control volume
170 */
171 ✗ DirichletValues internalDirichlet(const Element& element, const SubControlVolume& scv) const
172 705110 { return DirichletValues(analyticalSolution(scv.center())[Indices::pressureIdx]); }
173
174 /*!
175 * \brief Returns Dirichlet boundary values at a given position.
176 *
177 * \param globalPos The global position
178 */
179 DirichletValues dirichletAtPos(const GlobalPosition& globalPos) const
180 {
181 // use the values of the analytical solution
182 ✗ return analyticalSolution(globalPos, time_);
183 }
184
185 /*!
186 * \brief Evaluates the boundary conditions for a Neumann control volume.
187 *
188 * \param element The element for which the Neumann boundary condition is set
189 * \param fvGeometry The fvGeometry
190 * \param elemVolVars The element volume variables
191 * \param elemFaceVars The element face variables
192 * \param scvf The boundary sub control volume face
193 */
194 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
195 ✗ BoundaryFluxes neumann(const Element& element,
196 const FVElementGeometry& fvGeometry,
197 const ElementVolumeVariables& elemVolVars,
198 const ElementFluxVariablesCache& elemFluxVarsCache,
199 const SubControlVolumeFace& scvf) const
200 {
201 ✗ BoundaryFluxes values(0.0);
202
203 if constexpr (ParentType::isMomentumProblem())
204 {
205 ✗ const auto flux = [&](Scalar x, Scalar y)
206 {
207 ✗ Dune::FieldMatrix<Scalar, dimWorld, dimWorld> momentumFlux(0.0);
208 ✗ const Scalar t = time_;
209
210 ✗ momentumFlux[0][0] = -2*mu_*dxU_(x,y,t) + p_(x,y,t);
211 ✗ momentumFlux[0][1] = -mu_*(dyU_(x,y,t) + dxV_(x,y,t));
212 ✗ momentumFlux[1][0] = -mu_*(dyU_(x,y,t) + dxV_(x,y,t));
213 ✗ momentumFlux[1][1] = -2*mu_*dyV_(x,y,t) + p_(x,y,t);
214
215 ✗ if (this->enableInertiaTerms())
216 {
217 ✗ momentumFlux[0][0] += rho_*u_(x,y,t)*u_(x,y,t);
218 ✗ momentumFlux[0][1] += rho_*u_(x,y,t)*v_(x,y,t);
219 ✗ momentumFlux[1][0] += rho_*v_(x,y,t)*u_(x,y,t);
220 ✗ momentumFlux[1][1] += rho_*v_(x,y,t)*v_(x,y,t);
221 }
222
223 ✗ return momentumFlux;
224 };
225
226 ✗ flux(scvf.ipGlobal()[0], scvf.ipGlobal()[1]).mv(scvf.unitOuterNormal(), values);
227 }
228 else
229 {
230 ✗ const auto insideDensity = elemVolVars[scvf.insideScvIdx()].density();
231 ✗ values[Indices::conti0EqIdx] = insideDensity * (this->faceVelocity(element, fvGeometry, scvf) * scvf.unitOuterNormal());
232 }
233
234 ✗ return values;
235 }
236
237 /*!
238 * \brief Returns the analytical solution of the problem at a given time and position.
239 *
240 * \param globalPos The global position
241 * \param time The current simulation time
242 */
243 4397720 DirichletValues analyticalSolution(const GlobalPosition& globalPos, const Scalar time) const
244 {
245
2/2
✓ Branch 1 taken 2049 times.
✓ Branch 2 taken 430451 times.
 747500 const Scalar x = globalPos[0];
246
2/2
✓ Branch 1 taken 2049 times.
✓ Branch 2 taken 430451 times.
 1457610 const Scalar y = globalPos[1];
247 5855330 const Scalar t = time;
248 5855330 DirichletValues values;
249
250 if constexpr (ParentType::isMomentumProblem())
251 {
252 8795440 values[Indices::velocityXIdx] = u_(x,y,t);
253 8795440 values[Indices::velocityYIdx] = v_(x,y,t);
254 }
255 else
256
2/2
✓ Branch 1 taken 2049 times.
✓ Branch 2 taken 430451 times.
 752500 values[Indices::pressureIdx] = p_(x,y,t);
257
258 4397720 return values;
259 }
260
261 // \}
262
263 /*!
264 * \name Volume terms
265 */
266 // \{
267
268 /*!
269 * \brief Evaluates the initial value for a control volume.
270 *
271 * \param globalPos The global position
272 */
273 707200 InitialValues initialAtPos(const GlobalPosition &globalPos) const
274 {
275
2/2
✓ Branch 0 taken 338400 times.
✓ Branch 1 taken 15200 times.
 707200 if (isStationary_)
276 1128600 return InitialValues(0.0);
277 else
278 30400 return analyticalSolution(globalPos, 0.0);
279 }
280
281
282 /*!
283 * \brief Returns the analytical solution of the problem at a given position.
284 *
285 * \param globalPos The global position
286 */
287 DirichletValues analyticalSolution(const GlobalPosition& globalPos) const
288 {
289 1410220 return analyticalSolution(globalPos, time_);
290 }
291
292 /*!
293 * \brief Updates the time
294 */
295 ✗ void updateTime(const Scalar time)
296 {
297
2/4
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 34 times.
✗ Branch 5 not taken.
 68 time_ = time;
298 ✗ }
299
300 private:
301 ✗ Scalar f_(Scalar t) const
302 {
303 using std::sin;
304
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 619500850 if (isStationary_)
305 ✗ return 1.0;
306 else
307 222372750 return sin(2.0 * t);
308 }
309
310 ✗ Scalar df_(Scalar t) const
311 {
312 using std::cos;
313 ✗ if (isStationary_)
314 ✗ return 0.0;
315 else
316 24326800 return 2.0 * cos(2.0 * t);
317 }
318
319 ✗ Scalar f1_(Scalar x) const
320 2915220 { using std::cos; return -0.25 * cos(2.0 * x); }
321
322 ✗ Scalar df1_(Scalar x) const
323 ✗ { using std::sin; return 0.5 * sin(2.0 * x); }
324
325 ✗ Scalar f2_(Scalar x) const
326 414960640 { using std::cos; return -cos(x); }
327
328 ✗ Scalar df2_(Scalar x) const
329 685737440 { using std::sin; return sin(x); }
330
331 ✗ Scalar ddf2_(Scalar x) const
332 203082600 { using std::cos; return cos(x); }
333
334 ✗ Scalar dddf2_(Scalar x) const
335 67694200 { using std::sin; return -sin(x); }
336
337 2915220 Scalar p_(Scalar x, Scalar y, Scalar t) const
338
2/2
✓ Branch 0 taken 445102 times.
✓ Branch 1 taken 1012508 times.
 2915220 { return (f1_(x) + f1_(y)) * f_(t) * f_(t); }
339
340 ✗ Scalar dxP_ (Scalar x, Scalar y, Scalar t) const
341 ✗ { return df1_(x) * f_(t) * f_(t); }
342
343 ✗ Scalar dyP_ (Scalar x, Scalar y, Scalar t) const
344 ✗ { return df1_(y) * f_(t) * f_(t); }
345
346 105939020 Scalar u_(Scalar x, Scalar y, Scalar t) const
347
2/2
✓ Branch 0 taken 37983424 times.
✓ Branch 1 taken 67955596 times.
 105939020 { return f2_(x)*df2_(y) * f_(t); }
348
349 33847100 Scalar dtU_ (Scalar x, Scalar y, Scalar t) const
350
2/2
✓ Branch 0 taken 12163400 times.
✓ Branch 1 taken 21683700 times.
 33847100 { return f2_(x)*df2_(y) * df_(t); }
351
352 67694200 Scalar dxU_ (Scalar x, Scalar y, Scalar t) const
353
2/2
✓ Branch 0 taken 24326800 times.
✓ Branch 1 taken 43367400 times.
 67694200 { return df2_(x)*df2_(y) * f_(t); }
354
355 33847100 Scalar dyU_ (Scalar x, Scalar y, Scalar t) const
356
2/2
✓ Branch 0 taken 12163400 times.
✓ Branch 1 taken 21683700 times.
 33847100 { return f2_(x)*ddf2_(y) * f_(t); }
357
358 33847100 Scalar dxxU_ (Scalar x, Scalar y, Scalar t) const
359
2/2
✓ Branch 0 taken 12163400 times.
✓ Branch 1 taken 21683700 times.
 33847100 { return ddf2_(x)*df2_(y) * f_(t); }
360
361 33847100 Scalar dxyU_ (Scalar x, Scalar y, Scalar t) const
362
2/2
✓ Branch 0 taken 12163400 times.
✓ Branch 1 taken 21683700 times.
 33847100 { return df2_(x)*ddf2_(y) * f_(t); }
363
364 33847100 Scalar dyyU_ (Scalar x, Scalar y, Scalar t) const
365
2/2
✓ Branch 0 taken 12163400 times.
✓ Branch 1 taken 21683700 times.
 33847100 { return f2_(x)*dddf2_(y) * f_(t); }
366
367 105939020 Scalar v_(Scalar x, Scalar y, Scalar t) const
368
2/2
✓ Branch 0 taken 37983424 times.
✓ Branch 1 taken 67955596 times.
 105939020 { return -f2_(y)*df2_(x) * f_(t); }
369
370 33847100 Scalar dtV_ (Scalar x, Scalar y, Scalar t) const
371
2/2
✓ Branch 0 taken 12163400 times.
✓ Branch 1 taken 21683700 times.
 33847100 { return -f2_(y)*df2_(x) * df_(t); }
372
373 33847100 Scalar dxV_ (Scalar x, Scalar y, Scalar t) const
374
2/2
✓ Branch 0 taken 12163400 times.
✓ Branch 1 taken 21683700 times.
 33847100 { return -f2_(y)*ddf2_(x) * f_(t); }
375
376 67694200 Scalar dyV_ (Scalar x, Scalar y, Scalar t) const
377
2/2
✓ Branch 0 taken 24326800 times.
✓ Branch 1 taken 43367400 times.
 67694200 { return -df2_(y)*df2_(x) * f_(t); }
378
379 33847100 Scalar dyyV_ (Scalar x, Scalar y, Scalar t) const
380
2/2
✓ Branch 0 taken 12163400 times.
✓ Branch 1 taken 21683700 times.
 33847100 { return -ddf2_(y)*df2_(x) * f_(t); }
381
382 33847100 Scalar dxyV_ (Scalar x, Scalar y, Scalar t) const
383
2/2
✓ Branch 0 taken 12163400 times.
✓ Branch 1 taken 21683700 times.
 33847100 { return -df2_(y)*ddf2_(x) * f_(t); }
384
385 33847100 Scalar dxxV_ (Scalar x, Scalar y, Scalar t) const
386
2/2
✓ Branch 0 taken 12163400 times.
✓ Branch 1 taken 21683700 times.
 33847100 { return -f2_(y)*dddf2_(x) * f_(t); }
387
388 33847100 Scalar dxUU_ (Scalar x, Scalar y, Scalar t) const
389 33847100 { return 2.*u_(x,y,t)*dxU_(x,y,t); }
390
391 33847100 Scalar dyVV_ (Scalar x, Scalar y, Scalar t) const
392 33847100 { return 2.*v_(x,y,t)*dyV_(x,y,t); }
393
394 33847100 Scalar dxUV_ (Scalar x, Scalar y, Scalar t) const
395 33847100 { return v_(x,y,t)*dxU_(x,y,t) + u_(x,y,t)*dxV_(x,y,t); }
396
397 33847100 Scalar dyUV_ (Scalar x, Scalar y, Scalar t) const
398 33847100 { return v_(x,y,t)*dyU_(x,y,t) + u_(x,y,t)*dyV_(x,y,t); }
399
400 Scalar rho_;
401 Scalar mu_;
402 Scalar time_;
403
404 bool isStationary_;
405 bool useNeumann_;
406 };
407
408 } // end namespace Dumux
409
410 #endif
411