stokesfluxvariables.hh 9.81 KB
Newer Older
Andreas Lauser's avatar
Andreas Lauser committed
1
2
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
3
/*****************************************************************************
4
 *   See the file COPYING for full copying permissions.                      *
5
6
7
8
9
10
11
12
 *                                                                           *
 *   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          *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
14
15
16
17
18
19
20
21
 *   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
 * \brief This file contains the data which is required to calculate
Klaus Mosthaf's avatar
Klaus Mosthaf committed
22
 *        the fluxes of the Stokes model over a face of a finite volume.
23
 *
24
 * This means pressure gradients, phase densities at the integration point, etc.
25
26
27
28
29
 */
#ifndef DUMUX_STOKES_FLUX_VARIABLES_HH
#define DUMUX_STOKES_FLUX_VARIABLES_HH

#include <dumux/common/math.hh>
Andreas Lauser's avatar
Andreas Lauser committed
30
31
32
#include <dumux/common/valgrind.hh>

#include "stokesproperties.hh"
33
34
35
36
37

namespace Dumux
{

/*!
38
 * \ingroup BoxStokesModel
39
 * \ingroup ImplicitFluxVariables
40
 * \brief This template class contains the data which is required to
Klaus Mosthaf's avatar
Klaus Mosthaf committed
41
 *        calculate the mass and momentum fluxes over the face of a
42
 *        sub-control volume for the Stokes model.
43
 *
Klaus Mosthaf's avatar
Klaus Mosthaf committed
44
 * This means pressure gradients, phase densities, viscosities, etc.
45
 * at the integration point of the sub-control-volume face.
46
47
48
49
50
51
52
53
54
55
56
57
58
 */
template <class TypeTag>
class StokesFluxVariables
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;

    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;

    enum { dim = GridView::dimension };

    typedef typename GridView::template Codim<0>::Entity Element;
59
60
    typedef Dune::FieldVector<Scalar, dim> DimVector;
    typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
61
62
63
64
65
66
67

    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;

public:
    StokesFluxVariables(const Problem &problem,
                        const Element &element,
68
69
                        const FVElementGeometry &fvGeometry,
                        const int faceIdx,
70
                        const ElementVolumeVariables &elemVolVars,
71
72
                        const bool onBoundary = false)
        : fvGeometry_(fvGeometry), onBoundary_(onBoundary), faceIdx_(faceIdx)
73
74
75
76
77
78
79
80
81
82
83
    {
        calculateValues_(problem, element, elemVolVars);
        determineUpwindDirection_(elemVolVars);
    };

protected:
    void calculateValues_(const Problem &problem,
                          const Element &element,
                          const ElementVolumeVariables &elemVolVars)
    {
        // calculate gradients and secondary variables at IPs
84
85
86
        DimVector tmp(0.0);

        density_ = Scalar(0);
87
        dynamicViscosity_ = Scalar(0);
88
89
90
91
92
93
        pressure_ = Scalar(0);
        normalvelocity_ = Scalar(0);
        velocity_ = Scalar(0);
        pressureGrad_ = Scalar(0);
        velocityGrad_ = Scalar(0);
//        velocityDiv_ = Scalar(0);
94
95

        for (int idx = 0;
96
             idx < fvGeometry_.numScv;
97
98
99
             idx++) // loop over adjacent vertices
        {
            // phase density and viscosity at IP
100
            density_ += elemVolVars[idx].density() *
101
                face().shapeValue[idx];
102
            dynamicViscosity_ += elemVolVars[idx].dynamicViscosity() *
103
                face().shapeValue[idx];
104
            pressure_ += elemVolVars[idx].pressure() *
105
106
107
                face().shapeValue[idx];

            // velocity at the IP (fluxes)
108
            DimVector velocityTimesShapeValue = elemVolVars[idx].velocity();
109
            velocityTimesShapeValue *= face().shapeValue[idx];
110
            velocity_ += velocityTimesShapeValue;
111
112
113
114

            // the pressure gradient
            tmp = face().grad[idx];
            tmp *= elemVolVars[idx].pressure();
115
            pressureGrad_ += tmp;
116
117
            // take gravity into account
            tmp = problem.gravity();
118
            tmp *= density_;
119
            // pressure gradient including influence of gravity
120
            pressureGrad_ -= tmp;
121

Klaus Mosthaf's avatar
Klaus Mosthaf committed
122
            // the velocity gradients and divergence
123
124
125
126
            for (int dimIdx = 0; dimIdx<dim; ++dimIdx)
            {
                tmp = face().grad[idx];
                tmp *= elemVolVars[idx].velocity()[dimIdx];
127
                velocityGrad_[dimIdx] += tmp;
128

129
//                velocityDiv_ += face().grad[idx][dimIdx]*elemVolVars[idx].velocity()[dimIdx];
130
131
132
            }
        }

133
        normalvelocity_ = velocity_ * face().normal;
134

135
        Valgrind::CheckDefined(density_);
136
        Valgrind::CheckDefined(dynamicViscosity_);
137
138
139
140
141
        Valgrind::CheckDefined(normalvelocity_);
        Valgrind::CheckDefined(velocity_);
        Valgrind::CheckDefined(pressureGrad_);
        Valgrind::CheckDefined(velocityGrad_);
//        Valgrind::CheckDefined(velocityDiv_);
142
143
144
145
146
147
148
149
150
    };

