GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/freeflow/navierstokes/donea/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 111 139 79.9%
Functions: 35 140 25.0%
Branches: 96 151 63.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-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 (Donea 2003, \cite Donea2003).
11 */
12 #ifndef DUMUX_DONEA_TEST_PROBLEM_HH
13 #define DUMUX_DONEA_TEST_PROBLEM_HH
14
15 #include <dune/common/fmatrix.hh>
16 #include <dumux/common/math.hh>
17 #include <dumux/geometry/diameter.hh>
18
19 #include <dumux/common/properties.hh>
20 #include <dumux/common/parameters.hh>
21 #include <dumux/discretization/extrusion.hh>
22 #include <dumux/discretization/method.hh>
23 #include <dumux/freeflow/navierstokes/mass/1p/advectiveflux.hh>
24 #include <dumux/freeflow/navierstokes/mass/1p/localresidual.hh>
25
26 namespace Dumux {
27
28 /*!
29 * \ingroup NavierStokesTests
30 * \brief Test problem for the staggered grid (Donea 2003, \cite Donea2003).
31 *
32 * A two-dimensional Stokes flow in a square domain is considered.
33 * With the source terms as given in Donea 2003 \cite Donea2003, an analytical solution
34 * is available and can be compared to the numerical solution.
35 */
36 template <class TypeTag, class BaseProblem>
37 31 class DoneaTestProblem : public BaseProblem
38 {
39 using ParentType = BaseProblem;
40
41 using BoundaryTypes = typename ParentType::BoundaryTypes;
42 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
43 using FVElementGeometry = typename GridGeometry::LocalView;
44 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
45 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
46 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
47 using InitialValues = typename ParentType::InitialValues;
48 using Sources = typename ParentType::Sources;
49 using DirichletValues = typename ParentType::DirichletValues;
50 using BoundaryFluxes = typename ParentType::BoundaryFluxes;
51 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
52 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
53
54 static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
55 using Element = typename FVElementGeometry::Element;
56 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
57
58 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
59
60 public:
61 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
62
63 28 DoneaTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
64
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)
65 {
66
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 28 useNeumann_ = getParam<bool>("Problem.UseNeumann", false);
67
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 28 addBoxStabilization_ = getParam<bool>("Problem.AddBoxStabilization", false);
68
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 28 boxStabilizationParameter_ = getParam<Scalar>("Problem.StabilizationParameter", 0.1);
69
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 28 mu_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
70 28 }
71
72 17 DoneaTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
73
6/16
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 17 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 17 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 17 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 17 times.
✓ Branch 15 taken 17 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 51 : ParentType(gridGeometry)
74 {
75
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 17 useNeumann_ = getParam<bool>("Problem.UseNeumann", false);
76
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 17 addBoxStabilization_ = getParam<bool>("Problem.AddBoxStabilization", false);
77
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 17 boxStabilizationParameter_ = getParam<Scalar>("Problem.StabilizationParameter", 0.1);
78
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 17 mu_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
79 17 }
80
81 /*!
82 * \name Problem parameters
83 */
84 // \{
85
86 /*!
87 * \brief Return the sources within the domain.
88 *
89 * \param globalPos The global position
90 */
91 2249728 Sources sourceAtPos(const GlobalPosition &globalPos) const
92 {
93 if constexpr (ParentType::isMomentumProblem())
94 {
95 2249728 Sources source(0.0);
96 4499456 const Scalar x = globalPos[0];
97 4499456 const Scalar y = globalPos[1];
98
99 13498368 source[Indices::momentumXBalanceIdx] = -2.0*mu_*dxxU_(x,y) - mu_*dyyU_(x,y) - mu_*dxyV_(x,y) + dxP_(x,y);
100 11248640 source[Indices::momentumYBalanceIdx] = -2.0*mu_*dyyV_(x,y) - mu_*dxyU_(x,y) - mu_*dxxV_(x,y) + dyP_(x,y);
101 2249728 return source;
102 }
103 else
104 {
105
1/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 320320 times.
 320320 return Sources(0.0);
106 }
107 }
108 // \}
109
110 /*!
111 * \name Boundary conditions
112 */
113 // \{
114
115 /*!
116 * \brief Specifies which kind of boundary condition should be
117 * used for which equation on a given boundary control volume.
118 *
119 * \param globalPos The position of the center of the finite volume
120 */
121 21308 BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
122 {
123
3/4
✓ Branch 0 taken 948 times.
✓ Branch 1 taken 10240 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3328 times.
 35824 BoundaryTypes values;
124
125 // set Dirichlet values for the velocity and pressure everywhere
126 if constexpr (ParentType::isMomentumProblem())
127 {
128
2/2
✓ Branch 0 taken 6580 times.
✓ Branch 1 taken 8148 times.
 21308 if (useNeumann_)
129 {
130 static constexpr Scalar eps = 1e-8;
131
20/20
✓ Branch 0 taken 6106 times.
✓ Branch 1 taken 2042 times.
✓ Branch 2 taken 6106 times.
✓ Branch 3 taken 2042 times.
✓ Branch 4 taken 6106 times.
✓ Branch 5 taken 2042 times.
✓ Branch 6 taken 6106 times.
✓ Branch 7 taken 2042 times.
✓ Branch 8 taken 6106 times.
✓ Branch 9 taken 2042 times.
✓ Branch 10 taken 4070 times.
✓ Branch 11 taken 2036 times.
✓ Branch 12 taken 4070 times.
✓ Branch 13 taken 2036 times.
✓ Branch 14 taken 4070 times.
✓ Branch 15 taken 2036 times.
✓ Branch 16 taken 4070 times.
✓ Branch 17 taken 2036 times.
✓ Branch 18 taken 4070 times.
✓ Branch 19 taken 2036 times.
 45480 if ((globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps) || (globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps))
132 values.setAllNeumann();
133 else
134 values.setAllDirichlet();
135 }
136 else
137 values.setAllDirichlet();
138 }
139 else
140
6/8
✓ Branch 0 taken 948 times.
✓ Branch 1 taken 10240 times.
✓ Branch 2 taken 948 times.
✓ Branch 3 taken 10240 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3328 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 3328 times.
 29032 values.setNeumann(Indices::conti0EqIdx);
141
142 21308 return values;
143 }
144
145 /*!
146 * \brief Return dirichlet boundary values at a given position
147 *
148 * \param globalPos The global position
149 */
150 ✗ DirichletValues dirichletAtPos(const GlobalPosition& globalPos) const
151 7946 { return analyticalSolution(globalPos); }
152
153 /*!
154 * \brief Evaluates the boundary conditions for a Neumann control volume.
155 *
156 * \param element The element for which the Neumann boundary condition is set
157 * \param fvGeometry The fvGeometry
158 * \param elemVolVars The element volume variables
159 * \param elemFaceVars The element face variables
160 * \param scvf The boundary sub control volume face
161 */
162 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
163 65972 BoundaryFluxes neumann(const Element& element,
164 const FVElementGeometry& fvGeometry,
165 const ElementVolumeVariables& elemVolVars,
166 const ElementFluxVariablesCache& elemFluxVarsCache,
167 const SubControlVolumeFace& scvf) const
168 {
169 65972 BoundaryFluxes values(0.0);
170
171 if constexpr (ParentType::isMomentumProblem())
172 {
173 57704 const auto x = scvf.ipGlobal()[0];
174 57704 const auto y = scvf.ipGlobal()[1];
175
176 28852 Dune::FieldMatrix<Scalar, dimWorld, dimWorld> momentumFlux(0.0);
177 28852 momentumFlux[0][0] = -2.0*mu_*dxU_(x,y) + p_(x);
178 28852 momentumFlux[0][1] = -mu_*dyU_(x,y) - mu_*dxV_(x,y);
179 115408 momentumFlux[1][0] = momentumFlux[0][1];
180 28852 momentumFlux[1][1] = -2.0*mu_*dyV_(x,y) + p_(x);
181
182 28852 const auto normal = scvf.unitOuterNormal();
183 28852 momentumFlux.mv(normal, values);
184 }
185 else
186 {
187 90880 const auto insideDensity = elemVolVars[scvf.insideScvIdx()].density();
188 111360 values[Indices::conti0EqIdx] = this->faceVelocity(element, fvGeometry, scvf) * insideDensity * scvf.unitOuterNormal();
189
2/2
✓ Branch 0 taken 8320 times.
✓ Branch 1 taken 10240 times.
 37120 if (addBoxStabilization_)
190 16640 values[Indices::conti0EqIdx] += stabilizationFlux_(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
191 }
192
193 65972 return values;
194 }
195
196 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
197 ✗ Sources auxiliaryFlux(const Element& element,
198 const FVElementGeometry& fvGeometry,
199 const ElementVolumeVariables& elemVolVars,
200 const ElementFluxVariablesCache& elemFluxVarsCache,
201 const SubControlVolumeFace& scvf) const
202 {
203 // flux stabilization for the unstable Box-Box (P1-P1-CVFE) discretization scheme
204
1/4
✓ Branch 0 taken 249600 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 249600 Sources flux(0.0);
205
1/4
✓ Branch 0 taken 249600 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 249600 if (addBoxStabilization_)
206 {
207 249600 flux += stabilizationFlux_(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
208 using Extrusion = Extrusion_t<typename FVElementGeometry::GridGeometry>;
209 499200 flux *= Extrusion::area(fvGeometry, scvf);
210 }
211 ✗ return flux;
212 }
213
214 /*!
215 * \brief Return the analytical solution of the problem at a given position
216 *
217 * \param globalPos The global position
218 */
219 353014 DirichletValues analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
220 {
221
3/3
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 19258 times.
✓ Branch 2 taken 6996 times.
 392422 DirichletValues values(0.0);
222
223 if constexpr (ParentType::isMomentumProblem())
224 {
225 706028 const Scalar x = globalPos[0];
226 706028 const Scalar y = globalPos[1];
227 1412056 values[Indices::velocityXIdx] = f2_(x)*df2_(y);
228 1412056 values[Indices::velocityYIdx] = -f2_(y)*df2_(x);
229 }
230 else
231
3/3
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 19258 times.
✓ Branch 2 taken 6996 times.
 46480 values[Indices::pressureIdx] = p_(globalPos[0]);
232
233 353014 return values;
234 }
235
236 // \}
237
238 /*!
239 * \name Volume terms
240 */
241 // \{
242
243 //! TODO should these be spatial params?
244 ✗ Scalar pressureAtPos(const GlobalPosition& globalPos) const
245 1354944 { return p_(globalPos[0]); }
246
247 ✗ Scalar densityAtPos(const GlobalPosition& globalPos) const
248 ✗ { return 1.0; }
249
250 ✗ Scalar effectiveViscosityAtPos(const GlobalPosition& globalPos) const
251 ✗ { return 1.0; }
252
253 //! Enable internal Dirichlet constraints
254 static constexpr bool enableInternalDirichletConstraints()
255 { return !ParentType::isMomentumProblem(); }
256
257 /*!
258 * \brief Tag a degree of freedom to carry internal Dirichlet constraints.
259 * If true is returned for a dof, the equation for this dof is replaced
260 * by the constraint that its primary variable values must match the
261 * user-defined values obtained from the function internalDirichlet(),
262 * which must be defined in the problem.
263 *
264 * \param element The finite element
265 * \param scv The sub-control volume
266 */
267 197024 std::bitset<DirichletValues::dimension> hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
268 {
269 197024 std::bitset<DirichletValues::dimension> values;
270
271
2/2
✓ Branch 0 taken 62624 times.
✓ Branch 1 taken 134400 times.
 197024 if (!useNeumann_)
272 {
273
3/8
✓ Branch 1 taken 62624 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 62624 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 31312 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 219184 auto fvGeometry = localView(this->gridGeometry());
274
1/2
✓ Branch 1 taken 62624 times.
✗ Branch 2 not taken.
 62624 fvGeometry.bindElement(element);
275
276 62624 bool onBoundary = false;
277
5/6
✓ Branch 0 taken 250496 times.
✓ Branch 1 taken 62624 times.
✓ Branch 2 taken 250496 times.
✓ Branch 3 taken 62624 times.
✓ Branch 4 taken 125248 times.
✗ Branch 5 not taken.
 375744 for (const auto& scvf : scvfs(fvGeometry))
278
7/10
✓ Branch 0 taken 125248 times.
✓ Branch 1 taken 125248 times.
✓ Branch 2 taken 125248 times.
✓ Branch 3 taken 125248 times.
✓ Branch 4 taken 250496 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 250496 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 250496 times.
✗ Branch 9 not taken.
 876736 if (fvGeometry.scv(scvf.insideScvIdx()).dofIndex() == scv.dofIndex())
279
2/2
✓ Branch 0 taken 7184 times.
✓ Branch 1 taken 243312 times.
 257680 onBoundary = std::max(onBoundary, scvf.boundary());
280
281
2/2
✓ Branch 0 taken 7184 times.
✓ Branch 1 taken 55440 times.
 62624 if (onBoundary)
282 7184 values.set(0);
283
284 // TODO: only use one cell or pass fvGeometry to hasInternalDirichletConstraint
285
286 // if (scv.dofIndex() == 0)
287 // values.set(0);
288 // the pure Neumann problem is only defined up to a constant
289 // we create a well-posed problem by fixing the pressure at one dof
290 }
291
292 197024 return values;
293 }
294
295 /*!
296 * \brief Define the values of internal Dirichlet constraints for a degree of freedom.
297 * \param element The finite element
298 * \param scv The sub-control volume
299 */
300 ✗ DirichletValues internalDirichlet(const Element& element, const SubControlVolume& scv) const
301 14144 { return DirichletValues(analyticalSolution(scv.dofPosition())[Indices::pressureIdx]); }
302
303 private:
304 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
305 257920 Sources stabilizationFlux_(const Element& element,
306 const FVElementGeometry& fvGeometry,
307 const ElementVolumeVariables& elemVolVars,
308 const ElementFluxVariablesCache& elemFluxVarsCache,
309 const SubControlVolumeFace& scvf) const
310 {
311
1/2
✓ Branch 0 taken 257920 times.
✗ Branch 1 not taken.
 257920 Sources flux(0.0);
312
313 if constexpr (!ParentType::isMomentumProblem() && GridGeometry::discMethod == DiscretizationMethods::box)
314 {
315
1/2
✓ Branch 0 taken 257920 times.
✗ Branch 1 not taken.
 257920 if (addBoxStabilization_)
316 {
317 // add -rho*C*h^2/mu*(grad P - f) as stabilization, see for instance
318 // Quarteroni & Ruiz-Baier (2011) https://doi.org/10.1007/s00211-011-0373-4
319 // (note that most stabilization terms vanish due to the use of linear basis functions)
320 257920 const auto gradP = [&]
321 {
322
2/2
✓ Branch 0 taken 8320 times.
✓ Branch 1 taken 249600 times.
 257920 if (scvf.boundary())
323 {
324 8320 const auto& globalPos = scvf.ipGlobal();
325 41600 return Dune::FieldVector<Scalar, dimWorld>({
326 8320 dxP_(globalPos[0], globalPos[1]),
327 8320 dyP_(globalPos[0], globalPos[1])
328 16640 });
329 }
330 else
331 {
332 249600 const auto& fluxVarCache = elemFluxVarsCache[scvf];
333 249600 Dune::FieldVector<Scalar, dimWorld> gradP(0.0);
334
4/4
✓ Branch 0 taken 748800 times.
✓ Branch 1 taken 249600 times.
✓ Branch 2 taken 748800 times.
✓ Branch 3 taken 249600 times.
 1996800 for (const auto& scv : scvs(fvGeometry))
335 3744000 gradP.axpy(elemVolVars[scv].pressure(0), fluxVarCache.gradN(scv.indexInElement()));
336 249600 return gradP;
337 }
338 266240 }();
339
340 257920 const auto rho = densityAtPos(scvf.ipGlobal());
341 257920 const auto mu = effectiveViscosityAtPos(scvf.ipGlobal());
342 257920 const auto h = diameter(element.geometry());
343
344 515840 const auto f = this->sourceAtPos(scvf.ipGlobal());
345 1031680 flux[0] = -1.0*vtmv(scvf.unitOuterNormal(), boxStabilizationParameter_*h*h*rho/mu, gradP - f);
346 }
347 }
348
349 257920 return flux;
350 }
351
352 ✗ Scalar p_(Scalar x) const
353
3/3
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 19258 times.
✓ Branch 2 taken 6996 times.
 1401424 { return x*(1.0-x); }
354
355 ✗ Scalar dP_(Scalar x) const
356 1701376 { return 1.0 - 2.0*x; }
357
358 ✗ Scalar f2_(Scalar x) const
359 7787960 { return p_(x)*p_(x); /*=x^2*(1-2x+x^2)=x^2-2x^3+x^4*/ }
360
361 ✗ Scalar df2_(Scalar x) const
362 3906300 { return 2.0*x - 6.0*x*x + 4.0*x*x*x; }
363
364 ✗ Scalar ddf2_(Scalar x) const
365 3435392 { return 2.0 - 12.0*x + 12.0*x*x; }
366
367 ✗ Scalar dddf2_(Scalar x) const
368 3386112 { return - 12.0 + 24.0*x; }
369
370 ✗ Scalar dxP_ (Scalar x, Scalar y) const
371 2895232 { return dP_(x); }
372
373 ✗ Scalar dyP_ (Scalar x, Scalar y) const
374 ✗ { return 0.0; }
375
376 ✗ Scalar dxU_ (Scalar x, Scalar y) const
377 24640 { return df2_(x)*df2_(y); }
378
379 ✗ Scalar dyU_ (Scalar x, Scalar y) const
380 73920 { return f2_(x)*ddf2_(y); }
381
382 ✗ Scalar dxV_ (Scalar x, Scalar y) const
383 24640 { return -f2_(y)*ddf2_(x); }
384
385 ✗ Scalar dyV_ (Scalar x, Scalar y) const
386 24640 { return -df2_(y)*df2_(x); }
387
388 ✗ Scalar dxxU_ (Scalar x, Scalar y) const
389 5079168 { return ddf2_(x)*df2_(y); }
390
391 ✗ Scalar dxyU_ (Scalar x, Scalar y) const
392 1693056 { return df2_(x)*ddf2_(y); }
393
394 ✗ Scalar dyyU_ (Scalar x, Scalar y) const
395 5079168 { return f2_(x)*dddf2_(y); }
396
397 ✗ Scalar dyyV_ (Scalar x, Scalar y) const
398 5079168 { return -ddf2_(y)*df2_(x); }
399
400 ✗ Scalar dxyV_ (Scalar x, Scalar y) const
401 1693056 { return -df2_(y)*ddf2_(x); }
402
403 ✗ Scalar dxxV_ (Scalar x, Scalar y) const
404 5079168 { return -f2_(y)*dddf2_(x); }
405
406 bool useNeumann_;
407 bool addBoxStabilization_;
408 Scalar boxStabilizationParameter_;
409 Scalar mu_;
410 };
411
412 // enable auxiliary flux evaluation in the mass local residual
413 // this is used to implement the stabilization for the Box-Box discretization method
414 template <class TypeTag, class BaseProblem>
415 struct ImplementsAuxiliaryFluxNavierStokesMassOneP<DoneaTestProblem<TypeTag, BaseProblem>>
416 : public std::true_type
417 {};
418
419 } // end namespace Dumux
420
421 #endif
422