fluxvariablescachefiller.hh 35.5 KB
Newer Older
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
27
28
29
30
// -*- 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 3 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
 * \ingroup PorousmediumflowModels
 * \brief A helper class to fill the flux variables cache
 */
#ifndef DUMUX_POROUSMEDIUM_FLUXVARIABLESCACHE_FILLER_HH
#define DUMUX_POROUSMEDIUM_FLUXVARIABLESCACHE_FILLER_HH

#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>

#include <dumux/discretization/method.hh>
31
#include <dumux/discretization/extrusion.hh>
32
#include <dumux/flux/referencesystemformulation.hh>
33
#include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh>
34
35
36
37
38
39
40
41
42
43
44
45
46
47

namespace Dumux {

// forward declaration
template<class TypeTag, DiscretizationMethod discMethod>
class PorousMediumFluxVariablesCacheFillerImplementation;

/*!
 * \ingroup PorousmediumflowModels
 * \brief The flux variables cache filler class for porous media
 *
 * Helps filling the flux variables cache depending several policies
 */
template<class TypeTag>
48
using PorousMediumFluxVariablesCacheFiller = PorousMediumFluxVariablesCacheFillerImplementation<TypeTag, GetPropType<TypeTag, Properties::GridGeometry>::discMethod>;
49
50
51
52
53
54
55

//! Specialization of the flux variables cache filler for the cell centered tpfa method
template<class TypeTag>
class PorousMediumFluxVariablesCacheFillerImplementation<TypeTag, DiscretizationMethod::cctpfa>
{
    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
56
57
58
59
60
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using GridView = typename GridGeometry::GridView;
    using FVElementGeometry = typename GridGeometry::LocalView;
    using SubControlVolume = typename GridGeometry::SubControlVolume;
    using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
61
62
63
64
    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;

    using Element = typename GridView::template Codim<0>::Entity;

65
66
67
    static constexpr bool advectionEnabled = ModelTraits::enableAdvection();
    static constexpr bool diffusionEnabled = ModelTraits::enableMolecularDiffusion();
    static constexpr bool heatConductionEnabled = ModelTraits::enableEnergyBalance();
68

69
70
71
    static constexpr bool advectionIsSolDependent = getPropValue<TypeTag, Properties::SolutionDependentAdvection>();
    static constexpr bool diffusionIsSolDependent = getPropValue<TypeTag, Properties::SolutionDependentMolecularDiffusion>();
    static constexpr bool heatConductionIsSolDependent = getPropValue<TypeTag, Properties::SolutionDependentHeatConduction>();
72
73
74


public:
75
76
77
    static constexpr bool isSolDependent = (advectionEnabled && advectionIsSolDependent) ||
                                           (diffusionEnabled && diffusionIsSolDependent) ||
                                           (heatConductionEnabled && heatConductionIsSolDependent);
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