    void determineUpwindDirection_(const ElementVolumeVariables &elemVolVars)
    {

        // set the upstream and downstream vertices
        upstreamIdx_ = face().i;
        downstreamIdx_ = face().j;

151
        if (normalVelocity() < 0)
152
153
154
155
156
157
158
159
160
161
162
            std::swap(upstreamIdx_, downstreamIdx_);
    };

public:
    /*!
     * \brief The face of the current sub-control volume. This may be either
     *        an inner sub-control-volume face or a face on the boundary.
     */
    const SCVFace &face() const
    {
        if (onBoundary_)
163
            return fvGeometry_.boundaryFace[faceIdx_];
164
        else
165
            return fvGeometry_.subContVolFace[faceIdx_];
166
167
168
169
    }

    /*!
     * \brief Return the average volume of the upstream and the downstream sub-control volume;
Klaus Mosthaf's avatar
Klaus Mosthaf committed
170
     *        this is required for the stabilization.
171
172
173
     */
    const Scalar averageSCVVolume() const
    {
174
175
        return 0.5*(fvGeometry_.subContVol[upstreamIdx_].volume +
                fvGeometry_.subContVol[downstreamIdx_].volume);
176
177
178
    }

    /*!
Klaus Mosthaf's avatar
Klaus Mosthaf committed
179
     * \brief Return the pressure \f$\mathrm{[Pa]}\f$ at the integration
180
181
     *        point.
     */
182
183
184
185
186
187
188
189
190
    Scalar pressure() const
    { return pressure_; }

    /*!
     * \brief Return the mass density \f$ \mathrm{[kg/m^3]} \f$ at the integration
     *        point.
     */
    Scalar density() const
    { return density_; }
191

192
193
194
195
196
197
198
199
200
201
202
203
204
    /*!
     * \brief Return the dynamic viscosity \f$ \mathrm{[Pa\cdot s]} \f$ at the integration
     *        point.
     */
    Scalar dynamicViscosity() const
    { return dynamicViscosity_; }

    /*!
     * \brief Returns the kinematic viscosity \f$ \frac{m^2}{s} \f$ of the fluid in
     *        the sub-control volume.
     */
    Scalar kinematicViscosity() const
    { return dynamicViscosity_ / density_; }
205

206
    /*!
207
208
     * \brief Return the velocity \f$ \mathrm{[m/s]} \f$ at the integration
     *        point multiplied by the normal and the area.
209
     */
210
211
212
213
214
215
216
217
    Scalar normalVelocity() const
    { return normalvelocity_; }

    /*!
     * \brief Return the pressure gradient at the integration point.
     */
    const DimVector &pressureGrad() const
    { return pressureGrad_; }
218
219

    /*!
Klaus Mosthaf's avatar
Klaus Mosthaf committed
220
     * \brief Return the velocity vector at the integration point.
221
     */
222
223
224
225
226
227
228
229
230
    const DimVector &velocity() const
    { return velocity_; }

    /*!
     * \brief Return the velocity gradient at the integration
     *        point of a face.
     */
    const DimMatrix &velocityGrad() const
    { return velocityGrad_; }
231

232
233
234
235
236
237
238
239
240
241
242
243
    /*!
     * \brief Return the dynamic eddy viscosity 
     *        \f$\mathrm{[Pa \cdot s]} = \mathrm{[N \cdot s/m^2]}\f$ (if implemented).
     */
    const Scalar dynamicEddyViscosity() const
    { return kinematicEddyViscosity() * density(); }

    /*!
     * \brief Return the kinematic eddy viscosity 
     *        \f$\mathrm{[m^2/s]}\f$ (if implemented).
     */
    const Scalar kinematicEddyViscosity() const
244
245
    { return 0; }

246

Klaus Mosthaf's avatar
Klaus Mosthaf committed
247
248
249
250
//    /*!
//     * \brief Return the divergence of the normal velocity at the
//     *        integration point.
//     */
251
252
//    Scalar velocityDiv() const
//    { return velocityDiv_; }
253
254

    /*!
Klaus Mosthaf's avatar
Klaus Mosthaf committed
255
     * \brief Return the local index of the upstream sub-control volume.
256
257
258
259
260
     */
    int upstreamIdx() const
    { return upstreamIdx_; }

    /*!
Klaus Mosthaf's avatar
Klaus Mosthaf committed
261
     * \brief Return the local index of the downstream sub-control volume.
262
263
264
265
     */
    int downstreamIdx() const
    { return downstreamIdx_; }

Klaus Mosthaf's avatar
Klaus Mosthaf committed
266
267
268
269
    /*!
     * \brief Indicates if a face is on a boundary. Used for in the
     *        face() method (e.g. for outflow boundary conditions).
     */
270
271
272
273
    bool onBoundary() const
    { return onBoundary_; }

protected:
274
    const FVElementGeometry &fvGeometry_;
275
276
277
    const bool onBoundary_;

    // values at the integration point
278
    Scalar density_;
279
    Scalar dynamicViscosity_;
280
281
282
283
    Scalar pressure_;
    Scalar normalvelocity_;
//    Scalar velocityDiv_;
    DimVector velocity_;
284
285

    // gradients at the IPs
286
287
    DimVector pressureGrad_;
    DimMatrix velocityGrad_;
288
289
290
291
292
293
294
295
296

    // local index of the upwind vertex
    int upstreamIdx_;
    // local index of the downwind vertex
    int downstreamIdx_;
    // the index of the considered face
    int faceIdx_;
};

297
} // end namespace
298
299

#endif