localresidual.hh 17.8 KB
Newer Older
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       *
9
10
11
12
13
14
15
16
17
18
19
20
 *   (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
21
22
 * \ingroup NavierStokesModel
 * \copydoc Dumux::NavierStokesResidualImpl
23
24
25
26
27
 */
#ifndef DUMUX_STAGGERED_NAVIERSTOKES_LOCAL_RESIDUAL_HH
#define DUMUX_STAGGERED_NAVIERSTOKES_LOCAL_RESIDUAL_HH

#include <dumux/common/properties.hh>
28
#include <dumux/discretization/method.hh>
29
#include <dumux/assembly/staggeredlocalresidual.hh>
30
#include <dune/common/hybridutilities.hh>
31
#include <dumux/freeflow/nonisothermal/localresidual.hh>
32

33
namespace Dumux {
34

35
// forward declaration
36
template<class TypeTag, DiscretizationMethod discMethod>
37
class NavierStokesResidualImpl;
38

39
40
41
42
/*!
 * \ingroup NavierStokesModel
 * \brief Element-wise calculation of the Navier-Stokes residual for models using the staggered discretization
 */
43
template<class TypeTag>
44
class NavierStokesResidualImpl<TypeTag, DiscretizationMethod::staggered>
45
: public StaggeredLocalResidual<TypeTag>
46
47
48
{
    using ParentType = StaggeredLocalResidual<TypeTag>;
    friend class StaggeredLocalResidual<TypeTag>;
49

50
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
51
52
53
54
55
56
57
58
59
60

    using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
    using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
    using VolumeVariables = typename GridVolumeVariables::VolumeVariables;

    using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
    using ElementFluxVariablesCache = typename GridFluxVariablesCache::LocalView;

    using GridFaceVariables = typename GridVariables::GridFaceVariables;
    using ElementFaceVariables = typename GridFaceVariables::LocalView;
61

62
63
64
65
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
66
67
    using FVElementGeometry = typename FVGridGeometry::LocalView;
    using GridView = typename FVGridGeometry::GridView;
68
    using Element = typename GridView::template Codim<0>::Entity;
69
70
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
71
    using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
72
73
74
75
    using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
    using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
    using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
76

77
78
    using CellCenterResidual = CellCenterPrimaryVariables;
    using FaceResidual = FacePrimaryVariables;
79

80
    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
81

82
public:
83
    using EnergyLocalResidual = FreeFlowEnergyLocalResidual<FVGridGeometry, FluxVariables, ModelTraits::enableEnergyBalance(), (ModelTraits::numFluidComponents() > 1)>;
84

85
86
87
88
    // account for the offset of the cell center privars within the PrimaryVariables container
    static constexpr auto cellCenterOffset = ModelTraits::numEq() - CellCenterPrimaryVariables::dimension;
    static_assert(cellCenterOffset == ModelTraits::dim(), "cellCenterOffset must equal dim for staggered NavierStokes");

89
    //! Use the parent type's constructor
90
91
    using ParentType::ParentType;

92
    //! Evaluate fluxes entering or leaving the cell center control volume.
93
94
95
96
97
98
99
100
101
    CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
                                                        const Element &element,
                                                        const FVElementGeometry& fvGeometry,
                                                        const ElementVolumeVariables& elemVolVars,
                                                        const ElementFaceVariables& elemFaceVars,
                                                        const SubControlVolumeFace &scvf,
                                                        const ElementFluxVariablesCache& elemFluxVarsCache) const
    {
        FluxVariables fluxVars;
102
103
        CellCenterPrimaryVariables flux = fluxVars.computeMassFlux(problem, element, fvGeometry, elemVolVars,
                                                                   elemFaceVars, scvf, elemFluxVarsCache[scvf]);
104

105
        EnergyLocalResidual::heatFlux(flux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
106
107
108
109

        return flux;
    }

110
    //! Evaluate the source term for the cell center control volume.
111
112
113
114
115
116
117
    CellCenterPrimaryVariables computeSourceForCellCenter(const Problem& problem,
                                                          const Element &element,
                                                          const FVElementGeometry& fvGeometry,
                                                          const ElementVolumeVariables& elemVolVars,
                                                          const ElementFaceVariables& elemFaceVars,
                                                          const SubControlVolume &scv) const
    {
118
119
120
        CellCenterPrimaryVariables result(0.0);

        // get the values from the problem
121
        const auto sourceValues = problem.source(element, fvGeometry, elemVolVars, elemFaceVars, scv);
122
123

        // copy the respective cell center related values to the result
124
125
        for (int i = 0; i < result.size(); ++i)
            result[i] = sourceValues[i + cellCenterOffset];
126
127

        return result;
128
129
130
    }


131
    //! Evaluate the storage term for the cell center control volume.
132
133
134
135
136
    CellCenterPrimaryVariables computeStorageForCellCenter(const Problem& problem,
                                                           const SubControlVolume& scv,
                                                           const VolumeVariables& volVars) const
    {
        CellCenterPrimaryVariables storage;
137
        storage[Indices::conti0EqIdx - ModelTraits::dim()] = volVars.density();
138

139
        EnergyLocalResidual::fluidPhaseStorage(storage, volVars);
140

141
142
143
        return storage;
    }

144
    //! Evaluate the storage term for the face control volume.
145
146
147
    FacePrimaryVariables computeStorageForFace(const Problem& problem,
                                               const SubControlVolumeFace& scvf,
                                               const VolumeVariables& volVars,
148
                                               const ElementFaceVariables& elemFaceVars) const
149
150
    {
        FacePrimaryVariables storage(0.0);
151
        const Scalar velocity = elemFaceVars[scvf].velocitySelf();
152
153
154
155
        storage[0] = volVars.density() * velocity;
        return storage;
    }

156
    //! Evaluate the source term for the face control volume.
157
    FacePrimaryVariables computeSourceForFace(const Problem& problem,
158
159
                                              const Element& element,
                                              const FVElementGeometry& fvGeometry,
160
161
                                              const SubControlVolumeFace& scvf,
                                              const ElementVolumeVariables& elemVolVars,
162
                                              const ElementFaceVariables& elemFaceVars) const
163
164
    {
        FacePrimaryVariables source(0.0);
165
        const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
166
        source += problem.gravity()[scvf.directionIndex()] * insideVolVars.density();
167
        source += problem.source(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())];
168
169
170
171

        return source;
    }

172
    //! Evaluate the momentum flux for the face control volume.
173
174
175
176
177
    FacePrimaryVariables computeFluxForFace(const Problem& problem,
                                            const Element& element,
                                            const SubControlVolumeFace& scvf,
                                            const FVElementGeometry& fvGeometry,
                                            const ElementVolumeVariables& elemVolVars,
178
                                            const ElementFaceVariables& elemFaceVars,
179
180
181
                                            const ElementFluxVariablesCache& elemFluxVarsCache) const
    {
        FluxVariables fluxVars;
182
        return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache());
183
184
    }

185
    /*!
186
187
     * \brief Evaluate boundary conditions for a cell center dof
     */
188
189
190
191
192
193
194
195
    CellCenterResidual computeBoundaryFluxForCellCenter(const Problem& problem,
                                                        const Element& element,
                                                        const FVElementGeometry& fvGeometry,
                                                        const SubControlVolumeFace& scvf,
                                                        const ElementVolumeVariables& elemVolVars,
                                                        const ElementFaceVariables& elemFaceVars,
                                                        const ElementBoundaryTypes& elemBcTypes,
                                                        const ElementFluxVariablesCache& elemFluxVarsCache) const
196
    {
197
        CellCenterResidual result(0.0);
198

199
200
201
        if (scvf.boundary())
        {
            const auto bcTypes = problem.boundaryTypes(element, scvf);
202

203
            // no fluxes occur over symmetry boundaries
204
205
206
207
208
            if (bcTypes.isSymmetry())
                return result;

            // treat Dirichlet and outflow BCs
            result = computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
209

210
211
212
213
            // treat Neumann BCs, i.e. overwrite certain fluxes by user-specified values
            static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
            if (bcTypes.hasNeumann())
            {
214
215
216
217
                const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
                const auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf);

                for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
218
                {
219
220
                    if (bcTypes.isNeumann(eqIdx + cellCenterOffset))
                        result[eqIdx] = neumannFluxes[eqIdx + cellCenterOffset] * extrusionFactor * scvf.area();
221
                }
222
            }
223
224
225

            // account for wall functions, if used
            incorporateWallFunction_(result, problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars);
226
        }
227
        return result;
228
229
    }

