velocitygradients.hh 19.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
 * \ingroup NavierStokesModel
 * \copydoc Dumux::StaggeredVelocityGradients
 */
#ifndef DUMUX_NAVIERSTOKES_STAGGERED_VELOCITYGRADIENTS_HH
#define DUMUX_NAVIERSTOKES_STAGGERED_VELOCITYGRADIENTS_HH

#include <dune/common/std/optional.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/common/parameters.hh>

namespace Dumux {

/*!
 * \ingroup NavierStokesModel
 * \brief Helper class for calculating the velocity gradients for the Navier-Stokes model using the staggered grid discretization.
 */
37
template<class Scalar, class GridGeometry, class BoundaryTypes, class Indices>
38
39
class StaggeredVelocityGradients
{
40
41
    using FVElementGeometry = typename GridGeometry::LocalView;
    using GridView = typename GridGeometry::GridView;
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
    using Element = typename GridView::template Codim<0>::Entity;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;

public:

    /*!
     * \brief Returns the in-axis velocity gradient.
     *
     * \verbatim
     *              ---------=======                 == and # staggered half-control-volume
     *              |       #      | current scvf
     *              |       #      |                 # staggered face over which fluxes are calculated
     *   vel.Opp <~~|       O~~>   x~~~~> vel.Self
     *              |       #      |                 x dof position
     *        scvf  |       #      |
     *              --------========                 -- element
     *
     *                                               O position at which gradient is evaluated
     * \endverbatim
     */
    template<class FaceVariables>
    static Scalar velocityGradII(const SubControlVolumeFace& scvf,
                                 const FaceVariables& faceVars)
    {
        // The velocities of the dof at interest and the one of the opposite scvf.
        const Scalar velocitySelf = faceVars.velocitySelf();
        const Scalar velocityOpposite = faceVars.velocityOpposite();

        return ((velocityOpposite - velocitySelf) / scvf.selfToOppositeDistance()) * scvf.directionSign();
    }

    /*!
     * \brief Returns the velocity gradient perpendicular to the orientation of our current scvf.
     *
     * \verbatim
     *              ----------------
     *              |              |vel.
     *              |              |Parallel
     *              |              |~~~~>       ------->
     *              |              |             ------>     * gradient
     *              |              |              ----->
     *       scvf   ---------######O:::::::::      ---->     || and # staggered half-control-volume (own element)
     *              |      ||      | curr. ::       --->
     *              |      ||      | scvf  ::        -->     :: staggered half-control-volume (neighbor element)
     *              |      ||      x~~~~>  ::         ->
     *              |      ||      | vel.  ::                # lateral staggered faces over which fluxes are calculated
     *        scvf  |      ||      | Self  ::
     *              ---------#######:::::::::                x dof position
     *                 scvf
     *                                                       -- elements
     *
     *                                                       O position at which gradient is evaluated
     * \endverbatim
     */
    template<class Problem, class FaceVariables>
    static Scalar velocityGradIJ(const Problem& problem,
                                 const Element& element,
                                 const FVElementGeometry& fvGeometry,
                                 const SubControlVolumeFace& scvf,
                                 const FaceVariables& faceVars,
                                 const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
                                 const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
                                 const std::size_t localSubFaceIdx)
    {
        const auto eIdx = scvf.insideScvIdx();
        const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);

        // For the velocityGrad_ij derivative, get the velocities at the current (own) scvf
        // and at the parallel one at the neighboring scvf.
        const Scalar innerParallelVelocity = faceVars.velocitySelf();

