fickslaw.hh 15.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
 * \ingroup CCMpfaDiscretization
 * \brief Fick's law for cell-centered finite volume schemes with multi-point flux approximation
23
24
25
26
27
 */
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FICKS_LAW_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_FICKS_LAW_HH

#include <dumux/common/math.hh>
28
#include <dumux/common/properties.hh>
29
30
#include <dumux/discretization/methods.hh>

Dennis Gläser's avatar
Dennis Gläser committed
31
32
33
namespace Dumux
{
//! forward declaration of the method-specific implemetation
Timo Koch's avatar
Timo Koch committed
34
35
template<class TypeTag, DiscretizationMethods discMethod>
class FicksLawImplementation;
36

37
/*!
Dennis Gläser's avatar
Dennis Gläser committed
38
39
 * \ingroup CCMpfaDiscretization
 * \brief Fick's law for cell-centered finite volume schemes with multi-point flux approximation
40
41
42
43
44
45
46
 */
template <class TypeTag>
class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    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;
Dennis Gläser's avatar
Dennis Gläser committed
48

49
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
50
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
51
52
53
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
    using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
54
    using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
55
    using BalanceEqOpts = typename GET_PROP_TYPE(TypeTag, BalanceEqOpts);
56

57
58
    static constexpr int numComponents = GET_PROP_VALUE(TypeTag,NumComponents);
    using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>;
59

60
    //! Class that fills the cache corresponding to mpfa Fick's Law
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
    class MpfaFicksLawCacheFiller
    {
    public:
        //! Function to fill an MpfaFicksLawCache of a given scvf
        //! This interface has to be met by any diffusion-related cache filler class
        template<class FluxVariablesCacheFiller>
        static void fill(FluxVariablesCache& scvfFluxVarsCache,
                         unsigned int phaseIdx, unsigned int compIdx,
                         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
77
78
79
80
81
82
83
84
85
86
            // get interaction volume related data from the filler class & upate the cache
            if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
                scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.secondaryInteractionVolume(),
                                                  fluxVarsCacheFiller.secondaryIvLocalFaceData(),
                                                  fluxVarsCacheFiller.secondaryIvDataHandle(),
                                                  scvf, phaseIdx, compIdx);
            else
                scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.primaryInteractionVolume(),
                                                  fluxVarsCacheFiller.primaryIvLocalFaceData(),
                                                  fluxVarsCacheFiller.primaryIvDataHandle(),
                                                  scvf, phaseIdx, compIdx);
87
88
89
        }
    };

