GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/freeflow/navierstokes/staggered/localresidual.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 70 79 88.6%
Functions: 144 332 43.4%
Branches: 54 90 60.0%
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 NavierStokesModel
10 * \copydoc Dumux::NavierStokesResidualImpl
11 */
12 #ifndef DUMUX_STAGGERED_NAVIERSTOKES_LOCAL_RESIDUAL_HH
13 #define DUMUX_STAGGERED_NAVIERSTOKES_LOCAL_RESIDUAL_HH
14
15 #include <dune/common/hybridutilities.hh>
16
17 #include <dumux/common/properties.hh>
18 #include <dumux/discretization/method.hh>
19 #include <dumux/discretization/extrusion.hh>
20 #include <dumux/assembly/staggeredlocalresidual.hh>
21 #include <dumux/freeflow/nonisothermal/localresidual.hh>
22
23 namespace Dumux {
24
25 namespace Impl {
26 template<class T>
27 static constexpr bool isRotationalExtrusion = false;
28
29 template<int radialAxis>
30 static constexpr bool isRotationalExtrusion<RotationalExtrusion<radialAxis>> = true;
31 } // end namespace Impl
32
33 /*!
34 * \ingroup NavierStokesModel
35 * \brief Element-wise calculation of the Navier-Stokes residual for models using the staggered discretization
36 */
37
38 // forward declaration
39 template<class TypeTag, class DiscretizationMethod>
40 class NavierStokesResidualImpl;
41
42 template<class TypeTag>
43 class NavierStokesResidualImpl<TypeTag, DiscretizationMethods::Staggered>
44 : public StaggeredLocalResidual<TypeTag>
45 {
46 using ParentType = StaggeredLocalResidual<TypeTag>;
47 friend class StaggeredLocalResidual<TypeTag>;
48
49 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
50
51 using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
52 using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
53 using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
54
55 using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
56 using ElementFluxVariablesCache = typename GridFluxVariablesCache::LocalView;
57
58 using GridFaceVariables = typename GridVariables::GridFaceVariables;
59 using ElementFaceVariables = typename GridFaceVariables::LocalView;
60
61 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
62 using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;
63 using Problem = GetPropType<TypeTag, Properties::Problem>;
64 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
65 using FVElementGeometry = typename GridGeometry::LocalView;
66 using GridView = typename GridGeometry::GridView;
67 using Element = typename GridView::template Codim<0>::Entity;
68 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
69 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
70 using Extrusion = Extrusion_t<GridGeometry>;
71 using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
72 using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
73 using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
74 using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
75 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
76
77 static constexpr bool normalizePressure = getPropValue<TypeTag, Properties::NormalizePressure>();
78
79 using CellCenterResidual = CellCenterPrimaryVariables;
80 using FaceResidual = FacePrimaryVariables;
81
82 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
83
84 public:
85 using EnergyLocalResidual = FreeFlowEnergyLocalResidual<GridGeometry, FluxVariables, ModelTraits::enableEnergyBalance(), (ModelTraits::numFluidComponents() > 1)>;
86
87 // account for the offset of the cell center privars within the PrimaryVariables container
88 static constexpr auto cellCenterOffset = ModelTraits::numEq() - CellCenterPrimaryVariables::dimension;
89 static_assert(cellCenterOffset == ModelTraits::dim(), "cellCenterOffset must equal dim for staggered NavierStokes");
90
91 //! Use the parent type's constructor
92
6/8
✓ Branch 1 taken 1355444 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1355444 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2578 times.
✓ Branch 7 taken 50 times.
✓ Branch 8 taken 2578 times.
✓ Branch 9 taken 50 times.
 2716144 using ParentType::ParentType;
93
94 //! Evaluate fluxes entering or leaving the cell center control volume.
95 ✗ CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
96 const Element &element,
97 const FVElementGeometry& fvGeometry,
98 const ElementVolumeVariables& elemVolVars,
99 const ElementFaceVariables& elemFaceVars,
100 const SubControlVolumeFace &scvf,
101 const ElementFluxVariablesCache& elemFluxVarsCache) const
102 {
103 FluxVariables fluxVars;
104 118805706 CellCenterPrimaryVariables flux = fluxVars.computeMassFlux(problem, element, fvGeometry, elemVolVars,
105
1/8
✓ Branch 1 taken 138882 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 202698350 elemFaceVars, scvf, elemFluxVarsCache[scvf]);
106
107
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 99555458 EnergyLocalResidual::heatFlux(flux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
108
109 11358064 return flux;
110 }
111
112 //! Evaluate the source term for the cell center control volume.
113 ✗ CellCenterPrimaryVariables computeSourceForCellCenter(const Problem& problem,
114 const Element &element,
115 const FVElementGeometry& fvGeometry,
116 const ElementVolumeVariables& elemVolVars,
117 const ElementFaceVariables& elemFaceVars,
118 const SubControlVolume &scv) const
119 {
120
0/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 37620788 CellCenterPrimaryVariables result(0.0);
121
122 // get the values from the problem
123
0/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 37620788 const auto sourceValues = problem.source(element, fvGeometry, elemVolVars, elemFaceVars, scv);
124
125 // copy the respective cell center related values to the result
126
2/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 131014696 times.
✓ Branch 3 taken 37620788 times.
 168635484 for (int i = 0; i < result.size(); ++i)
127 393044088 result[i] = sourceValues[i + cellCenterOffset];
128
129 3746088 return result;
130 }
131
132
133 //! Evaluate the storage term for the cell center control volume.
134 ✗ CellCenterPrimaryVariables computeStorageForCellCenter(const Problem& problem,
135 const SubControlVolume& scv,
136 const VolumeVariables& volVars) const
137 {
138 5070312 CellCenterPrimaryVariables storage;
139 93959880 storage[Indices::conti0EqIdx - ModelTraits::dim()] = volVars.density();
140
141 65777128 EnergyLocalResidual::fluidPhaseStorage(storage, volVars);
142
143 ✗ return storage;
144 }
145
146 //! Evaluate the storage term for the face control volume.
147 ✗ FacePrimaryVariables computeStorageForFace(const Problem& problem,
148 const SubControlVolumeFace& scvf,
149 const VolumeVariables& volVars,
150 const ElementFaceVariables& elemFaceVars) const
151 {
152 135806428 FacePrimaryVariables storage(0.0);
153 135806428 const Scalar velocity = elemFaceVars[scvf].velocitySelf();
154 271612856 storage[0] = volVars.density() * velocity;
155 ✗ return storage;
156 }
157
158 //! Evaluate the source term for the face control volume.
159 150537748 FacePrimaryVariables computeSourceForFace(const Problem& problem,
160 const Element& element,
161 const FVElementGeometry& fvGeometry,
162 const SubControlVolumeFace& scvf,
163 const ElementVolumeVariables& elemVolVars,
164 const ElementFaceVariables& elemFaceVars) const
165 {
166 150537748 FacePrimaryVariables source(0.0);
167 301075496 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
168 602150992 source += problem.gravity()[scvf.directionIndex()] * insideVolVars.density();
169 150537748 source += problem.source(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())];
170
171 // Axisymmetric problems in 2D feature an extra source terms arising from the transformation to cylindrical coordinates.
172 // See Ferziger/Peric: Computational methods for fluid dynamics chapter 8.
173 // https://doi.org/10.1007/978-3-540-68228-8 (page 301)
174 if constexpr (ModelTraits::dim() == 2 && Impl::isRotationalExtrusion<Extrusion>)
175 {
176 if (scvf.directionIndex() == Extrusion::radialAxis)
177 {
178 // Velocity term
179 const auto r = scvf.center()[scvf.directionIndex()];
180 source -= -2.0*insideVolVars.effectiveViscosity() * elemFaceVars[scvf].velocitySelf() / (r*r);
181
182 // Pressure term (needed because we incorporate pressure in terms of a surface integral).
183 // grad(p) becomes div(pI) + (p/r)*n_r in cylindrical coordinates. The second term
184 // is new with respect to Cartesian coordinates and handled below as a source term.
185 const Scalar pressure =
186 normalizePressure ? insideVolVars.pressure() - problem.initial(scvf)[Indices::pressureIdx]
187 : insideVolVars.pressure();
188
189 source += pressure/r;
190 }
191 }
192
193 150537748 return source;
194 }
195
196 //! Evaluate the momentum flux for the face control volume.
197 ✗ FacePrimaryVariables computeFluxForFace(const Problem& problem,
198 const Element& element,
199 const SubControlVolumeFace& scvf,
200 const FVElementGeometry& fvGeometry,
201 const ElementVolumeVariables& elemVolVars,
202 const ElementFaceVariables& elemFaceVars,
203 const ElementFluxVariablesCache& elemFluxVarsCache) const
204 {
205 FluxVariables fluxVars;
206
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 141926782 return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache());
207 }
208
209 /*!
210 * \brief Evaluate boundary conditions for a cell center dof
211 */
212 9358728 CellCenterResidual computeBoundaryFluxForCellCenter(const Problem& problem,
213 const Element& element,
214 const FVElementGeometry& fvGeometry,
215 const SubControlVolumeFace& scvf,
216 const ElementVolumeVariables& elemVolVars,
217 const ElementFaceVariables& elemFaceVars,
218 const ElementBoundaryTypes& elemBcTypes,
219 const ElementFluxVariablesCache& elemFluxVarsCache) const
220 {
221
1/2
✓ Branch 0 taken 287554 times.
✗ Branch 1 not taken.
 9358728 CellCenterResidual result(0.0);
222
223
1/2
✓ Branch 0 taken 9358728 times.
✗ Branch 1 not taken.
 9358728 if (scvf.boundary())
224 {
225 9358728 const auto bcTypes = problem.boundaryTypes(element, scvf);
226
227 // no fluxes occur over symmetry boundaries
228
4/4
✓ Branch 0 taken 7676456 times.
✓ Branch 1 taken 1682272 times.
✓ Branch 2 taken 7676456 times.
✓ Branch 3 taken 1682272 times.
 18717456 if (bcTypes.isSymmetry())
229 551970 return result;
230
231 // treat Dirichlet and outflow BCs
232
1/2
✓ Branch 1 taken 138882 times.
✗ Branch 2 not taken.
 8806758 result = computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
233
234 // treat Neumann BCs, i.e. overwrite certain fluxes by user-specified values
235 static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
236
4/4
✓ Branch 0 taken 929888 times.
✓ Branch 1 taken 7876870 times.
✓ Branch 2 taken 929888 times.
✓ Branch 3 taken 7876870 times.
 17613516 if (bcTypes.hasNeumann())
237 {
238 1859776 const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
239
1/2
✓ Branch 1 taken 41192 times.
✗ Branch 2 not taken.
 929888 const auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
240
241
2/2
✓ Branch 0 taken 3159654 times.
✓ Branch 1 taken 929888 times.
 4089542 for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
242 {
243
2/2
✓ Branch 0 taken 1766956 times.
✓ Branch 1 taken 1392698 times.
 3159654 if (bcTypes.isNeumann(eqIdx + cellCenterOffset))
244 5529600 result[eqIdx] = neumannFluxes[eqIdx + cellCenterOffset] * extrusionFactor * Extrusion::area(fvGeometry, scvf);
245 }
246 }
247
248 // account for wall functions, if used
249 14593054 incorporateWallFunction_(result, problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars);
250 }
251 1150388 return result;
252 }
253
254 /*!
255 * \brief Evaluate Dirichlet (fixed value) boundary conditions for a face dof
256 */
257 150537748 void evalDirichletBoundariesForFace(FaceResidual& residual,
258 const Problem& problem,
259 const Element& element,
260 const FVElementGeometry& fvGeometry,
261 const SubControlVolumeFace& scvf,
262 const ElementVolumeVariables& elemVolVars,
263 const ElementFaceVariables& elemFaceVars,
264 const ElementBoundaryTypes& elemBcTypes,
265 const ElementFluxVariablesCache& elemFluxVarsCache) const
266 {
267
2/2
✓ Branch 0 taken 8610966 times.
✓ Branch 1 taken 141926782 times.
 150537748 if (scvf.boundary())
268 {
269 // handle the actual boundary conditions:
270 8610966 const auto bcTypes = problem.boundaryTypes(element, scvf);
271
272
2/2
✓ Branch 0 taken 3277000 times.
✓ Branch 1 taken 5333966 times.
 8610966 if(bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
273 {
274 // set a fixed value for the velocity for Dirichlet boundary conditions
275
2/2
✓ Branch 0 taken 5976 times.
✓ Branch 1 taken 1488 times.
 5333966 const Scalar velocity = elemFaceVars[scvf].velocitySelf();
276
2/2
✓ Branch 0 taken 5976 times.
✓ Branch 1 taken 1488 times.
 5441614 const Scalar dirichletValue = problem.dirichlet(element, scvf)[Indices::velocity(scvf.directionIndex())];
277 10667932 residual = velocity - dirichletValue;
278 }
279
4/4
✓ Branch 0 taken 489258 times.
✓ Branch 1 taken 2787742 times.
✓ Branch 2 taken 489258 times.
✓ Branch 3 taken 2787742 times.
 6554000 else if(bcTypes.isSymmetry())
280 {
281 // for symmetry boundary conditions, there is no flow across the boundary and
282 // we therefore treat it like a Dirichlet boundary conditions with zero velocity
283 489258 const Scalar velocity = elemFaceVars[scvf].velocitySelf();
284 489258 const Scalar fixedValue = 0.0;
285 978516 residual = velocity - fixedValue;
286 }
287 }
288 150537748 }
289
290 /*!
291 * \brief Evaluate boundary fluxes for a face dof
292 */
293 8610966 FaceResidual computeBoundaryFluxForFace(const Problem& problem,
294 const Element& element,
295 const FVElementGeometry& fvGeometry,
296 const SubControlVolumeFace& scvf,
297 const ElementVolumeVariables& elemVolVars,
298 const ElementFaceVariables& elemFaceVars,
299 const ElementBoundaryTypes& elemBcTypes,
300 const ElementFluxVariablesCache& elemFluxVarsCache) const
301 {
302
1/2
✓ Branch 0 taken 8610966 times.
✗ Branch 1 not taken.
 8610966 FaceResidual result(0.0);
303
304
1/2
✓ Branch 0 taken 8610966 times.
✗ Branch 1 not taken.
 8610966 if (scvf.boundary())
305 {
306 FluxVariables fluxVars;
307
308 // handle the actual boundary conditions:
309 8610966 const auto bcTypes = problem.boundaryTypes(element, scvf);
310
2/2
✓ Branch 0 taken 8437702 times.
✓ Branch 1 taken 173264 times.
 8610966 if (bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
311 {
312 // the source term has already been accounted for, here we
313 // add a given Neumann flux for the face on the boundary itself ...
314 346528 const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
315
2/4
✓ Branch 1 taken 173264 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 173264 times.
✗ Branch 5 not taken.
 173264 result = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
316
2/4
✓ Branch 1 taken 173264 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 173264 times.
✗ Branch 5 not taken.
 346528 * extrusionFactor * Extrusion::area(fvGeometry, scvf);
317
318 // ... and treat the fluxes of the remaining (frontal and lateral) faces of the staggered control volume
319
1/2
✓ Branch 1 taken 173264 times.
✗ Branch 2 not taken.
 173264 result += fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache());
320 }
321
2/2
✓ Branch 0 taken 5823224 times.
✓ Branch 1 taken 2614478 times.
 8437702 else if(bcTypes.isDirichlet(Indices::pressureIdx))
322 {
323 // we are at an "fixed pressure" boundary for which the resdiual of the momentum balance needs to be assembled
324 // as if it where inside the domain and not on the boundary (source term has already been acounted for)
325
1/2
✓ Branch 1 taken 2614478 times.
✗ Branch 2 not taken.
 2614478 result = fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache());
326
327 // incorporate the inflow or outflow contribution
328
1/2
✓ Branch 1 taken 79436 times.
✗ Branch 2 not taken.
 2614478 result += fluxVars.inflowOutflowBoundaryFlux(problem, scvf, fvGeometry, elemVolVars, elemFaceVars);
329 }
330 }
331 8610966 return result;
332 }
333
334 private:
335
336 //! do nothing if no turbulence model is used
337 template<class ...Args, bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<!turbulenceModel, int> = 0>
338 ✗ void incorporateWallFunction_(Args&&... args) const
339 ✗ {}
340
341 //! if a turbulence model is used, ask the problem is a wall function shall be employed and get the flux accordingly
342 template<bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<turbulenceModel, int> = 0>
343 2306480 void incorporateWallFunction_(CellCenterResidual& boundaryFlux,
344 const Problem& problem,
345 const Element& element,
346 const FVElementGeometry& fvGeometry,
347 const SubControlVolumeFace& scvf,
348 const ElementVolumeVariables& elemVolVars,
349 const ElementFaceVariables& elemFaceVars) const
350 {
351 static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
352 16185552 const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
353
354 // account for wall functions, if used
355
2/2
✓ Branch 0 taken 8304774 times.
✓ Branch 1 taken 2306480 times.
 16397550 for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
356 {
357 // use a wall function
358
2/2
✓ Branch 1 taken 662108 times.
✓ Branch 2 taken 7642666 times.
 8304774 if(problem.useWallFunction(element, scvf, eqIdx + cellCenterOffset))
359 {
360 1748276 boundaryFlux[eqIdx] = problem.wallFunction(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx]
361 1324216 * extrusionFactor * Extrusion::area(fvGeometry, scvf);
362 }
363 }
364 2306480 }
365
366 //! Returns the implementation of the problem (i.e. static polymorphism)
367 Implementation &asImp_()
368 { return *static_cast<Implementation *>(this); }
369
370 //! \copydoc asImp_()
371 const Implementation &asImp_() const
372 { return *static_cast<const Implementation *>(this); }
373 };
374
375 } // end namespace Dumux
376
377 #endif
378