stokesmodel.hh 8.69 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/>.   *
 *****************************************************************************/
#ifndef DUMUX_STOKES_MODEL_HH
#define DUMUX_STOKES_MODEL_HH

22
23
24
25
26
/*!
 * \file
 * \brief Base class for all models which use the Stokes box model.
 */

27
28
#include <dune/common/version.hh>

29
#include <dumux/implicit/common/implicitmodel.hh>
Klaus Mosthaf's avatar
Klaus Mosthaf committed
30

31
#include "stokeslocalresidual.hh"
32
#include "stokesnewtoncontroller.hh"
33
34
35
36
37
38
39
#include "stokeslocaljacobian.hh"
#include "stokesproblem.hh"
#include "stokesproperties.hh"

namespace Dumux
{
/*!
40
 * \ingroup BoxStokesModel
Thomas Fetzer's avatar
[docu]    
Thomas Fetzer committed
41
 * \brief Adaption of the box scheme to the Stokes model.
42
 *
43
 * This model implements laminar Stokes flow of a single fluid, solving the momentum balance equation
44
 * \f[
Thomas Fetzer's avatar
Thomas Fetzer committed
45
46
47
48
49
 *    \frac{\partial \left(\varrho_g {\boldsymbol{v}}_g\right)}{\partial t}
 *    + \text{div} \left( p_g {\bf {I}}
 *    - \mu_g \left( \textbf{grad}\, \boldsymbol{v}_g
 *                   + \textbf{grad}\, \boldsymbol{v}_g^T \right) \right)
 *    - \varrho_g {\bf g} = 0
50
 * \f]
51
 * By setting the property <code>EnableNavierStokes</code> to <code>true</code> the Navier-Stokes
52
 * equation can be solved. In this case an additional term
53
 * \f[
Thomas Fetzer's avatar
Thomas Fetzer committed
54
 *    + \text{div} \left( \varrho_g \boldsymbol{v}_g \boldsymbol{v}_g \right)
55
 * \f]
56
 * is added to the momentum balance equation.
57
 *
Thomas Fetzer's avatar
Thomas Fetzer committed
58
59
60
61
62
63
 * The mass balance equation:
 * \f[
 *    \frac{\partial \varrho_g}{\partial t}
 *    + \text{div} \left(\varrho_g {\boldsymbol{v}}_g\right) - q_g = 0
 * \f]
 *
Klaus Mosthaf's avatar
Klaus Mosthaf committed
64
 * This is discretized by a fully-coupled vertex-centered finite volume
65
 * (box) scheme in space and by the implicit Euler method in time.
66
 */
67
template<class TypeTag>
68
class StokesModel : public GET_PROP_TYPE(TypeTag, BaseModel)
69
70
71
72
73
74
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;

