localresidual.hh 16.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- 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 2 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
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/methods.hh>
29
#include <dumux/assembly/staggeredlocalresidual.hh>
30
#include <dune/common/hybridutilities.hh>
31
#include <dumux/freeflow/nonisothermal/localresidual.hh>
32
33
34
35

namespace Dumux
{

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

40
41
42
43
/*!
 * \ingroup NavierStokesModel
 * \brief Element-wise calculation of the Navier-Stokes residual for models using the staggered discretization
 */
44
template<class TypeTag>
45
class NavierStokesResidualImpl<TypeTag, DiscretizationMethod::staggered>
46
: public StaggeredLocalResidual<TypeTag>
47
48
49
{
    using ParentType = StaggeredLocalResidual<TypeTag>;
    friend class StaggeredLocalResidual<TypeTag>;
50
51
52
53
54
55
56
57
58
59
60
61

    using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);

    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;
62
63
64
65

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

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

80
81
    using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);

82
public:
83
    using EnergyLocalResidual = FreeFlowEnergyLocalResidual<FVGridGeometry, FluxVariables, ModelTraits::enableEnergyBalance(), (ModelTraits::numComponents() > 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);
183
184
    }

185
186
187
188
189
190
191
    /*!
     * \brief Sets a fixed Dirichlet value for a cell (such as pressure) at the boundary.
     *        This is a provisional alternative to setting the Dirichlet value on the boundary directly.
     */
    template<class BoundaryTypes>
    void setFixedCell(CellCenterResidual& residual,
                      const Problem& problem,
192
                      const Element& element,
193
194
195
196
197
                      const SubControlVolume& insideScv,
                      const ElementVolumeVariables& elemVolVars,
                      const BoundaryTypes& bcTypes) const
    {
        // set a fixed pressure for cells adjacent to a wall
198
        if(bcTypes.isDirichletCell(Indices::pressureIdx))
199
200
        {
            const auto& insideVolVars = elemVolVars[insideScv];
201
            residual[Indices::pressureIdx - cellCenterOffset] = insideVolVars.pressure() - problem.dirichlet(element, insideScv)[Indices::pressureIdx];
202
203
204
        }
    }

205
206
207
208
209
protected:

     /*!
     * \brief Evaluate boundary conditions for a cell center dof
     */
210
    template<class ElementBoundaryTypes>
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
    void evalBoundaryForCellCenter_(CellCenterResidual& residual,
                                    const Problem& problem,
                                    const Element& element,
                                    const FVElementGeometry& fvGeometry,
                                    const ElementVolumeVariables& elemVolVars,
                                    const ElementFaceVariables& elemFaceVars,
                                    const ElementBoundaryTypes& elemBcTypes,
                                    const ElementFluxVariablesCache& elemFluxVarsCache) const
    {
        for (auto&& scvf : scvfs(fvGeometry))
        {
            if (scvf.boundary())
            {
                const auto bcTypes = problem.boundaryTypes(element, scvf);

226
227
228
229
230
                // treat Dirichlet and outflow BCs
                FluxVariables fluxVars;
                auto boundaryFlux = fluxVars.computeMassFlux(problem, element, fvGeometry, elemVolVars,
                                                             elemFaceVars, scvf, elemFluxVarsCache[scvf]);

231
                EnergyLocalResidual::heatFlux(boundaryFlux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
232
233

                // treat Neumann BCs, i.e. overwrite certain fluxes by user-specified values
234
235
                if(bcTypes.hasNeumann())
                {
236
237
                    static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
                    for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
238
239
                    {
                        if(bcTypes.isNeumann(eqIdx + cellCenterOffset))
240
                        {
241
242
                            const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
                            boundaryFlux[eqIdx] = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx + cellCenterOffset]
243
                                                  * extrusionFactor * scvf.area();
244
                        }
245
                    }
246
247
248
249
                }

                residual += boundaryFlux;

250
251
                // if specified, set a fixed value at the center of a cell at the boundary
                const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
252
                asImp_().setFixedCell(residual, problem, element, scv, elemVolVars, bcTypes);
253
254
255
256
257
258
259
            }
        }
    }

     /*!
     * \brief Evaluate boundary conditions for a face dof
     */
260
    template<class ElementBoundaryTypes>
261
262
263
264
265
266
    void evalBoundaryForFace_(FaceResidual& residual,
                              const Problem& problem,
                              const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const SubControlVolumeFace& scvf,
                              const ElementVolumeVariables& elemVolVars,
267
                              const ElementFaceVariables& elemFaceVars,
268
269
270
271
272
273
274
275
                              const ElementBoundaryTypes& elemBcTypes,
                              const ElementFluxVariablesCache& elemFluxVarsCache) const
    {
        if (scvf.boundary())
        {
            // handle the actual boundary conditions:
            const auto bcTypes = problem.boundaryTypes(element, scvf);

276
            if(bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
277
            {
278
                // set a fixed value for the velocity for Dirichlet boundary conditions
279
                const Scalar velocity = elemFaceVars[scvf].velocitySelf();
280
                const Scalar dirichletValue = problem.dirichlet(element, scvf)[Indices::velocity(scvf.directionIndex())];
281
282
                residual = velocity - dirichletValue;
            }
283
            else if(bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
284
            {
285
                // set a given Neumann flux for the face on the boundary itself
286
                const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
287
                residual = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
288
289
290
291
                                           * extrusionFactor * scvf.area();

                // treat the remaining (normal) faces of the staggered control volume
                FluxVariables fluxVars;
292
                residual += fluxVars.computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars);
293
            }
294
            else if(bcTypes.isSymmetry())
295
            {
296
297
                // For symmetry boundary conditions, there is no flow accross the boundary and
                // we therefore treat it like a Dirichlet boundary conditions with zero velocity
298
                const Scalar velocity = elemFaceVars[scvf].velocitySelf();
299
300
301
                const Scalar fixedValue = 0.0;
                residual = velocity - fixedValue;
            }
302
            else if(bcTypes.isDirichlet(Indices::pressureIdx))
303
            {
304
305
                // If none of the above conditions apply, we are at an "fixed pressure" boundary for which the velocity needs to be assembled
                residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
306
            }
307
308
            else
                DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the boundary conditions for the momentum equations.");
309
310
311
312
313
314
315
316
317
318
319
320
321
        }
    }

private:

    //! 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); }
};
322
323

} // end namespace Dumux
324
325

#endif   // DUMUX_STAGGERED_NAVIERSTOKES_LOCAL_RESIDUAL_HH