fluxvariablescachefiller.hh 28.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// -*- 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/>.   *
 *****************************************************************************/
Dennis Gläser's avatar
Dennis Gläser committed
19
 /*!
20
 * \file
Dennis Gläser's avatar
Dennis Gläser committed
21
22
 * \ingroup CCMpfaDiscretization
 * \brief A helper class to fill the flux variable caches used in the flux constitutive laws
23
 */
24
25
#ifndef DUMUX_DISCRETIZATION_CCMPFA_FLUXVARSCACHE_FILLER_HH
#define DUMUX_DISCRETIZATION_CCMPFA_FLUXVARSCACHE_FILLER_HH
26

Dennis Gläser's avatar
Dennis Gläser committed
27
#include <dumux/common/properties.hh>
Dennis Gläser's avatar
Dennis Gläser committed
28
#include <dumux/common/parameters.hh>
Dennis Gläser's avatar
Dennis Gläser committed
29

30
#include <dumux/discretization/methods.hh>
31
#include <dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh>
Dennis Gläser's avatar
Dennis Gläser committed
32
#include <dumux/discretization/cellcentered/mpfa/localassembler.hh>
33
34
35

namespace Dumux
{
36

37
/*!
Dennis Gläser's avatar
Dennis Gläser committed
38
39
40
 * \ingroup CCMpfaDiscretization
 * \brief Helper class to fill the flux variables caches within
 *        the interaction volume around a given sub-control volume face.
41
42
 */
template<class TypeTag>
43
class CCMpfaFluxVariablesCacheFiller
44
45
46
{
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
Dennis Gläser's avatar
Dennis Gläser committed
47
48
    using Element = typename GridView::template Codim<0>::Entity;

49
50
51
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
52
    using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
Dennis Gläser's avatar
Dennis Gläser committed
53

Dennis Gläser's avatar
Dennis Gläser committed
54
    using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
Dennis Gläser's avatar
Dennis Gläser committed
55
56
    using PrimaryDataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
    using PrimaryLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData;
Dennis Gläser's avatar
Dennis Gläser committed
57
    using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
Dennis Gläser's avatar
Dennis Gläser committed
58
59
    using SecondaryDataHandle = typename SecondaryInteractionVolume::Traits::DataHandle;
    using SecondaryLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
60

Dennis Gläser's avatar
Dennis Gläser committed
61
62
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
63

64
65
66
67
68
69
70
71
    static constexpr bool doAdvection = GET_PROP_VALUE(TypeTag, EnableAdvection);
    static constexpr bool doDiffusion = GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion);
    static constexpr bool doHeatConduction = GET_PROP_VALUE(TypeTag, EnableEnergyBalance);

    static constexpr bool soldependentAdvection = GET_PROP_VALUE(TypeTag, SolutionDependentAdvection);
    static constexpr bool soldependentDiffusion = GET_PROP_VALUE(TypeTag, SolutionDependentMolecularDiffusion);
    static constexpr bool soldependentHeatConduction = GET_PROP_VALUE(TypeTag, SolutionDependentHeatConduction);

72
public:
Dennis Gläser's avatar
Dennis Gläser committed
73
74
75
76
    //! This cache filler is always solution-dependent, as it updates the
    //! vectors of cell unknowns with which the transmissibilities have to be
    //! multiplied in order to obtain the fluxes.
    static constexpr bool isSolDependent = true;
77

Dennis Gläser's avatar
Dennis Gläser committed
78
    //! The constructor. Sets problem pointer.
79
80
81
82
83
84
85
86
87
    CCMpfaFluxVariablesCacheFiller(const Problem& problem) : problemPtr_(&problem) {}