    //! The constructor. Sets the problem pointer
    PorousMediumFluxVariablesCacheFillerImplementation(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
     * \param elemVolVars The element volume variables
     * \param scvf The corresponding sub-control volume face
     * \param forceUpdateAll if true, forces all caches to be updated (even the solution-independent ones)
     */
94
    template<class FluxVariablesCacheContainer, class FluxVariablesCache>
95
96
97
98
99
100
101
102
103
104
105
    void fill(FluxVariablesCacheContainer& fluxVarsCacheContainer,
              FluxVariablesCache& scvfFluxVarsCache,
              const Element& element,
              const FVElementGeometry& fvGeometry,
              const ElementVolumeVariables& elemVolVars,
              const SubControlVolumeFace& scvf,
              bool forceUpdateAll = false)
    {
        // fill the physics-related quantities of the caches
        if (forceUpdateAll)
        {
106
107
108
109
110
111
            if constexpr (advectionEnabled)
                fillAdvection_(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
            if constexpr (diffusionEnabled)
                fillDiffusion_(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
            if constexpr (heatConductionEnabled)
                fillHeatConduction_(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
112
113
114
        }
        else
        {
115
            if constexpr (advectionEnabled && advectionIsSolDependent)
116
                fillAdvection_(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
117
            if constexpr (diffusionEnabled && diffusionIsSolDependent)
118
                fillDiffusion_(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
119
            if constexpr (heatConductionEnabled && heatConductionIsSolDependent)
120
121
122
123
124
125
126
127
128
129
                fillHeatConduction_(scvfFluxVarsCache, element, fvGeometry, elemVolVars, scvf);
        }
    }

private:

    const Problem& problem() const
    { return *problemPtr_; }

    //! method to fill the advective quantities
130
131
132
133
134
135
    template<class FluxVariablesCache>
    void fillAdvection_(FluxVariablesCache& scvfFluxVarsCache,
                        const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
                        const SubControlVolumeFace& scvf)
136
137
138
139
140
141
142
143
144
    {
        using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
        using AdvectionFiller = typename AdvectionType::Cache::Filler;

        // forward to the filler for the advective quantities
        AdvectionFiller::fill(scvfFluxVarsCache, problem(), element, fvGeometry, elemVolVars, scvf, *this);
    }

    //! method to fill the diffusive quantities
145
146
147
148
149
150
    template<class FluxVariablesCache>
    void fillDiffusion_(FluxVariablesCache& scvfFluxVarsCache,
                        const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
                        const SubControlVolumeFace& scvf)
151
152
153
154
155
156
157
158
159
    {
        using DiffusionType = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
        using DiffusionFiller = typename DiffusionType::Cache::Filler;
        using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;

        static constexpr int numPhases = ModelTraits::numFluidPhases();
        static constexpr int numComponents = ModelTraits::numFluidComponents();

        // forward to the filler of the diffusive quantities
160
161
162
        if constexpr (FluidSystem::isTracerFluidSystem())
            for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
                for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
163
                    DiffusionFiller::fill(scvfFluxVarsCache, phaseIdx, compIdx, problem(), element, fvGeometry, elemVolVars, scvf, *this);
164
165
166
167
168
        else
            for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
                for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
                    if (compIdx != FluidSystem::getMainComponent(phaseIdx))
                        DiffusionFiller::fill(scvfFluxVarsCache, phaseIdx, compIdx, problem(), element, fvGeometry, elemVolVars, scvf, *this);
169
170
171
    }

    //! method to fill the quantities related to heat conduction
172
173
174
175
176
177
    template<class FluxVariablesCache>
    void fillHeatConduction_(FluxVariablesCache& scvfFluxVarsCache,
                             const Element& element,
                             const FVElementGeometry& fvGeometry,
                             const ElementVolumeVariables& elemVolVars,
                             const SubControlVolumeFace& scvf)
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
    {
        using HeatConductionType = GetPropType<TypeTag, Properties::HeatConductionType>;
        using HeatConductionFiller = typename HeatConductionType::Cache::Filler;

        // forward to the filler of the diffusive quantities
        HeatConductionFiller::fill(scvfFluxVarsCache, problem(), element, fvGeometry, elemVolVars, scvf, *this);
    }

    const Problem* problemPtr_;
};

//! Specialization of the flux variables cache filler for the cell centered mpfa method
template<class TypeTag>
class PorousMediumFluxVariablesCacheFillerImplementation<TypeTag, DiscretizationMethod::ccmpfa>
{
    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
195
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
196
    using Element = typename GridView::template Codim<0>::Entity;
197
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
198

199
200
201
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GridGeometry::LocalView;
    using MpfaHelper = typename GridGeometry::MpfaHelper;
202
203
    using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
    using Extrusion = Extrusion_t<GridGeometry>;
204
205
206
207
208
209
210
211
212
213
214
215
216
    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
    using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;

    using PrimaryInteractionVolume = GetPropType<TypeTag, Properties::PrimaryInteractionVolume>;
    using PrimaryDataHandle = typename ElementFluxVariablesCache::PrimaryIvDataHandle;
    using PrimaryLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData;
    using SecondaryInteractionVolume = GetPropType<TypeTag, Properties::SecondaryInteractionVolume>;
    using SecondaryDataHandle = typename ElementFluxVariablesCache::SecondaryIvDataHandle;
    using SecondaryLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;

    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;

217
218
219
    static constexpr bool advectionEnabled = ModelTraits::enableAdvection();
    static constexpr bool diffusionEnabled = ModelTraits::enableMolecularDiffusion();
    static constexpr bool heatConductionEnabled = ModelTraits::enableEnergyBalance();
220

221
222
223
    static constexpr bool advectionIsSolDependent = getPropValue<TypeTag, Properties::SolutionDependentAdvection>();
    static constexpr bool diffusionIsSolDependent = getPropValue<TypeTag, Properties::SolutionDependentMolecularDiffusion>();
    static constexpr bool heatConductionIsSolDependent = getPropValue<TypeTag, Properties::SolutionDependentHeatConduction>();
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246

public:
    //! 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;

    //! The constructor. Sets problem pointer.
    PorousMediumFluxVariablesCacheFillerImplementation(const Problem& problem)
    : problemPtr_(&problem) {}

    /*!
     * \brief function to fill the flux variables caches
     *
     * \param fluxVarsCacheStorage Class that holds the scvf flux vars caches
     * \param scvfFluxVarsCache The flux var cache to be updated corresponding to the given scvf
     * \param ivDataStorage Class that stores the interaction volumes & handles
     * \param element The finite element
     * \param fvGeometry The finite volume geometry
     * \param elemVolVars The element volume variables (primary/secondary variables)
     * \param scvf The corresponding sub-control volume face
     * \param forceUpdateAll if true, forces all caches to be updated (even the solution-independent ones)
     */
247
    template<class FluxVarsCacheStorage, class FluxVariablesCache, class IVDataStorage>
248
249
250
251
252
253
254
255
256
257
258
    void fill(FluxVarsCacheStorage& fluxVarsCacheStorage,
              FluxVariablesCache& scvfFluxVarsCache,
              IVDataStorage& ivDataStorage,
              const FVElementGeometry& fvGeometry,
              const ElementVolumeVariables& elemVolVars,
              const SubControlVolumeFace& scvf,
              bool forceUpdateAll = false)
    {
        // Set pointers
        fvGeometryPtr_ = &fvGeometry;
        elemVolVarsPtr_ = &elemVolVars;
259
        const auto& gridGeometry = fvGeometry.gridGeometry();
260
261
262
263

        // 1. prepare interaction volume (iv)
        // 2. solve for all transmissibilities and store them in data handles
        // 3. set pointers to transmissibilities in caches of all the scvfs of the iv
264
        if (gridGeometry.vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
265
266
267
        {
            if (forceUpdateAll)
            {
268
                // create new interaction volume
269
                const auto ivIndexInContainer = ivDataStorage.secondaryInteractionVolumes.size();
270
                const auto& indexSet = gridGeometry.gridInteractionVolumeIndexSets().secondaryIndexSet(scvf);
271
272
273
274
                ivDataStorage.secondaryInteractionVolumes.emplace_back();
                secondaryIv_ = &ivDataStorage.secondaryInteractionVolumes.back();
                secondaryIv_->bind(indexSet, problem(), fvGeometry);

275
                // create the corresponding data handle
276
277
                ivDataStorage.secondaryDataHandles.emplace_back();
                secondaryIvDataHandle_ = &ivDataStorage.secondaryDataHandles.back();
278
                prepareDataHandle_(*secondaryIv_, *secondaryIvDataHandle_, forceUpdateAll);
279
280

                // fill the caches for all the scvfs in the interaction volume
281
                fillCachesInInteractionVolume_<FluxVariablesCache>(fluxVarsCacheStorage, *secondaryIv_, ivIndexInContainer);
282
283
284
            }
            else
            {
285
                // get previously created interaction volume/handle
286
287
288
                const auto ivIndexInContainer = scvfFluxVarsCache.ivIndexInContainer();
                secondaryIv_ = &ivDataStorage.secondaryInteractionVolumes[ivIndexInContainer];
                secondaryIvDataHandle_ = &ivDataStorage.secondaryDataHandles[ivIndexInContainer];
289
                prepareDataHandle_(*secondaryIv_, *secondaryIvDataHandle_, forceUpdateAll);
290
291

                // fill the caches for all the scvfs in the interaction volume
292
                fillCachesInInteractionVolume_<FluxVariablesCache>(fluxVarsCacheStorage, *secondaryIv_, ivIndexInContainer);
293
294
            }
        }
295
296

        // primary interaction volume type
297
298
299
300
        else
        {
            if (forceUpdateAll)
            {
301
                // create new interaction volume
302
                const auto ivIndexInContainer = ivDataStorage.primaryInteractionVolumes.size();
303
                const auto& indexSet = gridGeometry.gridInteractionVolumeIndexSets().primaryIndexSet(scvf);
304
305
306
307
                ivDataStorage.primaryInteractionVolumes.emplace_back();
                primaryIv_ = &ivDataStorage.primaryInteractionVolumes.back();
                primaryIv_->bind(indexSet, problem(), fvGeometry);

308
                // create the corresponding data handle
309
310
                ivDataStorage.primaryDataHandles.emplace_back();
                primaryIvDataHandle_ = &ivDataStorage.primaryDataHandles.back();
311
                prepareDataHandle_(*primaryIv_, *primaryIvDataHandle_, forceUpdateAll);
312
313

                // fill the caches for all the scvfs in the interaction volume
314
                fillCachesInInteractionVolume_<FluxVariablesCache>(fluxVarsCacheStorage, *primaryIv_, ivIndexInContainer);
315
316
317
            }
            else
            {
318
                // get previously created interaction volume/handle
319
320
321
                const auto ivIndexInContainer = scvfFluxVarsCache.ivIndexInContainer();
                primaryIv_ = &ivDataStorage.primaryInteractionVolumes[ivIndexInContainer];
                primaryIvDataHandle_ = &ivDataStorage.primaryDataHandles[ivIndexInContainer];
322
                prepareDataHandle_(*primaryIv_, *primaryIvDataHandle_, forceUpdateAll);
323
324

                // fill the caches for all the scvfs in the interaction volume
325
                fillCachesInInteractionVolume_<FluxVariablesCache>(fluxVarsCacheStorage, *primaryIv_, ivIndexInContainer);
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
            }
        }
    }

    //! returns the stored interaction volume pointer
    const PrimaryInteractionVolume& primaryInteractionVolume() const
    { return *primaryIv_; }

    //! returns the stored interaction volume pointer
    const SecondaryInteractionVolume& secondaryInteractionVolume() const
    { return *secondaryIv_; }

    //! returns the stored data handle pointer
    const PrimaryDataHandle& primaryIvDataHandle() const
    { return *primaryIvDataHandle_; }

    //! returns the stored data handle pointer
    const SecondaryDataHandle& secondaryIvDataHandle() const
    { return *secondaryIvDataHandle_; }

    //! returns the currently stored iv-local face data object
    const PrimaryLocalFaceData& primaryIvLocalFaceData() const
    { return *primaryLocalFaceData_; }

    //! returns the currently stored iv-local face data object
    const SecondaryLocalFaceData& secondaryIvLocalFaceData() const
    { return *secondaryLocalFaceData_; }

private:

    const Problem& problem() const { return *problemPtr_; }
    const FVElementGeometry& fvGeometry() const { return *fvGeometryPtr_; }
    const ElementVolumeVariables& elemVolVars() const { return *elemVolVarsPtr_; }

    //! Method to fill the flux var caches within an interaction volume
361
    template<class FluxVariablesCache, class FluxVarsCacheStorage, class InteractionVolume>
362
363
    void fillCachesInInteractionVolume_(FluxVarsCacheStorage& fluxVarsCacheStorage,
                                        InteractionVolume& iv,
364
                                        unsigned int ivIndexInContainer)
365
366
367
    {
        // determine if secondary interaction volumes are used here
        static constexpr bool isSecondary = MpfaHelper::considerSecondaryIVs()
368
                                            && std::is_same_v<InteractionVolume, SecondaryInteractionVolume>;
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393

        // 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.
        const auto numGlobalScvfs = iv.localFaceData().size();
        std::vector<const SubControlVolumeFace*> ivScvfs(numGlobalScvfs);
        std::vector<FluxVariablesCache*> ivFluxVarCaches(numGlobalScvfs);

        unsigned int i = 0;
        for (const auto& d : iv.localFaceData())
        {
            // obtain the scvf
            const auto& scvfJ = fvGeometry().scvf(d.gridScvfIndex());
            ivScvfs[i] = &scvfJ;
            ivFluxVarCaches[i] = &fluxVarsCacheStorage[scvfJ];
            ivFluxVarCaches[i]->setIvIndexInContainer(ivIndexInContainer);
            ivFluxVarCaches[i]->setUpdateStatus(true);
            ivFluxVarCaches[i]->setSecondaryIvUsage(isSecondary);
            ivFluxVarCaches[i]->setIvLocalFaceIndex(d.ivLocalScvfIndex());
            if (dim < dimWorld)
                if (d.isOutsideFace())
                    ivFluxVarCaches[i]->setIndexInOutsideFaces(d.scvfLocalOutsideScvfIndex());
            i++;
        }

394
        if constexpr (advectionEnabled)
395
            fillAdvection_(iv, ivScvfs, ivFluxVarCaches);
396
        if constexpr (diffusionEnabled)
397
            fillDiffusion_(iv, ivScvfs, ivFluxVarCaches);
398
        if constexpr (heatConductionEnabled)
399
            fillHeatConduction_(iv, ivScvfs, ivFluxVarCaches);
400
401
402
    }

    //! fills the advective quantities (enabled advection)
403
    template<class InteractionVolume, class FluxVariablesCache>
404
405
    void fillAdvection_(InteractionVolume& iv,
                        const std::vector<const SubControlVolumeFace*>& ivScvfs,
406
                        const std::vector<FluxVariablesCache*>& ivFluxVarCaches)
407
408
409
410
411
412
413
414
415
    {
        using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
        using AdvectionFiller = typename AdvectionType::Cache::Filler;

        // fill advection caches
        for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
        {
            // set pointer to current local face data object
            // ifs are evaluated at compile time and are optimized away
416
            if (std::is_same_v<PrimaryInteractionVolume, SecondaryInteractionVolume>)
417
418
419
420
421
            {
                // we cannot make a disctinction, thus we set both pointers
                primaryLocalFaceData_ = &(iv.localFaceData()[i]);
                secondaryLocalFaceData_ = &(iv.localFaceData()[i]);
            }
422
            else if (std::is_same_v<InteractionVolume, PrimaryInteractionVolume>)
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
                primaryLocalFaceData_ = &(iv.localFaceData()[i]);
            else
                secondaryLocalFaceData_ = &(iv.localFaceData()[i]);

            // fill this scvfs cache
            AdvectionFiller::fill(*ivFluxVarCaches[i],
                                  problem(),
                                  iv.element(iv.localFaceData()[i].ivLocalInsideScvIndex()),
                                  fvGeometry(),
                                  elemVolVars(),
                                  *ivScvfs[i],
                                  *this);
        }
    }

    //! fills the diffusive quantities (diffusion enabled)
439
    template<class InteractionVolume, class FluxVariablesCache>
440
441
    void fillDiffusion_(InteractionVolume& iv,
                        const std::vector<const SubControlVolumeFace*>& ivScvfs,
442
                        const std::vector<FluxVariablesCache*>& ivFluxVarCaches)
443
444
445
446
447
448
449
450
451
452
453
454
    {
        using DiffusionType = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
        using DiffusionFiller = typename DiffusionType::Cache::Filler;

        static constexpr int numPhases = ModelTraits::numFluidPhases();
        static constexpr int numComponents = ModelTraits::numFluidComponents();

        for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
        {
            for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
            {
                using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
455
456
457
                if constexpr (!FluidSystem::isTracerFluidSystem())
                    if (compIdx == FluidSystem::getMainComponent(phaseIdx))
                        continue;
458
459
460
461
462
463

                // fill diffusion caches
                for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
                {
                    // set pointer to current local face data object
                    // ifs are evaluated at compile time and are optimized away
464
                    if constexpr (std::is_same_v<PrimaryInteractionVolume, SecondaryInteractionVolume>)
465
466
467
468
469
                    {
                        // we cannot make a disctinction, thus we set both pointers
                        primaryLocalFaceData_ = &(iv.localFaceData()[i]);
                        secondaryLocalFaceData_ = &(iv.localFaceData()[i]);
                    }
470
                    else if constexpr (std::is_same_v<InteractionVolume, PrimaryInteractionVolume>)
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
                        primaryLocalFaceData_ = &(iv.localFaceData()[i]);
                    else
                        secondaryLocalFaceData_ = &(iv.localFaceData()[i]);

                    // fill this scvfs cache
                    DiffusionFiller::fill(*ivFluxVarCaches[i],
                                          phaseIdx,
                                          compIdx,
                                          problem(),
                                          iv.element(iv.localFaceData()[i].ivLocalInsideScvIndex()),
                                          fvGeometry(),
                                          elemVolVars(),
                                          *ivScvfs[i],
                                          *this);
                }
            }
        }
    }

    //! fills the quantities related to heat conduction (heat conduction enabled)
491
    template<class InteractionVolume, class FluxVariablesCache>
492
493
    void fillHeatConduction_(InteractionVolume& iv,
                             const std::vector<const SubControlVolumeFace*>& ivScvfs,
494
                             const std::vector<FluxVariablesCache*>& ivFluxVarCaches)
495
496
497
498
499
500
501
502
503
    {
        using HeatConductionType = GetPropType<TypeTag, Properties::HeatConductionType>;
        using HeatConductionFiller = typename HeatConductionType::Cache::Filler;

        // fill heat conduction caches
        for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
        {
            // set pointer to current local face data object
            // ifs are evaluated at compile time and are optimized away
504
            if constexpr (std::is_same_v<PrimaryInteractionVolume, SecondaryInteractionVolume>)
505
506
507
508
509
            {
                // we cannot make a disctinction, thus we set both pointers
                primaryLocalFaceData_ = &(iv.localFaceData()[i]);
                secondaryLocalFaceData_ = &(iv.localFaceData()[i]);
            }
510
            else if constexpr (std::is_same_v<InteractionVolume, PrimaryInteractionVolume>)
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
                primaryLocalFaceData_ = &(iv.localFaceData()[i]);
            else
                secondaryLocalFaceData_ = &(iv.localFaceData()[i]);

            // fill this scvfs cache
            HeatConductionFiller::fill(*ivFluxVarCaches[i],
                                       problem(),
                                       iv.element(iv.localFaceData()[i].ivLocalInsideScvIndex()),
                                       fvGeometry(),
                                       elemVolVars(),
                                       *ivScvfs[i],
                                       *this);
        }
    }

526
527
    //! Solves the local systems and stores the result in the handles
    template< class InteractionVolume, class DataHandle>
528
    void prepareDataHandle_([[maybe_unused]] InteractionVolume& iv, [[maybe_unused]] DataHandle& handle, [[maybe_unused]] bool forceUpdate)
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
    {
        // (maybe) solve system subject to intrinsic permeability
        if constexpr (advectionEnabled)
        {
            using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
            if constexpr (AdvectionType::discMethod == DiscretizationMethod::ccmpfa)
                prepareAdvectionHandle_(iv, handle, forceUpdate);
        }

        // (maybe) solve system subject to diffusion tensors
        if constexpr (diffusionEnabled)
        {
            using DiffusionType = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
            if constexpr (DiffusionType::discMethod == DiscretizationMethod::ccmpfa)
                prepareDiffusionHandles_(iv, handle, forceUpdate);
        }

        // (maybe) solve system subject to thermal conductivity
        if constexpr (heatConductionEnabled)
        {
            using HeatConductionType = GetPropType<TypeTag, Properties::HeatConductionType>;
            if constexpr (HeatConductionType::discMethod == DiscretizationMethod::ccmpfa)
                prepareHeatConductionHandle_(iv, handle, forceUpdate);
        }
    }

555
    //! prepares the quantities necessary for advective fluxes in the handle
556
    template<class InteractionVolume, class DataHandle>
557
    void prepareAdvectionHandle_(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll)
558
559
560
561
562
563
    {
        // get instance of the interaction volume-local assembler
        using Traits = typename InteractionVolume::Traits;
        using IvLocalAssembler = typename Traits::template LocalAssembler<Problem, FVElementGeometry, ElementVolumeVariables>;
        IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());

564
565
566
        // lambda to obtain the permeability tensor
        auto getK = [] (const auto& volVars) { return volVars.permeability(); };

567
        // Assemble T only if permeability is sol-dependent or if update is forced
568
        if (forceUpdateAll || advectionIsSolDependent)
569
            localAssembler.assembleMatrices(handle.advectionHandle(), iv, getK);
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589

        // assemble pressure vectors
        for (unsigned int pIdx = 0; pIdx < ModelTraits::numFluidPhases(); ++pIdx)
        {
            // set context in handle
            handle.advectionHandle().setPhaseIndex(pIdx);

            // maybe (re-)assemble gravity contribution vector
            auto getRho = [pIdx] (const auto& volVars) { return volVars.density(pIdx); };
            static const bool enableGravity = getParamFromGroup<bool>(problem().paramGroup(), "Problem.EnableGravity");
            if (enableGravity)
                localAssembler.assembleGravity(handle.advectionHandle(), iv, getRho);

            // reassemble pressure vector
            auto getPressure = [pIdx] (const auto& volVars) { return volVars.pressure(pIdx); };
            localAssembler.assembleU(handle.advectionHandle(), iv, getPressure);
        }
    }

    //! prepares the quantities necessary for diffusive fluxes in the handle
590
    template<class InteractionVolume, class DataHandle>
591
592
593
    void prepareDiffusionHandles_(InteractionVolume& iv,
                                  DataHandle& handle,
                                  bool forceUpdateAll)
594
    {
595
596
597
598
599
600
        for (unsigned int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
        {
            for (unsigned int compIdx = 0; compIdx < ModelTraits::numFluidComponents(); ++compIdx)
            {
                // skip main component
                using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
601
602
603
                if constexpr (!FluidSystem::isTracerFluidSystem())
                    if (compIdx == FluidSystem::getMainComponent(phaseIdx))
                        continue;
604

605
606
607
608
609
610
611
612
613
614
                // fill data in the handle
                handle.diffusionHandle().setPhaseIndex(phaseIdx);
                handle.diffusionHandle().setComponentIndex(compIdx);

                using DiffusionType = GetPropType<TypeTag, Properties::MolecularDiffusionType>;

                // get instance of the interaction volume-local assembler
                using Traits = typename InteractionVolume::Traits;
                using IvLocalAssembler = typename Traits::template LocalAssembler<Problem, FVElementGeometry, ElementVolumeVariables>;
                IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
615

616
617
                // maybe (re-)assemble matrices
                if (forceUpdateAll || diffusionIsSolDependent)
618
                {
619
                    // lambda to obtain diffusion coefficient
620
                    const auto getD = [&](const auto& volVars)
621
622
                    {
                        if constexpr (FluidSystem::isTracerFluidSystem())
623
                            return volVars.effectiveDiffusionCoefficient(0, 0, compIdx);
624
                        else
625
                            return volVars.effectiveDiffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
626
627
                    };

628
                    // Effective diffusion coefficients might be zero if saturation = 0.
629
                    // Compute epsilon to detect obsolete rows in the iv-local matrices during assembly
630
631
632
633
634
635
636
637
638
639
                    static const auto zeroD = getParamFromGroup<Scalar>(
                        problem().paramGroup(),
                        "Mpfa.ZeroEffectiveDiffusionCoefficientThreshold",
                        1e-16
                    );

                    // compute a representative transmissibility for this interaction volume, using
                    // the threshold for zero diffusion coefficients, and use this as epsilon
                    const auto& scv = fvGeometry().scv(iv.localScv(0).gridScvIndex());
                    const auto& scvf = fvGeometry().scvf(iv.localScvf(0).gridScvfIndex());
640
                    const auto& vv = elemVolVars()[scv];
641
642
643
                    const auto eps = Extrusion::area(scvf)*computeTpfaTransmissibility(
                        scvf, scv, zeroD, vv.extrusionFactor()
                    );
644
645

                    localAssembler.assembleMatrices(handle.diffusionHandle(), iv, getD, eps);
646
                }
647

648
                // assemble vector of mole fractions
649
650
651
652
653
654
                auto getMassOrMoleFraction = [phaseIdx, compIdx] (const auto& volVars)
                {
                    return (DiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged) ? volVars.massFraction(phaseIdx, compIdx) :
                                                                                                                       volVars.moleFraction(phaseIdx, compIdx);
                };

655
656
657
                localAssembler.assembleU(handle.diffusionHandle(), iv, getMassOrMoleFraction);
            }
        }
658
659
660
    }

    //! prepares the quantities necessary for conductive fluxes in the handle
661
    template<class InteractionVolume, class DataHandle>
662
    void prepareHeatConductionHandle_(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll)
663
664
665
666
667
668
    {
        // get instance of the interaction volume-local assembler
        using Traits = typename InteractionVolume::Traits;
        using IvLocalAssembler = typename Traits::template LocalAssembler<Problem, FVElementGeometry, ElementVolumeVariables>;
        IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());

669
670
671
        // lambda to obtain the effective thermal conductivity
        auto getLambda = [] (const auto& volVars) { return volVars.effectiveThermalConductivity(); };

672
        // maybe (re-)assemble matrices
673
        if (forceUpdateAll || heatConductionIsSolDependent)
674
            localAssembler.assembleMatrices(handle.heatConductionHandle(), iv, getLambda);
675
676

        // assemble vector of temperatures
677
678
        auto getTemperature = [] (const auto& volVars) { return volVars.temperature(); };
        localAssembler.assembleU(handle.heatConductionHandle(), iv, getTemperature);
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
    }

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

    // We store pointers to an inner and a boundary interaction volume.
    // These are updated during the filling of the caches and the
    // physics-related caches have access to them
    PrimaryInteractionVolume* primaryIv_;
    SecondaryInteractionVolume* secondaryIv_;

    // pointer to the current interaction volume data handle
    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_;
};

} // end namespace Dumux

#endif