GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/freeflow/navierstokesnc/channel/problem.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 48 55 87.3%
Functions: 11 20 55.0%
Branches: 64 90 71.1%
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 NavierStokesNCTests
10 * \brief Channel flow test for the multi-component staggered grid (Navier-)Stokes model.
11 */
12
13 #ifndef DUMUX_CHANNEL_NC_TEST_PROBLEM_HH
14 #define DUMUX_CHANNEL_NC_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/1pnc/advectiveflux.hh>
25
26 namespace Dumux {
27
28 /*!
29 * \ingroup NavierStokesNCTests
30 * \brief Test problem for the one-phase (Navier-)Stokes model.
31 *
32 * Flow from left to right in a channel is considered. A parabolic velocity
33 * profile is set at the left boundary, while the pressure is set to
34 * a fixed value on the right boundary. The top and bottom boundaries
35 * represent solid walls with no-slip/no-flow conditions.
36 */
37 template <class TypeTag, class BaseProblem>
38 class ChannelNCTestProblem : public BaseProblem
39 {
40 using ParentType = BaseProblem;
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 Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
47 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
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 SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
54 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
55
56 static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
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 using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
64
65 static constexpr bool useMoles = ModelTraits::useMoles();
66 static constexpr bool isNonisothermal = ModelTraits::enableEnergyBalance();
67 static constexpr auto compIdx = 1;
68
69 public:
70 24 ChannelNCTestProblem(std::shared_ptr<const GridGeometry> gridGeometry,
71 std::shared_ptr<CouplingManager> couplingManager)
72 : ParentType(gridGeometry, couplingManager)
73
8/24
✓ 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 21 taken 12 times.
✗ 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 31 not taken.
✗ Branch 32 not taken.
 96 , eps_(1e-6)
74 {
75
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 24 FluidSystem::init();
76
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 24 inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
77 24 }
78
79 /*!
80 * \brief Returns a reference pressure at a given sub control volume face.
81 * This pressure is subtracted from the actual pressure for the momentum balance
82 * which potentially helps to improve numerical accuracy by avoiding issues related do floating point arithmetic.
83 */
84 ✗ Scalar referencePressure(const Element& element,
85 const FVElementGeometry& fvGeometry,
86 const SubControlVolumeFace& scvf) const
87 ✗ { return 1.1e5; }
88
89 /*!
90 * \brief Specifies which kind of boundary condition should be
91 * used for which equation on a given boundary control volume.
92 *
93 * \param globalPos The position of the center of the finite volume
94 */
95 1399076 BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
96 {
97 1399076 BoundaryTypes values;
98
99 if constexpr (ParentType::isMomentumProblem())
100 {
101
2/2
✓ Branch 0 taken 152696 times.
✓ Branch 1 taken 130502 times.
 566396 values.setDirichlet(Indices::velocityXIdx);
102
2/2
✓ Branch 0 taken 152696 times.
✓ Branch 1 taken 130502 times.
 566396 values.setDirichlet(Indices::velocityYIdx);
103
104 1132792 if (isOutlet_(globalPos))
105 values.setAllNeumann();
106 }
107 else
108 {
109 832680 values.setAllNeumann();
110
111 1665360 if (isInlet_(globalPos))
112 {
113 97040 values.setDirichlet(Indices::pressureIdx);
114 97040 values.setDirichlet(Indices::conti0EqIdx + compIdx);
115 if constexpr (isNonisothermal)
116 39290 values.setDirichlet(Indices::temperatureIdx);
117 }
118 }
119
120 1399076 return values;
121 }
122
123 /*!
124 * \brief Evaluates the boundary conditions for a Dirichlet control volume.
125 *
126 * \param globalPos The center of the finite volume which ought to be set.
127 */
128 21440 DirichletValues dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
129 {
130 541286 const auto& globalPos = scvf.ipGlobal();
131 562726 DirichletValues values = initialAtPos(globalPos);
132
133 if constexpr (!ParentType::isMomentumProblem())
134 {
135 // give the system some time so that the pressure can equilibrate, then start the injection of the tracer
136 42880 if (isInlet_(globalPos))
137 {
138
6/6
✓ Branch 0 taken 1420 times.
✓ Branch 1 taken 9300 times.
✓ Branch 2 taken 1420 times.
✓ Branch 3 taken 9300 times.
✓ Branch 4 taken 1420 times.
✓ Branch 5 taken 9300 times.
 64320 values[Indices::pressureIdx] = this->couplingManager().cellPressure(element, scvf);
139
140
2/2
✓ Branch 0 taken 420 times.
✓ Branch 1 taken 1000 times.
 2840 if (time() >= 10.0 || inletVelocity_ < eps_)
141 {
142 19440 Scalar moleFracTransportedComp = 1e-3;
143 if constexpr (useMoles)
144 31060 values[Indices::conti0EqIdx+compIdx] = moleFracTransportedComp;
145 else
146 {
147 Scalar averageMolarMassPhase = moleFracTransportedComp * FluidSystem::molarMass(compIdx)
148 + (1. - moleFracTransportedComp) * FluidSystem::molarMass(1-compIdx);
149 values[Indices::conti0EqIdx+compIdx] = moleFracTransportedComp * FluidSystem::molarMass(compIdx)
150 / averageMolarMassPhase;
151 }
152
153 if constexpr (isNonisothermal)
154 15640 values[Indices::temperatureIdx] = 293.15;
155 }
156 }
157 }
158
159 541286 return values;
160 }
161
162 /*!
163 * \brief Evaluates the boundary conditions for a Neumann control volume.
164 *
165 * \param element The element for which the Neumann boundary condition is set
166 * \param fvGeometry The fvGeometry
167 * \param elemVolVars The element volume variables
168 * \param elemFaceVars The element face variables
169 * \param scvf The boundary sub control volume face
170 */
171 template<class ElementVolumeVariables, class ElementFluxVariablesCache>
172 571200 BoundaryFluxes neumann(const Element& element,
173 const FVElementGeometry& fvGeometry,
174 const ElementVolumeVariables& elemVolVars,
175 const ElementFluxVariablesCache& elemFluxVarsCache,
176 const SubControlVolumeFace& scvf) const
177 {
178 571200 BoundaryFluxes values(0.0);
179
180 if constexpr (ParentType::isMomentumProblem())
181 {
182 using FluxHelper = NavierStokesMomentumBoundaryFluxHelper;
183 ✗ values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf,
184 elemVolVars, elemFluxVarsCache,
185 referencePressure(element, fvGeometry, scvf), true /*zeroNormalVelocityGradient*/);
186 }
187 else
188 {
189 1142400 const auto insideDensity = elemVolVars[scvf.insideScvIdx()].density();
190
4/4
✓ Branch 2 taken 37800 times.
✓ Branch 3 taken 247800 times.
✓ Branch 4 taken 37800 times.
✓ Branch 5 taken 247800 times.
 2856000 values[Indices::conti0EqIdx] = this->faceVelocity(element, fvGeometry, scvf) * insideDensity * scvf.unitOuterNormal();
191
192 using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
193
2/2
✓ Branch 0 taken 37800 times.
✓ Branch 1 taken 247800 times.
 571200 if (isOutlet_(scvf.ipGlobal()))
194 75600 values = FluxHelper::scalarOutflowFlux( *this, element, fvGeometry, scvf, elemVolVars);
195 }
196
197 571200 return values;
198 }
199
200 /*!
201 * \brief Evaluates the initial value for a control volume.
202 *
203 * \param globalPos The global position
204 */
205 1061502 InitialValues initialAtPos(const GlobalPosition& globalPos) const
206 {
207 1073972 InitialValues values(0.0);
208
209 if constexpr (ParentType::isMomentumProblem())
210 {
211 // parabolic velocity profile
212 4246008 values[Indices::velocityXIdx] = parabolicProfile(globalPos[1], inletVelocity_);
213 }
214 else
215 {
216
1/2
✓ Branch 0 taken 10720 times.
✗ Branch 1 not taken.
 12470 values[Indices::pressureIdx] = 1.1e+5;
217 if constexpr (isNonisothermal)
218
1/2
✓ Branch 0 taken 4285 times.
✗ Branch 1 not taken.
 7785 values[Indices::temperatureIdx] = 283.15;
219 }
220
221 1067937 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 2123004 const Scalar yMin = this->gridGeometry().bBoxMin()[1];
233 2123004 const Scalar yMax = this->gridGeometry().bBoxMax()[1];
234 530751 return vMax * (y - yMin)*(yMax - y) / (0.25*(yMax - yMin)*(yMax - yMin));
235 }
236
237 void setTimeLoop(TimeLoopPtr timeLoop)
238
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 12 { timeLoop_ = timeLoop; }
239
240 Scalar time() const
241
6/6
✓ Branch 0 taken 1420 times.
✓ Branch 1 taken 9300 times.
✓ Branch 2 taken 1420 times.
✓ Branch 3 taken 9300 times.
✓ Branch 4 taken 1420 times.
✓ Branch 5 taken 9300 times.
 32160 { return timeLoop_->time(); }
242
243 //! Enable internal Dirichlet constraints
244 static constexpr bool enableInternalDirichletConstraints()
245 { return !ParentType::isMomentumProblem(); }
246
247 /*!
248 * \brief Tag a degree of freedom to carry internal Dirichlet constraints.
249 * If true is returned for a dof, the equation for this dof is replaced
250 * by the constraint that its primary variable values must match the
251 * user-defined values obtained from the function internalDirichlet(),
252 * which must be defined in the problem.
253 *
254 * \param element The finite element
255 * \param scv The sub-control volume
256 */
257 ✗ std::bitset<DirichletValues::dimension> hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
258 {
259 1273400 std::bitset<DirichletValues::dimension> values;
260 ✗ return values;
261 }
262
263 /*!
264 * \brief Define the values of internal Dirichlet constraints for a degree of freedom.
265 * \param element The finite element
266 * \param scv The sub-control volume
267 */
268 ✗ DirichletValues internalDirichlet(const Element& element, const SubControlVolume& scv) const
269 600000 { return DirichletValues(1.1e5); }
270
271 private:
272 ✗ bool isInlet_(const GlobalPosition& globalPos) const
273
6/8
✓ Branch 0 taken 10720 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10720 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 48520 times.
✓ Branch 5 taken 367820 times.
✓ Branch 6 taken 48520 times.
✓ Branch 7 taken 367820 times.
 854120 { return globalPos[0] < eps_; }
274
275 bool isOutlet_(const GlobalPosition& globalPos) const
276
20/20
✓ Branch 0 taken 37800 times.
✓ Branch 1 taken 247800 times.
✓ Branch 2 taken 37800 times.
✓ Branch 3 taken 247800 times.
✓ Branch 4 taken 37800 times.
✓ Branch 5 taken 247800 times.
✓ Branch 6 taken 37800 times.
✓ Branch 7 taken 247800 times.
✓ Branch 8 taken 37800 times.
✓ Branch 9 taken 247800 times.
✓ Branch 10 taken 152696 times.
✓ Branch 11 taken 130502 times.
✓ Branch 12 taken 152696 times.
✓ Branch 13 taken 130502 times.
✓ Branch 14 taken 152696 times.
✓ Branch 15 taken 130502 times.
✓ Branch 16 taken 152696 times.
✓ Branch 17 taken 130502 times.
✓ Branch 18 taken 152696 times.
✓ Branch 19 taken 130502 times.
 2843990 { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
277
278 const Scalar eps_;
279 Scalar inletVelocity_;
280 TimeLoopPtr timeLoop_;
281 };
282 } // end namespace Dumux
283
284 #endif
285