        const auto outerParallelVelocity = [&]()
        {
            if (!lateralScvf.boundary())
                return faceVars.velocityParallel(localSubFaceIdx, 0);
            else if (lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex())))
            {
                // Sample the value of the Dirichlet BC at the center of the staggered lateral face.
                const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
                return problem.dirichlet(element, lateralScvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(scvf.directionIndex())];
            }
            else if (lateralFaceBoundaryTypes->isBJS(Indices::velocity(scvf.directionIndex())))
            {
126
127
                return beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf,  faceVars,
                                                          currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
            }
            else
                DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary type at " << lateralScvf.center());
        }();

        // The velocity gradient already accounts for the orientation
        // of the staggered face's outer normal vector. This also correctly accounts for the reduced
        // distance used in the gradient if the lateral scvf lies on a boundary.
        return (outerParallelVelocity - innerParallelVelocity)
               / scvf.parallelDofsDistance(localSubFaceIdx, 0) * lateralScvf.directionSign();
    }

    /*!
     * \brief Returns the velocity gradient in line with our current scvf.
     *
     * \verbatim
     *                      ^       gradient
     *                      |  ^
     *                      |  |  ^
     *                      |  |  |  ^
     *                      |  |  |  |  ^
     *                      |  |  |  |  |  ^
     *                      |  |  |  |  |  |
     *
     *              ----------------
     *              |              |
     *              |    in.norm.  |
     *              |       vel.   |
     *              |       ^      |        ^ out.norm.vel.
     *              |       |      |        |
     *       scvf   ---------######O:::::::::       || and # staggered half-control-volume (own element)
     *              |      ||      | curr. ::
     *              |      ||      | scvf  ::       :: staggered half-control-volume (neighbor element)
     *              |      ||      x~~~~>  ::
     *              |      ||      | vel.  ::       # lateral staggered faces over which fluxes are calculated
     *        scvf  |      ||      | Self  ::
     *              ---------#######:::::::::       x dof position
     *                 scvf
     *                                              -- elements
     *
     *                                              O position at which gradient is evaluated
     * \endverbatim
     */
    template<class Problem, class FaceVariables>
    static Scalar velocityGradJI(const Problem& problem,
                                 const Element& element,
                                 const FVElementGeometry& fvGeometry,
                                 const SubControlVolumeFace& scvf,
                                 const FaceVariables& faceVars,
                                 const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
                                 const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
                                 const std::size_t localSubFaceIdx)
    {
        const auto eIdx = scvf.insideScvIdx();
        const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);

        // Assume a zero velocity gradient for pressure boundary conditions.
        if (currentScvfBoundaryTypes && currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx))
            return 0.0;

        // For the velocityGrad_ji gradient, get the velocities perpendicular to the velocity at the current scvf.
        // The inner one is located at staggered face within the own element,
        // the outer one at the respective staggered face of the element on the other side of the
        // current scvf.
        const Scalar innerLateralVelocity = faceVars.velocityLateralInside(localSubFaceIdx);
        const Scalar outerLateralVelocity = [&]()
        {
            if (!scvf.boundary())
                return faceVars.velocityLateralOutside(localSubFaceIdx);
            else if (currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralScvf.directionIndex())))
            {
                // Sample the value of the Dirichlet BC at the center of the lateral face intersecting with the boundary.
                const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
                return problem.dirichlet(element, scvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(lateralScvf.directionIndex())];
            }
            else if (currentScvfBoundaryTypes->isBJS(Indices::velocity(lateralScvf.directionIndex())))
            {
205
206
                return beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf,  faceVars,
                                                          currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
207
208
209
210
211
212
213
214
215
216
217
218
219
            }
            else
                DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary types at " << lateralScvf.center());
        }();

        // Calculate the velocity gradient in positive coordinate direction.
        const Scalar lateralDeltaV = scvf.normalInPosCoordDir()
                                    ? (outerLateralVelocity - innerLateralVelocity)
                                    : (innerLateralVelocity - outerLateralVelocity);

        return lateralDeltaV / scvf.pairData(localSubFaceIdx).lateralDistance;
    }