    /*!
     * \brief function to fill the flux variables caches
     *
     * \param fluxVarsCacheContainer Either the element or global flux variables cache
     * \param scvfFluxVarsCache The flux var cache to be updated corresponding to the given scvf
     * \param element The finite element
     * \param fvGeometry The finite volume geometry
Dennis Gläser's avatar
Dennis Gläser committed
88
     * \param elemVolVars The element volume variables (primary/secondary variables)
89
     * \param scvf The corresponding sub-control volume face
90
     * \param forceUpdateAll if true, forces all caches to be updated (even the solution-independent ones)
91
92
93
94
95
96
97
98
     */
    template<class FluxVariablesCacheContainer>
    void fill(FluxVariablesCacheContainer& fluxVarsCacheContainer,
              FluxVariablesCache& scvfFluxVarsCache,
              const Element& element,
              const FVElementGeometry& fvGeometry,
              const ElementVolumeVariables& elemVolVars,
              const SubControlVolumeFace& scvf,
99
              bool forceUpdateAll = false)
100
    {
101
102
103
104
105
106
        // Set pointers
        elementPtr_ = &element;
        fvGeometryPtr_ = &fvGeometry;
        elemVolVarsPtr_ = &elemVolVars;

        // prepare interaction volume and fill caches of all the scvfs connected to it
Dennis Gläser's avatar
Dennis Gläser committed
107
108
        const auto& fvGridGeometry = fvGeometry.fvGridGeometry();
        if (fvGridGeometry.vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
109
        {
110
            if (forceUpdateAll)
Dennis Gläser's avatar
Dennis Gläser committed
111
112
            {
                // the local index of the interaction volume to be created in its container
Dennis Gläser's avatar
Dennis Gläser committed
113
                const auto ivIndexInContainer = fluxVarsCacheContainer.secondaryInteractionVolumes().size();
Dennis Gläser's avatar
Dennis Gläser committed
114
115

                // prepare the locally cached boundary interaction volume
116
                const auto& indexSet = fvGridGeometry.gridInteractionVolumeIndexSets().secondaryIndexSet(scvf);
Dennis Gläser's avatar
Dennis Gläser committed
117
118
                fluxVarsCacheContainer.secondaryInteractionVolumes().emplace_back();
                secondaryIv_ = &fluxVarsCacheContainer.secondaryInteractionVolumes().back();
119
                secondaryIv_->setUpLocalScope(indexSet, problem(), fvGeometry);
Dennis Gläser's avatar
Dennis Gläser committed
120
121

                // prepare the corresponding data handle
Dennis Gläser's avatar
Dennis Gläser committed
122
123
124
                fluxVarsCacheContainer.secondaryDataHandles().emplace_back();
                secondaryIvDataHandle_ = &fluxVarsCacheContainer.secondaryDataHandles().back();
                secondaryIvDataHandle_->resize(*secondaryIv_);
125

Dennis Gläser's avatar
Dennis Gläser committed
126
                // fill the caches for all the scvfs in the interaction volume
Dennis Gläser's avatar
Dennis Gläser committed
127
                fillCachesInInteractionVolume_(fluxVarsCacheContainer, *secondaryIv_, *secondaryIvDataHandle_, ivIndexInContainer, true);
Dennis Gläser's avatar
Dennis Gläser committed
128
129
130
131
            }
            else
            {
                const auto ivIndexInContainer = scvfFluxVarsCache.ivIndexInContainer();
Dennis Gläser's avatar
Dennis Gläser committed
132
133
                secondaryIv_ = &fluxVarsCacheContainer.secondaryInteractionVolumes()[ivIndexInContainer];
                secondaryIvDataHandle_ = &fluxVarsCacheContainer.secondaryDataHandles()[ivIndexInContainer];
Dennis Gläser's avatar
Dennis Gläser committed
134
135

                // fill the caches for all the scvfs in the interaction volume
Dennis Gläser's avatar
Dennis Gläser committed
136
                fillCachesInInteractionVolume_(fluxVarsCacheContainer, *secondaryIv_, *secondaryIvDataHandle_, ivIndexInContainer);
Dennis Gläser's avatar
Dennis Gläser committed
137
            }
138
139
140
        }
        else
        {
141
            if (forceUpdateAll)
Dennis Gläser's avatar
Dennis Gläser committed
142
143
            {
                // the local index of the interaction volume to be created in its container
Dennis Gläser's avatar
Dennis Gläser committed
144
                const auto ivIndexInContainer = fluxVarsCacheContainer.primaryInteractionVolumes().size();
Dennis Gläser's avatar
Dennis Gläser committed
145
146

                // prepare the locally cached boundary interaction volume
147
                const auto& indexSet = fvGridGeometry.gridInteractionVolumeIndexSets().primaryIndexSet(scvf);
Dennis Gläser's avatar
Dennis Gläser committed
148
149
                fluxVarsCacheContainer.primaryInteractionVolumes().emplace_back();
                primaryIv_ = &fluxVarsCacheContainer.primaryInteractionVolumes().back();
150
                primaryIv_->setUpLocalScope(indexSet, problem(), fvGeometry);
151

Dennis Gläser's avatar
Dennis Gläser committed
152
                // prepare the corresponding data handle
Dennis Gläser's avatar
Dennis Gläser committed
153
154
155
                fluxVarsCacheContainer.primaryDataHandles().emplace_back();
                primaryIvDataHandle_ = &fluxVarsCacheContainer.primaryDataHandles().back();
                primaryIvDataHandle_->resize(*primaryIv_);
Dennis Gläser's avatar
Dennis Gläser committed
156
157

                // fill the caches for all the scvfs in the interaction volume
Dennis Gläser's avatar
Dennis Gläser committed
158
                fillCachesInInteractionVolume_(fluxVarsCacheContainer, *primaryIv_, *primaryIvDataHandle_, ivIndexInContainer, true);
Dennis Gläser's avatar
Dennis Gläser committed
159
160
161
162
            }
            else
            {
                const auto ivIndexInContainer = scvfFluxVarsCache.ivIndexInContainer();
Dennis Gläser's avatar
Dennis Gläser committed
163
164
                primaryIv_ = &fluxVarsCacheContainer.primaryInteractionVolumes()[ivIndexInContainer];
                primaryIvDataHandle_ = &fluxVarsCacheContainer.primaryDataHandles()[ivIndexInContainer];
Dennis Gläser's avatar
Dennis Gläser committed
165
166

                // fill the caches for all the scvfs in the interaction volume
Dennis Gläser's avatar
Dennis Gläser committed
167
                fillCachesInInteractionVolume_(fluxVarsCacheContainer, *primaryIv_, *primaryIvDataHandle_, ivIndexInContainer);
Dennis Gläser's avatar
Dennis Gläser committed
168
            }
169
        }
170
    }
171

Dennis Gläser's avatar
Dennis Gläser committed
172
    //! returns the stored interaction volume pointer
Dennis Gläser's avatar
Dennis Gläser committed
173
174
175
    const PrimaryInteractionVolume& primaryInteractionVolume() const
    { return *primaryIv_; }

Dennis Gläser's avatar
Dennis Gläser committed
176
    //! returns the stored interaction volume pointer
Dennis Gläser's avatar
Dennis Gläser committed
177
178
    const SecondaryInteractionVolume& secondaryInteractionVolume() const
    { return *secondaryIv_; }
179

Dennis Gläser's avatar
Dennis Gläser committed
180
181
182
    //! returns the stored data handle pointer
    const PrimaryDataHandle& primaryIvDataHandle() const
    { return *primaryIvDataHandle_; }
183

Dennis Gläser's avatar
Dennis Gläser committed
184
185
186
    //! returns the stored data handle pointer
    const SecondaryDataHandle& secondaryIvDataHandle() const
    { return *secondaryIvDataHandle_; }
187

Dennis Gläser's avatar
Dennis Gläser committed
188
189
190
    //! returns the currently stored iv-local face data object
    const PrimaryLocalFaceData& primaryIvLocalFaceData() const
    { return *primaryLocalFaceData_; }
191

Dennis Gläser's avatar
Dennis Gläser committed
192
193
194
    //! returns the currently stored iv-local face data object
    const SecondaryLocalFaceData& secondaryIvLocalFaceData() const
    { return *secondaryLocalFaceData_; }
195

Dennis Gläser's avatar
Dennis Gläser committed
196
private:
197

Dennis Gläser's avatar
Dennis Gläser committed
198
199
200
201
    const Problem& problem() const { return *problemPtr_; }
    const Element& element() const { return *elementPtr_; }
    const FVElementGeometry& fvGeometry() const { return *fvGeometryPtr_; }
    const ElementVolumeVariables& elemVolVars() const { return *elemVolVarsPtr_; }
202

203
    //! Method to fill the flux var caches within an interaction volume
Dennis Gläser's avatar
Dennis Gläser committed
204
205
206
    template<class FluxVarsCacheContainer, class InteractionVolume, class DataHandle>
    void fillCachesInInteractionVolume_(FluxVarsCacheContainer& fluxVarsCacheContainer,
                                        InteractionVolume& iv,
Dennis Gläser's avatar
Dennis Gläser committed
207
208
                                        DataHandle& handle,
                                        unsigned int ivIndexInContainer,
209
                                        bool forceUpdateAll = false)
210
    {
Dennis Gläser's avatar
Dennis Gläser committed
211
212
213
        // First we upate data which are not dependent on the physical processes.
        // We store pointers to the other flux var caches, so that we have to obtain
        // this data only once and can use it again in the sub-cache fillers.
214
215
216
        const auto numGlobalScvfs = iv.localFaceData().size();
        std::vector<const SubControlVolumeFace*> ivScvfs(numGlobalScvfs);
        std::vector<FluxVariablesCache*> ivFluxVarCaches(numGlobalScvfs);
217

218
219
220
221
222
223
224
225
226
227
228
        unsigned int i = 0;
        for (const auto& d : iv.localFaceData())
        {
            // obtain the scvf
            const auto& scvfJ = fvGeometry().scvf(d.globalScvfIndex());
            ivScvfs[i] = &scvfJ;
            ivFluxVarCaches[i] = &fluxVarsCacheContainer[scvfJ];
            ivFluxVarCaches[i]->setIvIndexInContainer(ivIndexInContainer);
            ivFluxVarCaches[i]->setUpdateStatus(true);
            i++;
        }
229

Dennis Gläser's avatar
Dennis Gläser committed
230
231
232
        fillAdvection(iv, handle, ivScvfs, ivFluxVarCaches, forceUpdateAll);
        fillDiffusion(iv, handle, ivScvfs, ivFluxVarCaches, forceUpdateAll);
        fillHeatConduction(iv, handle, ivScvfs, ivFluxVarCaches, forceUpdateAll);
233
234
    }

Dennis Gläser's avatar
Dennis Gläser committed
235
236
237
238
239
240
241
242
243
244
    //! fills the advective quantities (enabled advection)
    template< class InteractionVolume,
              class DataHandle,
              bool enableAdvection = doAdvection,
              typename std::enable_if_t<enableAdvection, int> = 0 >
    void fillAdvection(InteractionVolume& iv,
                       DataHandle& handle,
                       const std::vector<const SubControlVolumeFace*>& ivScvfs,
                       const std::vector<FluxVariablesCache*>& ivFluxVarCaches,
                       bool forceUpdateAll = false)
245
246
    {
        using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
247
        using AdvectionFiller = typename AdvectionType::Cache::Filler;
248

249
250
        static constexpr auto AdvectionMethod = AdvectionType::myDiscretizationMethod;
        using LambdaFactory = TensorLambdaFactory<TypeTag, AdvectionMethod>;
251

Dennis Gläser's avatar
Dennis Gläser committed
252
        // skip the following if advection doesn't use mpfa
253
        if (AdvectionMethod == DiscretizationMethods::CCMpfa)
Dennis Gläser's avatar
Dennis Gläser committed
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
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
312
313
314
315
316
317
318
319
320
321
322
323
324
325
        {
            // get instance of the interaction volume-local assembler
            using IvLocalAssembler = InteractionVolumeAssembler< InteractionVolume >;
            IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());

            // Use different assembly if gravity is enabled
            static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");

            // Assemble T only if permeability is sol-dependent or if update is forced
            if (forceUpdateAll || soldependentAdvection)
            {
                // distinguish between normal/surface grids (optimized away by compiler)
                if (dim < dimWorld)
                {
                    if (enableGravity)
                        localAssembler.assembleWithGravity( handle.advectionTout(),
                                                            handle.advectionT(),
                                                            handle.gravityOutside(),
                                                            handle.gravity(),
                                                            handle.advectionCA(),
                                                            handle.advectionA(),
                                                            iv,
                                                            LambdaFactory::getAdvectionLambda() );

                    else
                        localAssembler.assemble( handle.advectionTout(),
                                                 handle.advectionT(),
                                                 iv,
                                                 LambdaFactory::getAdvectionLambda() );
                }

                // normal grids
                else
                {
                    if (enableGravity)
                        localAssembler.assembleWithGravity( handle.advectionT(),
                                                            handle.gravity(),
                                                            handle.advectionCA(),
                                                            iv,
                                                            LambdaFactory::getAdvectionLambda() );
                    else
                        localAssembler.assemble( handle.advectionT(),
                                                 iv,
                                                 LambdaFactory::getAdvectionLambda() );
                }
            }

            // (maybe) only reassemble gravity vector
            else if (enableGravity)
            {
                if (dim == dimWorld)
                    localAssembler.assembleGravity( handle.gravity(),
                                                    iv,
                                                    handle.advectionCA(),
                                                    LambdaFactory::getAdvectionLambda() );
                else
                    localAssembler.assembleGravity( handle.gravity(),
                                                    handle.gravityOutside(),
                                                    iv,
                                                    handle.advectionCA(),
                                                    handle.advectionA(),
                                                    LambdaFactory::getAdvectionLambda() );
            }

            // assemble pressure vectors
            for (unsigned int pIdx = 0; pIdx < GET_PROP_VALUE(TypeTag, NumPhases); ++pIdx)
            {
                const auto& evv = &elemVolVars();
                auto getPressure = [&evv, pIdx] (auto volVarIdx) { return (evv->operator[](volVarIdx)).pressure(pIdx); };
                localAssembler.assemble(handle.pressures(pIdx), iv, getPressure);
            }
        }
326

Dennis Gläser's avatar
Dennis Gläser committed
327
        // fill advection caches
328
        for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
Dennis Gläser's avatar
Dennis Gläser committed
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
        {
            // set pointer to current local face data object
            // ifs are evaluated at compile time and are optimized away
            if (std::is_same<PrimaryInteractionVolume, SecondaryInteractionVolume>::value)
            {
                // we cannot make a disctinction, thus we set both pointers
                primaryLocalFaceData_ = &(iv.localFaceData()[i]);
                secondaryLocalFaceData_ = &(iv.localFaceData()[i]);
            }
            else if (std::is_same<InteractionVolume, PrimaryInteractionVolume>::value)
                primaryLocalFaceData_ = &(iv.localFaceData()[i]);
            else
                secondaryLocalFaceData_ = &(iv.localFaceData()[i]);

            // fill this scvfs cache
344
            AdvectionFiller::fill(*ivFluxVarCaches[i],
345
                                  problem(),
346
                                  iv.element(iv.localFaceData()[i].ivLocalInsideScvIndex()),
347
348
                                  fvGeometry(),
                                  elemVolVars(),
349
                                  *ivScvfs[i],
350
                                  *this);
Dennis Gläser's avatar
Dennis Gläser committed
351
        }
352
    }
353

354
    //! do nothing if advection is not enabled
Dennis Gläser's avatar
Dennis Gläser committed
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
    template< class InteractionVolume,
              class DataHandle,
              bool enableAdvection = doAdvection,
              typename std::enable_if_t<!enableAdvection, int> = 0 >
    void fillAdvection(InteractionVolume& iv,
                       DataHandle& handle,
                       const std::vector<const SubControlVolumeFace*>& ivScvfs,
                       const std::vector<FluxVariablesCache*>& ivFluxVarCaches,
                       bool forceUpdateAll = false)
    {}

