problem.hh 26.7 KB
Newer Older
Thomas Fetzer's avatar
Thomas Fetzer committed
1
2
3
4
5
6
7
// -*- 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    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
Thomas Fetzer's avatar
Thomas Fetzer committed
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 *   (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 RANSModel
 * \copydoc Dumux::RANSProblem
 */
#ifndef DUMUX_RANS_PROBLEM_HH
#define DUMUX_RANS_PROBLEM_HH

Kilian Weishaupt's avatar
Kilian Weishaupt committed
27
#include <dune/common/fmatrix.hh>
Thomas Fetzer's avatar
Thomas Fetzer committed
28
29
#include <dumux/common/properties.hh>
#include <dumux/common/staggeredfvproblem.hh>
30
#include <dumux/discretization/localview.hh>
31
#include <dumux/discretization/staggered/elementsolution.hh>
32
#include <dumux/discretization/method.hh>
Thomas Fetzer's avatar
Thomas Fetzer committed
33
34
35
36
#include <dumux/freeflow/navierstokes/problem.hh>

#include "model.hh"

37
namespace Dumux {
Thomas Fetzer's avatar
Thomas Fetzer committed
38

39
40
41
42
43
44
//! forward declare
template<class TypeTag, TurbulenceModel turbulenceModel>
class RANSProblemImpl;

//! the turbulence-model-specfic RANS problem
template<class TypeTag>
45
using RANSProblem = RANSProblemImpl<TypeTag, GetPropType<TypeTag, Properties::ModelTraits>::turbulenceModel()>;
46

Thomas Fetzer's avatar
Thomas Fetzer committed
47
48
49
50
/*!
 * \ingroup RANSModel
 * \brief Reynolds-Averaged Navier-Stokes problem base class.
 *
51
52
53
 * This implements some base functionality for RANS models.
 * Especially vectors containing all wall-relevant properties, which are accessed
 * by the volumevariables.
Thomas Fetzer's avatar
Thomas Fetzer committed
54
55
 */
template<class TypeTag>
56
class RANSProblemBase : public NavierStokesProblem<TypeTag>
Thomas Fetzer's avatar
Thomas Fetzer committed
57
{
58
    using ParentType = NavierStokesProblem<TypeTag>;
59
    using Implementation = GetPropType<TypeTag, Properties::Problem>;
Thomas Fetzer's avatar
Thomas Fetzer committed
60

61
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
Thomas Fetzer's avatar
Thomas Fetzer committed
62

63
64
65
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GridGeometry::LocalView;
    using GridView = typename GridGeometry::GridView;
66
    using Element = typename GridView::template Codim<0>::Entity;
67
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
68
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
69
70
    using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
Kilian Weishaupt's avatar
Kilian Weishaupt committed
71
    using PrimaryVariables = typename VolumeVariables::PrimaryVariables;
72
73
74
    using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
    using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
Thomas Fetzer's avatar
Thomas Fetzer committed
75

76
    using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
Kilian Weishaupt's avatar
Kilian Weishaupt committed
77
78

    static constexpr auto dim = GridView::dimension;
79
    static constexpr int numCorners = SubControlVolumeFace::numCornersPerFace;
Kilian Weishaupt's avatar
Kilian Weishaupt committed
80
    using DimVector = GlobalPosition;
81
    using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
Thomas Fetzer's avatar
Thomas Fetzer committed
82
83

public:
84
85
    /*!
     * \brief The constructor
86
     * \param gridGeometry The finite volume grid geometry
87
88
     * \param paramGroup The parameter group in which to look for runtime parameters first (default is "")
     */
89
90
    RANSProblemBase(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "")
    : ParentType(gridGeometry, paramGroup)
91
    { }
92
93

    /*!
94
95
96
97
     * \brief Update the static (solution independent) relations to the walls
     *
     * This function determines all element with a wall intersection,
     * the wall distances and the relation to the neighboring elements.
98
     */
99
    void updateStaticWallProperties()
100
    {
101
        using std::abs;
102
        std::cout << "Update static wall properties. ";
103
        calledUpdateStaticWallProperties = true;
104

105
        // update size and initial values of the global vectors
106
107
108
109
110
111
112
113
114
        wallElementIdx_.resize(this->gridGeometry().elementMapper().size());
        wallDistance_.resize(this->gridGeometry().elementMapper().size(), std::numeric_limits<Scalar>::max());
        neighborIdx_.resize(this->gridGeometry().elementMapper().size());
        cellCenter_.resize(this->gridGeometry().elementMapper().size(), GlobalPosition(0.0));
        velocity_.resize(this->gridGeometry().elementMapper().size(), DimVector(0.0));
        velocityMaximum_.resize(this->gridGeometry().elementMapper().size(), DimVector(0.0));
        velocityGradients_.resize(this->gridGeometry().elementMapper().size(), DimMatrix(0.0));
        stressTensorScalarProduct_.resize(this->gridGeometry().elementMapper().size(), 0.0);
        vorticityTensorScalarProduct_.resize(this->gridGeometry().elementMapper().size(), 0.0);
115
116
        flowNormalAxis_.resize(this->gridGeometry().elementMapper().size(), fixedFlowNormalAxis_);
        wallNormalAxis_.resize(this->gridGeometry().elementMapper().size(), fixedWallNormalAxis_);
117
118
        kinematicViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
        sandGrainRoughness_.resize(this->gridGeometry().elementMapper().size(), 0.0);
119

120
121
122
123
124
125
126
127
        // store the element indicies for all elements with an intersection on the wall
        std::vector<unsigned int> wallElementIndicies;

        // for each wall element, store the location of the face center and each corner.
        std::vector<std::array<GlobalPosition, numCorners+1>> wallPositions;

        // for each wall element, store the faces normal axis
        std::vector<unsigned int> wallFaceNormalAxis;
128

129
130
        const auto gridView = this->gridGeometry().gridView();
        auto fvGeometry = localView(this->gridGeometry());
131

132
133
        for (const auto& element : elements(gridView))
        {
134
135
            fvGeometry.bindElement(element);
            for (const auto& scvf : scvfs(fvGeometry))
136
            {
137
                // only search for walls at a global boundary
138
                if (!scvf.boundary())
139
140
                    continue;

141
                if (asImp_().isOnWall(scvf))
142
                {
143
144
145
146
147
148
149
150
151
152
153
154
                    // element has an scvf on the wall, store element index
                    wallElementIndicies.push_back(this->gridGeometry().elementMapper().index(element));

                    // store the location of the wall adjacent face's center and all corners
                    std::array<GlobalPosition, numCorners+1> wallElementPosition;
                    wallElementPosition[0] = scvf.center();
                    for (int i = 1; i <= numCorners; i++)
                        wallElementPosition[i] = scvf.corner(i-1);
                    wallPositions.push_back(wallElementPosition);

                    // Store the wall adjacent face's normal direction
                    wallFaceNormalAxis.push_back(scvf.directionIndex());
155
156
157
                }
            }
        }
158
        // output the number of wall adjacent faces. Check that this is non-zero.
159
        std::cout << "NumWallIntersections=" << wallPositions.size() << std::endl;
160
161
162
163
        if (wallPositions.size() == 0)
            DUNE_THROW(Dune::InvalidStateException,
                       "No wall intersections have been found. Make sure that the isOnWall(globalPos) is working properly.");

164

165
        // search for shortest distance to the wall for each element
166
167
        for (const auto& element : elements(gridView))
        {
168
            unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
169
170

            // Store the cell center position for each element
171
            cellCenter_[elementIdx] = element.geometry().center();
172

173
174
            for (unsigned int i = 0; i < wallPositions.size(); ++i)
            {
175
176
177
178
179
                // Find the minimum distance from the cell center to the wall face (center and corners)
                std::array<Scalar,numCorners+1> cellToWallDistances;
                for (unsigned int j = 0; j < wallPositions[i].size(); j++)
                    cellToWallDistances[j] = (cellCenter(elementIdx) - wallPositions[i][j]).two_norm();
                Scalar distanceToWall = *std::min_element(cellToWallDistances.begin(), cellToWallDistances.end());
180

181
                if (distanceToWall < wallDistance_[elementIdx])
182
                {
183
184
185
186
187
                    wallDistance_[elementIdx] = distanceToWall;
                    wallElementIdx_[elementIdx] = wallElementIndicies[i];
                    if ( !(hasParam("RANS.WallNormalAxis")) )
                        wallNormalAxis_[elementIdx] = wallFaceNormalAxis[i];
                    sandGrainRoughness_[elementIdx] = asImp_().sandGrainRoughnessAtPos(wallPositions[i][0]);
188
189
190
                }
            }
        }
191

192
        // search for neighbor Idxs
193
194
        for (const auto& element : elements(gridView))
        {
195
            unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
196
197
            for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
            {
198
199
                neighborIdx_[elementIdx][dimIdx][0] = elementIdx;
                neighborIdx_[elementIdx][dimIdx][1] = elementIdx;
200
            }
201
202

            for (const auto& intersection : intersections(gridView, element))
203
            {
204
                if (intersection.boundary())
205
206
                    continue;

207
                unsigned int neighborIdx = this->gridGeometry().elementMapper().index(intersection.outside());
208
209
                for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
                {
210
                    if (abs(cellCenter_[elementIdx][dimIdx] - cellCenter_[neighborIdx][dimIdx]) > 1e-8)
211
                    {
212
                        if (cellCenter_[elementIdx][dimIdx] > cellCenter_[neighborIdx][dimIdx])
213
                        {
214
                            neighborIdx_[elementIdx][dimIdx][0] = neighborIdx;
215
                        }
216
                        if (cellCenter_[elementIdx][dimIdx] < cellCenter_[neighborIdx][dimIdx])
217
                        {
218
                            neighborIdx_[elementIdx][dimIdx][1] = neighborIdx;
219
220
221
222
223
                        }
                    }
                }
            }
        }
224
225
226
    }

    /*!
227
228
229
230
231
232
233
     * \brief Update the dynamic (solution dependent) relations to the walls
     *
     * The basic function calcuates the cell-centered velocities and
     * the respective gradients.
     * Further, the kinematic viscosity at the wall is stored.
     *
     * \param curSol The solution vector.
234
     */
235
    void updateDynamicWallProperties(const SolutionVector& curSol)
236
    {
237
        using std::abs;
238
239
        using std::max;
        using std::min;
240
        std::cout << "Update dynamic wall properties." << std::endl;
241
242
243
        if (!calledUpdateStaticWallProperties)
            DUNE_THROW(Dune::InvalidStateException,
                       "You have to call updateStaticWallProperties() once before you call updateDynamicWallProperties().");
244

245
        // re-initialize min and max values
246
247
        velocityMaximum_.assign(this->gridGeometry().elementMapper().size(), DimVector(1e-16));
        velocityMinimum_.assign(this->gridGeometry().elementMapper().size(), DimVector(std::numeric_limits<Scalar>::max()));
248

249
        // calculate cell-center-averaged velocities
250
        for (const auto& element : elements(this->gridGeometry().gridView()))
251
        {
252
            auto fvGeometry = localView(this->gridGeometry());
253
            fvGeometry.bindElement(element);
254
            unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
255
256
257
258
259
260

            // calculate velocities
            DimVector velocityTemp(0.0);
            for (auto&& scvf : scvfs(fvGeometry))
            {
                const int dofIdxFace = scvf.dofIndex();
261
                const auto numericalSolutionFace = curSol[GridGeometry::faceIdx()][dofIdxFace][Indices::velocity(scvf.directionIndex())];
262
263
264
                velocityTemp[scvf.directionIndex()] += numericalSolutionFace;
            }
            for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
265
                velocity_[elementIdx][dimIdx] = velocityTemp[dimIdx] * 0.5; // faces are equidistant to cell center
266
        }
267

268
        // calculate cell-center-averaged velocity gradients, maximum, and minimum values
269
        for (const auto& element : elements(this->gridGeometry().gridView()))
270
        {
271
            unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
272
            unsigned int wallElementIdx = wallElementIdx_[elementIdx];
273

274
            Scalar maxVelocity = 0.0;
275
276
277
278
            for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
            {
                for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
                {
279
280
281
282
283
                    velocityGradients_[elementIdx][velIdx][dimIdx]
                        = (velocity_[neighborIdx_[elementIdx][dimIdx][1]][velIdx]
                              - velocity_[neighborIdx_[elementIdx][dimIdx][0]][velIdx])
                          / (cellCenter_[neighborIdx_[elementIdx][dimIdx][1]][dimIdx]
                              - cellCenter_[neighborIdx_[elementIdx][dimIdx][0]][dimIdx]);
284
285
286
                    if (abs(cellCenter_[neighborIdx_[elementIdx][dimIdx][1]][dimIdx]
                            - cellCenter_[neighborIdx_[elementIdx][dimIdx][0]][dimIdx]) < 1e-8)
                        velocityGradients_[elementIdx][velIdx][dimIdx] = 0.0;
287
                }
288

289
                if (abs(velocity_[elementIdx][dimIdx]) > abs(velocityMaximum_[wallElementIdx][dimIdx]))
290
                {
291
                    velocityMaximum_[wallElementIdx][dimIdx] = velocity_[elementIdx][dimIdx];
292
                }
293
                if (abs(velocity_[elementIdx][dimIdx]) < abs(velocityMinimum_[wallElementIdx][dimIdx]))
294
                {
295
                    velocityMinimum_[wallElementIdx][dimIdx] = velocity_[elementIdx][dimIdx];
296
                }
297
298
299

                if (0 <= flowNormalAxis && flowNormalAxis < dim)
                {
300
                    flowNormalAxis_[elementIdx] = flowNormalAxis;
301
                }
302
                else if (abs(maxVelocity) < abs(velocity_[elementIdx][dimIdx]))
303
                {
304
305
                    maxVelocity = abs(velocity_[elementIdx][dimIdx]);
                    flowNormalAxis_[elementIdx] = dimIdx;
306
                }
307
            }
308

309
            auto fvGeometry = localView(this->gridGeometry());
310
311
312
            fvGeometry.bindElement(element);
            for (auto&& scvf : scvfs(fvGeometry))
            {
313
314
315
                // adapt calculations for Dirichlet condition
                unsigned int scvfNormDim = scvf.directionIndex();
                if (scvf.boundary())
316
317
318
                {
                    for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
                    {
319
320
321
                        if (!asImp_().boundaryTypes(element, scvf).isDirichlet(Indices::velocity(velIdx)))
                            continue;

322
323
                        Scalar dirichletVelocity = asImp_().dirichlet(element, scvf)[Indices::velocity(velIdx)];

324
325
326
                        unsigned int neighborIdx = neighborIdx_[elementIdx][scvfNormDim][0];
                        if (scvf.center()[scvfNormDim] < cellCenter_[elementIdx][scvfNormDim])
                            neighborIdx = neighborIdx_[elementIdx][scvfNormDim][1];
327

328
329
330
                        velocityGradients_[elementIdx][velIdx][scvfNormDim]
                            = (velocity_[neighborIdx][velIdx] - dirichletVelocity)
                              / (cellCenter_[neighborIdx][scvfNormDim] - scvf.center()[scvfNormDim]);
331
332
333
334
335
336
337
338
339
340
341
342
                    }
                }

                // Calculate the BJS-velocity by accounting for all sub faces.
                std::vector<int> bjsNumFaces(dim, 0);
                std::vector<unsigned int> bjsNeighbor(dim, 0);
                DimVector bjsVelocityAverage(0.0);
                DimVector normalNormCoordinate(0.0);
                unsigned int velIdx = Indices::velocity(scvfNormDim);
                const int numSubFaces = scvf.pairData().size();
                for(int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
                {
343
                    const auto& lateralFace = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData()[localSubFaceIdx].localLateralFaceIdx);
344
345

                    // adapt calculations for Beavers-Joseph-Saffman condition
346
                    unsigned int normalNormDim = lateralFace.directionIndex();
347
                    if (lateralFace.boundary() && (asImp_().boundaryTypes(element, lateralFace).isBeaversJoseph(Indices::velocity(velIdx))))
348
                    {
349
                        unsigned int neighborIdx = neighborIdx_[elementIdx][normalNormDim][0];
350
                        if (lateralFace.center()[normalNormDim] < cellCenter_[elementIdx][normalNormDim])
351
352
                            neighborIdx = neighborIdx_[elementIdx][normalNormDim][1];

353
                        const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx());
354
                        bjsVelocityAverage[normalNormDim] += ParentType::beaversJosephVelocity(element, scv, scvf, lateralFace, velocity_[elementIdx][velIdx], 0.0);
355
356
357
                        if (bjsNumFaces[normalNormDim] > 0 && neighborIdx != bjsNeighbor[normalNormDim])
                            DUNE_THROW(Dune::InvalidStateException, "Two different neighborIdx should not occur");
                        bjsNeighbor[normalNormDim] = neighborIdx;
358
                        normalNormCoordinate[normalNormDim] = lateralFace.center()[normalNormDim];
359
                        bjsNumFaces[normalNormDim]++;
360
361
                    }
                }
362
363
364
365
366
                for (unsigned dirIdx = 0; dirIdx < dim; ++dirIdx)
                {
                    if (bjsNumFaces[dirIdx] == 0)
                        continue;

367
                    unsigned int neighborIdx = bjsNeighbor[dirIdx];
368
369
                    bjsVelocityAverage[dirIdx] /= bjsNumFaces[dirIdx];

370
371
372
                    velocityGradients_[elementIdx][velIdx][dirIdx]
                        = (velocity_[neighborIdx][velIdx] - bjsVelocityAverage[dirIdx])
                          / (cellCenter_[neighborIdx][dirIdx] - normalNormCoordinate[dirIdx]);
373
374

                }
375
            }
376
        }
377

378
        // calculate or call all secondary variables
379
        for (const auto& element : elements(this->gridGeometry().gridView()))
380
        {
381
            unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
382

383
384
385
386
387
            Dune::FieldMatrix<Scalar, GridView::dimension, GridView::dimension> stressTensor(0.0);
            for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
            {
                for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
                {
388
389
                    stressTensor[dimIdx][velIdx] = 0.5 * velocityGradients_[elementIdx][dimIdx][velIdx]
                                                   + 0.5 * velocityGradients_[elementIdx][velIdx][dimIdx];
390
391
              }
            }
392
            stressTensorScalarProduct_[elementIdx] = 0.0;
393
394
395
396
            for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
            {
                for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
                {
397
                    stressTensorScalarProduct_[elementIdx] += stressTensor[dimIdx][velIdx] * stressTensor[dimIdx][velIdx];
398
399
400
                }
            }

401
402
403
404
405
            Dune::FieldMatrix<Scalar, GridView::dimension, GridView::dimension> vorticityTensor(0.0);
            for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
            {
                for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
                {
406
407
                    vorticityTensor[dimIdx][velIdx] = 0.5 * velocityGradients_[elementIdx][dimIdx][velIdx]
                                                      - 0.5 * velocityGradients_[elementIdx][velIdx][dimIdx];
408
409
              }
            }
410
            vorticityTensorScalarProduct_[elementIdx] = 0.0;
411
412
413
414
            for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
            {
                for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
                {
415
                    vorticityTensorScalarProduct_[elementIdx] += vorticityTensor[dimIdx][velIdx] * vorticityTensor[dimIdx][velIdx];
416
417
418
                }
            }

419
            auto fvGeometry = localView(this->gridGeometry());
420
421
422
423
            fvGeometry.bindElement(element);
            for (auto&& scv : scvs(fvGeometry))
            {
                const int dofIdx = scv.dofIndex();
424
425

                // construct a privars object from the cell center solution vector
426
                const auto& cellCenterPriVars = curSol[GridGeometry::cellCenterIdx()][dofIdx];
427
                PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars);
428
                auto elemSol = elementSolution<typename GridGeometry::LocalView>(std::move(priVars));
429

430
                VolumeVariables volVars;
431
                volVars.update(elemSol, asImp_(), element, scv);
432
                kinematicViscosity_[elementIdx] = volVars.viscosity() / volVars.density();
433
            }
434
        }
435
436
    }

437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
    /*!
     * \brief Returns whether a wall function should be used at a given face
     *
     * \param element The element.
     * \param scvf The sub control volume face.
     * \param eqIdx The equation index.
     */
    bool useWallFunction(const Element& element,
                         const SubControlVolumeFace& scvf,
                         const int& eqIdx) const
    { return false; }

    /*!
     * \brief Returns an additional wall function momentum flux
     */
    template<class ElementVolumeVariables, class ElementFaceVariables>
    FacePrimaryVariables wallFunction(const Element& element,
                                      const FVElementGeometry& fvGeometry,
                                      const ElementVolumeVariables& elemVolVars,
                                      const ElementFaceVariables& elemFaceVars,
                                      const SubControlVolumeFace& scvf,
458
                                      const SubControlVolumeFace& lateralBoundaryFace) const
459
460
461
462
463
464
465
466
467
468
469
470
471
    { return FacePrimaryVariables(0.0); }

    /*!
     * \brief  Returns an additional wall function flux for cell-centered quantities
     */
    template<class ElementVolumeVariables, class ElementFaceVariables>
    CellCenterPrimaryVariables wallFunction(const Element& element,
                                            const FVElementGeometry& fvGeometry,
                                            const ElementVolumeVariables& elemVolVars,
                                            const ElementFaceVariables& elemFaceVars,
                                            const SubControlVolumeFace& scvf) const
    { return CellCenterPrimaryVariables(0.0); }

472
473
474
475
476
477
478
479
480
481
    /*!
     * \brief Returns whether a given sub control volume face is on a wall
     *
     * \param scvf The sub control volume face.
     */
    bool isOnWall(const SubControlVolumeFace& scvf) const
    {
        return asImp_().isOnWallAtPos(scvf.center());
    }

482
483
    /*!
     * \brief Returns whether a given point is on a wall
484
     *
485
     * \param globalPos The position in global coordinates.
486
     */
487
    bool isOnWallAtPos(const GlobalPosition &globalPos) const
488
489
490
491
    {
        // Throw an exception if no walls are implemented
        DUNE_THROW(Dune::InvalidStateException,
                   "The problem does not provide an isOnWall() method.");
Thomas Fetzer's avatar
Thomas Fetzer committed
492
493
    }

494
495
496
497
498
499
    bool hasChannelGeometry() const
    {
        static const bool hasChannelGeometry = getParamFromGroup<bool>(this->paramGroup(), "RANS.HasChannelGeometry");
        return hasChannelGeometry;
    }

500
501
502
503
504
505
506
507
508
509
    /*!
     * \brief Returns the sand-grain roughness \f$\mathrm{[m]}\f$ at a given position
     *
     * \param globalPos The position in global coordinates.
     */
    Scalar sandGrainRoughnessAtPos(const GlobalPosition &globalPos) const
    {
        return 0.0;
    }

510
511
512
513
514
515
    /*!
     * \brief Returns the Karman constant
     */
    const Scalar karmanConstant() const
    { return 0.41; }

516
    //! \brief Returns the \f$ \beta_{\omega} \f$ constant
517
518
519
520
521
    const Scalar betaOmega() const
    {
        return 0.0708;
    }

522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
    /*!
     * \brief Return the turbulent Prandtl number \f$ [-] \f$ which is used to convert
     *        the eddy viscosity to an eddy thermal conductivity
     */
    Scalar turbulentPrandtlNumber() const
    {
        static const Scalar turbulentPrandtlNumber
            = getParamFromGroup<Scalar>(this->paramGroup(), "RANS.TurbulentPrandtlNumber", 1.0);
        return turbulentPrandtlNumber;
    }

    /*!
     * \brief Return the turbulent Schmidt number \f$ [-] \f$ which is used to convert
     *        the eddy viscosity to an eddy diffusivity
     */
    Scalar turbulentSchmidtNumber() const
    {
        static const Scalar turbulentSchmidtNumber
            = getParamFromGroup<Scalar>(this->paramGroup(), "RANS.TurbulentSchmidtNumber", 1.0);
        return turbulentSchmidtNumber;
    }

544
public:
545
546
547
548
549
550
551
552
553
554
555
556
557
    int wallNormalAxis(const int elementIdx) const
    { return wallNormalAxis_[elementIdx]; }

    int flowNormalAxis(const int elementIdx) const
    { return flowNormalAxis_[elementIdx]; }

    unsigned int wallElementIndex(const int elementIdx) const
    { return wallElementIdx_[elementIdx]; }

    Scalar wallDistance(const int elementIdx) const
    { return wallDistance_[elementIdx]; }


558
    bool calledUpdateStaticWallProperties = false;
559
560
561
562
563
564

    const int fixedFlowNormalAxis_ = getParam<int>("RANS.FlowNormalAxis", 0);
    const int fixedWallNormalAxis_ = getParam<int>("RANS.WallNormalAxis", 1);

    std::vector<unsigned int> wallNormalAxis_;
    std::vector<unsigned int> flowNormalAxis_;
565
    std::vector<Scalar> wallDistance_;
566
    std::vector<unsigned int> wallElementIdx_;
567
    std::vector<std::array<std::array<unsigned int, 2>, dim>> neighborIdx_;
568
    std::vector<GlobalPosition> cellCenter_;
569
570
571
572
    std::vector<DimVector> velocity_;
    std::vector<DimVector> velocityMaximum_;
    std::vector<DimVector> velocityMinimum_;
    std::vector<DimMatrix> velocityGradients_;
573

574
    std::vector<Scalar> stressTensorScalarProduct_;
575
    std::vector<Scalar> vorticityTensorScalarProduct_;
576

577
    std::vector<Scalar> kinematicViscosity_;
578
    std::vector<Scalar> sandGrainRoughness_;
Thomas Fetzer's avatar
Thomas Fetzer committed
579

580
private:
Thomas Fetzer's avatar
Thomas Fetzer committed
581
582
583
584
585
586
587
588
589
    //! Returns the implementation of the problem (i.e. static polymorphism)
    Implementation &asImp_()
    { return *static_cast<Implementation *>(this); }

    //! \copydoc asImp_()
    const Implementation &asImp_() const
    { return *static_cast<const Implementation *>(this); }
};

590
} // end namespace Dumux
Thomas Fetzer's avatar
Thomas Fetzer committed
591
592

#endif