interactionvolumedatahandle.hh 15 KB
Newer Older
Dennis Gläser's avatar
Dennis Gläser committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
// -*- 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 Data handle class for interaction volumes of mpfa methods.
 *        This class is passed to interaction volumes to store the necessary data in it.
 */
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEDATAHANDLE_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_INTERACTIONVOLUMEDATAHANDLE_HH

Dennis Gläser's avatar
Dennis Gläser committed
27
28
#include <vector>

Timo Koch's avatar
Timo Koch committed
29
30
#include <dumux/common/properties.hh>

Dennis Gläser's avatar
Dennis Gläser committed
31
32
33
namespace Dumux
{

Dennis Gläser's avatar
Dennis Gläser committed
34
35
36
37
38
39
40
//! Empty data handle class
template<class InteractionVolume>
class EmptyDataHandle
{
public:
    void resize(const InteractionVolume& iv) {}
};
41

Dennis Gläser's avatar
Dennis Gläser committed
42
43
44
45
//! Data handle for quantities related to advection
template<class TypeTag, class InteractionVolume, bool EnableAdvection>
class AdvectionDataHandle
{
46
    // export matrix & vector types from interaction volume
Dennis Gläser's avatar
Dennis Gläser committed
47
48
49
50
    using Matrix = typename InteractionVolume::Traits::Matrix;
    using Vector = typename InteractionVolume::Traits::Vector;
    using Scalar = typename Vector::value_type;
    using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType;
51

Dennis Gläser's avatar
Dennis Gläser committed
52
    using OutsideDataContainer = std::vector< std::vector<Vector> >;
53

Dennis Gläser's avatar
Dennis Gläser committed
54
55
56
    static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
    static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
    static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
57

Dennis Gläser's avatar
Dennis Gläser committed
58
public:
59

Dennis Gläser's avatar
Dennis Gläser committed
60
61
62
    //! prepare data handle for subsequent fill (normal grids)
    template< int d = dim, int dw = dimWorld, std::enable_if_t<(d==dw), int> = 0>
    void resize(const InteractionVolume& iv)
Dennis Gläser's avatar
Dennis Gläser committed
63
    {
Dennis Gläser's avatar
Dennis Gläser committed
64
65
66
67
68
69
70
71
        // resize transmissibility matrix & pressure vectors
        T_.resize(iv.numFaces(), iv.numKnowns());
        for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
            p_[pIdx].resize(iv.numKnowns());

        // maybe resize gravity container
        static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
        if (enableGravity)
Dennis Gläser's avatar
Dennis Gläser committed
72
        {
Dennis Gläser's avatar
Dennis Gläser committed
73
74
75
            CA_.resize(iv.numFaces(), iv.numUnknowns());
            for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
                g_[pIdx].resize(iv.numFaces());
Dennis Gläser's avatar
Dennis Gläser committed
76
        }
Dennis Gläser's avatar
Dennis Gläser committed
77
    }
Dennis Gläser's avatar
Dennis Gläser committed
78
79


Dennis Gläser's avatar
Dennis Gläser committed
80
81
82
83
84
    //! prepare data handle for subsequent fill (surface grids)
    template< int d = dim, int dw = dimWorld, std::enable_if_t<(d<dw), int> = 0>
    void resize(const InteractionVolume& iv)
    {
        static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
Dennis Gläser's avatar
Dennis Gläser committed
85

Dennis Gläser's avatar
Dennis Gläser committed
86
        if (!enableGravity)
Dennis Gläser's avatar
Dennis Gläser committed
87
        {
Dennis Gläser's avatar
Dennis Gläser committed
88
89
90
91
92
93
94
95
96
97
98
99
            T_.resize(iv.numFaces(), iv.numKnowns());
            outsideT_.resize(iv.numFaces());
            for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
                p_[pIdx].resize(iv.numKnowns());

            for (LocalIndexType i = 0; i < iv.numFaces(); ++i)
            {
                const auto numNeighbors = iv.localScvf(i).neighboringLocalScvIndices().size() - 1;
                outsideT_[i].resize(numNeighbors);
                for (LocalIndexType j = 0; j < numNeighbors; ++j)
                    outsideT_[i][j].resize(iv.numKnowns());
            }
Dennis Gläser's avatar
Dennis Gläser committed
100
101
        }

Dennis Gläser's avatar
Dennis Gläser committed
102
        else
Dennis Gläser's avatar
Dennis Gläser committed
103
        {
Dennis Gläser's avatar
Dennis Gläser committed
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
            T_.resize(iv.numFaces(), iv.numKnowns());
            CA_.resize(iv.numFaces(), iv.numUnknowns());
            A_.resize(iv.numUnknowns(), iv.numUnknowns());
            outsideT_.resize(iv.numFaces());

            for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
            {
                p_[pIdx].resize(iv.numKnowns());
                g_[pIdx].resize(iv.numFaces());
                outsideG_[pIdx].resize(iv.numFaces());
            }

            for (LocalIndexType i = 0; i < iv.numFaces(); ++i)
            {
                const auto numNeighbors = iv.localScvf(i).neighboringLocalScvIndices().size() - 1;
                outsideT_[i].resize(numNeighbors);

                for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
                    outsideG_[pIdx][i].resize(numNeighbors);

                for (LocalIndexType j = 0; j < numNeighbors; ++j)
                    outsideT_[i][j].resize(iv.numKnowns());
            }
Dennis Gläser's avatar
Dennis Gläser committed
127
128

        }
Dennis Gläser's avatar
Dennis Gläser committed
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
    }