    //! fills the diffusive quantities (diffusion enabled)
    template< class InteractionVolume,
              class DataHandle,
              bool enableDiffusion = doDiffusion,
              typename std::enable_if_t<enableDiffusion, int> = 0 >
    void fillDiffusion(InteractionVolume& iv,
                       DataHandle& handle,
                       const std::vector<const SubControlVolumeFace*>& ivScvfs,
                       const std::vector<FluxVariablesCache*>& ivFluxVarCaches,
                       bool forceUpdateAll = false)
376
377
    {
        using DiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
378
        using DiffusionFiller = typename DiffusionType::Cache::Filler;
379

380
381
382
        static constexpr auto DiffusionMethod = DiffusionType::myDiscretizationMethod;
        using LambdaFactory = TensorLambdaFactory<TypeTag, DiffusionMethod>;

383
384
        static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
        static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
385

Dennis Gläser's avatar
Dennis Gläser committed
386
387
388
389
        // get instance of the interaction volume-local assembler
        using IvLocalAssembler = InteractionVolumeAssembler< InteractionVolume >;
        IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());

390
        for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
391
        {
392
            for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
393
            {
394
395
396
397
                if (phaseIdx == compIdx)
                    continue;

                // solve the local system subject to the diffusion tensor (if uses mpfa)
398
                if (DiffusionMethod == DiscretizationMethods::CCMpfa)
Dennis Gläser's avatar
Dennis Gläser committed
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
                {
                    // update the context in handle
                    handle.setPhaseIndex(phaseIdx);
                    handle.setComponentIndex(compIdx);

                    // assemble T
                    if (forceUpdateAll || soldependentDiffusion)
                    {
                        if (dim < dimWorld)
                            localAssembler.assemble( handle.diffusionTout(),
                                                     handle.diffusionT(),
                                                     iv,
                                                     LambdaFactory::getDiffusionLambda(phaseIdx, compIdx) );
                        else
                            localAssembler. assemble( handle.diffusionT(),
                                                      iv,
                                                      LambdaFactory::getDiffusionLambda(phaseIdx, compIdx) );
                    }

                    // assemble vector of mole fractions
                    const auto& evv = &elemVolVars();
                    auto getMoleFraction = [&evv, phaseIdx, compIdx] (auto volVarIdx)
                                           { return (evv->operator[](volVarIdx)).moleFraction(phaseIdx, compIdx); };
                    localAssembler.assemble(handle.moleFractions(phaseIdx, compIdx), iv, getMoleFraction);
                }
424

Dennis Gläser's avatar
Dennis Gläser committed
425
                // fill diffusion caches
426
                for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
Dennis Gläser's avatar
Dennis Gläser committed
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
                {
                    // set pointer to current local face data object
                    // ifs are evaluated at compile time and are optimized away
                    if (std::is_same<PrimaryInteractionVolume, SecondaryInteractionVolume>::value)
                    {
                        // we cannot make a disctinction, thus we set both pointers
                        primaryLocalFaceData_ = &(iv.localFaceData()[i]);
                        secondaryLocalFaceData_ = &(iv.localFaceData()[i]);
                    }
                    else if (std::is_same<InteractionVolume, PrimaryInteractionVolume>::value)
                        primaryLocalFaceData_ = &(iv.localFaceData()[i]);
                    else
                        secondaryLocalFaceData_ = &(iv.localFaceData()[i]);

                    // fill this scvfs cache
442
                    DiffusionFiller::fill(*ivFluxVarCaches[i],
443
444
445
                                          phaseIdx,
                                          compIdx,
                                          problem(),
446
                                          iv.element(iv.localFaceData()[i].ivLocalInsideScvIndex()),
447
448
                                          fvGeometry(),
                                          elemVolVars(),
449
                                          *ivScvfs[i],
450
                                          *this);
Dennis Gläser's avatar
Dennis Gläser committed
451
                }
452
453
            }
        }
454
455
    }

