GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 139 168 82.7%
Functions: 338 463 73.0%
Branches: 164 299 54.8%
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::NavierStokesFluxVariablesImpl
11 */
12 #ifndef DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH
13 #define DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH
14
15 #include <array>
16 #include <optional>
17 #include <type_traits>
18
19 #include <dumux/common/math.hh>
20 #include <dumux/common/exceptions.hh>
21 #include <dumux/common/parameters.hh>
22 #include <dumux/common/properties.hh>
23 #include <dumux/common/typetraits/problem.hh>
24
25 #include <dumux/flux/fluxvariablesbase.hh>
26 #include <dumux/discretization/method.hh>
27 #include <dumux/discretization/extrusion.hh>
28
29 #include "staggeredupwindhelper.hh"
30 #include "velocitygradients.hh"
31
32 namespace Dumux {
33
34 /*!
35 * \ingroup NavierStokesModel
36 * \brief The flux variables class for the Navier-Stokes model using the staggered grid discretization.
37 */
38
39 // forward declaration
40 template<class TypeTag, class DiscretizationMethod>
41 class NavierStokesFluxVariablesImpl;
42
43 template<class TypeTag>
44 class NavierStokesFluxVariablesImpl<TypeTag, DiscretizationMethods::Staggered>
45 : public FluxVariablesBase<GetPropType<TypeTag, Properties::Problem>,
46 typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView,
47 typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView,
48 typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView>
49 {
50 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
51
52 using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
53 using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
54 using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
55
56 using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
57 using FluxVariablesCache = typename GridFluxVariablesCache::FluxVariablesCache;
58
59 using GridFaceVariables = typename GridVariables::GridFaceVariables;
60 using ElementFaceVariables = typename GridFaceVariables::LocalView;
61 using FaceVariables = typename GridFaceVariables::FaceVariables;
62
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 Scalar = GetPropType<TypeTag, Properties::Scalar>;
69 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
70 using Indices = typename ModelTraits::Indices;
71 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
72 using FaceFrontalSubControlVolumeFace = typename GridGeometry::Traits::FaceFrontalSubControlVolumeFace;
73 using FaceLateralSubControlVolumeFace = typename GridGeometry::Traits::FaceLateralSubControlVolumeFace;
74 using Extrusion = Extrusion_t<GridGeometry>;
75 using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
76 using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
77 using BoundaryTypes = typename ProblemTraits<Problem>::BoundaryTypes;
78 using VelocityGradients = StaggeredVelocityGradients<Scalar, GridGeometry, BoundaryTypes, Indices>;
79
80 static constexpr bool normalizePressure = getPropValue<TypeTag, Properties::NormalizePressure>();
81
82 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
83
84 static constexpr auto cellCenterIdx = GridGeometry::cellCenterIdx();
85 static constexpr auto faceIdx = GridGeometry::faceIdx();
86
87 static constexpr int upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
88 static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
89
90 public:
91
92 using HeatConductionType = GetPropType<TypeTag, Properties::HeatConductionType>;
93
94
95 /*!
96 * \brief Returns the advective flux over a sub control volume face.
97 * \param problem The object specifying the problem which ought to be simulated
98 * \param fvGeometry The element-centered finite volume geometry view
99 * \param elemVolVars All volume variables for the element
100 * \param elemFaceVars The face variables
101 * \param scvf The sub control volume face
102 * \param upwindTerm The upwind term (i.e. the advectively transported quantity)
103 */
104 template<class UpwindTerm>
105 1037732436 static Scalar advectiveFluxForCellCenter(const Problem& problem,
106 const FVElementGeometry& fvGeometry,
107 const ElementVolumeVariables& elemVolVars,
108 const ElementFaceVariables& elemFaceVars,
109 const SubControlVolumeFace &scvf,
110 UpwindTerm upwindTerm)
111 {
112
2/2
✓ Branch 0 taken 111 times.
✓ Branch 1 taken 536756455 times.
 1037732436 const Scalar velocity = elemFaceVars[scvf].velocitySelf();
113
2/2
✓ Branch 0 taken 111 times.
✓ Branch 1 taken 536756455 times.
 1037732436 const bool insideIsUpstream = scvf.directionSign() == sign(velocity);
114
5/8
✓ Branch 0 taken 111 times.
✓ Branch 1 taken 536756455 times.
✓ Branch 3 taken 111 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 111 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 111 times.
✗ Branch 10 not taken.
 1037732436 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
115
116 2075464872 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
117 2075464872 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
118
119
2/2
✓ Branch 0 taken 275445746 times.
✓ Branch 1 taken 261310820 times.
 1037732436 const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
120
2/2
✓ Branch 0 taken 275445746 times.
✓ Branch 1 taken 261310820 times.
 1037732436 const auto& downstreamVolVars = insideIsUpstream ? outsideVolVars : insideVolVars;
121
122 1637180256 const Scalar flux = (upwindWeight * upwindTerm(upstreamVolVars) +
123 1237548376 (1.0 - upwindWeight) * upwindTerm(downstreamVolVars))
124 1237548376 * velocity * Extrusion::area(fvGeometry, scvf) * scvf.directionSign();
125
126 1037732436 return flux * extrusionFactor_(elemVolVars, scvf);
127 }
128
129 /*!
130 * \brief Computes the flux for the cell center residual (mass balance).
131 *
132 * \verbatim
133 * scvf
134 * ----------------
135 * | # current scvf # scvf over which fluxes are calculated
136 * | #
137 * | x #~~~~> vel.Self x dof position
138 * | #
139 * scvf | # -- element
140 * ----------------
141 * scvf
142 * \endverbatim
143 */
144 ✗ CellCenterPrimaryVariables computeMassFlux(const Problem& problem,
145 const Element& element,
146 const FVElementGeometry& fvGeometry,
147 const ElementVolumeVariables& elemVolVars,
148 const ElementFaceVariables& elemFaceVars,
149 const SubControlVolumeFace& scvf,
150 const FluxVariablesCache& fluxVarsCache)
151 {
152 // The advectively transported quantity (i.e density for a single-phase system).
153 322239468 auto upwindTerm = [](const auto& volVars) { return volVars.density(); };
154
155 // Call the generic flux function.
156
1/8
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 138882 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 100386724 CellCenterPrimaryVariables result(0.0);
157
5/10
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 138882 times.
✓ Branch 5 taken 140546 times.
✓ Branch 6 taken 41192 times.
✓ Branch 7 taken 97690 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 8126 times.
✗ Branch 10 not taken.
 100386724 result[Indices::conti0EqIdx - ModelTraits::dim()] = advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTerm);
158
159 ✗ return result;
160 }
161
162 /*!
163 * \brief Returns the momentum flux over all staggered faces.
164 */
165 47252532 FacePrimaryVariables computeMomentumFlux(const Problem& problem,
166 const Element& element,
167 const SubControlVolumeFace& scvf,
168 const FVElementGeometry& fvGeometry,
169 const ElementVolumeVariables& elemVolVars,
170 const ElementFaceVariables& elemFaceVars,
171 const GridFluxVariablesCache& gridFluxVarsCache)
172 {
173 47252532 return computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) +
174 47252532 computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache);
175 }
176
177 /*!
178 * \brief Returns the frontal part of the momentum flux.
179 * This treats the flux over the staggered face at the center of an element,
180 * parallel to the current scvf where the velocity dof of interest lives.
181 *
182 * \verbatim
183 * scvf
184 * ---------======= == and # staggered half-control-volume
185 * | # | current scvf
186 * | # | # staggered face over which fluxes are calculated
187 * vel.Opp <~~| #~~> x~~~~> vel.Self
188 * | # | x dof position
189 * scvf | # |
190 * --------======== -- element
191 * scvf
192 * \endverbatim
193 */
194 144714524 FacePrimaryVariables computeFrontalMomentumFlux(const Problem& problem,
195 const Element& element,
196 const SubControlVolumeFace& scvf,
197 const FVElementGeometry& fvGeometry,
198 const ElementVolumeVariables& elemVolVars,
199 const ElementFaceVariables& elemFaceVars,
200 const GridFluxVariablesCache& gridFluxVarsCache)
201 {
202
2/2
✓ Branch 0 taken 134869060 times.
✓ Branch 1 taken 9845464 times.
 144714524 FacePrimaryVariables frontalFlux(0.0);
203
204 // The velocities of the dof at interest and the one of the opposite scvf.
205
2/2
✓ Branch 0 taken 134869060 times.
✓ Branch 1 taken 9845464 times.
 144714524 const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf();
206
2/2
✓ Branch 0 taken 134869060 times.
✓ Branch 1 taken 9845464 times.
 144714524 const Scalar velocityOpposite = elemFaceVars[scvf].velocityOpposite();
207
2/2
✓ Branch 0 taken 134869060 times.
✓ Branch 1 taken 9845464 times.
 144714524 const auto& faceVars = elemFaceVars[scvf];
208
209 // Advective flux.
210
2/2
✓ Branch 0 taken 134869060 times.
✓ Branch 1 taken 9845464 times.
 144714524 if (problem.enableInertiaTerms())
211 {
212 // Get the average velocity at the center of the element (i.e. the location of the staggered face).
213 134869060 const Scalar transportingVelocity = (velocitySelf + velocityOpposite) * 0.5;
214 134869060 const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
215
216 269738120 StaggeredUpwindHelper<TypeTag, upwindSchemeOrder> upwindHelper(element, fvGeometry, scvf, elemFaceVars, elemVolVars, gridFluxVarsCache.staggeredUpwindMethods());
217 269738120 frontalFlux += upwindHelper.computeUpwindFrontalMomentum(selfIsUpstream)
218 134869060 * transportingVelocity * -1.0 * scvf.directionSign();
219 }
220
221 // The volume variables within the current element. We only require those (and none of neighboring elements)
222 // because the fluxes are calculated over the staggered face at the center of the element.
223 289429048 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
224
225 // Diffusive flux.
226
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 144714474 times.
 144714524 const Scalar velocityGrad_ii = VelocityGradients::velocityGradII(scvf, faceVars) * scvf.directionSign();
227
228
3/4
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 144714474 times.
✓ Branch 3 taken 50 times.
✗ Branch 4 not taken.
 144714574 static const bool enableUnsymmetrizedVelocityGradient
229
2/4
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
 150 = getParamFromGroup<bool>(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
230
1/2
✓ Branch 0 taken 144714524 times.
✗ Branch 1 not taken.
 144714524 const Scalar factor = enableUnsymmetrizedVelocityGradient ? 1.0 : 2.0;
231 289429048 frontalFlux -= factor * insideVolVars.effectiveViscosity() * velocityGrad_ii;
232
233 // The pressure term.
234 // If specified, the pressure can be normalized using the initial value on the scfv of interest.
235 // The scvf is used to normalize by the same value from the left and right side.
236 // Can potentially help to improve the condition number of the system matrix.
237 144714524 const Scalar pressure = normalizePressure ?
238 347738314 insideVolVars.pressure() - problem.initial(scvf)[Indices::pressureIdx]
239 : insideVolVars.pressure();
240
241 // Account for the orientation of the staggered face's normal outer normal vector
242 // (pointing in opposite direction of the scvf's one).
243 144714524 frontalFlux += pressure * -1.0 * scvf.directionSign();
244
245 // Account for the staggered face's area. For rectangular elements, this equals the area of the scvf
246 // our velocity dof of interest lives on but with adjusted centroid
247 289429048 const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
248 289429048 FaceFrontalSubControlVolumeFace frontalFace(scv.center(), scvf.area());
249 578858096 return frontalFlux * Extrusion::area(fvGeometry, frontalFace) * insideVolVars.extrusionFactor();
250 }
251
252 /*!
253 * \brief Returns the momentum flux over the staggered faces
254 * perpendicular to the scvf where the velocity dof of interest
255 * lives (coinciding with the element's scvfs).
256 *
257 * \verbatim
258 * scvf
259 * ---------####### || and # staggered half-control-volume
260 * | || | current scvf
261 * | || | # normal staggered sub faces over which fluxes are calculated
262 * | || x~~~~> vel.Self
263 * | || | x dof position
264 * scvf | || |
265 * --------######## -- element
266 * scvf
267 * \endverbatim
268 */
269 144714524 FacePrimaryVariables computeLateralMomentumFlux(const Problem& problem,
270 const Element& element,
271 const SubControlVolumeFace& scvf,
272 const FVElementGeometry& fvGeometry,
273 const ElementVolumeVariables& elemVolVars,
274 const ElementFaceVariables& elemFaceVars,
275 const GridFluxVariablesCache& gridFluxVarsCache)
276 {
277
2/2
✓ Branch 0 taken 2787742 times.
✓ Branch 1 taken 141926782 times.
 144714524 FacePrimaryVariables lateralFlux(0.0);
278
2/2
✓ Branch 0 taken 2787742 times.
✓ Branch 1 taken 141926782 times.
 144714524 const auto& faceVars = elemFaceVars[scvf];
279
2/2
✓ Branch 0 taken 2787742 times.
✓ Branch 1 taken 141926782 times.
 144714524 const std::size_t numSubFaces = scvf.pairData().size();
280
281 // If the current scvf is on a boundary, check if there is a Neumann BC for the stress in tangential direction.
282 // Create a boundaryTypes object (will be empty if not at a boundary).
283 144714524 std::optional<BoundaryTypes> currentScvfBoundaryTypes;
284
2/2
✓ Branch 0 taken 2787742 times.
✓ Branch 1 taken 141926782 times.
 144714524 if (scvf.boundary())
285 2787742 currentScvfBoundaryTypes.emplace(problem.boundaryTypes(element, scvf));
286
287 // Account for all sub faces.
288
2/2
✓ Branch 0 taken 289429048 times.
✓ Branch 1 taken 144714524 times.
 434143572 for (int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
289 {
290
2/2
✓ Branch 0 taken 14551668 times.
✓ Branch 1 taken 274877380 times.
 289429048 const auto eIdx = scvf.insideScvIdx();
291 // Get the face normal to the face the dof lives on. The staggered sub face coincides with half of this lateral face.
292
4/4
✓ Branch 0 taken 14551668 times.
✓ Branch 1 taken 274877380 times.
✓ Branch 2 taken 14551668 times.
✓ Branch 3 taken 274877380 times.
 578858096 const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
293
294 // Create a boundaryTypes object (will be empty if not at a boundary).
295 289429048 std::optional<BoundaryTypes> lateralFaceBoundaryTypes;
296
297 // Check if there is face/element parallel to our face of interest where the dof lives on. If there is no parallel neighbor,
298 // we are on a boundary where we have to check for boundary conditions.
299
2/2
✓ Branch 0 taken 14551668 times.
✓ Branch 1 taken 274877380 times.
 289429048 if (lateralFace.boundary())
300 {
301 // Retrieve the boundary types that correspond to the center of the lateral scvf. As a convention, we always query
302 // the type of BCs at the center of the element's "actual" lateral scvf (not the face of the staggered control volume).
303 // The value of the BC will be evaluated at the center of the staggered face.
304 // --------T######V || frontal face of staggered half-control-volume
305 // | || | current scvf # lateral staggered face of interest (may lie on a boundary)
306 // | || | x dof position
307 // | || x~~~~> vel.Self -- element boundaries
308 // | || | T position at which the type of boundary conditions will be evaluated
309 // | || | (center of lateral scvf)
310 // ---------------- V position at which the value of the boundary conditions will be evaluated
311 // (center of the staggered lateral face)
312 14551668 lateralFaceBoundaryTypes.emplace(problem.boundaryTypes(element, lateralFace));
313 }
314
315 // If the current scvf is on a boundary and if there is a Neumann or Beavers-Joseph-(Saffmann) BC for the stress in tangential direction,
316 // assign this value for the lateral momentum flux. No further calculations are required. We assume that all lateral faces
317 // have the same type of BC (Neumann or Beavers-Joseph-(Saffmann)), but we sample the value at their actual positions.
318
4/4
✓ Branch 0 taken 5575484 times.
✓ Branch 1 taken 283853564 times.
✓ Branch 2 taken 5575484 times.
✓ Branch 3 taken 283853564 times.
 578858096 if (currentScvfBoundaryTypes)
319 {
320 // Handle Neumann BCs.
321
2/4
✓ Branch 0 taken 5575484 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5575484 times.
✗ Branch 3 not taken.
 11150968 if (currentScvfBoundaryTypes->isNeumann(Indices::velocity(lateralFace.directionIndex())))
322 {
323 // Get the location of the lateral staggered face's center.
324 ✗ FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area());
325 ✗ const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
326 ✗ const auto lateralBoundaryFace = makeStaggeredBoundaryFace(scvf, lateralStaggeredFaceCenter);
327 ✗ lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, lateralBoundaryFace)[Indices::velocity(lateralFace.directionIndex())]
328 ✗ * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(fvGeometry, lateralScvf) * lateralFace.directionSign();
329
330 continue;
331 }
332
333 // Treat the edge case when our current scvf has a BJ condition while the lateral face has a Neumann BC.
334 // It is not clear how to evaluate the BJ condition here.
335 // For symmetry reasons, our own scvf should then have the same Neumann flux as the lateral face.
336 // TODO: We should clarify if this is the correct approach.
337
10/10
✓ Branch 0 taken 346528 times.
✓ Branch 1 taken 5228956 times.
✓ Branch 2 taken 346528 times.
✓ Branch 3 taken 5228956 times.
✓ Branch 4 taken 346528 times.
✓ Branch 5 taken 5228956 times.
✓ Branch 6 taken 18242 times.
✓ Branch 7 taken 328286 times.
✓ Branch 8 taken 18242 times.
✓ Branch 9 taken 328286 times.
 16726452 if (currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralFace.directionIndex())) && lateralFaceBoundaryTypes &&
338
2/4
✓ Branch 0 taken 18242 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 18242 times.
✗ Branch 3 not taken.
 36484 lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex())))
