GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/freeflow/navierstokes/staggered/velocitygradients.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 64 70 91.4%
Functions: 156 264 59.1%
Branches: 69 150 46.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::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 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 144714524 const Scalar velocitySelf = faceVars.velocitySelf();
57 144714524 const Scalar velocityOpposite = faceVars.velocityOpposite();
58
59
4/4
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 144714474 times.
✓ Branch 2 taken 50 times.
✓ Branch 3 taken 144714474 times.
289429048 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 284152147 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 284152147 const auto eIdx = scvf.insideScvIdx();
116 568304294 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 284152147 const Scalar innerParallelVelocity = faceVars.velocitySelf();
121
122 284152147 const auto outerParallelVelocity = [&]()
123 {
124
2/2
✓ Branch 0 taken 275283830 times.
✓ Branch 1 taken 8868317 times.
284152147 if (!lateralScvf.boundary())
125 275283830 return faceVars.velocityParallel(localSubFaceIdx, 0);
126
4/4
✓ Branch 0 taken 255244 times.
✓ Branch 1 taken 8613073 times.
✓ Branch 2 taken 255244 times.
✓ Branch 3 taken 8613073 times.
17736634 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 25839219 const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
130
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
8613073 const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralScvf, lateralBoundaryFacePos);
131
6/7
✓ Branch 0 taken 198168 times.
✓ Branch 1 taken 8414905 times.
✓ Branch 2 taken 201108 times.
✓ Branch 3 taken 8411965 times.
✓ Branch 4 taken 14868 times.
✓ Branch 5 taken 289112 times.
✗ Branch 6 not taken.
9009409 return problem.dirichlet(element, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())];
132 }
133
3/6
✓ Branch 0 taken 255244 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 255244 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 255244 times.
✗ Branch 5 not taken.
765732 else if (lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())))
134 {
135 765732 return beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars,
136 255244 currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
137 }
138 else
139 DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary type at " << lateralScvf.center());
140 559691221 }();
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 284152147 return (outerParallelVelocity - innerParallelVelocity)
146 568304294 / 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 283228418 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
2/2
✓ Branch 0 taken 346528 times.
✓ Branch 1 taken 282881890 times.
283228418 const auto eIdx = scvf.insideScvIdx();
191
4/4
✓ Branch 0 taken 346528 times.
✓ Branch 1 taken 282881890 times.
✓ Branch 2 taken 346528 times.
✓ Branch 3 taken 282881890 times.
566456836 const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
192
193 // Assume a zero velocity gradient for pressure boundary conditions.
194
6/8
✓ Branch 0 taken 346528 times.
✓ Branch 1 taken 282881890 times.
✓ Branch 2 taken 346528 times.
✓ Branch 3 taken 282881890 times.
✓ Branch 4 taken 346528 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 346528 times.
✗ Branch 7 not taken.
566456836 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 283228418 const Scalar innerLateralVelocity = faceVars.velocityLateralInside(localSubFaceIdx);
202 283228418 const Scalar outerLateralVelocity = [&]()
203 {
204
2/2
✓ Branch 0 taken 282881890 times.
✓ Branch 1 taken 346528 times.
283228418 if (!scvf.boundary())
205 282881890 return faceVars.velocityLateralOutside(localSubFaceIdx);
206
2/4
✓ Branch 0 taken 346528 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 346528 times.
✗ Branch 3 not taken.
693056 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
3/6
✓ Branch 0 taken 346528 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 346528 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 346528 times.
✗ Branch 5 not taken.
1039584 else if (currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())))
214 {
215 1039584 return beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars,
216 346528 currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
217 }
218 else
219 DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary types at " << lateralScvf.center());
220 566456836 }();
221
222 // Calculate the velocity gradient in positive coordinate direction.
223 283228418 const Scalar lateralDeltaV = scvf.normalInPosCoordDir()
224
4/4
✓ Branch 0 taken 141440847 times.
✓ Branch 1 taken 141787571 times.
✓ Branch 2 taken 141440847 times.
✓ Branch 3 taken 141787571 times.
566456836 ? (outerLateralVelocity - innerLateralVelocity)
225 : (innerLateralVelocity - outerLateralVelocity);
226
227 566456836 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 496136 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 496136 const auto eIdx = scvf.insideScvIdx();
261 992272 const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
262 496136 const Scalar innerLateralVelocity = faceVars.velocityLateralInside(localSubFaceIdx);
263
264 496136 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
4/6
✓ Branch 0 taken 17 times.
✓ Branch 1 taken 496119 times.
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 17 times.
✗ Branch 7 not taken.
496136 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 429176 times.
✓ Branch 1 taken 66960 times.
496136 if (unsymmetrizedGradientForBJ)
273 return 0.0;
274
275
2/2
✓ Branch 0 taken 22726 times.
✓ Branch 1 taken 406450 times.
429176 if (lateralScvf.boundary())
276 {
277
6/8
✓ Branch 0 taken 13488 times.
✓ Branch 1 taken 9238 times.
✓ Branch 2 taken 13488 times.
✓ Branch 3 taken 9238 times.
✓ Branch 4 taken 13488 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 13488 times.
✗ Branch 7 not taken.
72428 if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
278
1/4
✓ Branch 0 taken 13488 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
13488 lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())))
279 return 0.0;
280 }
281
282 1679752 return velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
283 916074 }();
284
285 return problem.beaversJosephVelocity(element,
286 fvGeometry.scv(scvf.insideScvIdx()),
287 lateralScvf,
288 scvf, /*on boundary*/
289 innerLateralVelocity,
290 1488408 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 367152 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 367152 const auto eIdx = scvf.insideScvIdx();
321 734304 const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
322 367152 const Scalar innerParallelVelocity = faceVars.velocitySelf();
323
324 367152 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
4/6
✓ Branch 0 taken 17 times.
✓ Branch 1 taken 367135 times.
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 17 times.
✗ Branch 7 not taken.
367152 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 311752 times.
✓ Branch 1 taken 55400 times.
367152 if (unsymmetrizedGradientForBJ)
333 return 0.0;
334
335
2/2
✓ Branch 0 taken 7028 times.
✓ Branch 1 taken 304724 times.
311752 if (scvf.boundary())
336 {
337
2/8
✗ Branch 0 not taken.
✓ Branch 1 taken 7028 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7028 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
14056 if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
338 currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())))
339 return 0.0;
340 }
341
342 914172 return velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
343 671876 }();
344
345 return problem.beaversJosephVelocity(element,
346 fvGeometry.scv(scvf.insideScvIdx()),
347 scvf,
348 lateralScvf, /*on boundary*/
349 innerParallelVelocity,
350 1101456 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 static const GlobalPosition& lateralStaggeredFaceCenter_(const SubControlVolumeFace& scvf, const int localSubFaceIdx)
370 {
371 17226146 return scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter;
372 };
373 };
374
375 } // end namespace Dumux
376
377 #endif
378