456
    //! do nothing if diffusion is not enabled
Dennis Gläser's avatar
Dennis Gläser committed
457
458
459
460
461
    template< class InteractionVolume,
              class DataHandle,
              bool enableDiffusion = doDiffusion,
              typename std::enable_if_t<!enableDiffusion, int> = 0 >
    void fillDiffusion(InteractionVolume& iv,
Dennis Gläser's avatar
Dennis Gläser committed
462
                       DataHandle& handle,
463
                       const std::vector<const SubControlVolumeFace*>& ivScvfs,
Dennis Gläser's avatar
Dennis Gläser committed
464
465
466
467
468
469
470
471
472
473
474
475
476
477
                       const std::vector<FluxVariablesCache*>& ivFluxVarCaches,
                       bool forceUpdateAll = false)
    {}

    //! fills the quantities related to heat conduction (heat conduction enabled)
    template< class InteractionVolume,
              class DataHandle,
              bool enableHeatConduction = doHeatConduction,
              typename std::enable_if_t<enableHeatConduction, int> = 0 >
    void fillHeatConduction(InteractionVolume& iv,
                            DataHandle& handle,
                            const std::vector<const SubControlVolumeFace*>& ivScvfs,
                            const std::vector<FluxVariablesCache*>& ivFluxVarCaches,
                            bool forceUpdateAll = false)