339 {
340 ✗ FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area());
341 ✗ const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
342 ✗ const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralFace, lateralStaggeredFaceCenter);
343 ✗ lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())]
344 ✗ * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(fvGeometry, lateralScvf) * lateralFace.directionSign();
345
346 continue;
347 }
348 }
349
350 // Check if the lateral face (perpendicular to our current scvf) lies on a boundary. If yes, boundary conditions might need to be treated
351 // and further calculations can be skipped.
352
2/2
✓ Branch 0 taken 14551668 times.
✓ Branch 1 taken 274877380 times.
 289429048 if (lateralFace.boundary())
353 {
354 // Check if we have a symmetry boundary condition. If yes, the tangential part of the momentum flux can be neglected
355 // and we may skip any further calculations for the given sub face.
356
6/6
✓ Branch 0 taken 13735545 times.
✓ Branch 1 taken 816123 times.
✓ Branch 2 taken 13735545 times.
✓ Branch 3 taken 816123 times.
✓ Branch 4 taken 13735545 times.
✓ Branch 5 taken 816123 times.
 43655004 if (lateralFaceBoundaryTypes->isSymmetry())
357 continue;
358
359 // Handle Neumann boundary conditions. No further calculations are then required for the given sub face.
360
2/4
✓ Branch 0 taken 13735545 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 13735545 times.
✗ Branch 3 not taken.
 27471090 if (lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex())))
