darcyslaw.hh 12.4 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
Dennis Gläser's avatar
Dennis Gläser committed
21
22
23
 * \ingroup CCMpfaDiscretization
 * \brief Darcy's Law for cell-centered finite volume schemes
 *        with multi-point flux approximation.
24
25
26
27
 */
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH

Dennis Gläser's avatar
Dennis Gläser committed
28
29
30
#include <dune/common/densevector.hh>
#include <dune/common/densematrix.hh>

31
#include <dumux/common/properties.hh>
32
33
#include <dumux/common/parameters.hh>
#include <dumux/discretization/methods.hh>
34

Dennis Gläser's avatar
Dennis Gläser committed
35
36
#include <dumux/discretization/methods.hh>

37
namespace Dumux
38
{
Dennis Gläser's avatar
Dennis Gläser committed
39
//! forward declaration of the method-specific implementation
Timo Koch's avatar
Timo Koch committed
40
41
template<class TypeTag, DiscretizationMethods discMethod>
class DarcysLawImplementation;
42
43

/*!
Dennis Gläser's avatar
Dennis Gläser committed
44
45
 * \ingroup CCMpfaDiscretization
 * \brief Darcy's law for cell-centered finite volume schemes with multi-point flux approximation.
46
47
48
49
 */
template <class TypeTag>
class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
50
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
51
52
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
53
    using Element = typename GridView::template Codim<0>::Entity;
Dennis Gläser's avatar
Dennis Gläser committed
54

55
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
56
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
57
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
58
59
    using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
    using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
60

61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
    //! 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)
        {
Dennis Gläser's avatar
Dennis Gläser committed
76
            // get interaction volume related data from the filler class & upate the cache
77
78
            if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
                scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.secondaryInteractionVolume(),
Dennis Gläser's avatar
Dennis Gläser committed
79
80
                                                  fluxVarsCacheFiller.secondaryIvLocalFaceData(),
                                                  fluxVarsCacheFiller.secondaryIvDataHandle(),
81
82
83
                                                  scvf);
            else
                scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.primaryInteractionVolume(),
Dennis Gläser's avatar
Dennis Gläser committed
84
85
                                                  fluxVarsCacheFiller.primaryIvLocalFaceData(),
                                                  fluxVarsCacheFiller.primaryIvDataHandle(),
86
87
88
89
                                                  scvf);
        }
    };

90
91
92
    //! The cache used in conjunction with the mpfa Darcy's Law
    class MpfaDarcysLawCache
    {
Dennis Gläser's avatar
Dennis Gläser committed
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
        static constexpr int dim = GridView::dimension;
        static constexpr int dimWorld = GridView::dimensionworld;
        static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);

        // In the current implementation of the flux variables cache we cannot
        // make a disctinction between dynamic (mpfa-o) and static (mpfa-l)
        // matrix and vector types, as currently the cache class can only be templated
        // by a type tag (and there can only be one). We use a dynamic vector here to
        // make sure it works in case one of the two used interaction volume types uses
        // dynamic types performance is thus lowered for schemes using static types.
        // TODO: this has to be checked thoroughly as soon as a scheme using static types
        //       is implemented. One idea to overcome the performance drop could be only
        //       storing the iv-local index here and obtain tij always from the datahandle
        //       of the fluxVarsCacheContainer
        using Vector = Dune::DynamicVector< Scalar >;
        using Matrix = Dune::DynamicMatrix< Scalar >;

        using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
        using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
        using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;

        static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value,
                       "The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" );
        static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value,
                       "The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" );

        using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
        using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
        using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;

        static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value,
                       "The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" );
        static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value,
                       "The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" );
127
128

    public:
129
130
131
        // export the filler type
        using Filler = MpfaDarcysLawCacheFiller;

Dennis Gläser's avatar
Dennis Gläser committed
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
        /*!
         * \brief Update cached objects (transmissibilities and gravity)
         *
         * \tparam InteractionVolume The (mpfa scheme-specific) interaction volume
         * \tparam LocalFaceData The object used to store iv-local info on an scvf
         * \tparam DataHandle The object used to store transmissibility matrices etc.
         *
         * \param iv The interaction volume this scvf is embedded in
         * \param localFaceData iv-local info on this scvf
         * \param dataHandle Transmissibility matrix & gravity data of this iv
         * \param scvf The sub-control volume face
         */
        template< class InteractionVolume, class LocalFaceData, class DataHandle >
        void updateAdvection(const InteractionVolume& iv,
                             const LocalFaceData& localFaceData,
                             const DataHandle& dataHandle,
                             const SubControlVolumeFace &scvf)