    enum {
        dim = GridView::dimension,
75
        dimWorld = GridView::dimensionworld
76
    };
Klaus Mosthaf's avatar
Klaus Mosthaf committed
77
    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;

    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;

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

public:
    /*!
     * \brief Calculate the fluxes across a certain layer in the domain.
     * The layer is situated perpendicular to the coordinate axis "coord" and cuts
Klaus Mosthaf's avatar
Klaus Mosthaf committed
94
     * the axis at the value "coordVal".
95
     *
Klaus Mosthaf's avatar
Klaus Mosthaf committed
96
97
98
99
     * \param globalSol The global solution vector
     * \param flux A vector to store the flux
     * \param axis The dimension, perpendicular to which the layer is situated
     * \param coordVal The (Scalar) coordinate on the axis, at which the layer is situated
100
     */
Klaus Mosthaf's avatar
Klaus Mosthaf committed
101
    void calculateFluxAcrossLayer(const SolutionVector &globalSol, Dune::FieldVector<Scalar, numEq> &flux, int axis, Scalar coordVal)
102
103
104
105
    {
        GlobalPosition globalI, globalJ;
        PrimaryVariables tmpFlux(0.0);

106
        FVElementGeometry fvGeometry;
107
108
109
        ElementVolumeVariables elemVolVars;

        // Loop over elements
110
111
112
        ElementIterator eIt = this->problem_.gridView().template begin<0>();
        ElementIterator eEndIt = this->problem_.gridView().template end<0>();
        for (; eIt != eEndIt; ++eIt)
113
        {
114
            if (eIt->partitionType() != Dune::InteriorEntity)
115
116
                continue;

117
118
119
            fvGeometry.update(this->gridView_(), *eIt);
            elemVolVars.update(this->problem_(), *eIt, fvGeometry);
            this->localResidual().evalFluxes(*eIt, elemVolVars);
120
121
122

            bool hasLeft = false;
            bool hasRight = false;
123
            for (int i = 0; i < fvGeometry.numScv; i++) {
124
                const GlobalPosition &globalPos = fvGeometry.subContVol[i].global;
Klaus Mosthaf's avatar
Klaus Mosthaf committed
125
                if (globalI[axis] < coordVal)
126
                    hasLeft = true;
Klaus Mosthaf's avatar
Klaus Mosthaf committed
127
                else if (globalI[axis] >= coordVal)
128
129
130
131
132
                    hasRight = true;
            }
            if (!hasLeft || !hasRight)
                continue;

133
            for (int i = 0; i < fvGeometry.numScv; i++) {
134
                const GlobalPosition &globalPos = fvGeometry.subContVol[i].global;
Klaus Mosthaf's avatar
Klaus Mosthaf committed
135
                if (globalI[axis] < coordVal)
136
137
138
139
140
141
142
                    flux += this->localResidual().residual(i);
            }
        }

        flux = this->problem_.gridView().comm().sum(flux);
    }

Bernd Flemisch's avatar
Bernd Flemisch committed
143
    //! \copydoc ImplicitModel::addOutputVtkFields
144
    template <class MultiWriter>
145
146
    void addOutputVtkFields(const SolutionVector &sol,
                            MultiWriter &writer)
147
148
149
150
151
152
    {
        typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
        typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > VelocityField;

        // create the required scalar fields
        unsigned numVertices = this->gridView_().size(dim);
153
        ScalarField &pn = *writer.allocateManagedBuffer(numVertices);
154
155
156
157
158
159
160
161
        ScalarField &delP = *writer.allocateManagedBuffer(numVertices);
        ScalarField &rho = *writer.allocateManagedBuffer(numVertices);
        ScalarField &mu = *writer.allocateManagedBuffer(numVertices);
        VelocityField &velocity = *writer.template allocateManagedBuffer<Scalar, dim> (numVertices);

        unsigned numElements = this->gridView_().size(0);
        ScalarField &rank = *writer.allocateManagedBuffer(numElements);

162
        FVElementGeometry fvGeometry;
163
164
165
        VolumeVariables volVars;
        ElementBoundaryTypes elemBcTypes;

166
167
168
        ElementIterator eIt = this->gridView_().template begin<0>();
        ElementIterator eEndIt = this->gridView_().template end<0>();
        for (; eIt != eEndIt; ++eIt)
169
        {
170
171
172
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
            int eIdx = this->elementMapper().index(*eIt);
#else
173
            int eIdx = this->elementMapper().map(*eIt);
174
#endif
175
            rank[eIdx] = this->gridView_().comm().rank();
176

177
178
            fvGeometry.update(this->gridView_(), *eIt);
            elemBcTypes.update(this->problem_(), *eIt);
179

180
181
182
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
            int numLocalVerts = eIt->subEntities(dim);
#else
183
            int numLocalVerts = eIt->template count<dim>();
184
#endif
185
186
            for (int i = 0; i < numLocalVerts; ++i)
            {
187
188
189
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
                int vIdxGlobal = this->vertexMapper().subIndex(*eIt, i, dim);
#else
190
                int vIdxGlobal = this->vertexMapper().map(*eIt, i, dim);
191
#endif
192
                volVars.update(sol[vIdxGlobal],
193
                               this->problem_(),
194
                               *eIt,
195
                               fvGeometry,
196
197
198
                               i,
                               false);

199
200
201
202
203
                pn[vIdxGlobal] = volVars.pressure();
                delP[vIdxGlobal] = volVars.pressure() - 1e5;
                rho[vIdxGlobal] = volVars.density();
                mu[vIdxGlobal] = volVars.dynamicViscosity();
                velocity[vIdxGlobal] = volVars.velocity();
204
            }
205
        }
206
        writer.attachVertexData(pn, "P");
207
208
209
210
211
212
213
214
        writer.attachVertexData(delP, "delP");
        writer.attachVertexData(rho, "rho");
        writer.attachVertexData(mu, "mu");
        writer.attachVertexData(velocity, "v", dim);
    }
};
}

Andreas Lauser's avatar
Andreas Lauser committed
215
216
#include "stokespropertydefaults.hh"

217
#endif