darcyslaw.hh 10.7 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
{
Timo Koch's avatar
Timo Koch committed
33
34
35
// forward declarations
template<class TypeTag, DiscretizationMethods discMethod>
class DarcysLawImplementation;
36
37

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

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

59
    static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
60
61
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
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
87
88
89
    //! 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);
        }
    };

90
91
92
93
    //! 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
94
95
        using Stencil = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer;
        using DirichletDataContainer = typename PrimaryInteractionVolume::DirichletDataContainer;
96
97

    public:
98
99
100
        // export the filler type
        using Filler = MpfaDarcysLawCacheFiller;

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

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

            // 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()];
119
120
        }

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

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

Dennis Gläser's avatar
Dennis Gläser committed
128
129
        //! 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
130
        bool advectionSwitchFluxSign() const { return advectionSwitchFluxSign_; }
Dennis Gläser's avatar
Dennis Gläser committed
131
132
133

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

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

143
public:
144
145
146
    // state the discretization method this implementation belongs to
    static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;

147
    // export the type for the corresponding cache
148
149
    using Cache = MpfaDarcysLawCache;

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

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

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

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

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

Dennis Gläser's avatar
Dennis Gläser committed
184
            scvfFlux += tij[i++]*h;
185
186
        }

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

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

Dennis Gläser's avatar
Dennis Gläser committed
201
202
            scvfFlux += tij[i++]*h;
        }
203

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

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

        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);
223
                for (const auto outsideIdx : scvf.outsideScvIndices())
224
225
226
227
228
229
230
                    rho += elemVolVars[outsideIdx].density(phaseIdx);
                return rho/(scvf.outsideScvIndices().size()+1);
            }
            else
                return elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
        }
    }
231
232
233
234
235
};

} // end namespace

#endif