fourierslaw.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
Dennis Gläser's avatar
Dennis Gläser committed
21
22
 * \ingroup CCMpfaDiscretization
 * \brief Fourier's law for cell-centered finite volume schemes with multi-point flux approximation
23
24
25
26
 */
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FOURIERS_LAW_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_FOURIERS_LAW_HH

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

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

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

/*!
Dennis Gläser's avatar
Dennis Gläser committed
39
40
41
* \ingroup CCMpfaDiscretization
* \brief Fourier's law for cell-centered finite volume schemes with two-point flux approximation
*/
42
43
44
45
46
template <class TypeTag>
class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
47
48
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using Element = typename GridView::template Codim<0>::Entity;
49
50
51
52
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
    using ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
53
    using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
54
55
    using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);

56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
    //! Class that fills the cache corresponding to mpfa Darcy's Law
    class MpfaFouriersLawCacheFiller
    {
    public:
        //! Function to fill an MpfaDarcysLawCache of a given scvf
        //! This interface has to be met by any cache filler class for heat conduction quantities
        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)
        {
71
72
73
          // get interaction volume from the flux vars cache filler & upate the cache
          if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
              scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.secondaryInteractionVolume(),
Dennis Gläser's avatar
Dennis Gläser committed
74
75
                                                     fluxVarsCacheFiller.secondaryIvLocalFaceData(),
                                                     fluxVarsCacheFiller.secondaryIvDataHandle(),
76
77
78
                                                     scvf);
          else
              scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.primaryInteractionVolume(),
Dennis Gläser's avatar
Dennis Gläser committed
79
80
                                                     fluxVarsCacheFiller.primaryIvLocalFaceData(),
                                                     fluxVarsCacheFiller.primaryIvDataHandle(),
81
                                                     scvf);
82
83
84
        }
    };

85
86
87
    //! The cache used in conjunction with the mpfa Fourier's Law
    class MpfaFouriersLawCache
    {
Dennis Gläser's avatar
Dennis Gläser committed
88
89
90
91
92
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
        // 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;

122
    public:
123
124
125
        // export filler type
        using Filler = MpfaFouriersLawCacheFiller;

Dennis Gläser's avatar
Dennis Gläser committed
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
        /*!
         * \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>
        void updateHeatConduction(const InteractionVolume& iv,
                                  const LocalFaceData& localFaceData,
                                  const DataHandle& dataHandle,
                                  const SubControlVolumeFace &scvf)
143
        {
Dennis Gläser's avatar
Dennis Gläser committed
144
            switchFluxSign_ = localFaceData.isOutside();
145

146
            // store pointer to the temperature vector of this iv
Dennis Gläser's avatar
Dennis Gläser committed
147
            Tj_ = &dataHandle.temperatures();
148

Dennis Gläser's avatar
Dennis Gläser committed
149
            const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
150
            if (dim == dimWorld)
Dennis Gläser's avatar
Dennis Gläser committed
151
                Tij_ = &dataHandle.heatConductionT()[ivLocalIdx];
152
            else
Dennis Gläser's avatar
Dennis Gläser committed
153
154
                Tij_ = localFaceData.isOutside() ? &dataHandle.heatConductionTout()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()]
                                                 : &dataHandle.heatConductionT()[ivLocalIdx];
155
156
        }

Dennis Gläser's avatar
Dennis Gläser committed
157
158
        //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions
        const Vector& heatConductionTij() const { return *Tij_; }
159

Dennis Gläser's avatar
Dennis Gläser committed
160
161
162
163
164
        //! 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 heatConductionSwitchFluxSign() const { return switchFluxSign_; }
165

Dennis Gläser's avatar
Dennis Gläser committed
166
167
        //! The cell (& Dirichlet) temperatures within this interaction volume
        const Vector& temperatures() const { return *Tj_; }
168
169

    private:
Dennis Gläser's avatar
Dennis Gläser committed
170
171
172
        bool switchFluxSign_;
        const Vector* Tij_;   //!< The transmissibilities such that f = Tij*Tj
        const Vector* Tj_;    //!< The interaction-volume wide temperature Tj
173
174
    };

175
176
177
178
public:
    // state the discretization method this implementation belongs to
    static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;

179
180
181
    // state the type for the corresponding cache and its filler
    using Cache = MpfaFouriersLawCache;

Dennis Gläser's avatar
Dennis Gläser committed
182
    //! Compute the conductive flux across an scvf
183
184
185
186
187
188
189
190
191
    static Scalar flux(const Problem& problem,
                       const Element& element,
                       const FVElementGeometry& fvGeometry,
                       const ElementVolumeVariables& elemVolVars,
                       const SubControlVolumeFace& scvf,
                       const ElementFluxVarsCache& elemFluxVarsCache)
    {
        const auto& fluxVarsCache = elemFluxVarsCache[scvf];
        const auto& tij = fluxVarsCache.heatConductionTij();
Dennis Gläser's avatar
Dennis Gläser committed
192
        const auto& Tj = fluxVarsCache.temperatures();
193

194
        // compute Tij*tj
Dennis Gläser's avatar
Dennis Gläser committed
195
196
197
        Scalar flux = tij*Tj;
        if (fluxVarsCache.heatConductionSwitchFluxSign())
            flux *= -1.0;
198

Dennis Gläser's avatar
Dennis Gläser committed
199
        return flux;
200
201
202
203
204
205
    }
};

} // end namespace Dumux

#endif