    //! Access to the iv-wide pressure of one phase
    const Vector& pressures(unsigned int pIdx) const { return p_[pIdx]; }
    Vector& pressures(unsigned int pIdx) { return p_[pIdx]; }

    //! Access to the gravitational flux contributions for one phase
    const Vector& gravity(unsigned int pIdx) const { return g_[pIdx]; }
    Vector& gravity(unsigned int pIdx) { return g_[pIdx]; }

    //! Access to the gravitational flux contributions for all phases
    const std::array< Vector, numPhases >& gravity() const { return g_; }
    std::array< Vector, numPhases >& gravity() { return g_; }

    //! Projection matrix for gravitational acceleration
    const Matrix& advectionCA() const { return CA_; }
    Matrix& advectionCA() { return CA_; }

    //! Additional projection matrix needed on surface grids
    const Matrix& advectionA() const { return A_; }
    Matrix& advectionA() { return A_; }

    //! The transmissibilities associated with advective fluxes
    const Matrix& advectionT() const { return T_; }
    Matrix& advectionT() { return T_; }

    //! The transmissibilities for "outside" faces (used on surface grids)
    const std::vector< std::vector<Vector> >& advectionTout() const { return outsideT_; }
    std::vector< std::vector<Vector> >& advectionTout() { return outsideT_; }

    //! The gravitational acceleration for "outside" faces (used on surface grids)
    const std::array< std::vector<Vector>, numPhases >& gravityOutside() const { return outsideG_; }
    std::array< std::vector<Vector>, numPhases >& gravityOutside() { return outsideG_; }

