1p2clocalresidual.hh 12.5 KB
Newer Older
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
 *   Copyright (C) 2011 by Katherina Baber                                   *
5
6
7
 *   Copyright (C) 2009 by Karin Erbertseder                                 *
 *   Copyright (C) 2009 by Andreas Lauser                                    *
 *   Copyright (C) 2008 by Bernd Flemisch                                    *
Andreas Lauser's avatar
Andreas Lauser committed
8
 *   Institute for Modelling Hydraulic and Environmental Systems             *
9
10
11
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
12
 *   This program is free software: you can redistribute it and/or modify    *
13
 *   it under the terms of the GNU General Public License as published by    *
14
15
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
16
 *                                                                           *
17
18
19
20
21
22
23
 *   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/>.   *
24
 *****************************************************************************/
Karin Erbertseder's avatar
Karin Erbertseder committed
25
26
27
/*!
 * \file
 *
Karin Erbertseder's avatar
Karin Erbertseder committed
28
 * \brief Element-wise calculation the local Jacobian for the single-phase,
Karin Erbertseder's avatar
Karin Erbertseder committed
29
30
31
 *        two-component model in the BOX scheme.
 */

32
33
#ifndef DUMUX_ONEP_TWOC_LOCAL_RESIDUAL_HH
#define DUMUX_ONEP_TWOC_LOCAL_RESIDUAL_HH
Katherina Baber's avatar
Katherina Baber committed
34
#define VELOCITY_OUTPUT 1 //1 turns velocity output on, 0 turns it off
35

36
#include <dumux/boxmodels/common/boxmodel.hh>
37
38

#include <dumux/boxmodels/1p2c/1p2cproperties.hh>
39
40
#include <dumux/boxmodels/1p2c/1p2cvolumevariables.hh>
#include <dumux/boxmodels/1p2c/1p2cfluxvariables.hh>
41
42
43
44
45
46
47
48

#include <dune/common/collectivecommunication.hh>
#include <vector>
#include <iostream>