478
    {
479
        using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
480
        using HeatConductionFiller = typename HeatConductionType::Cache::Filler;
481

482
483
        static constexpr auto HeatConductionMethod = HeatConductionType::myDiscretizationMethod;
        using LambdaFactory = TensorLambdaFactory<TypeTag, HeatConductionMethod>;
484

485
486
        // maybe solve the local system subject to fourier coefficient
        if (HeatConductionMethod == DiscretizationMethods::CCMpfa)
Dennis Gläser's avatar
Dennis Gläser committed
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
        {
            // get instance of the interaction volume-local assembler
            using IvLocalAssembler = InteractionVolumeAssembler< InteractionVolume >;
            IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());

            if (forceUpdateAll || soldependentAdvection)
            {
                if (dim < dimWorld)
                    localAssembler.assemble( handle.heatConductionTout(),
                                             handle.heatConductionT(),
                                             iv,
                                             LambdaFactory::getHeatConductionLambda() );
                else
                    localAssembler.assemble( handle.heatConductionT(),
                                             iv,
                                             LambdaFactory::getHeatConductionLambda() );
            }

            // assemble vector of temperatures
            const auto& evv = &elemVolVars();
            auto getMoleFraction = [&evv] (auto volVarIdx) { return (evv->operator[](volVarIdx)).temperature(); };
            localAssembler.assemble(handle.temperatures(), iv, getMoleFraction);
        }
