richardslocalresidual.hh 6.34 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// $Id: richardslocalresidual.hh 3738 2010-06-15 14:01:09Z lauser $
/*****************************************************************************
 *   Copyright (C) 2007 by Peter Bastian                                     *
 *   Institute of Parallel and Distributed System                            *
 *   Department Simulation of Large Systems                                  *
 *   University of Stuttgart, Germany                                        *
 *                                                                           *
 *   Copyright (C) 2009 by Onur Dogan                                        *
 *   Copyright (C) 2009 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.                       *
 *****************************************************************************/
Bernd Flemisch's avatar
Bernd Flemisch committed
22
23
24
25
26
/*!
 * \file
 *
 * \brief Element-wise calculation of the residual for the Richards box model.
 */
27
28
#ifndef DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
#define DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
29
30
31

#include <dumux/boxmodels/common/boxlocalresidual.hh>

32
#include "richardsvolumevariables.hh"
33

34
#include "richardsfluxvariables.hh"
35
36
37
38
39

namespace Dumux
{
/*!
 * \ingroup RichardsModel
40
 * \brief Element-wise calculation of the residual for the Richards box model.
41
42
43
44
45
 */
template<class TypeTag>
class RichardsLocalResidual : public BoxLocalResidual<TypeTag>
{
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
46
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
47
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
48
49
50
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
    
51
52
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(RichardsIndices)) Indices;
    enum {
53
        dimWorld = GridView::dimensionworld,
54

55
        contiEqIdx = Indices::contiEqIdx,
56

57
58
59
        wPhaseIdx = Indices::wPhaseIdx,
        nPhaseIdx = Indices::nPhaseIdx,
    };
60

61
    typedef Dune::FieldVector<Scalar, dimWorld> Vector;
62
    static const Scalar mobilityUpwindAlpha = GET_PROP_VALUE(TypeTag, PTAG(MobilityUpwindAlpha));
63
    
64
65
66
67
68
69
70
71
public:
    /*!
     * \brief Evaluate the rate of change of all conservation
     *        quantites (e.g. phase mass) within a sub control
     *        volume of a finite volume element for the Richards
     *        model.
     *
     * This function should not include the source and sink terms.
72
73
74
75
76
     *
     * \param result Stores the average mass per unit volume for each phase [kg/m^3]
     * \param scvIdx The sub control volume index of the current element
     * \param usePrevSol Calculate the storage term of the previous solution
     *                   instead of the model's current solution.
77
     */
78
    void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
79
80
81
82
83
84
    {
        // 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.
85
86
87
88
        const VolumeVariables &volVars = 
            usePrevSol ? 
            this->prevVolVars_(scvIdx) :
            this->curVolVars_(scvIdx);
89
90

        // partial time derivative of the wetting phase mass
91
92
93
94
        result[contiEqIdx] =
            volVars.density(wPhaseIdx)
            * volVars.saturation(wPhaseIdx)
            * volVars.porosity();
95
96
97
98
99
100
    }


    /*!
     * \brief Evaluates the mass flux over a face of a subcontrol
     *        volume.
101
102
103
104
105
106
     *
     * 
     * \param flux Stores the total mass fluxes over a sub-control volume face
     *             of the current element [kg/s]
     * \param scvfIdx The sub control volume face index inside the current 
     *                element
107
     */
108
    void computeFlux(PrimaryVariables &flux, int scvfIdx) const
109
    {
110
111
112
        FluxVariables fluxVars(this->problem_(),
                               this->elem_(),
                               this->fvElemGeom_(),
113
                               scvfIdx,
114
115
116
117
118
119
120
121
122
                               this->curVolVars_());

        // calculate the flux in the normal direction of the
        // current sub control volume face
        Vector tmpVec;
        fluxVars.intrinsicPermeability().mv(fluxVars.potentialGradW(),
                                            tmpVec);
        Scalar normalFlux = - (tmpVec*fluxVars.face().normal);
        
123
        // data attached to upstream and the downstream vertices
124
125
126
127
128
129
130
131
132
133
        // of the current phase
        const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(normalFlux));
        const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(normalFlux));
        
        flux[contiEqIdx] =
            normalFlux
            *
            ((    mobilityUpwindAlpha)*up.density(wPhaseIdx)*up.mobility(wPhaseIdx) 
             +
             (1 - mobilityUpwindAlpha)*dn.density(wPhaseIdx)*dn.mobility(wPhaseIdx));
134
135
136
137
    }

    /*!
     * \brief Calculate the source term of the equation
138
     *
139
140
     * \param q Stores the average source term of all phases inside a 
     *          sub-control volume of the current element [kg/(m^3 s)]
141
     * \param scvIdx The sub control volume index inside the current 
142
     *               element
143
     */
144
    void computeSource(PrimaryVariables &q, int scvIdx)
145
146
147
148
    {
        this->problem_().source(q,
                                this->elem_(),
                                this->fvElemGeom_(),
149
                                scvIdx);
150
151
152
153
154
155
    }
};

};

#endif