361 {
362 // Construct a temporary scvf which corresponds to the staggered sub face, featuring the location
363 // the staggered faces's center.
364 ✗ FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area());
365 ✗ const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
366 ✗ const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralFace, lateralStaggeredFaceCenter);
367 ✗ lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())]
368 ✗ * elemVolVars[lateralFace.insideScvIdx()].extrusionFactor() * Extrusion::area(fvGeometry, lateralScvf);
369
370 continue;
371 }
372
373 // Handle wall-function fluxes (only required for RANS models)
374
2/2
✓ Branch 1 taken 3007923 times.
✓ Branch 2 taken 508848 times.
 3516771 if (incorporateWallFunction_(lateralFlux, problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars, localSubFaceIdx))
375 continue;
376 }
377
378 // Check the consistency of the boundary conditions, exactly one of the following must be set
379
4/4
✓ Branch 0 taken 13226697 times.
✓ Branch 1 taken 274877380 times.
✓ Branch 2 taken 13226697 times.
✓ Branch 3 taken 274877380 times.
 576208154 if (lateralFaceBoundaryTypes)
380 {
381 13226697 std::bitset<3> admittableBcTypes;
382
4/4
✓ Branch 0 taken 8854829 times.
✓ Branch 1 taken 4371868 times.
✓ Branch 2 taken 8854829 times.
✓ Branch 3 taken 4371868 times.
 26453394 admittableBcTypes.set(0, lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx));
