GCC Code Coverage Report


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