230
    /*!
231
     * \brief Evaluate Dirichlet (fixed value) boundary conditions for a face dof
232
     */
233
234
235
236
237
238
239
240
241
    void evalDirichletBoundariesForFace(FaceResidual& residual,
                                        const Problem& problem,
                                        const Element& element,
                                        const FVElementGeometry& fvGeometry,
                                        const SubControlVolumeFace& scvf,
                                        const ElementVolumeVariables& elemVolVars,
                                        const ElementFaceVariables& elemFaceVars,
                                        const ElementBoundaryTypes& elemBcTypes,
                                        const ElementFluxVariablesCache& elemFluxVarsCache) const
242
243
244
245
246
247
    {
        if (scvf.boundary())
        {
            // handle the actual boundary conditions:
            const auto bcTypes = problem.boundaryTypes(element, scvf);

248
            if(bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
249
            {
250
                // set a fixed value for the velocity for Dirichlet boundary conditions
251
                const Scalar velocity = elemFaceVars[scvf].velocitySelf();
252
                const Scalar dirichletValue = problem.dirichlet(element, scvf)[Indices::velocity(scvf.directionIndex())];
253
254
                residual = velocity - dirichletValue;
            }
255
            else if(bcTypes.isSymmetry())
256
            {
257
                // for symmetry boundary conditions, there is no flow accross the boundary and
258
                // we therefore treat it like a Dirichlet boundary conditions with zero velocity
259
                const Scalar velocity = elemFaceVars[scvf].velocitySelf();
260
261
262
                const Scalar fixedValue = 0.0;
                residual = velocity - fixedValue;
            }
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
        }
    }

    /*!
     * \brief Evaluate boundary boundary fluxes for a face dof
     */
    FaceResidual computeBoundaryFluxForFace(const Problem& problem,
                                            const Element& element,
                                            const FVElementGeometry& fvGeometry,
                                            const SubControlVolumeFace& scvf,
                                            const ElementVolumeVariables& elemVolVars,
                                            const ElementFaceVariables& elemFaceVars,
                                            const ElementBoundaryTypes& elemBcTypes,
                                            const ElementFluxVariablesCache& elemFluxVarsCache) const
    {
        FaceResidual result(0.0);

        if (scvf.boundary())
        {
            // handle the actual boundary conditions:
            const auto bcTypes = problem.boundaryTypes(element, scvf);
            if (bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
            {
                // the source term has already been accounted for, here we
                // add a given Neumann flux for the face on the boundary itself ...
                const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
                result = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
                                         * extrusionFactor * scvf.area();

                // ... and treat the fluxes of the remaining (frontal and lateral) faces of the staggered control volume
                result += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
            }
295
            else if(bcTypes.isDirichlet(Indices::pressureIdx))
296
            {
297
298
                // if none of the above conditions apply, we are at an "fixed pressure" boundary for which the resdiual of the momentum balance needs to be assembled
                // as if it where inside the domain and not on the boundary (source term has already been acounted for)
299
                result = computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
300
301
            }
        }
302
        return result;
303
304
305
306
    }

private:

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
335
336
    //! do nothing if no turbulence model is used
    template<class ...Args, bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<!turbulenceModel, int> = 0>
    void incorporateWallFunction_(Args&&... args) const
    {};

    //! if a turbulence model is used, ask the problem is a wall function shall be employed and get the flux accordingly
    template<bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<turbulenceModel, int> = 0>
    void incorporateWallFunction_(CellCenterResidual& boundaryFlux,
                                  const Problem& problem,
                                  const Element& element,
                                  const FVElementGeometry& fvGeometry,
                                  const SubControlVolumeFace& scvf,
                                  const ElementVolumeVariables& elemVolVars,
                                  const ElementFaceVariables& elemFaceVars) const
    {
        static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
        const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();

        // account for wall functions, if used
        for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
        {
            // use a wall function
            if(problem.useWallFunction(element, scvf, eqIdx + cellCenterOffset))
            {
                boundaryFlux[eqIdx] = problem.wallFunction(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx]
                                                           * extrusionFactor * scvf.area();
            }
        }
    };

337
338
339
340
341
342
343
344
    //! 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); }
};
345
346

} // end namespace Dumux
347
348

#endif   // DUMUX_STAGGERED_NAVIERSTOKES_LOCAL_RESIDUAL_HH