220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
    /*!
     * \brief Returns the Beavers-Jospeh slip velocity for a scvf which lies on the boundary itself.
     *
     * \verbatim
     *                  in.norm.  B-J slip
     *                     vel.   vel.
     *                     ^       ^
     *                     |       |
     *       scvf   ---------######|*               * boundary
     *              |      ||      |* curr.
     *              |      ||      |* scvf          || and # staggered half-control-volume (own element)
     *              |      ||      x~~~~>
     *              |      ||      |* vel.          # lateral staggered faces
     *        scvf  |      ||      |* Self
     *              ---------#######*                x dof position
     *                 scvf
     *                                              -- element
     * \endverbatim
     *
     */
    template<class Problem, class FaceVariables>
    static Scalar beaversJosephVelocityAtCurrentScvf(const Problem& problem,
                                                     const Element& element,
                                                     const FVElementGeometry& fvGeometry,
                                                     const SubControlVolumeFace& scvf,
                                                     const FaceVariables& faceVars,
                                                     const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
                                                     const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
                                                     const std::size_t localSubFaceIdx)
    {
        const auto eIdx = scvf.insideScvIdx();
        const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
        const Scalar innerLateralVelocity = faceVars.velocityLateralInside(localSubFaceIdx);

        const auto tangentialVelocityGradient = [&]()
        {
            // If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for
            // the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face
            // (towards the current scvf).
            static const bool unsymmetrizedGradientForBJ = getParamFromGroup<bool>(problem.paramGroup(),
                                                           "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);

            if (unsymmetrizedGradientForBJ)
                return 0.0;

            if (lateralScvf.boundary())
            {
                if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
                    lateralFaceBoundaryTypes->isBJS(Indices::velocity(scvf.directionIndex())))
                    return 0.0;
            }

            return velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
        }();

275
276
277
278
279
        return problem.beaversJosephVelocity(element,
                                             fvGeometry.scv(scvf.insideScvIdx()),
                                             lateralScvf,
                                             scvf, /*on boundary*/
                                             innerLateralVelocity,
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
                                             tangentialVelocityGradient);
    }

    /*!
     * \brief Returns the Beavers-Jospeh slip velocity for a lateral scvf which lies on the boundary.
     *
     * \verbatim
     *                             B-J slip                  * boundary
     *              ************** vel. *****
     *       scvf   ---------##### ~~~~> ::::                || and # staggered half-control-volume (own element)
     *              |      ||      | curr. ::
     *              |      ||      | scvf  ::                :: staggered half-control-volume (neighbor element)
     *              |      ||      x~~~~>  ::
     *              |      ||      | vel.  ::                # lateral staggered faces
     *        scvf  |      ||      | Self  ::
     *              ---------#######:::::::::                x dof position
     *                 scvf
     *                                                       -- elements
     * \endverbatim
     */
    template<class Problem, class FaceVariables>
    static Scalar beaversJosephVelocityAtLateralScvf(const Problem& problem,
                                                     const Element& element,
                                                     const FVElementGeometry& fvGeometry,
                                                     const SubControlVolumeFace& scvf,
                                                     const FaceVariables& faceVars,
                                                     const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
                                                     const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
                                                     const std::size_t localSubFaceIdx)
    {
        const auto eIdx = scvf.insideScvIdx();
        const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
        const Scalar innerParallelVelocity = faceVars.velocitySelf();

        const auto tangentialVelocityGradient = [&]()
        {
            // If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for
            // the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face
            // (towards the current scvf).
            static const bool unsymmetrizedGradientForBJ = getParamFromGroup<bool>(problem.paramGroup(),
                                                           "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);

            if (unsymmetrizedGradientForBJ)
                return 0.0;

            if (scvf.boundary())
            {
                if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
                    currentScvfBoundaryTypes->isBJS(Indices::velocity(lateralScvf.directionIndex())))
                    return 0.0;
            }

            return velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
        }();

335
336
337
338
339
        return problem.beaversJosephVelocity(element,
                                             fvGeometry.scv(scvf.insideScvIdx()),
                                             scvf,
                                             lateralScvf, /*on boundary*/
                                             innerParallelVelocity,
340
341
342
                                             tangentialVelocityGradient);
    }

343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
private:

    /*!
     * \brief Get the location of the lateral staggered face's center.
     *        Only needed for boundary conditions if the current scvf or the lateral one is on a bounary.
     *
     * \verbatim
     *      --------#######o                 || frontal face of staggered half-control-volume
     *      |      ||      | current scvf    #  lateral staggered face of interest (may lie on a boundary)
     *      |      ||      |                 x  dof position
     *      |      ||      x~~~~> vel.Self   -- element boundaries, current scvf may lie on a boundary
     *      |      ||      |                 o  position at which the boundary conditions will be evaluated
     *      |      ||      |                    (lateralStaggeredFaceCenter)
     *      ----------------
     * \endverbatim
     */
    static const GlobalPosition& lateralStaggeredFaceCenter_(const SubControlVolumeFace& scvf, const int localSubFaceIdx)
    {
        return scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter;
    };
};

} // end namespace Dumux

#endif // DUMUX_NAVIERSTOKES_STAGGERED_VELOCITYGRADIENTS_HH