510

Dennis Gläser's avatar
Dennis Gläser committed
511
        // fill heat conduction caches
512
        for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
Dennis Gläser's avatar
Dennis Gläser committed
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
        {
            // set pointer to current local face data object
            // ifs are evaluated at compile time and are optimized away
            if (std::is_same<PrimaryInteractionVolume, SecondaryInteractionVolume>::value)
            {
                // we cannot make a disctinction, thus we set both pointers
                primaryLocalFaceData_ = &(iv.localFaceData()[i]);
                secondaryLocalFaceData_ = &(iv.localFaceData()[i]);
            }
            else if (std::is_same<InteractionVolume, PrimaryInteractionVolume>::value)
                primaryLocalFaceData_ = &(iv.localFaceData()[i]);
            else
                secondaryLocalFaceData_ = &(iv.localFaceData()[i]);

            // fill this scvfs cache
528
            HeatConductionFiller::fill(*ivFluxVarCaches[i],
529
                                       problem(),
530
                                       iv.element(iv.localFaceData()[i].ivLocalInsideScvIndex()),
531
532
                                       fvGeometry(),
                                       elemVolVars(),
533
                                       *ivScvfs[i],
534
                                       *this);
Dennis Gläser's avatar
Dennis Gläser committed
535
        }
536
537
    }