    //! The gravitational acceleration for one phase on "outside" faces (used on surface grids)
    const std::vector<Vector>& gravityOutside(unsigned int pIdx) const { return outsideG_[pIdx]; }
    std::vector<Vector>& gravityOutside(unsigned int pIdx) { return outsideG_[pIdx]; }

private:
    //! advection-related variables
    Matrix T_;                                               //!< The transmissibilities such that f_i = T_ij*p_j
    Matrix CA_;                                              //!< Matrix to project gravitational acceleration to all scvfs
    Matrix A_;                                               //!< Matrix additionally needed for the projection on surface grids
    std::array< Vector, numPhases > p_;                      //!< The interaction volume-wide phase pressures
    std::array< Vector, numPhases > g_;                      //!< The gravitational acceleration at each scvf (only for enabled gravity)
    std::vector< std::vector<Vector> > outsideT_;            //!< The transmissibilities for "outside" faces (only on surface grids)
175
    std::array< std::vector<Vector>, numPhases > outsideG_;  //!< The gravitational acceleration on "outside" faces (only on surface grids)
Dennis Gläser's avatar
Dennis Gläser committed
176
177
178
179
180
181
};

//! Data handle for quantities related to diffusion
template<class TypeTag, class InteractionVolume, bool EnableDiffusion>
class DiffusionDataHandle
{
182
    // export matrix & vector types from interaction volume
Dennis Gläser's avatar
Dennis Gläser committed
183
184
185
    using Matrix = typename InteractionVolume::Traits::Matrix;
    using Vector = typename InteractionVolume::Traits::Vector;
    using OutsideTContainer = std::vector< std::vector<Vector> >;
Dennis Gläser's avatar
Dennis Gläser committed
186

Dennis Gläser's avatar
Dennis Gläser committed
187
188
    static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
    static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
189

Dennis Gläser's avatar
Dennis Gläser committed
190
191
    static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
    static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
Dennis Gläser's avatar
Dennis Gläser committed
192

Dennis Gläser's avatar
Dennis Gläser committed
193
194
195
196
public:
    //! diffusion caches need to set phase and component index
    void setPhaseIndex(unsigned int phaseIdx) { phaseIdx_ = phaseIdx; }
    void setComponentIndex(unsigned int compIdx) { compIdx_ = compIdx; }
Dennis Gläser's avatar
Dennis Gläser committed
197

Dennis Gläser's avatar
Dennis Gläser committed
198
199
200
201
    //! prepare data handle for subsequent fill
    void resize(const InteractionVolume& iv)
    {
        for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
Dennis Gläser's avatar
Dennis Gläser committed
202
        {
Dennis Gläser's avatar
Dennis Gläser committed
203
204
205
206
207
            for (unsigned int cIdx = 0; cIdx < numComponents; ++cIdx)
            {
                if (pIdx == cIdx)
                    continue;

208
                // resize transmissibility matrix & mole fraction vector
Dennis Gläser's avatar
Dennis Gläser committed
209
210
211
                T_[pIdx][cIdx].resize(iv.numFaces(), iv.numKnowns());
                xj_[pIdx][cIdx].resize(iv.numKnowns());

212
                // resize outsideTij on surface grids
Dennis Gläser's avatar
Dennis Gläser committed
213
214
215
216
217
218
219
220
221
                if (dim < dimWorld)
                {
                    outsideT_[pIdx][cIdx].resize(iv.numFaces());

                    using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType;
                    for (LocalIndexType i = 0; i < iv.numFaces(); ++i)
                      outsideT_[pIdx][cIdx][i].resize(iv.localScvf(i).neighboringLocalScvIndices().size()-1);
                }
            }
Dennis Gläser's avatar
Dennis Gläser committed
222
        }
Dennis Gläser's avatar
Dennis Gläser committed
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
    }

    //! Access to the iv-wide mole fractions of a component in one phase
    const Vector& moleFractions(unsigned int pIdx, unsigned int compIdx) const { return xj_[pIdx][compIdx]; }
    Vector& moleFractions(unsigned int pIdx, unsigned int compIdx) { return xj_[pIdx][compIdx]; }

    //! The transmissibilities associated with diffusive fluxes
    const Matrix& diffusionT() const { return T_[phaseIdx_][compIdx_]; }
    Matrix& diffusionT() { return T_[phaseIdx_][compIdx_]; }