383
4/4
✓ Branch 0 taken 4627112 times.
✓ Branch 1 taken 8599585 times.
✓ Branch 2 taken 4627112 times.
✓ Branch 3 taken 8599585 times.
 31080506 admittableBcTypes.set(1, lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex())));
384
6/6
✓ Branch 0 taken 255244 times.
✓ Branch 1 taken 12971453 times.
✓ Branch 2 taken 255244 times.
✓ Branch 3 taken 12971453 times.
✓ Branch 4 taken 255244 times.
✓ Branch 5 taken 12971453 times.
 39680091 admittableBcTypes.set(2, lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())));
385
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 13226697 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13226697 times.
 26453394 if (admittableBcTypes.count() != 1)
386 {
387 ✗ DUNE_THROW(Dune::InvalidStateException, "Invalid boundary conditions for lateral scvf "
388 "for the momentum equations at global position " << lateralStaggeredFaceCenter_(scvf, localSubFaceIdx)
389 << ", current scvf global position " << scvf.center());
390 }
391 }
392
393 // If none of the above boundary conditions apply for the given sub face, proceed to calculate the tangential momentum flux.
394
2/2
✓ Branch 0 taken 268413149 times.
✓ Branch 1 taken 19690928 times.
 288104077 if (problem.enableInertiaTerms())