namespace Dumux
{
/*!
49
 *
50
 * \ingroup OnePTwoCBoxModel
51
 * \ingroup BoxLocalResidual
52
53
 * \brief Calculate the local Jacobian for the single-phase,
 *        two-component model in the BOX scheme.
54
 *
Karin Erbertseder's avatar
Karin Erbertseder committed
55
 *  This class is used to fill the gaps in BoxLocalResidual for the 1P-2C flow.
56
57
 */
template<class TypeTag>
58
class OnePTwoCLocalResidual : public BoxLocalResidual<TypeTag>
59
60
{
protected:
61
    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
62

63
64
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
Katherina Baber's avatar
Katherina Baber committed
65
66
    typedef typename GridView::IntersectionIterator IntersectionIterator;

67
68
69
70
71
    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
72

73
    typedef typename GET_PROP_TYPE(TypeTag, OnePTwoCIndices) Indices;
Katherina Baber's avatar
Katherina Baber committed
74

75
    enum
76
77
        {
            numEq = GET_PROP_VALUE(TypeTag, NumEq),
78

79
            // indices of the primary variables
80

81
            comp1Idx = Indices::comp1Idx,
82

83
84
85
86
            // indices of the equations
            contiEqIdx = Indices::contiEqIdx,
            transEqIdx = Indices::transEqIdx
        };
87

88
    static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
Andreas Lauser's avatar
Andreas Lauser committed
89

90
91
92


public:
93
94
95
96
97
    /*!
     * \brief Constructor. Sets the upwind weight.
     */
    OnePTwoCLocalResidual()
    {
98
        // retrieve the upwind weight for the mass conservation equations. Use the value
99
100
        // specified via the property system as default, and overwrite
        // it by the run-time parameter from the Dune::ParameterTree
101
        upwindWeight_ = GET_PARAM(TypeTag, Scalar, UpwindWeight);
102
103
    };

104
    /*!
Karin Erbertseder's avatar
Karin Erbertseder committed
105
     * \brief Evaluate the amount of all conservation quantities
106
     *        (e.g. phase mass) within a finite volume.
Karin Erbertseder's avatar
Karin Erbertseder committed
107
     *
Karin Erbertseder's avatar
Karin Erbertseder committed
108
     *        \param result The mass of the component within the sub-control volume
Karin Erbertseder's avatar
Karin Erbertseder committed
109
     *        \param scvIdx The index of the considered face of the sub control volume
Karin Erbertseder's avatar
Karin Erbertseder committed
110
     *        \param usePrevSol Evaluate function with solution of current or previous time step
111
     */
112
    void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
113
114
115
116
117
118
    {
        // 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.
119
120
121
122
        const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
        const VolumeVariables &volVars = elemVolVars[scvIdx];

        result = 0;
123
        if(!useMoles)
124
125
126
        {
            // storage term of continuity equation - massfractions
            result[contiEqIdx] +=
Andreas Lauser's avatar
Andreas Lauser committed
127
                volVars.fluidState().density(/*phaseIdx*/0)*volVars.porosity();
128
129
            //storage term of the transport equation - massfractions
            result[transEqIdx] +=
Andreas Lauser's avatar
Andreas Lauser committed
130
                volVars.fluidState().density(/*phaseIdx*/0) * volVars.fluidState().massFraction(/*phaseIdx*/0, comp1Idx) * volVars.porosity();
131
132
133
134
        }
        else
        {
            // storage term of continuity equation- molefractions
Andreas Lauser's avatar
Andreas Lauser committed
135
            //careful: molarDensity changes with moleFrac!
136
            result[contiEqIdx] += volVars.molarDensity()*volVars.porosity();
137
138
            // storage term of the transport equation - molefractions
            result[transEqIdx] +=
Andreas Lauser's avatar
Andreas Lauser committed
139
                volVars.fluidState().molarDensity(/*phaseIdx*/0)*volVars.fluidState().moleFraction(/*phaseIdx*/0, comp1Idx) *
140
141
                volVars.porosity();
        }
Katherina Baber's avatar
Katherina Baber committed
142

143
144
145
146
147
    }

    /*!
     * \brief Evaluates the mass flux over a face of a subcontrol
     *        volume.
Karin Erbertseder's avatar
Karin Erbertseder committed
148
     *
Karin Erbertseder's avatar
Karin Erbertseder committed
149
     *        \param flux The flux over the SCV (sub-control-volume) face for each component
Karin Erbertseder's avatar
Karin Erbertseder committed
150
     *        \param faceId The index of the considered face of the sub control volume
151
     */
152
    void computeFlux(PrimaryVariables &flux, int faceId, bool onBoundary=false) const
153
154
    {
        flux = 0;
Andreas Lauser's avatar
Andreas Lauser committed
155
156
157
158
        FluxVariables fluxVars(this->problem_(),
                               this->elem_(),
                               this->fvElemGeom_(),
                               faceId,
159
160
                               this->curVolVars_(),
                               onBoundary);
161

162
163
164
165
166
        asImp_()->computeAdvectiveFlux(flux, fluxVars);
        asImp_()->computeDiffusiveFlux(flux, fluxVars);
    }

    /*!
167
168
169
170
     * \brief Evaluates the advective mass flux of all components over
     *        a face of a subcontrol volume.
     *
     * \param flux The advective flux over the sub-control-volume face for each component
171
     * \param fluxVars The flux variables at the current SCV
172
     */
173
174
175
176
177
    void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
    {
        ////////
        // advective fluxes of all components in all phases
        ////////
Karin Erbertseder's avatar
Karin Erbertseder committed
178

Katherina Baber's avatar
Katherina Baber committed
179
        // data attached to upstream and the downstream vertices
180
181
        // of the current phase
        const VolumeVariables &up =
182
            this->curVolVars_(fluxVars.upstreamIdx());
183
        const VolumeVariables &dn =
184
            this->curVolVars_(fluxVars.downstreamIdx());
185

186
        if(!useMoles)
187
188
189
190
        {
            // total mass flux - massfraction
            //KmvpNormal is the Darcy velocity multiplied with the normal vector, calculated in 1p2cfluxvariables.hh
            flux[contiEqIdx] +=
191
192
193
194
                fluxVars.KmvpNormal() *
                ((     upwindWeight_)*up.density()/up.viscosity()
                 +
                 ((1 - upwindWeight_)*dn.density()/dn.viscosity()));
195
196
197

            // advective flux of the second component - massfraction
            flux[transEqIdx] +=
198
                fluxVars.KmvpNormal() *
Andreas Lauser's avatar
Andreas Lauser committed
199
                ((    upwindWeight_)*up.fluidState().density(/*phaseIdx=*/0) * up.fluidState().massFraction(/*phaseIdx=*/0, comp1Idx)/up.viscosity()
200
                 +
Andreas Lauser's avatar
Andreas Lauser committed
201
                 (1 - upwindWeight_)*dn.fluidState().density(/*phaseIdx=*/0)*dn.fluidState().massFraction(/*phaseIdx=*/0, comp1Idx)/dn.viscosity());
202
203
204
205
206
207
        }
        else
        {
            // total mass flux - molefraction
            //KmvpNormal is the Darcy velocity multiplied with the normal vector, calculated in 1p2cfluxvariables.hh
            flux[contiEqIdx] +=
208
209
210
211
                fluxVars.KmvpNormal() *
                ((     upwindWeight_)*up.molarDensity()/up.viscosity()
                 +
                 ((1 - upwindWeight_)*dn.molarDensity()/dn.viscosity()));
212
213
214

            // advective flux of the second component -molefraction
            flux[transEqIdx] +=
215
                fluxVars.KmvpNormal() *
Andreas Lauser's avatar
Andreas Lauser committed
216
217
218
                ((    upwindWeight_)*up.molarDensity() * up.fluidState().moleFraction(/*phaseIdx=*/0, comp1Idx)/up.viscosity()
                 +
                 (1 - upwindWeight_)*dn.molarDensity() * dn.fluidState().moleFraction(/*phaseIdx=*/0, comp1Idx)/dn.viscosity());
219
        }
220

221
222
    }

223
    /*!
224
225
226
227
     * \brief Adds the diffusive mass flux of all components over
     *        a face of a subcontrol volume.
     *
     * \param flux The diffusive flux over the sub-control-volume face for each component
228
     * \param fluxVars The flux variables at the current SCV
229
     */
230
231
    void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
    {
Katherina Baber's avatar
Katherina Baber committed
232
        Scalar tmp(0);
233

234
        // diffusive flux of second component
235
        if(!useMoles)
236
237
        {
            // diffusive flux of the second component - massfraction
238
            tmp = -(fluxVars.massFracGrad(comp1Idx)*fluxVars.face().normal);
239
            tmp *= fluxVars.porousDiffCoeff() * fluxVars.densityAtIP();
240
241
242
243
244
245

            flux[transEqIdx] += tmp;// * FluidSystem::molarMass(comp1Idx);
        }
        else
        {
            // diffusive flux of the second component - molefraction
246
            tmp = -(fluxVars.moleFracGrad(comp1Idx)*fluxVars.face().normal);
247
            tmp *= fluxVars.porousDiffCoeff() * fluxVars.molarDensityAtIP();
248

249
250
251
252
253
            // dispersive flux of second component - molefraction
            //            Vector normalDisp;
            //            fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
            //            tmp -= fluxVars.molarDensityAtIP()*
            //                (normalDisp * fluxVars.moleFracGrad(comp1Idx));
254

255
            flux[transEqIdx] += tmp;
256
        }
257
    }
258

259
260
    /*!
     * \brief Calculate the source term of the equation
Karin Erbertseder's avatar
Karin Erbertseder committed
261
     *        \param q The source/sink in the SCV for each component
Karin Erbertseder's avatar
Karin Erbertseder committed
262
263
     *        \param localVertexIdx The index of the vertex of the sub control volume
     *
264
     */
265
    void computeSource(PrimaryVariables &q, int localVertexIdx)
266
    {
267
268
269
270
271
        this->problem_().boxSDSource(q,
                                     this->elem_(),
                                     this->fvElemGeom_(),
                                     localVertexIdx,
                                     this->curVolVars_());
272
    }
Andreas Lauser's avatar
Andreas Lauser committed
273

274
    /*!
275
276
     * \brief Add Outflow boundary conditions for a single sub-control
     *        volume face to the local residual.
277
     */
278
    void evalOutflowSegment(const IntersectionIterator &isIt,
279
280
281
282
                            int scvIdx,
                            int boundaryFaceIdx)
    {
        const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
283
        // deal with outflow boundaries
284
        if (bcTypes.hasOutflow())
285
        {
286
287
            //calculate outflow fluxes
            PrimaryVariables values(0.0);
288
            asImp_()->computeFlux(values, boundaryFaceIdx, true);
289
            Valgrind::CheckDefined(values);
290

291
            for (int equationIdx = 0; equationIdx < numEq; ++equationIdx)
292
            {
293
294
295
296
                if (!bcTypes.isOutflow(equationIdx) )
                    continue;
                // deduce outflow
                this->residual_[scvIdx][equationIdx] += values[equationIdx];
297
            }
298
299
        }
    }
300

301
302
303
304
    Implementation *asImp_()
    { return static_cast<Implementation *> (this); }
    const Implementation *asImp_() const
    { return static_cast<const Implementation *> (this); }
305
306

private:
307
    Scalar upwindWeight_;
308
309
310
311
312
};

}

#endif