darcyslaw.hh 10.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   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          *
 *   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/>.   *
 *****************************************************************************/
/*!
 * \file
21
 * \brief This file contains the class which is required to calculate
22
 *        volume and mass fluxes of fluid phases over a face of a finite volume by means
23
24
 *        of the Darcy approximation. This specializations is for cell-centered schemes
 *        using multi-point flux approximation.
25
26
27
28
 */
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH

29
#include <dumux/common/properties.hh>
30

31
namespace Dumux
32
33
34
{

/*!
35
 * \ingroup Mpfa
Dennis Gläser's avatar
Dennis Gläser committed
36
 * \brief Specialization of Darcy's Law for the CCMpfa method.
37
38
39
40
 */
template <class TypeTag>
class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
41
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
42
43
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
44
    using Element = typename GridView::template Codim<0>::Entity;
45
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
46
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
47
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
48
49
    using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
    using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
50
51

    // Always use the dynamic type for vectors (compatibility with the boundary)
Dennis Gläser's avatar
Dennis Gläser committed
52
53
54
    using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
    using CoefficientVector = typename PrimaryInteractionVolume::Traits::DynamicVector;
    using DataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
55

56
    static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
57
58
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
59

60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
    //! Class that fills the cache corresponding to mpfa Darcy's Law
    class MpfaDarcysLawCacheFiller
    {
    public:
        //! Function to fill an MpfaDarcysLawCache of a given scvf
        //! This interface has to be met by any advection-related cache filler class
        template<class FluxVariablesCacheFiller>
        static void fill(FluxVariablesCache& scvfFluxVarsCache,
                         const Problem& problem,
                         const Element& element,
                         const FVElementGeometry& fvGeometry,
                         const ElementVolumeVariables& elemVolVars,
                         const SubControlVolumeFace& scvf,
                         const FluxVariablesCacheFiller& fluxVarsCacheFiller)
        {
            // get interaction volume from the flux vars cache filler & upate the cache
            if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
                scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.secondaryInteractionVolume(),
                                                  fluxVarsCacheFiller.dataHandle(),
                                                  scvf);
            else
                scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.primaryInteractionVolume(),
                                                  fluxVarsCacheFiller.dataHandle(),
                                                  scvf);
        }
    };

87
88
89
90
    //! The cache used in conjunction with the mpfa Darcy's Law
    class MpfaDarcysLawCache
    {
        // We always use the dynamic types here to be compatible on the boundary
Dennis Gläser's avatar
Dennis Gläser committed
91
92
        using Stencil = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer;
        using DirichletDataContainer = typename PrimaryInteractionVolume::DirichletDataContainer;
93
94

    public:
95
96
97
        // export the filler type
        using Filler = MpfaDarcysLawCacheFiller;

98
99
        //! update cached objects
        template<class InteractionVolume>
Dennis Gläser's avatar
Dennis Gläser committed
100
        void updateAdvection(const InteractionVolume& iv, const DataHandle& dataHandle, const SubControlVolumeFace &scvf)
101
102
103
        {
            const auto& localFaceData = iv.getLocalFaceData(scvf);

Dennis Gläser's avatar
Dennis Gläser committed
104
            // update the quantities that are equal for all phases
105
            advectionSwitchFluxSign_ = localFaceData.isOutside();
Dennis Gläser's avatar
Dennis Gläser committed
106
107
            advectionVolVarsStencil_ = &dataHandle.volVarsStencil();
            advectionDirichletData_ = &dataHandle.dirichletData();
108
109
110
111
112
113
114
115

            // the transmissibilities on surface grids have to be obtained from the outside
            if (dim == dimWorld)
                advectionTij_ = &dataHandle.T()[localFaceData.ivLocalScvfIndex()];
            else
                advectionTij_ = localFaceData.isOutside() ?
                                &dataHandle.outsideTij()[localFaceData.ivLocalOutsideScvfIndex()] :
                                &dataHandle.T()[localFaceData.ivLocalScvfIndex()];
116
117
        }

Dennis Gläser's avatar
Dennis Gläser committed
118
        //! Returns the stencil for advective scvf flux computation
119
        const Stencil& advectionVolVarsStencil() const { return *advectionVolVarsStencil_; }
120
121
122

        //! Returns the transmissibilities associated with the volume variables
        //! All phases flow through the same rock, thus, tij are equal for all phases
123
        const CoefficientVector& advectionTij() const { return *advectionTij_; }
124

Dennis Gläser's avatar
Dennis Gläser committed
125
126
        //! On faces that are "outside" w.r.t. a face in the interaction volume,
        //! we have to take the negative value of the fluxes, i.e. multiply by -1.0
127
        bool advectionSwitchFluxSign() const { return advectionSwitchFluxSign_; }
Dennis Gläser's avatar
Dennis Gläser committed
128
129
130

        //! Returns the data on dirichlet boundary conditions affecting
        //! the flux computation on this face
131
        const DirichletDataContainer& advectionDirichletData() const { return *advectionDirichletData_; }
132
133

    private:
Dennis Gläser's avatar
Dennis Gläser committed
134
135
136
137
        bool advectionSwitchFluxSign_;
        const Stencil* advectionVolVarsStencil_;
        const CoefficientVector* advectionTij_;
        const DirichletDataContainer* advectionDirichletData_;
138
139
    };