90
91
92
    //! The cache used in conjunction with the mpfa Fick's Law
    class MpfaFicksLawCache
    {
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
        // 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!" );

        static constexpr int dim = GridView::dimension;
        static constexpr int dimWorld = GridView::dimensionworld;
        static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
127
128

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

Dennis Gläser's avatar
Dennis Gläser committed
132
133
134
135
136
137
138
139
140
141
142
143
144
        /*!
         * \brief Update cached objects (transmissibilities)
         *
         * \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>
145
        void updateDiffusion(const InteractionVolume& iv,
Dennis Gläser's avatar
Dennis Gläser committed
146
                             const LocalFaceData& localFaceData,
147
148
                             const DataHandle& dataHandle,
                             const SubControlVolumeFace &scvf,
149
150
                             unsigned int phaseIdx, unsigned int compIdx)
        {
Dennis Gläser's avatar
Dennis Gläser committed
151
            switchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutside();
152

153
            // store pointer to the mole fraction vector of this iv
Dennis Gläser's avatar
Dennis Gläser committed
154
            xj_[phaseIdx][compIdx] = &dataHandle.moleFractions(phaseIdx, compIdx);
155

Dennis Gläser's avatar
Dennis Gläser committed
156
            const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
157
            if (dim == dimWorld)
Dennis Gläser's avatar
Dennis Gläser committed
158
                Tij_[phaseIdx][compIdx] = &dataHandle.diffusionT()[ivLocalIdx];
159
            else
Dennis Gläser's avatar
Dennis Gläser committed
160
161
                Tij_[phaseIdx][compIdx] = localFaceData.isOutside() ? &dataHandle.diffusionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()]
                                                                    : &dataHandle.diffusionT()[ivLocalIdx];
162
163
        }

Dennis Gläser's avatar
Dennis Gläser committed
164
165
166
167
        //! 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.
168
        bool diffusionSwitchFluxSign(unsigned int phaseIdx, unsigned int compIdx) const
Dennis Gläser's avatar
Dennis Gläser committed
169
        { return switchFluxSign_[phaseIdx][compIdx]; }
170

Dennis Gläser's avatar
Dennis Gläser committed
171
172
173
        //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions
        const Vector& diffusionTij(unsigned int phaseIdx, unsigned int compIdx) const
        { return *Tij_[phaseIdx][compIdx]; }
174

Dennis Gläser's avatar
Dennis Gläser committed
175
176
177
        //! The cell (& Dirichlet) mole fractions within this interaction volume
        const Vector& moleFractions(unsigned int phaseIdx, unsigned int compIdx) const
        { return *xj_[phaseIdx][compIdx]; }
178
179

    private:
Dennis Gläser's avatar
Dennis Gläser committed
180
181
182
        std::array< std::array<bool, numComponents>, numPhases > switchFluxSign_;
        std::array< std::array<const Vector*, numComponents>, numPhases > Tij_; //!< The transmissibilities such that f = Tij*xj
        std::array< std::array<const Vector*, numComponents>, numPhases > xj_;  //!< The interaction-volume wide mole fractions xj
183
184
    };

185
public:
186
187
188
    // state the discretization method this implementation belongs to
    static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;

189
190
191
    // state the type for the corresponding cache and its filler
    using Cache = MpfaFicksLawCache;

Dennis Gläser's avatar
Dennis Gläser committed
192
193
194
195
196
197
198
199
    //! Compute the diffusive flux across an scvf
    static ComponentFluxVector flux(const Problem& problem,
                                    const Element& element,
                                    const FVElementGeometry& fvGeometry,
                                    const ElementVolumeVariables&  elemVolVars,
                                    const SubControlVolumeFace& scvf,
                                    const int phaseIdx,
                                    const ElementFluxVariablesCache& elemFluxVarsCache)
200
    {
201
        ComponentFluxVector componentFlux(0.0);
202
203
        for (int compIdx = 0; compIdx < numComponents; compIdx++)
        {
204
            if(compIdx == FluidSystem::getMainComponent(phaseIdx))
205
                continue;
206
207

            // get the scaling factor for the effective diffusive fluxes
208
            const auto effFactor = computeEffectivityFactor(elemVolVars, scvf, phaseIdx);
209
210
211

            // if factor is zero, the flux will end up zero anyway
            if (effFactor == 0.0)
212
                continue;
213
214

            // calculate the density at the interface
215
            const auto rho = interpolateDensity(elemVolVars, scvf, phaseIdx);
216

217
218
219
            // prepare computations
            const auto& fluxVarsCache = elemFluxVarsCache[scvf];
            const auto& tij = fluxVarsCache.diffusionTij(phaseIdx, compIdx);
Dennis Gläser's avatar
Dennis Gläser committed
220
            const auto& xj = fluxVarsCache.moleFractions(phaseIdx, compIdx);
221

222
            // calculate Tij*xj
Dennis Gläser's avatar
Dennis Gläser committed
223
224
225
            Scalar flux = tij*xj;
            if (fluxVarsCache.diffusionSwitchFluxSign(phaseIdx, compIdx))
                flux *= -1.0;
226

Dennis Gläser's avatar
Dennis Gläser committed
227
            componentFlux[compIdx] = flux*rho*effFactor;
228
        }
229
230
231

        // accumulate the phase component flux
        for(int compIdx = 0; compIdx < numComponents; compIdx++)
232
            if(compIdx != FluidSystem::getMainComponent(phaseIdx) && BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem())
233
                componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
234

235
        return componentFlux;
236
237
    }

238
private:
Dennis Gläser's avatar
Dennis Gläser committed
239
    //! compute the density at branching facets for network grids as arithmetic mean
240
    static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars,
241
                                     const SubControlVolumeFace& scvf,
242
                                     const unsigned int phaseIdx)
243
244
245
246
    {
        // use arithmetic mean of the densities around the scvf
        if (!scvf.boundary())
        {
247
248
249
            Scalar rho = elemVolVars[scvf.insideScvIdx()].molarDensity(phaseIdx);
            for (const auto outsideIdx : scvf.outsideScvIndices())
                rho += elemVolVars[outsideIdx].molarDensity(phaseIdx);
250
251
252
            return rho/(scvf.outsideScvIndices().size()+1);
        }
        else
253
            return elemVolVars[scvf.outsideScvIdx()].molarDensity(phaseIdx);
254
255
256
257
258
259
    }

    //! Here we want to calculate the factors with which the diffusion coefficient has to be
    //! scaled to get the effective diffusivity. For this we use the effective diffusivity with
    //! a diffusion coefficient of 1.0 as input. Then we scale the transmissibilites during flux
    //! calculation (above) with the harmonic average of the two factors
260
    static Scalar computeEffectivityFactor(const ElementVolumeVariables& elemVolVars,
261
                                           const SubControlVolumeFace& scvf,
262
                                           const unsigned int phaseIdx)
263
    {
264
265
        using EffDiffModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);

266
        // use the harmonic mean between inside and outside
267
        const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
268
        const auto factor = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(),
269
270
                                                               insideVolVars.saturation(phaseIdx),
                                                               /*Diffusion coefficient*/ 1.0);
271
272
273

        if (!scvf.boundary())
        {
274
275
            // interpret outside factor as arithmetic mean
            Scalar outsideFactor = 0.0;
276
            for (const auto outsideIdx : scvf.outsideScvIndices())
277
278
279
            {
                const auto& outsideVolVars = elemVolVars[outsideIdx];
                outsideFactor += EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(),
280
281
                                                                    outsideVolVars.saturation(phaseIdx),
                                                                    /*Diffusion coefficient*/ 1.0);
282
283
            }
            outsideFactor /= scvf.outsideScvIndices().size();
284
285

            // use the harmonic mean of the two
286
            return harmonicMean(factor, outsideFactor);
287
288
289
290
291
292
293
294
295
        }

        return factor;
    }
};

} // end namespace

#endif