395
1/2
✓ Branch 1 taken 29634048 times.
✗ Branch 2 not taken.
 268413149 lateralFlux += computeAdvectivePartOfLateralMomentumFlux_(problem, fvGeometry, element,
396 scvf, elemVolVars, elemFaceVars,
397 gridFluxVarsCache,
398 currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
399 localSubFaceIdx);
400
401
1/2
✓ Branch 1 taken 49324976 times.
✗ Branch 2 not taken.
 288104077 lateralFlux += computeDiffusivePartOfLateralMomentumFlux_(problem, fvGeometry, element,
402 scvf, elemVolVars, faceVars,
403 currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
404 localSubFaceIdx);
405 }
406 144714524 return lateralFlux;
407 }
408
409
410 /*!
411 * \brief Returns the momentum flux over an inflow or outflow boundary face.
412 *
413 * \verbatim
414 * scvf //
415 * ---------=======// == and # staggered half-control-volume
416 * | || #// current scvf
417 * | || #// # staggered boundary face over which fluxes are calculated
418 * | || x~~~~> vel.Self
419 * | || #// x dof position
420 * scvf | || #//
421 * --------========// -- element
422 * scvf //
423 * // boundary
424 * \endverbatim
425 */
426 2614478 FacePrimaryVariables inflowOutflowBoundaryFlux(const Problem& problem,
427 const SubControlVolumeFace& scvf,
428 const FVElementGeometry& fvGeometry,
429 const ElementVolumeVariables& elemVolVars,
430 const ElementFaceVariables& elemFaceVars) const
431 {
432 2614478 FacePrimaryVariables inOrOutflow(0.0);
433 2614478 const auto& element = fvGeometry.element();
434 5228956 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
435
436 // Advective momentum flux.
437
2/2
✓ Branch 0 taken 2543714 times.
✓ Branch 1 taken 70764 times.
 2614478 if (problem.enableInertiaTerms())
438 {
439 2543714 const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf();
440 5087428 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
441
4/4
✓ Branch 0 taken 18215 times.
✓ Branch 1 taken 2525499 times.
✓ Branch 2 taken 18215 times.
✓ Branch 3 taken 2525499 times.
 5087428 const auto& upVolVars = (scvf.directionSign() == sign(velocitySelf)) ?
442 insideVolVars : outsideVolVars;
443
444 5087428 inOrOutflow += velocitySelf * velocitySelf * upVolVars.density();
445 }
446
447 // Apply a pressure at the boundary.
448 2614478 const Scalar boundaryPressure = normalizePressure
449 2688922 ? (problem.dirichlet(element, scvf)[Indices::pressureIdx] -
450 2617454 problem.initial(scvf)[Indices::pressureIdx])
451 : problem.dirichlet(element, scvf)[Indices::pressureIdx];
452 5228956 inOrOutflow += boundaryPressure;
453
454 // Account for the orientation of the face at the boundary,
455 5228956 return inOrOutflow * scvf.directionSign() * Extrusion::area(fvGeometry, scvf) * insideVolVars.extrusionFactor();
456 }
457
458 private:
459
460 /*!
461 * \brief Returns the advective momentum flux over the staggered face perpendicular to the scvf
462 * where the velocity dof of interest lives (coinciding with the element's scvfs).
463 *
464 * \verbatim
465 * ----------------
466 * | |
467 * | transp. |
468 * | vel. |~~~~> vel.Parallel
469 * | ^ |
470 * | | |
471 * scvf ---------####### || and # staggered half-control-volume
472 * | || | current scvf
473 * | || | # normal staggered faces over which fluxes are calculated
474 * | || x~~~~> vel.Self
475 * | || | x dof position
476 * scvf | || |
477 * ---------####### -- elements
478 * scvf
479 * \endverbatim
480 */
481 268413149 FacePrimaryVariables computeAdvectivePartOfLateralMomentumFlux_(const Problem& problem,
482 const FVElementGeometry& fvGeometry,
483 const Element& element,
484 const SubControlVolumeFace& scvf,
485 const ElementVolumeVariables& elemVolVars,
486 const ElementFaceVariables& elemFaceVars,
487 const GridFluxVariablesCache& gridFluxVarsCache,
488 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
489 const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
490 const int localSubFaceIdx)
491 {
492 268413149 const auto eIdx = scvf.insideScvIdx();
493 536826298 const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
494
495 // Get the transporting velocity, located at the scvf perpendicular to the current scvf where the dof
496 // of interest is located.
497 536826298 const Scalar transportingVelocity = [&]()
498 {
499
2/2
✓ Branch 0 taken 263224686 times.
✓ Branch 1 taken 5188463 times.
 268413149 const auto& faceVars = elemFaceVars[scvf];
500
2/2
✓ Branch 0 taken 263224686 times.
✓ Branch 1 taken 5188463 times.
 268413149 if (!scvf.boundary())
501 263224686 return faceVars.velocityLateralInside(localSubFaceIdx);
502 else
503 {
504 // Create a boundaryTypes object. Get the boundary conditions. We sample the type of BC at the center of the current scvf.
505 5188463 const auto bcTypes = problem.boundaryTypes(element, scvf);
506
507
1/2
✓ Branch 0 taken 5188463 times.
✗ Branch 1 not taken.
 5188463 if (bcTypes.isDirichlet(Indices::velocity(lateralFace.directionIndex())))
508 {
509 // Construct a temporary scvf which corresponds to the staggered sub face, featuring the location
510 // the staggered faces's center.
511 ✗ const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
512 ✗ const auto lateralBoundaryFace = makeStaggeredBoundaryFace(scvf, lateralBoundaryFacePos);
513 ✗ return problem.dirichlet(element, lateralBoundaryFace)[Indices::velocity(lateralFace.directionIndex())];
514 }
515
4/4
✓ Branch 0 taken 149608 times.
✓ Branch 1 taken 5038855 times.
✓ Branch 2 taken 149608 times.
✓ Branch 3 taken 5038855 times.
 10376926 else if (bcTypes.isBeaversJoseph(Indices::velocity(lateralFace.directionIndex())))
516 {
517 448824 return VelocityGradients::beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars,
518 149608 currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
519 }
520 else
521 5038855 return faceVars.velocityLateralInside(localSubFaceIdx);
522 }
523 536826298 }();
524
525 268413149 const bool selfIsUpstream = lateralFace.directionSign() == sign(transportingVelocity);
526 536826298 StaggeredUpwindHelper<TypeTag, upwindSchemeOrder> upwindHelper(element, fvGeometry, scvf, elemFaceVars, elemVolVars, gridFluxVarsCache.staggeredUpwindMethods());
527 268413149 FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area());
528 return upwindHelper.computeUpwindLateralMomentum(selfIsUpstream, lateralFace, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes)
529 268413149 * transportingVelocity * lateralFace.directionSign() * Extrusion::area(fvGeometry, lateralScvf) * extrusionFactor_(elemVolVars, lateralFace);
530 }
531
532 /*!
533 * \brief Returns the diffusive momentum flux over the staggered face perpendicular to the scvf
534 * where the velocity dof of interest lives (coinciding with the element's scvfs).
535 *
536 * \verbatim
537 * ----------------
538 * | |vel.
539 * | in.norm. |Parallel
540 * | vel. |~~~~>
541 * | ^ | ^ out.norm.vel.
542 * | | | |
543 * scvf ---------#######::::::::: || and # staggered half-control-volume (own element)
544 * | || | curr. ::
545 * | || | scvf :: :: staggered half-control-volume (neighbor element)
546 * | || x~~~~> ::
547 * | || | vel. :: # lateral staggered faces over which fluxes are calculated
548 * scvf | || | Self ::
549 * ---------#######::::::::: x dof position
550 * scvf
551 * -- elements
552 * \endverbatim
553 */
554 288104077 FacePrimaryVariables computeDiffusivePartOfLateralMomentumFlux_(const Problem& problem,
555 const FVElementGeometry& fvGeometry,
556 const Element& element,
557 const SubControlVolumeFace& scvf,
558 const ElementVolumeVariables& elemVolVars,
559 const FaceVariables& faceVars,
560 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
561 const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
562 const int localSubFaceIdx)
563 {
564
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 288104027 times.
 288104077 const auto eIdx = scvf.insideScvIdx();
565
4/4
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 288104027 times.
✓ Branch 2 taken 50 times.
✓ Branch 3 taken 288104027 times.
 576208154 const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
566
567
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 288104027 times.
 288104077 FacePrimaryVariables lateralDiffusiveFlux(0.0);
568
569
3/4
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 288104027 times.
✓ Branch 3 taken 50 times.
✗ Branch 4 not taken.
 288104127 static const bool enableUnsymmetrizedVelocityGradient
570
2/4
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
 150 = getParamFromGroup<bool>(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
571
572 // Get the volume variables of the own and the neighboring element. The neighboring
573 // element is adjacent to the staggered face normal to the current scvf
574 // where the dof of interest is located.
575 576208154 const auto& insideVolVars = elemVolVars[lateralFace.insideScvIdx()];
576 576208154 const auto& outsideVolVars = elemVolVars[lateralFace.outsideScvIdx()];
577
578 // Get the averaged viscosity at the staggered face normal to the current scvf.
579 288104077 const Scalar muAvg = lateralFace.boundary()
580
2/2
✓ Branch 0 taken 13226697 times.
✓ Branch 1 taken 274877380 times.
 288104077 ? insideVolVars.effectiveViscosity()
581 824632140 : (insideVolVars.effectiveViscosity() + outsideVolVars.effectiveViscosity()) * 0.5;
582
583 // Consider the shear stress caused by the gradient of the velocities normal to our face of interest.
584
1/2
✓ Branch 0 taken 288104077 times.
✗ Branch 1 not taken.
 288104077 if (!enableUnsymmetrizedVelocityGradient)
585 {
586 304684810 if (!scvf.boundary() ||
587
8/10
✓ Branch 0 taken 5526911 times.
✓ Branch 1 taken 282577166 times.
✓ Branch 2 taken 5526911 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 5526911 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 346528 times.
✓ Branch 7 taken 5180383 times.
✓ Branch 8 taken 346528 times.
✓ Branch 9 taken 234936 times.
 294212452 currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralFace.directionIndex())) ||
588
3/4
✓ Branch 0 taken 346528 times.
✓ Branch 1 taken 5180383 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4945447 times.
 10472358 currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralFace.directionIndex())))
