2p2clocalresidual.hh 12.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// $Id: 2p2clocalresidual.hh 3795 2010-06-25 16:08:04Z melanie $
/*****************************************************************************
 *   Copyright (C) 2008-2009 by Klaus Mosthaf                                *
 *   Copyright (C) 2008-2009 by Bernd Flemisch                               *
 *   Copyright (C) 2009-2010 by Andreas Lauser                               *
 *   Institute of Hydraulic Engineering                                      *
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
 *   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, as long as this copyright notice    *
 *   is included in its original form.                                       *
 *                                                                           *
 *   This program is distributed WITHOUT ANY WARRANTY.                       *
 *****************************************************************************/
Melanie Darcis's avatar
doc    
Melanie Darcis committed
18
19
20
21
22
23
24
/*!
 * \file
 *
 * \brief Element-wise calculation of the Jacobian matrix for problems
 *        using the two-phase two-component box model.
 */

25
26
#ifndef DUMUX_NEW_2P2C_LOCAL_RESIDUAL_BASE_HH
#define DUMUX_NEW_2P2C_LOCAL_RESIDUAL_BASE_HH
27
28
29
30
31
32

#include <dumux/boxmodels/common/boxmodel.hh>
#include <dumux/common/math.hh>

#include "2p2cproperties.hh"

33
#include "2p2cvolumevariables.hh"
34

35
#include "2p2cfluxvariables.hh"
36
37
38
39
40
41
42
43
44
45
46
47

#include "2p2cnewtoncontroller.hh"

#include <iostream>
#include <vector>

//#define VELOCITY_OUTPUT 1 // uncomment this line if an output of the velocity is needed

namespace Dumux
{
/*!
 * \ingroup TwoPTwoCModel
Melanie Darcis's avatar
doc    
Melanie Darcis committed
48
49
 * \brief Element-wise calculation of the Jacobian matrix for problems
 *        using the two-phase two-component box model.
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
 *
 * This class is used to fill the gaps in BoxLocalResidual for the 2P-2C flow.
 */
template<class TypeTag>
class TwoPTwoCLocalResidual: public BoxLocalResidual<TypeTag>
{
protected:
    typedef TwoPTwoCLocalResidual<TypeTag> ThisType;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) Implementation;
    typedef BoxLocalResidual<TypeTag> ParentType;

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
69
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;

    typedef TwoPTwoCFluidState<TypeTag> FluidState;

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;

    enum
    {
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld,

        numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
        numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
        numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),

        pressureIdx = Indices::pressureIdx,
        switchIdx = Indices::switchIdx,

        contiLEqIdx = Indices::contiLEqIdx,
        contiGEqIdx = Indices::contiGEqIdx,

        lPhaseIdx = Indices::lPhaseIdx,
        gPhaseIdx = Indices::gPhaseIdx,

        lCompIdx = Indices::lCompIdx,
        gCompIdx = Indices::gCompIdx,

        lPhaseOnly = Indices::lPhaseOnly,
        gPhaseOnly = Indices::gPhaseOnly,
        bothPhases = Indices::bothPhases,

        plSg = TwoPTwoCFormulation::plSg,
        pgSl = TwoPTwoCFormulation::pgSl,
        formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation))
    };

Andreas Lauser's avatar
Andreas Lauser committed
106
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
107
108
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
Andreas Lauser's avatar
Andreas Lauser committed
109
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
110
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
111

Andreas Lauser's avatar
Andreas Lauser committed
112
113
114
    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::ctype CoordScalar;

115
116
117
118
119
120
121
122
123
    typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
    typedef Dune::FieldVector<Scalar, dim> LocalPosition;
    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;

    static const Scalar mobilityUpwindAlpha =
            GET_PROP_VALUE(TypeTag, PTAG(MobilityUpwindAlpha));

