localresidual.hh 18.8 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
32
33
34

namespace Dumux
{

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
51
52
53
54
55
56
57
58
59
60

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

    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
65
66
67
    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
    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
72
73
    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
    using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
    using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
74
    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
75

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

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

81
82
public:

83
84
85
86
    // 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");

87
    //! Use the parent type's constructor
88
89
    using ParentType::ParentType;

90
    //! Evaluate fluxes entering or leaving the cell center control volume.
91
92
93
94
95
96
97
98
99
100
101
102
    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;
        CellCenterPrimaryVariables flux = fluxVars.computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars,
                                                 elemFaceVars, scvf, elemFluxVarsCache[scvf]);

103
        computeFluxForCellCenterNonIsothermal_(std::integral_constant<bool, ModelTraits::enableEnergyBalance()>(),
104
                                               problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache, flux);
105
106
107
108

        return flux;
    }

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

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

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

        return result;
127
128
129
    }


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

138
        computeStorageForCellCenterNonIsothermal_(std::integral_constant<bool, ModelTraits::enableEnergyBalance() >(),
139
                                                  problem, scv, volVars, storage);
140

141
142
143
        return storage;
    }

144
    //! Evaluate the storage term for the face control volume.
145
146
147
148
149
150
151
152
153
154
155
    FacePrimaryVariables computeStorageForFace(const Problem& problem,
                                               const SubControlVolumeFace& scvf,
                                               const VolumeVariables& volVars,
                                               const ElementFaceVariables& elementFaceVars) const
    {
        FacePrimaryVariables storage(0.0);
        const Scalar velocity = elementFaceVars[scvf].velocitySelf();
        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
162
163
164
165
166
167
                                              const SubControlVolumeFace& scvf,
                                              const ElementVolumeVariables& elemVolVars,
                                              const ElementFaceVariables& elementFaceVars) const
    {
        FacePrimaryVariables source(0.0);
        const auto insideScvIdx = scvf.insideScvIdx();
        const auto& insideVolVars = elemVolVars[insideScvIdx];
        source += problem.gravity()[scvf.directionIndex()] * insideVolVars.density();
168
        source += problem.source(element, fvGeometry, elemVolVars, elementFaceVars, scvf)[Indices::velocity(scvf.directionIndex())];
169
170
171
172

        return source;
    }

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

protected:

188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
    //  /*!
    //  * \brief Evaluate boundary conditions
    //  */
    // template<class ElementBoundaryTypes>
    // void evalBoundary_(const Element& element,
    //                    const FVElementGeometry& fvGeometry,
    //                    const ElementVolumeVariables& elemVolVars,
    //                    const ElementFaceVariables& elemFaceVars,
    //                    const ElementBoundaryTypes& elemBcTypes,
    //                    const ElementFluxVariablesCache& elemFluxVarsCache)
    // {
    //     evalBoundaryForCellCenter_(element, fvGeometry, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
    //     for (auto&& scvf : scvfs(fvGeometry))
    //     {
    //         evalBoundaryForFace_(element, fvGeometry, scvf, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
    //     }
    // }
205
206
207
208

     /*!
     * \brief Evaluate boundary conditions for a cell center dof
     */
209
    template<class ElementBoundaryTypes>
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
    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())
            {
                auto boundaryFlux = computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);

                // handle the actual boundary conditions:
                const auto bcTypes = problem.boundaryTypes(element, scvf);

                if(bcTypes.hasNeumann())
                {
230
                    static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
231
                    // handle Neumann BCs, i.e. overwrite certain fluxes by user-specified values
232
                    for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
233
234
                    {
                        if(bcTypes.isNeumann(eqIdx + cellCenterOffset))
235
236
                        {
                            const auto extrusionFactor = 1.0; //TODO: get correct extrusion factor
237
238
                            boundaryFlux[eqIdx] = problem.neumann(element, fvGeometry, elemVolVars, scvf)[eqIdx + cellCenterOffset]
                                                  * extrusionFactor * scvf.area();
239
                        }
240
                    }
241
242
243
244
245
246
247
248
249
250
251
252
253
                }

                residual += boundaryFlux;

                asImp_().setFixedCell_(residual, problem, fvGeometry.scv(scvf.insideScvIdx()), elemVolVars, bcTypes);
            }
        }
    }

    /*!
     * \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.
     */
254
    template<class BoundaryTypes>
255
256
257
258
259
260
261
    void setFixedCell_(CellCenterResidual& residual,
                       const Problem& problem,
                       const SubControlVolume& insideScv,
                       const ElementVolumeVariables& elemVolVars,
                       const BoundaryTypes& bcTypes) const
    {
        // set a fixed pressure for cells adjacent to a wall
262
        if(bcTypes.isDirichletCell(Indices::conti0EqIdx))
263
264
        {
            const auto& insideVolVars = elemVolVars[insideScv];
265
            residual[Indices::conti0EqIdx - cellCenterOffset] = insideVolVars.pressure() - problem.dirichletAtPos(insideScv.center())[Indices::pressureIdx];
266
267
268
269
270
271
        }
    }

     /*!
     * \brief Evaluate boundary conditions for a face dof
     */