589 {
590 282923694 const Scalar velocityGrad_ji = VelocityGradients::velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
591 // Account for the orientation of the staggered normal face's outer normal vector.
592 282923694 lateralDiffusiveFlux -= muAvg * velocityGrad_ji * lateralFace.directionSign();
593 }
594 }
595
596 // Consider the shear stress caused by the gradient of the velocities parallel to our face of interest.
597 // If we have a Dirichlet condition for the pressure at the lateral face we assume to have a zero velocityGrad_ij velocity gradient
598 // so we can skip the computation.
599
6/6
✓ Branch 0 taken 13226697 times.
✓ Branch 1 taken 274877380 times.
✓ Branch 2 taken 8854829 times.
✓ Branch 3 taken 4371868 times.
✓ Branch 4 taken 8854829 times.
✓ Branch 5 taken 4371868 times.
 288104077 if (!lateralFace.boundary() || !lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx))
600 {
601 283732209 const Scalar velocityGrad_ij = VelocityGradients::velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
602 283732209 lateralDiffusiveFlux -= muAvg * velocityGrad_ij * lateralFace.directionSign();
603 }
604
605 // Account for the area of the staggered lateral face (0.5 of the coinciding scfv).
606 288104077 FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area());
607 288104077 return lateralDiffusiveFlux * Extrusion::area(fvGeometry, lateralScvf) * extrusionFactor_(elemVolVars, lateralFace);
608 }
609
610 /*!
611 * \brief Get the location of the lateral staggered face's center.
612 * Only needed for boundary conditions if the current scvf or the lateral one is on a boundary.
613 *
614 * \verbatim
615 * --------#######o || frontal face of staggered half-control-volume
616 * | || | current scvf # lateral staggered face of interest (may lie on a boundary)
617 * | || | x dof position
618 * | || x~~~~> vel.Self -- element boundaries, current scvf may lie on a boundary
619 * | || | o position at which the boundary conditions will be evaluated
620 * | || | (lateralStaggeredFaceCenter)
621 * ----------------
622 * \endverbatim
623 */
624 ✗ const GlobalPosition& lateralStaggeredFaceCenter_(const SubControlVolumeFace& scvf, const int localSubFaceIdx) const
625 {
626 ✗ return scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter;
627 };
628
629 /*!
630 * \brief Get the location of the lateral staggered sub control volume face center.
631 * Only needed for boundary conditions if the current scvf or the lateral one is on a boundary.
632 *
633 * \verbatim
634 * --------###o#### || frontal face of staggered half-control-volume
635 * | || | current scvf # lateral staggered face of interest (may lie on a boundary)
636 * | || | x dof position
637 * | || x~~~~> vel.Self -- element boundaries, current scvf may lie on a boundary
638 * | || | o center of the lateral staggered scvf
639 * | || |
640 * ----------------
641 * \endverbatim
642 */
643 ✗ GlobalPosition lateralStaggeredSCVFCenter_(const SubControlVolumeFace& lateralFace,
644 const SubControlVolumeFace& currentFace,
645 const int localSubFaceIdx) const
646 {
647 ✗ auto pos = lateralStaggeredFaceCenter_(currentFace, localSubFaceIdx) + lateralFace.center();
648 ✗ pos *= 0.5;
649 ✗ return pos;
650 }
651
652 //! helper function to get the averaged extrusion factor for a face
653 1093782640 static Scalar extrusionFactor_(const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf)
654 {
655 2187565280 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
656 2187565280 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
657
1/2
✓ Branch 0 taken 1093782640 times.
✗ Branch 1 not taken.
 1093782640 return harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor());