public:
Andreas Lauser's avatar
Andreas Lauser committed
124
125
126
    /*!
     * \brief Evaluate the storage term of the current solution in a
     *        single phase.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
127
128
129
     *
     * \param element The element
     * \param phaseIdx The index of the fluid phase
Andreas Lauser's avatar
Andreas Lauser committed
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
     */
    void evalPhaseStorage(const Element &element, int phaseIdx)
    {
        FVElementGeometry fvGeom;
        fvGeom.update(this->gridView_(), element);
        ElementBoundaryTypes bcTypes;
        bcTypes.update(this->problem_(), element, fvGeom);
        ElementVolumeVariables volVars;
        volVars.update(this->problem_(), element, fvGeom, false);
        
        this->residual_.resize(fvGeom.numVertices);
        this->residual_ = 0;

        this->elemPtr_ = &element;
        this->fvElemGeomPtr_ = &fvGeom;
        this->bcTypesPtr_ = &bcTypes;
        this->prevVolVarsPtr_ = 0;
        this->curVolVarsPtr_ = &volVars;
        evalPhaseStorage_(phaseIdx);
    }
    
151
    /*!
Melanie Darcis's avatar
doc    
Melanie Darcis committed
152
     * \brief Evaluate the amount all conservation quantities
153
154
155
156
     *        (e.g. phase mass) within a sub-control volume.
     *
     * The result should be averaged over the volume (e.g. phase mass
     * inside a sub control volume divided by the volume)
Melanie Darcis's avatar
doc    
Melanie Darcis committed
157
158
159
160
     *
     *  \param result The mass of the component within the sub-control volume
     *  \param scvIdx The SCV (sub-control-volume) index
     *  \param usePrevSol Evaluate function with solution of current or previous time step
161
     */
162
    void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
163
164
165
166
167
168
    {
        // if flag usePrevSol is set, the solution from the previous
        // time step is used, otherwise the current solution is
        // used. The secondary variables are used accordingly.  This
        // is required to compute the derivative of the storage term
        // using the implicit euler method.
Andreas Lauser's avatar
Andreas Lauser committed
169
        const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_()
170
                : this->curVolVars_();
Andreas Lauser's avatar
Andreas Lauser committed
171
        const VolumeVariables &volVars = elemVolVars[scvIdx];
172
173
174
175
176
177
178
179

        // compute storage term of all components within all phases
        result = 0;
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
        {
            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
            {
                int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
Andreas Lauser's avatar
Andreas Lauser committed
180
181
182
                result[eqIdx] += volVars.density(phaseIdx)
                        * volVars.saturation(phaseIdx)
                        * volVars.fluidState().massFrac(phaseIdx, compIdx);
183
184
            }
        }
Andreas Lauser's avatar
Andreas Lauser committed
185
        result *= volVars.porosity();
186
187
188
189
    }

    /*!
     * \brief Evaluates the total flux of all conservation quantities
Melanie Darcis's avatar
doc    
Melanie Darcis committed
190
191
192
193
     *        over a face of a sub-control volume.
     *
     * \param flux The flux over the SCV (sub-control-volume) face for each component
     * \param faceIdx The index of the SCV face
194
     */
195
    void computeFlux(PrimaryVariables &flux, int faceIdx) const
196
    {
197
        FluxVariables vars(this->problem_(),
198
                      this->elem_(),
199
200
201
                      this->fvElemGeom_(),
                      faceIdx,
                      this->curVolVars_());
202
203
204
205
206
207
208
209
210

        flux = 0;
        asImp_()->computeAdvectiveFlux(flux, vars);
        asImp_()->computeDiffusiveFlux(flux, vars);
    }

    /*!
     * \brief Evaluates the advective mass flux of all components over
     *        a face of a subcontrol volume.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
211
212
     *
     * \param flux The advective flux over the sub-control-volume face for each component
Bernd Flemisch's avatar
Bernd Flemisch committed
213
     * \param vars The flux variables at the current SCV
214
     */
215
    void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
