GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/freeflow/navierstokes/channel/1d/problem.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 51 71 71.8%
Functions: 4 18 22.2%
Branches: 59 106 55.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 1-D Navier-Stokes model with an analytical solution.
11 *
12 * \copydoc Dumux::NavierStokesAnalyticProblem
13 */
14 #ifndef DUMUX_DONEA_TEST_PROBLEM_HH
15 #define DUMUX_DONEA_TEST_PROBLEM_HH
16
17 #include <dune/common/fmatrix.hh>
18
19 #include <dumux/common/properties.hh>
20 #include <dumux/common/parameters.hh>
21 #include <dumux/freeflow/navierstokes/boundarytypes.hh>
22
23 namespace Dumux {
24
25 /*!
26 * \ingroup NavierStokesTests
27 * \brief Test for the 1-D Navier-Stokes model with an analytical solution.
28 */
29 template <class TypeTag, class BaseProblem>
30 2 class NavierStokesAnalyticProblem : public BaseProblem
31 {
32 using ParentType = BaseProblem;
33
34 using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
35 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
36 using InitialValues = typename ParentType::InitialValues;
37 using Sources = typename ParentType::Sources;
38 using DirichletValues = typename ParentType::DirichletValues;
39 using BoundaryFluxes = typename ParentType::BoundaryFluxes;
40 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
41 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
42 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
43 using SubControlVolume = typename GridGeometry::SubControlVolume;
44 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
45 using FVElementGeometry = typename GridGeometry::LocalView;
46
47 static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
48 using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
49 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
50 using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
51 using DimVector = GlobalPosition;
52 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
53
54 public:
55 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
56
57 4 NavierStokesAnalyticProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
58
7/20
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 2 times.
✓ Branch 18 taken 2 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.
 16 : ParentType(gridGeometry, couplingManager)
59 {
60
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 density_ = getParam<Scalar>("Component.LiquidDensity");
61
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity");
62 4 }
63
64 /*!
65 * \brief Returns the sources within the domain.
66 *
67 * \param globalPos The global position
68 */
69 3200 Sources sourceAtPos(const GlobalPosition &globalPos) const
70 {
71 3200 Sources source(0.0);
72
73 if constexpr (!ParentType::isMomentumProblem())
74 {
75 // mass balance - term div(rho*v)
76 ✗ for (unsigned int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
77 {
78 ✗ source[Indices::conti0EqIdx] += dvdx(globalPos)[dimIdx][dimIdx];
79 }
80 ✗ source[Indices::conti0EqIdx] *= density_;
81 }
82 else
83 {
84 // momentum balance
85
2/2
✓ Branch 0 taken 1600 times.
✓ Branch 1 taken 1600 times.
 6400 for (unsigned int velIdx = 0; velIdx < dimWorld; ++velIdx)
86 {
87
2/2
✓ Branch 0 taken 1600 times.
✓ Branch 1 taken 1600 times.
 6400 for (unsigned int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
88 {
89 // inertia term
90
1/2
✓ Branch 0 taken 1600 times.
✗ Branch 1 not taken.
 3200 if (this->enableInertiaTerms())
91 6400 source[Indices::velocity(velIdx)] += density_ * dv2dx(globalPos)[velIdx][dimIdx];
92
93 // viscous term (molecular)
94 9600 source[Indices::velocity(velIdx)] -= density_ * kinematicViscosity_* dvdx2(globalPos)[velIdx][dimIdx];
95
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1599 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 3200 static const bool enableUnsymmetrizedVelocityGradient = getParam<bool>("FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
96
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1600 times.
 3200 if (!enableUnsymmetrizedVelocityGradient)
97 3200 source[Indices::velocity(velIdx)] -= density_ * kinematicViscosity_* dvdx2(globalPos)[dimIdx][velIdx];
98 }
99 // pressure term
100 6400 source[Indices::velocity(velIdx)] += dpdx(globalPos)[velIdx];
101
102 // gravity term
103
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 static const bool enableGravity = getParam<bool>("Problem.EnableGravity");
104
1/2
✓ Branch 0 taken 1600 times.
✗ Branch 1 not taken.
 3200 if (enableGravity)
105 {
106 ✗ source[Indices::velocity(velIdx)] -= density_ * this->gravity()[velIdx];
107 }
108 }
109 }
110
111 3200 return source;
112 }
113 // \}
114 /*!
115 * \name Boundary conditions
116 */
117 // \{
118
119 /*!
120 * \brief Specifies which kind of boundary condition should be
121 * used for which equation on a given boundary control volume.
122 *
123 * \param globalPos The position of the center of the finite volume
124 */
125 ✗ BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
126 {
127
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 60 times.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 26 times.
 96 BoundaryTypes values;
128
129 if constexpr (ParentType::isMomentumProblem())
130 {
131 // set Dirichlet values for the velocity everywhere
132
2/8
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 20 values.setDirichlet(Indices::momentumXBalanceIdx);
133 }
134 else
135
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 60 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 26 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 26 times.
 172 values.setNeumann(Indices::pressureIdx);
136
137 ✗ return values;
138 }
139
140 /*!
141 * \brief Returns Dirichlet boundary values at a given position
142 *
143 * \param globalPos The global position
144 */
145 ✗ DirichletValues dirichletAtPos(const GlobalPosition& globalPos) const
146 {
147 // use the values of the analytical solution
148 20 return analyticalSolution(globalPos);
149 }
150
151 /*!
152 * \brief Returns the analytical solution of the problem at a given position
153 *
154 * \param globalPos The global position
155 * \param time A parameter for consistent signatures. It is ignored here as this is a stationary test.
156 */
157 ✗ DirichletValues analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
158 {
159 554 DirichletValues values;
160
161 if constexpr (ParentType::isMomentumProblem())
162 596 values[Indices::velocityXIdx] = v(globalPos);
163 else
164 576 values[Indices::pressureIdx] = p(globalPos);
165
166 ✗ return values;
167 }
168
169 //! \brief The velocity
170 ✗ const DimVector v(const DimVector& globalPos) const
171 {
172
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 189 times.
 298 DimVector v(0.0);
173
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 189 times.
 298 v[0] = 2.0 * globalPos[0] * globalPos[0] * globalPos[0];
174 ✗ return v;
175 }
176
177 //! \brief The velocity gradient
178 ✗ const DimMatrix dvdx(const DimVector& globalPos) const
179 {
180 ✗ DimMatrix dvdx(0.0);
181
0/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 3200 dvdx[0][0] = 6.0 * globalPos[0] * globalPos[0];
182 ✗ return dvdx;
183 }
184
185 //! \brief The gradient of the velocity squared (using product rule -> nothing to do here)
186 const DimMatrix dv2dx(const DimVector& globalPos) const
187 {
188 1600 DimMatrix dv2dx;
189
2/2
✓ Branch 0 taken 1600 times.
✓ Branch 1 taken 1600 times.
 3200 for (unsigned int velIdx = 0; velIdx < dimWorld; ++velIdx)
190 {
191 3200 for (unsigned int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
192 {
193 6400 dv2dx[velIdx][dimIdx] = dvdx(globalPos)[velIdx][dimIdx] * v(globalPos)[dimIdx]
194 6400 + dvdx(globalPos)[dimIdx][dimIdx] * v(globalPos)[velIdx];
195 }
196 }
197 return dv2dx;
198 }
199
200 //! \brief The gradient of the velocity gradient
201 ✗ const DimMatrix dvdx2(const DimVector& globalPos) const
202 {
203
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1599 times.
 3200 DimMatrix dvdx2(0.0);
204
4/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1599 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 1599 times.
 6400 dvdx2[0][0] = 12.0 * globalPos[0];
205 ✗ return dvdx2;
206 }
207
208 //! \brief The pressure
209 ✗ const Scalar p(const DimVector& globalPos) const
210
2/2
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 84 times.
 128 { return 2.0 - 2.0 * globalPos[0]; }
211
212 //! \brief The pressure gradient
213 ✗ const DimVector dpdx(const DimVector& globalPos) const
214 {
215
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1599 times.
 1600 DimVector dpdx(0.0);
216
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1599 times.
 1600 dpdx[0] = -2.0;
217 ✗ return dpdx;
218 }
219
220 //! Enable internal Dirichlet constraints
221 static constexpr bool enableInternalDirichletConstraints()
222 { return !ParentType::isMomentumProblem(); }
223
224 /*!
225 * \brief Tag a degree of freedom to carry internal Dirichlet constraints.
226 * If true is returned for a dof, the equation for this dof is replaced
227 * by the constraint that its primary variable values must match the
228 * user-defined values obtained from the function internalDirichlet(),
229 * which must be defined in the problem.
230 *
231 * \param element The finite element
232 * \param scv The sub-control volume
233 */
234 790 std::bitset<DirichletValues::dimension> hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
235 {
236 790 std::bitset<DirichletValues::dimension> values;
237
238 1580 auto fvGeometry = localView(this->gridGeometry());
239 790 fvGeometry.bindElement(element);
240
241 790 bool onBoundary = false;
242
6/6
✓ Branch 0 taken 1580 times.
✓ Branch 1 taken 790 times.
✓ Branch 2 taken 1580 times.
✓ Branch 3 taken 790 times.
✓ Branch 4 taken 40 times.
✓ Branch 5 taken 1540 times.
 4740 for (const auto& scvf : scvfs(fvGeometry))
243
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 1540 times.
 1620 onBoundary = std::max(onBoundary, scvf.boundary());
244
245
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 750 times.
 790 if (onBoundary)
246 40 values.set(0);
247
248 // TODO: only use one cell or pass fvGeometry to hasInternalDirichletConstraint
249
250 // if (scv.dofIndex() == 0)
251 // values.set(0);
252 // the pure Neumann problem is only defined up to a constant
253 // we create a well-posed problem by fixing the pressure at one dof
254 790 return values;
255 }
256
257 /*!
258 * \brief Define the values of internal Dirichlet constraints for a degree of freedom.
259 * \param element The finite element
260 * \param scv The sub-control volume
261 */
262 ✗ DirichletValues internalDirichlet(const Element& element, const SubControlVolume& scv) const
263 320 { return DirichletValues(analyticalSolution(scv.center())[Indices::pressureIdx]); }
264
265 // \}
266
267 /*!
268 * \name Volume terms
269 */
270 // \{
271
272 /*!
273 * \brief Evaluates the initial value for a control volume.
274 *
275 * \param globalPos The global position
276 */
277 InitialValues initialAtPos(const GlobalPosition& globalPos) const
278 {
279 return analyticalSolution(globalPos);
280 }
281
282 private:
283 Scalar density_;
284 Scalar kinematicViscosity_;
285 };
286 } // end namespace Dumux
287
288 #endif
289