658 }
659
660 //! do nothing if no turbulence model is used
661 template<class ...Args, bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<!turbulenceModel, int> = 0>
662 ✗ bool incorporateWallFunction_(Args&&... args) const
663 ✗ { return false; }
664
665 //! if a turbulence model is used, ask the problem is a wall function shall be employed and get the flux accordingly
666 template<bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<turbulenceModel, int> = 0>
667 3516771 bool incorporateWallFunction_(FacePrimaryVariables& lateralFlux,
668 const Problem& problem,
669 const Element& element,
670 const FVElementGeometry& fvGeometry,
671 const SubControlVolumeFace& scvf,
672 const ElementVolumeVariables& elemVolVars,
673 const ElementFaceVariables& elemFaceVars,
674 const std::size_t localSubFaceIdx) const
675 {
676 3516771 const auto eIdx = scvf.insideScvIdx();
677 7033542 const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
678
679
2/2
✓ Branch 1 taken 508848 times.
✓ Branch 2 taken 3007923 times.
 3516771 if (problem.useWallFunction(element, lateralFace, Indices::velocity(scvf.directionIndex())))
680 {
681 508848 FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area());
682
0/2
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1526544 const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralFace, lateralStaggeredFaceCenter_(scvf, localSubFaceIdx));
683
1/2
✓ Branch 1 taken 508848 times.
✗ Branch 2 not taken.
 508848 lateralFlux += problem.wallFunction(element, fvGeometry, elemVolVars, elemFaceVars, scvf, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())]
684
2/4
✓ Branch 1 taken 508848 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 508848 times.
✗ Branch 4 not taken.
 508848 * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(fvGeometry, lateralScvf);
685
1/2
✓ Branch 0 taken 508848 times.
✗ Branch 1 not taken.
 508848 return true;
686 }
687 else
688 return false;
689 }
690 };
691
692 } // end namespace Dumux
693
694 #endif
695