216
217
218
219
220
221
222
223
    {
        ////////
        // advective fluxes of all components in all phases
        ////////
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
        {
            // data attached to upstream and the downstream vertices
            // of the current phase
224
            const VolumeVariables &up =
225
226
227
                this->curVolVars_(vars.upstreamIdx(phaseIdx));
            const VolumeVariables &dn =
                this->curVolVars_(vars.downstreamIdx(phaseIdx));
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252

            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
            {
                int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
                // add advective flux of current component in current
                // phase
                if (mobilityUpwindAlpha > 0.0)
                    // upstream vertex
                    flux[eqIdx] += vars.KmvpNormal(phaseIdx)
                            * mobilityUpwindAlpha * (up.density(phaseIdx)
                            * up.mobility(phaseIdx) * up.fluidState().massFrac(
                            phaseIdx, compIdx));
                if (mobilityUpwindAlpha < 1.0)
                    // downstream vertex
                    flux[eqIdx] += vars.KmvpNormal(phaseIdx) * (1
                            - mobilityUpwindAlpha) * (dn.density(phaseIdx)
                            * dn.mobility(phaseIdx) * dn.fluidState().massFrac(
                            phaseIdx, compIdx));
            }
        }
    }

    /*!
     * \brief Adds the diffusive mass flux of all components over
     *        a face of a subcontrol volume.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
253
254
     *
     * \param flux The diffusive flux over the sub-control-volume face for each component
Bernd Flemisch's avatar
Bernd Flemisch committed
255
     * \param vars The flux variables at the current SCV
256
     */
257
    void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
258
259
    {
        // add diffusive flux of gas component in liquid phase
260
261
        Scalar tmp =
            - vars.porousDiffCoeff(lPhaseIdx) *
262
263
264
265
266
267
            vars.molarDensityAtIP(lPhaseIdx) *
            (vars.molarConcGrad(lPhaseIdx) * vars.face().normal);
        flux[contiGEqIdx] += tmp * FluidSystem::molarMass(gCompIdx);
        flux[contiLEqIdx] -= tmp * FluidSystem::molarMass(lCompIdx);

        // add diffusive flux of liquid component in gas phase
268
269
        tmp =
            - vars.porousDiffCoeff(gPhaseIdx) *
270
271
272
273
274
275
276
277
            vars.molarDensityAtIP(gPhaseIdx) *
            (vars.molarConcGrad(gPhaseIdx) * vars.face().normal);
        flux[contiLEqIdx] += tmp * FluidSystem::molarMass(lCompIdx);
        flux[contiGEqIdx] -= tmp * FluidSystem::molarMass(gCompIdx);
    }

    /*!
     * \brief Calculate the source term of the equation
Melanie Darcis's avatar
doc    
Melanie Darcis committed
278
279
280
     *
     * \param q The source/sink in the SCV for each component
     * \param localVertexIdx The index of the SCV
281
     */
282
    void computeSource(PrimaryVariables &q, int localVertexIdx)
283
284
285
286
287
288
289
290
    {
        this->problem_().source(q,
                                this->elem_(),
                                this->fvElemGeom_(),
                                localVertexIdx);
    }

protected:
Andreas Lauser's avatar
Andreas Lauser committed
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
    void evalPhaseStorage_(int phaseIdx)
    {
        // evaluate the storage terms of a single phase
        for (int i=0; i < this->fvElemGeom_().numVertices; i++) {
            PrimaryVariables &result = this->residual_[i];
            const ElementVolumeVariables &elemVolVars = this->curVolVars_();
            const VolumeVariables &volVars = elemVolVars[i];
            
            // compute storage term of all components within all phases
            result = 0;
            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
            {
                int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
                result[eqIdx] += volVars.density(phaseIdx)
                    * volVars.saturation(phaseIdx)
                    * volVars.fluidState().massFrac(phaseIdx, compIdx);
            }

            result *= volVars.porosity();
            result *= this->fvElemGeom_().subContVol[i].volume;
        }
    }


315
316
317
318
319
320
321
322
323
324
325
326
327
328
    Implementation *asImp_()
    {
        return static_cast<Implementation *> (this);
    }
    const Implementation *asImp_() const
    {
        return static_cast<const Implementation *> (this);
    }

};

} // end namepace

#endif