149
        {
Dennis Gläser's avatar
Dennis Gläser committed
150
151
152
153
            switchFluxSign_ = localFaceData.isOutside();

            for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
                pj_[pIdx] = &dataHandle.pressures(pIdx);
154

Dennis Gläser's avatar
Dennis Gläser committed
155
156
            static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
            const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
157

Dennis Gläser's avatar
Dennis Gläser committed
158
            // standard grids
159
            if (dim == dimWorld)
Dennis Gläser's avatar
Dennis Gläser committed
160
161
162
163
164
165
166
167
168
            {
                Tij_ = &dataHandle.advectionT()[ivLocalIdx];

                if (enableGravity)
                    for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
                        g_[phaseIdx] = dataHandle.gravity(phaseIdx)[ivLocalIdx];
            }

            // surface grids
169
            else
Dennis Gläser's avatar
Dennis Gläser committed
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
            {
                if (!localFaceData.isOutside())
                {
                    Tij_ = &dataHandle.advectionT()[ivLocalIdx];

                    if (enableGravity)
                        for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
                            g_[phaseIdx] = dataHandle.gravity(phaseIdx)[ivLocalIdx];
                }
                else
                {
                    const auto idxInOutsideFaces = localFaceData.scvfLocalOutsideScvfIndex();
                    Tij_ = &dataHandle.advectionTout()[ivLocalIdx][idxInOutsideFaces];

                    if (enableGravity)
                        for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
                            g_[phaseIdx] = dataHandle.gravityOutside(phaseIdx)[ivLocalIdx][idxInOutsideFaces];
                }
            }
189
190
        }

Dennis Gläser's avatar
Dennis Gläser committed
191
192
        //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions
        const Vector& advectionTij() const { return *Tij_; }
193

Dennis Gläser's avatar
Dennis Gläser committed
194
195
        //! The cell (& Dirichlet) pressures within this interaction volume
        const Vector& pressures(unsigned int phaseIdx) const { return *pj_[phaseIdx]; }
196

Dennis Gläser's avatar
Dennis Gläser committed
197
198
        //! The gravitational acceleration for a phase on this scvf
        Scalar gravity(unsigned int phaseIdx) const { return g_[phaseIdx]; }
Dennis Gläser's avatar
Dennis Gläser committed
199

Dennis Gläser's avatar
Dennis Gläser committed
200
201
202
203
204
        //! In the interaction volume-local system of eq we have one unknown per face.
        //! On scvfs on this face, but in "outside" (neighbor) elements of it, we have
        //! to take the negative value of the fluxes due to the flipped normal vector.
        //! This function returns whether or not this scvf is an "outside" face in the iv.
        bool advectionSwitchFluxSign() const { return switchFluxSign_; }
205
206

    private:
Dennis Gläser's avatar
Dennis Gläser committed
207
208
209
210
        bool switchFluxSign_;
        const Vector* Tij_;                       //!< The transmissibilities such that f = Tij*pj
        std::array< Scalar, numPhases > g_;       //!< Gravitational flux contribution on this face
        std::array<const Vector*, numPhases> pj_; //!< The interaction-volume wide phase pressures pj
211
212
    };

213
public:
214
215
216
    // state the discretization method this implementation belongs to
    static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;

217
    // export the type for the corresponding cache
218
219
    using Cache = MpfaDarcysLawCache;

Dennis Gläser's avatar
Dennis Gläser committed
220
    //! Compute the advective flux across an scvf
221
222
223
224
225
226
    static Scalar flux(const Problem& problem,
                       const Element& element,
                       const FVElementGeometry& fvGeometry,
                       const ElementVolumeVariables& elemVolVars,
                       const SubControlVolumeFace& scvf,
                       const unsigned int phaseIdx,
227
                       const ElementFluxVariablesCache& elemFluxVarsCache)
228
    {
229
230
        const auto& fluxVarsCache = elemFluxVarsCache[scvf];
        const auto& tij = fluxVarsCache.advectionTij();
Dennis Gläser's avatar
Dennis Gläser committed
231
        const auto& pj = fluxVarsCache.pressures(phaseIdx);
Dennis Gläser's avatar
Dennis Gläser committed
232

233
        // compute t_ij*p_j
Dennis Gläser's avatar
Dennis Gläser committed
234
        Scalar scvfFlux = tij*pj;
235

236
        // maybe add gravitational acceleration
Dennis Gläser's avatar
Dennis Gläser committed
237
238
239
        static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
        if (enableGravity)
            scvfFlux += fluxVarsCache.gravity(phaseIdx);
240

241
        // switch the sign if necessary
Dennis Gläser's avatar
Dennis Gläser committed
242
243
        if (fluxVarsCache.advectionSwitchFluxSign())
            scvfFlux *= -1.0;
244

Dennis Gläser's avatar
Dennis Gläser committed
245
        return scvfFlux;
246
    }
247
248
249
250
251
};

} // end namespace

#endif