538
    //! do nothing if heat conduction is disabled
Dennis Gläser's avatar
Dennis Gläser committed
539
540
541
542
543
544
545
546
547
548
    template< class InteractionVolume,
              class DataHandle,
              bool enableHeatConduction = doHeatConduction,
              typename std::enable_if_t<!enableHeatConduction, int> = 0 >
    void fillHeatConduction(InteractionVolume& iv,
                            DataHandle& handle,
                            const std::vector<const SubControlVolumeFace*>& ivScvfs,
                            const std::vector<FluxVariablesCache*>& ivFluxVarCaches,
                            bool forceUpdateAll = false)
    {}
549
550
551
552
553
554

    const Problem* problemPtr_;
    const Element* elementPtr_;
    const FVElementGeometry* fvGeometryPtr_;
    const ElementVolumeVariables* elemVolVarsPtr_;

555
556
    // We store pointers to an inner and a boundary interaction volume.
    // These are updated during the filling of the caches and the
557
    // physics-related caches have access to them
Dennis Gläser's avatar
Dennis Gläser committed
558
559
560
561
    PrimaryInteractionVolume* primaryIv_;
    SecondaryInteractionVolume* secondaryIv_;

    // pointer to the current interaction volume data handle
Dennis Gläser's avatar
Dennis Gläser committed
562
563
564
565
566
567
568
569
570
    PrimaryDataHandle* primaryIvDataHandle_;
    SecondaryDataHandle* secondaryIvDataHandle_;

    // We do an interaction volume-wise filling of the caches
    // While filling, we store a pointer to the current localScvf
    // face data object of the IV so that the individual caches
    // can access it and don't have to retrieve it again
    const PrimaryLocalFaceData* primaryLocalFaceData_;
    const SecondaryLocalFaceData* secondaryLocalFaceData_;
571
572
};

Dennis Gläser's avatar
Dennis Gläser committed
573
} // end namespace Dumux
574
575

#endif