fourierslaw.hh 9.51 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// -*- 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
 * \brief This file contains the data which is required to calculate
Dennis Gläser's avatar
Dennis Gläser committed
22
 *        heat conduction fluxes with Fourier's law for cell-centered MPFA models.
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
 */
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FOURIERS_LAW_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_FOURIERS_LAW_HH

#include <dune/common/float_cmp.hh>

#include <dumux/common/math.hh>
#include <dumux/common/parameters.hh>
#include <dumux/implicit/properties.hh>
#include <dumux/discretization/methods.hh>

namespace Dumux
{

/*!
 * \ingroup FouriersLaw
 * \brief Specialization of Fourier's Law for the CCMpfa method.
 */
template <class TypeTag>
class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
44
45
    using Implementation = typename GET_PROP_TYPE(TypeTag, HeatConductionType);

46
47
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
48
    using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
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
56
57
    using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);

    // Always use the dynamic type for vectors (compatibility with the boundary)
    using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
58
    using CoefficientVector = typename BoundaryInteractionVolume::Vector;
59
60
61
62
63

    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using Element = typename GridView::template Codim<0>::Entity;
    using IndexType = typename GridView::IndexSet::IndexType;

64
65
66
67
68
    static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
    static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);

    static constexpr int energyEqIdx = GET_PROP_TYPE(TypeTag, Indices)::energyEqIdx;

69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
    //! 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)
        {
            // get interaction volume from the flux vars cache filler & upate the cache
            if (problem.model().fvGridGeometry().isInBoundaryInteractionVolume(scvf))
                scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.boundaryInteractionVolume(), scvf);
            else
                scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.interactionVolume(), scvf);
        }
    };

92
93
94
95
96
    //! The cache used in conjunction with the mpfa Fourier's Law
    class MpfaFouriersLawCache
    {
        using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
    public:
97
98
99
        // export filler type
        using Filler = MpfaFouriersLawCacheFiller;

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
127
128
129
130
131
132
        // update cached objects for heat conduction
        template<typename InteractionVolume>
        void updateHeatConduction(const InteractionVolume& iv, const SubControlVolumeFace &scvf)
        {
            const auto& localFaceData = iv.getLocalFaceData(scvf);
            heatConductionVolVarsStencil_ = iv.volVarsStencil();
            heatConductionTij_ = iv.getTransmissibilities(localFaceData);
            heatNeumannFlux_ = iv.getNeumannFlux(localFaceData, energyEqIdx);
        }

        //! Returns the volume variables indices necessary for heat conduction flux
        //! computation. This includes all participating boundary volume variables
        //! and it can be different for the phases & components.
        const Stencil& heatConductionVolVarsStencil() const
        { return heatConductionVolVarsStencil_; }

        //! Returns the transmissibilities associated with the volume variables
        //! This can be different for the phases & components.
        const CoefficientVector& heatConductionTij() const
        { return heatConductionTij_; }

        //! If the useTpfaBoundary property is set to false, the boundary conditions
        //! are put into the local systems leading to possible contributions on all faces
        Scalar heatNeumannFlux() const
        { return heatNeumannFlux_; }

    private:
        // Quantities associated with heat conduction
        Stencil heatConductionVolVarsStencil_;
        CoefficientVector heatConductionTij_;
        Scalar heatNeumannFlux_;
    };

133
134
135
136
public:
    // state the discretization method this implementation belongs to
    static const DiscretizationMethods myDiscretizationMethod = DiscretizationMethods::CCMpfa;

137
138
139
    // state the type for the corresponding cache and its filler
    using Cache = MpfaFouriersLawCache;

140
141
142
143
144
145
146
147
148
149
150
    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& volVarsStencil = fluxVarsCache.heatConductionVolVarsStencil();
        const auto& tij = fluxVarsCache.heatConductionTij();

151
        const bool isInteriorBoundary = enableInteriorBoundaries && fluxVarsCache.isInteriorBoundary();
152
        // For interior Neumann boundaries when using Tpfa on boundaries, return the user-specified flux
153
154
155
156
157
        if (isInteriorBoundary
            && useTpfaBoundary
            && fluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorNeumann)
            return scvf.area()*
                   elemVolVars[scvf.insideScvIdx()].extrusionFactor()*
158
                   problem.neumann(element, fvGeometry, elemVolVars, scvf)[energyEqIdx];
159

160
161
162
163
164
165
        // calculate Tij*tj
        Scalar flux(0.0);
        unsigned int localIdx = 0;
        for (const auto volVarIdx : volVarsStencil)
            flux += tij[localIdx++]*elemVolVars[volVarIdx].temperature();

166
167
168
169
170
        // if no interior boundaries are present, return heat conduction flux
        if (!enableInteriorBoundaries)
            return useTpfaBoundary ? flux : flux + fluxVarsCache.heatNeumannFlux();

        // Handle interior boundaries
171
        flux += Implementation::computeInteriorBoundaryContribution(fvGeometry, elemVolVars, fluxVarsCache);
172
173
174
175

        // return overall resulting flux
        return useTpfaBoundary ? flux : flux + fluxVarsCache.heatNeumannFlux();
    }
176

177
    static Scalar computeInteriorBoundaryContribution(const FVElementGeometry& fvGeometry,
178
                                                      const ElementVolumeVariables& elemVolVars,
179
180
181
182
183
184
185
186
                                                      const FluxVariablesCache& fluxVarsCache)
    {
        // obtain the transmissibilites associated with all pressures
        const auto& tij = fluxVarsCache.heatConductionTij();

        // the interior dirichlet boundaries local indices start after
        // the cell and the domain Dirichlet boundary pressures
        const auto startIdx = fluxVarsCache.heatConductionVolVarsStencil().size();
187

188
189
        // add interior Dirichlet boundary contributions
        Scalar flux = 0.0;
190
191
        for (auto&& data : fluxVarsCache.interiorBoundaryData())
            if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
192
                flux += tij[startIdx + data.localIndexInInteractionVolume()]*data.facetVolVars(fvGeometry).temperature();
193

194
        return flux;
195
196
197
198
199
200
    }
};

} // end namespace Dumux

#endif