140
public:
141
142
143
    // state the discretization method this implementation belongs to
    static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;

144
    // export the type for the corresponding cache
145
146
    using Cache = MpfaDarcysLawCache;

147
148
149
150
151
152
    static Scalar flux(const Problem& problem,
                       const Element& element,
                       const FVElementGeometry& fvGeometry,
                       const ElementVolumeVariables& elemVolVars,
                       const SubControlVolumeFace& scvf,
                       const unsigned int phaseIdx,
153
                       const ElementFluxVariablesCache& elemFluxVarsCache)
154
    {
155
        static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
156

157
        // Calculate the interface density for gravity evaluation
158
        const auto rho = interpolateDensity(elemVolVars, scvf, phaseIdx);
Dennis Gläser's avatar
Dennis Gläser committed
159

160
        // prepare computations
Dennis Gläser's avatar
Dennis Gläser committed
161
        unsigned int i = 0;
162
163
164
        Scalar scvfFlux(0.0);
        const auto& fluxVarsCache = elemFluxVarsCache[scvf];
        const auto& tij = fluxVarsCache.advectionTij();
Dennis Gläser's avatar
Dennis Gläser committed
165
166
167

        // add contributions from cell-centered unknowns
        for (const auto volVarIdx : fluxVarsCache.advectionVolVarsStencil())
168
169
170
171
        {
            const auto& volVars = elemVolVars[volVarIdx];
            Scalar h = volVars.pressure(phaseIdx);

172
            // if gravity is enabled, add gravitational acceleration
173
            if (gravity)
174
175
            {
                // gravitational acceleration in the center of the actual element
Dennis Gläser's avatar
Dennis Gläser committed
176
                const auto x = fvGeometry.scv(volVarIdx).center();
177
178
179
                const auto g = problem.gravityAtPos(x);
                h -= rho*(g*x);
            }
180

Dennis Gläser's avatar
Dennis Gläser committed
181
            scvfFlux += tij[i++]*h;
182
183
        }

Dennis Gläser's avatar
Dennis Gläser committed
184
185
186
        // add contributions from possible dirichlet boundary conditions
        for (const auto& d : fluxVarsCache.advectionDirichletData())
        {
187
            const auto& volVars = elemVolVars[d.volVarIndex()];
Dennis Gläser's avatar
Dennis Gläser committed
188
189
190
191
192
            Scalar h = volVars.pressure(phaseIdx);

            // maybe add gravitational acceleration
            if (gravity)
            {
193
                const auto x = d.ipGlobal();
Dennis Gläser's avatar
Dennis Gläser committed
194
195
196
                const auto g = problem.gravityAtPos(x);
                h -= rho*(g*x);
            }
197

Dennis Gläser's avatar
Dennis Gläser committed
198
199
            scvfFlux += tij[i++]*h;
        }
200

Dennis Gläser's avatar
Dennis Gläser committed
201
202
        // return the flux (maybe adjust the sign)
        return fluxVarsCache.advectionSwitchFluxSign() ? -scvfFlux : scvfFlux;
203
    }
204

Dennis Gläser's avatar
Dennis Gläser committed
205
private:
206
    static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars,
207
                                     const SubControlVolumeFace& scvf,
Dennis Gläser's avatar
Dennis Gläser committed
208
                                     const unsigned int phaseIdx)
209
    {
210
        static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
211
212
213
214
215
216
217
218
219

        if (!gravity)
            return Scalar(0.0);
        else
        {
            // use arithmetic mean of the densities around the scvf
            if (!scvf.boundary())
            {
                Scalar rho = elemVolVars[scvf.insideScvIdx()].density(phaseIdx);
220
                for (const auto outsideIdx : scvf.outsideScvIndices())
221
222
223
224
225
226
227
                    rho += elemVolVars[outsideIdx].density(phaseIdx);
                return rho/(scvf.outsideScvIndices().size()+1);
            }
            else
                return elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
        }
    }
228
229
230
231
232
};

} // end namespace

#endif