GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/freeflow/navierstokes/staggered/velocitygradients.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 71 78 91.0%
Functions: 156 264 59.1%
Branches: 53 122 43.4%

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::StaggeredVelocityGradients
11 */
12 #ifndef DUMUX_NAVIERSTOKES_STAGGERED_VELOCITYGRADIENTS_HH
13 #define DUMUX_NAVIERSTOKES_STAGGERED_VELOCITYGRADIENTS_HH
14
15 #include <optional>
16 #include <dumux/common/exceptions.hh>
17 #include <dumux/common/parameters.hh>
18
19 namespace Dumux {
20
21 /*!
22 * \ingroup NavierStokesModel
23 * \brief Helper class for calculating the velocity gradients for the Navier-Stokes model using the staggered grid discretization.
24 */
25 template<class Scalar, class GridGeometry, class BoundaryTypes, class Indices>
26 class StaggeredVelocityGradients
27 {
28 using FVElementGeometry = typename GridGeometry::LocalView;
29 using GridView = typename GridGeometry::GridView;
30 using Element = typename GridView::template Codim<0>::Entity;
31 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
32 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
33
34 public:
35
36 /*!
37 * \brief Returns the in-axis velocity gradient.
38 *
39 * \verbatim
40 * ---------======= == and # staggered half-control-volume
41 * | # | current scvf
42 * | # | # staggered face over which fluxes are calculated
43 * vel.Opp <~~| O~~> x~~~~> vel.Self
44 * | # | x dof position
45 * scvf | # |
46 * --------======== -- element
47 *
48 * O position at which gradient is evaluated
49 * \endverbatim
50 */
51 template<class FaceVariables>
52
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1947594 times.
145578344 static Scalar velocityGradII(const SubControlVolumeFace& scvf,
53 const FaceVariables& faceVars)
54 {
55 // The velocities of the dof at interest and the one of the opposite scvf.
56
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1947594 times.
145578344 const Scalar velocitySelf = faceVars.velocitySelf();
57 145578344 const Scalar velocityOpposite = faceVars.velocityOpposite();
58
59
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 145578293 times.
145578344 return ((velocityOpposite - velocitySelf) / scvf.selfToOppositeDistance()) * scvf.directionSign();
60 }
61
62 /*!
63 * \brief Returns the velocity gradient perpendicular to the orientation of our current scvf.
64 *
65 * \verbatim
66 * ----------------
67 * | |vel.
68 * | |Parallel
69 * | |~~~~> ------->
70 * | | ------> * gradient
71 * | | ----->
72 * scvf ---------######O::::::::: ----> || and # staggered half-control-volume (own element)
73 * | || | curr. :: --->
74 * | || | scvf :: --> :: staggered half-control-volume (neighbor element)
75 * | || x~~~~> :: ->
76 * | || | vel. :: # lateral staggered faces over which fluxes are calculated
77 * scvf | || | Self ::
78 * ---------#######::::::::: x dof position
79 * scvf
80 * -- elements
81 *
82 * O position at which gradient is evaluated
83 * \endverbatim
84 *
85 * ------------
86 * | xxxx s
87 * | xxxx a
88 * | xxxx s
89 * -----------O-----------
90 * | yyyy s zzzz |
91 * | yyyy b zzzz |
92 * | yyyy s zzzz |
93 * -----------------------
94 *
95 * In a corner geometry (scvf is sas or sbs), we calculate the velocity gradient at O, by
96 * (velocity(a)-velocity(b))/distance(a,b) for the half-control volumes x and y, but by
97 * (velocity(O)-velocity(b))/distance(O,b) for z. This does not harm flux continuity (x and y use the same
98 * formulation). We do this different formulation for y (approximate gradient by central differncing) and
99 * z (approximate gradient by forward/backward differencing), because it is the natural way of implementing
100 * it and it is not clear which gradient is the better approximation in this case anyway.
101 * Particularly, for the frequent case of no-slip, no-flow boundaries, the velocity would be zero at O and a
102 * and thus, the gradient within the flow domain might be better approximated by velocity(b)/distanc(O,b)
103 * than by velocity(b)/distance(a,b).
104 */
105 template<class Problem, class FaceVariables>
106 285829870 static Scalar velocityGradIJ(const Problem& problem,
107 const Element& element,
108 const FVElementGeometry& fvGeometry,
109 const SubControlVolumeFace& scvf,
110 const FaceVariables& faceVars,
111 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
112 const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
113 const std::size_t localSubFaceIdx)
114 {
115 285829870 const auto eIdx = scvf.insideScvIdx();
116 285829870 const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
117
118 // For the velocityGrad_ij derivative, get the velocities at the current (own) scvf
119 // and at the parallel one at the neighboring scvf.
120 285829870 const Scalar innerParallelVelocity = faceVars.velocitySelf();
121
122 857489610 const auto outerParallelVelocity = [&]()
123 {
124
2/2
✓ Branch 0 taken 276911036 times.
✓ Branch 1 taken 8918834 times.
285829870 if (!lateralScvf.boundary())
125 276911036 return faceVars.velocityParallel(localSubFaceIdx, 0);
126
3/4
✓ Branch 0 taken 255422 times.
✓ Branch 1 taken 8663412 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 255422 times.
8918834 else if (lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex())))
127 {
128 // Sample the value of the Dirichlet BC at the center of the staggered lateral face.
129 8663412 const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
130 17326824 const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralScvf, lateralBoundaryFacePos);
131
2/3
✓ Branch 1 taken 478856 times.
✗ Branch 2 not taken.
✓ Branch 0 taken 11928 times.
8982176 return problem.dirichlet(element, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())];
132 8663412 }
133
1/2
✓ Branch 0 taken 255422 times.
✗ Branch 1 not taken.
255422 else if (lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())))
134 {
135 255422 return beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars,
136 255422 currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
137 }
138 else
139 DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary type at " << lateralScvf.center());
140 285829870 }();
141
142 // The velocity gradient already accounts for the orientation
143 // of the staggered face's outer normal vector. This also correctly accounts for the reduced
144 // distance used in the gradient if the lateral scvf lies on a boundary.
145 285829870 return (outerParallelVelocity - innerParallelVelocity)
146 285829870 / scvf.parallelDofsDistance(localSubFaceIdx, 0) * lateralScvf.directionSign();
147 }
148
149 /*!
150 * \brief Returns the velocity gradient in line with our current scvf.
151 *
152 * \verbatim
153 * ^ gradient
154 * | ^
155 * | | ^
156 * | | | ^
157 * | | | | ^
158 * | | | | | ^
159 * | | | | | |
160 *
161 * ----------------
162 * | |
163 * | in.norm. |
164 * | vel. |
165 * | ^ | ^ out.norm.vel.
166 * | | | |
167 * scvf ---------######O::::::::: || and # staggered half-control-volume (own element)
168 * | || | curr. ::
169 * | || | scvf :: :: staggered half-control-volume (neighbor element)
170 * | || x~~~~> ::
171 * | || | vel. :: # lateral staggered faces over which fluxes are calculated
172 * scvf | || | Self ::
173 * ---------#######::::::::: x dof position
174 * scvf
175 * -- elements
176 *
177 * O position at which gradient is evaluated
178 * \endverbatim
179 */
180 template<class Problem, class FaceVariables>
181
2/2
✓ Branch 0 taken 346772 times.
✓ Branch 1 taken 284551500 times.
284898272 static Scalar velocityGradJI(const Problem& problem,
182 const Element& element,
183 const FVElementGeometry& fvGeometry,
184 const SubControlVolumeFace& scvf,
185 const FaceVariables& faceVars,
186 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
187 const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
188 const std::size_t localSubFaceIdx)
189 {
190 284898272 const auto eIdx = scvf.insideScvIdx();
191
2/2
✓ Branch 0 taken 346772 times.
✓ Branch 1 taken 284551500 times.
284898272 const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
192
193 // Assume a zero velocity gradient for pressure boundary conditions.
194
4/6
✓ Branch 0 taken 346772 times.
✓ Branch 1 taken 284551500 times.
✓ Branch 2 taken 346772 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 346772 times.
✗ Branch 5 not taken.
284898272 if (currentScvfBoundaryTypes && currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx))
195 return 0.0;
196
197 // For the velocityGrad_ji gradient, get the velocities perpendicular to the velocity at the current scvf.
198 // The inner one is located at staggered face within the own element,
199 // the outer one at the respective staggered face of the element on the other side of the
200 // current scvf.
201 284898272 const Scalar innerLateralVelocity = faceVars.velocityLateralInside(localSubFaceIdx);
202 854694816 const Scalar outerLateralVelocity = [&]()
203 {
204
2/2
✓ Branch 0 taken 284551500 times.
✓ Branch 1 taken 346772 times.
284898272 if (!scvf.boundary())
205 284551500 return faceVars.velocityLateralOutside(localSubFaceIdx);
206
2/4
✓ Branch 0 taken 346772 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 346772 times.
346772 else if (currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralScvf.directionIndex())))
207 {
208 // Sample the value of the Dirichlet BC at the center of the lateral face intersecting with the boundary.
209 const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
210 const auto lateralBoundaryFace = makeStaggeredBoundaryFace(scvf, lateralBoundaryFacePos);
211 return problem.dirichlet(element, lateralBoundaryFace)[Indices::velocity(lateralScvf.directionIndex())];
212 }
213
1/2
✓ Branch 0 taken 346772 times.
✗ Branch 1 not taken.
346772 else if (currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())))
214 {
215 346772 return beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars,
216 346772 currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
217 }
218 else
219 DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary types at " << lateralScvf.center());
220
2/2
✓ Branch 1 taken 142275652 times.
✓ Branch 2 taken 142622620 times.
284898272 }();
221
222 // Calculate the velocity gradient in positive coordinate direction.
223 569796544 const Scalar lateralDeltaV = scvf.normalInPosCoordDir()
224
2/2
✓ Branch 0 taken 142275652 times.
✓ Branch 1 taken 142622620 times.
284898272 ? (outerLateralVelocity - innerLateralVelocity)
225 : (innerLateralVelocity - outerLateralVelocity);
226
227 284898272 return lateralDeltaV / scvf.pairData(localSubFaceIdx).lateralDistance;
228 }
229
230 /*!
231 * \brief Returns the Beavers-Jospeh slip velocity for a scvf which lies on the boundary itself.
232 *
233 * \verbatim
234 * in.norm. B-J slip
235 * vel. vel.
236 * ^ ^
237 * | |
238 * scvf ---------######|* * boundary
239 * | || |* curr.
240 * | || |* scvf || and # staggered half-control-volume (own element)
241 * | || x~~~~>
242 * | || |* vel. # lateral staggered faces
243 * scvf | || |* Self
244 * ---------#######* x dof position
245 * scvf
246 * -- element
247 * \endverbatim
248 *
249 */
250 template<class Problem, class FaceVariables>
251 496380 static Scalar beaversJosephVelocityAtCurrentScvf(const Problem& problem,
252 const Element& element,
253 const FVElementGeometry& fvGeometry,
254 const SubControlVolumeFace& scvf,
255 const FaceVariables& faceVars,
256 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
257 const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
258 const std::size_t localSubFaceIdx)
259 {
260 496380 const auto eIdx = scvf.insideScvIdx();
261 496380 const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
262 496380 const Scalar innerLateralVelocity = faceVars.velocityLateralInside(localSubFaceIdx);
263
264 1489140 const auto tangentialVelocityGradient = [&]()
265 {
266 // If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for
267 // the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face
268 // (towards the current scvf).
269
3/4
✓ Branch 0 taken 17 times.
✓ Branch 1 taken 496363 times.
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
496397 static const bool unsymmetrizedGradientForBJ = getParamFromGroup<bool>(problem.paramGroup(),
270
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
34 "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
271
272
2/2
✓ Branch 0 taken 429420 times.
✓ Branch 1 taken 66960 times.
496380 if (unsymmetrizedGradientForBJ)
273 return 0.0;
274
275
2/2
✓ Branch 0 taken 22752 times.
✓ Branch 1 taken 406668 times.
429420 if (lateralScvf.boundary())
276 {
277
4/6
✓ Branch 0 taken 13501 times.
✓ Branch 1 taken 9251 times.
✓ Branch 2 taken 13501 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 13501 times.
✗ Branch 5 not taken.
509881 if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
278
1/2
✓ Branch 0 taken 13501 times.
✗ Branch 1 not taken.
13501 lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())))
279 return 0.0;
280 }
281
282 420169 return velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
283 496380 }();
284
285 496380 return problem.beaversJosephVelocity(element,
286 496380 fvGeometry.scv(scvf.insideScvIdx()),
287 lateralScvf,
288 scvf, /*on boundary*/
289 innerLateralVelocity,
290 496380 tangentialVelocityGradient);
291 }
292
293 /*!
294 * \brief Returns the Beavers-Jospeh slip velocity for a lateral scvf which lies on the boundary.
295 *
296 * \verbatim
297 * B-J slip * boundary
298 * ************** vel. *****
299 * scvf ---------##### ~~~~> :::: || and # staggered half-control-volume (own element)
300 * | || | curr. ::
301 * | || | scvf :: :: staggered half-control-volume (neighbor element)
302 * | || x~~~~> ::
303 * | || | vel. :: # lateral staggered faces
304 * scvf | || | Self ::
305 * ---------#######::::::::: x dof position
306 * scvf
307 * -- elements
308 * \endverbatim
309 */
310 template<class Problem, class FaceVariables>
311 367330 static Scalar beaversJosephVelocityAtLateralScvf(const Problem& problem,
312 const Element& element,
313 const FVElementGeometry& fvGeometry,
314 const SubControlVolumeFace& scvf,
315 const FaceVariables& faceVars,
316 const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
317 const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
318 const std::size_t localSubFaceIdx)
319 {
320 367330 const auto eIdx = scvf.insideScvIdx();
321 367330 const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
322 367330 const Scalar innerParallelVelocity = faceVars.velocitySelf();
323
324 1101990 const auto tangentialVelocityGradient = [&]()
325 {
326 // If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for
327 // the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face
328 // (towards the current scvf).
329
3/4
✓ Branch 0 taken 17 times.
✓ Branch 1 taken 367313 times.
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
367347 static const bool unsymmetrizedGradientForBJ = getParamFromGroup<bool>(problem.paramGroup(),
330
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
34 "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
331
332
2/2
✓ Branch 0 taken 311930 times.
✓ Branch 1 taken 55400 times.
367330 if (unsymmetrizedGradientForBJ)
333 return 0.0;
334
335
2/2
✓ Branch 0 taken 7038 times.
✓ Branch 1 taken 304892 times.
311930 if (scvf.boundary())
336 {
337
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 7038 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
367330 if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
338 currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())))
339 return 0.0;
340 }
341
342 304892 return velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
343 367330 }();
344
345 367330 return problem.beaversJosephVelocity(element,
346 367330 fvGeometry.scv(scvf.insideScvIdx()),
347 scvf,
348 lateralScvf, /*on boundary*/
349 innerParallelVelocity,
350 367330 tangentialVelocityGradient);
351 }
352
353 private:
354
355 /*!
356 * \brief Get the location of the lateral staggered face's center.
357 * Only needed for boundary conditions if the current scvf or the lateral one is on a boundary.
358 *
359 * \verbatim
360 * --------#######o || frontal face of staggered half-control-volume
361 * | || | current scvf # lateral staggered face of interest (may lie on a boundary)
362 * | || | x dof position
363 * | || x~~~~> vel.Self -- element boundaries, current scvf may lie on a boundary
364 * | || | o position at which the boundary conditions will be evaluated
365 * | || | (lateralStaggeredFaceCenter)
366 * ----------------
367 * \endverbatim
368 */
369 8663412 static const GlobalPosition& lateralStaggeredFaceCenter_(const SubControlVolumeFace& scvf, const int localSubFaceIdx)
370 {
371 8663412 return scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter;
372 };
373 };
374
375 } // end namespace Dumux
376
377 #endif
378