    //! The transmissibilities for "outside" faces (used on surface grids)
    const OutsideTContainer& diffusionTout() const { return outsideT_[phaseIdx_][compIdx_]; }
    OutsideTContainer& diffusionTout() { return outsideT_[phaseIdx_][compIdx_]; }

private:
    //! diffusion-related variables
239
240
    unsigned int phaseIdx_;                                         //!< The phase index set for the context
    unsigned int compIdx_;                                          //!< The component index set for the context
Dennis Gläser's avatar
Dennis Gläser committed
241
242
243
244
245
246
247
248
249
250
251
252
    std::array< std::array<Matrix, numComponents>, numPhases > T_;  //!< The transmissibilities such that f_i = T_ij*x_j
    std::array< std::array<Vector, numComponents>, numPhases > xj_; //!< The interaction volume-wide mole fractions
    std::array< std::array<OutsideTContainer, numComponents>, numPhases> outsideT_;
};

//! Data handle for quantities related to heat conduction
template<class TypeTag, class InteractionVolume, bool EnableHeatConduction>
class HeatConductionDataHandle
{
    //! export matrix & vector types from interaction volume
    using Matrix = typename InteractionVolume::Traits::Matrix;
    using Vector = typename InteractionVolume::Traits::Vector;
Dennis Gläser's avatar
Dennis Gläser committed
253

Dennis Gläser's avatar
Dennis Gläser committed
254
255
    static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
    static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
Dennis Gläser's avatar
Dennis Gläser committed
256

Dennis Gläser's avatar
Dennis Gläser committed
257
258
259
260
261
262
263
public:
    //! prepare data handle for subsequent fill
    void resize(const InteractionVolume& iv)
    {
        //! resize transmissibility matrix & temperature vector
        T_.resize(iv.numFaces(), iv.numKnowns());
        Tj_.resize(iv.numKnowns());
Dennis Gläser's avatar
Dennis Gläser committed
264

Dennis Gläser's avatar
Dennis Gläser committed
265
266
        //! resize outsideTij on surface grids
        if (dim < dimWorld)
Dennis Gläser's avatar
Dennis Gläser committed
267
        {
Dennis Gläser's avatar
Dennis Gläser committed
268
            outsideT_.resize(iv.numFaces());
Dennis Gläser's avatar
Dennis Gläser committed
269

Dennis Gläser's avatar
Dennis Gläser committed
270
271
272
            using LocalIndexType = typename InteractionVolume::Traits::LocalIndexType;
            for (LocalIndexType i = 0; i < iv.numFaces(); ++i)
              outsideT_[i].resize(iv.localScvf(i).neighboringLocalScvIndices().size()-1);
Dennis Gläser's avatar
Dennis Gläser committed
273
        }
Dennis Gläser's avatar
Dennis Gläser committed
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
    }

    //! Access to the iv-wide temperatures
    const Vector& temperatures() const { return Tj_; }
    Vector& temperatures() { return Tj_; }

    //! The transmissibilities associated with conductive fluxes
    const Matrix& heatConductionT() const { return T_; }
    Matrix& heatConductionT() { return T_; }

    //! The transmissibilities for "outside" faces (used on surface grids)
    const std::vector< std::vector<Vector> >& heatConductionTout() const { return outsideT_; }
    std::vector< std::vector<Vector> >& heatConductionTout() { return outsideT_; }

private:
    // heat conduction-related variables
    Matrix T_;                                    //!< The transmissibilities such that f_i = T_ij*T_j
    Vector Tj_;                                   //!< The interaction volume-wide temperatures
    std::vector< std::vector<Vector> > outsideT_; //!< The transmissibilities for "outside" faces (only necessary on surface grids)
};

//! Process-dependent data handles when related process is disabled
template<class TypeTag, class InteractionVolume>
class AdvectionDataHandle<TypeTag, InteractionVolume, false> : public EmptyDataHandle<InteractionVolume> {};
template<class TypeTag, class InteractionVolume>
class DiffusionDataHandle<TypeTag, InteractionVolume, false> : public EmptyDataHandle<InteractionVolume> {};
template<class TypeTag, class InteractionVolume>
class HeatConductionDataHandle<TypeTag, InteractionVolume, false> : public EmptyDataHandle<InteractionVolume> {};

//! Interaction volume data handle class
template<class TypeTag, class InteractionVolume>
class InteractionVolumeDataHandle : public AdvectionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableAdvection)>,
                                    public DiffusionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion)>,
                                    public HeatConductionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>
{
    using AdvectionHandle = AdvectionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableAdvection)>;
    using DiffusionHandle = DiffusionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion)>;
    using HeatConductionHandle = HeatConductionDataHandle<TypeTag, InteractionVolume, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>;
Dennis Gläser's avatar
Dennis Gläser committed
312

Dennis Gläser's avatar
Dennis Gläser committed
313
314
315
316
317
318
319
320
321
public:
    //! prepare data handles for subsequent fills
    void resize(const InteractionVolume& iv)
    {
        AdvectionHandle::resize(iv);
        DiffusionHandle::resize(iv);
        HeatConductionHandle::resize(iv);
    }
};
Dennis Gläser's avatar
Dennis Gläser committed
322

323
} // end namespace Dumux
Dennis Gläser's avatar
Dennis Gläser committed
324
325

#endif