272
    template<class ElementBoundaryTypes>
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
    void evalBoundaryForFace_(FaceResidual& residual,
                              const Problem& problem,
                              const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const SubControlVolumeFace& scvf,
                              const ElementVolumeVariables& elemVolVars,
                              const ElementFaceVariables& elementFaceVars,
                              const ElementBoundaryTypes& elemBcTypes,
                              const ElementFluxVariablesCache& elemFluxVarsCache) const
    {
        if (scvf.boundary())
        {
            // handle the actual boundary conditions:
            const auto bcTypes = problem.boundaryTypes(element, scvf);

            // set a fixed value for the velocity for Dirichlet boundary conditions
289
            if(bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
290
291
            {
                const Scalar velocity = elementFaceVars[scvf].velocitySelf();
292
                const Scalar dirichletValue = problem.dirichlet(element, scvf)[Indices::velocity(scvf.directionIndex())];
293
294
295
296
297
298
299
300
301
302
303
304
305
306
                residual = velocity - dirichletValue;
            }

            // For symmetry boundary conditions, there is no flow accross the boundary and
            // we therefore treat it like a Dirichlet boundary conditions with zero velocity
            if(bcTypes.isSymmetry())
            {
                // const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
                const Scalar velocity = elementFaceVars[scvf].velocitySelf();
                const Scalar fixedValue = 0.0;
                residual = velocity - fixedValue;
            }

            // outflow condition for the momentum balance equation
307
            if(bcTypes.isOutflow(Indices::velocity(scvf.directionIndex())))
308
            {
309
                if(bcTypes.isDirichlet(Indices::conti0EqIdx))
310
311
312
313
314
315
316
                    residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elementFaceVars, elemFluxVarsCache);
                else
                    DUNE_THROW(Dune::InvalidStateException, "Face at " << scvf.center()  << " has an outflow BC for the momentum balance but no Dirichlet BC for the pressure!");
            }
        }
    }

317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
    //! Evaluate energy fluxes entering or leaving the cell center control volume for non isothermal models
    void computeFluxForCellCenterNonIsothermal_(std::true_type,
                                                const Problem& problem,
                                                const Element &element,
                                                const FVElementGeometry& fvGeometry,
                                                const ElementVolumeVariables& elemVolVars,
                                                const ElementFaceVariables& elemFaceVars,
                                                const SubControlVolumeFace &scvf,
                                                const ElementFluxVariablesCache& elemFluxVarsCache,
                                                CellCenterPrimaryVariables& flux) const
    {
        // if we are on an inflow/outflow boundary, use the volVars of the element itself
        // TODO: catch neumann and outflow in localResidual's evalBoundary_()
        bool isOutflow = false;
        if(scvf.boundary())
        {
            const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
                if(bcTypes.isOutflow(Indices::energyBalanceIdx))
                    isOutflow = true;
        }

        auto upwindTerm = [](const auto& volVars) { return volVars.density() * volVars.enthalpy(); };
        using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);

341
342
        flux[Indices::energyBalanceIdx - cellCenterOffset] = FluxVariables::advectiveFluxForCellCenter(elemVolVars, elemFaceVars, scvf, upwindTerm, isOutflow);
        flux[Indices::energyBalanceIdx - cellCenterOffset] += HeatConductionType::diffusiveFluxForCellCenter(problem, element, fvGeometry, elemVolVars, scvf);
343
344
345
346
347
348
349
350
351
352
353
354
355
    }

    //! Evaluate energy fluxes entering or leaving the cell center control volume for non isothermal models
    template <typename... Args>
    void computeFluxForCellCenterNonIsothermal_(std::false_type, Args&&... args) const {}

    //! Evaluate energy storage for non isothermal models
    void computeStorageForCellCenterNonIsothermal_(std::true_type,
                                                   const Problem& problem,
                                                   const SubControlVolume& scv,
                                                   const VolumeVariables& volVars,
                                                   CellCenterPrimaryVariables& storage) const
    {
356
        storage[Indices::energyBalanceIdx - cellCenterOffset] = volVars.density() * volVars.internalEnergy();
357
358
359
360
361
    }

    //! Evaluate energy storage for isothermal models
    template <typename... Args>
    void computeStorageForCellCenterNonIsothermal_(std::false_type, Args&&... args) const {}
362
363
364
365
366
367
368
369
370
371
372

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); }
};
373
374

} // end namespace Dumux
375
376

#endif   // DUMUX_STAGGERED_NAVIERSTOKES_LOCAL_RESIDUAL_HH