diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
index 5a90cc053a209e2b0dd5f93f714733c626088d04..b0ced7e5123746dafb0f6d0b5776cd8f5a5eafd2 100644
--- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
+++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
@@ -252,84 +252,8 @@ private:
using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
using AdvectionFiller = typename AdvectionType::Cache::Filler;
- static constexpr auto AdvectionMethod = AdvectionType::myDiscretizationMethod;
- using LambdaFactory = TensorLambdaFactory<TypeTag, AdvectionMethod>;
-
- // skip the following if advection doesn't use mpfa
- if (AdvectionMethod == DiscretizationMethods::CCMpfa)
- {
- // get instance of the interaction volume-local assembler
- using IVTraits = typename InteractionVolume::Traits;
- using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
- 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);
- }
- }
+ // fill data in the handle
+ fillAdvectionHandle(iv, handle, forceUpdateAll);
// fill advection caches
for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
@@ -384,17 +308,9 @@ private:
using DiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
using DiffusionFiller = typename DiffusionType::Cache::Filler;
- static constexpr auto DiffusionMethod = DiffusionType::myDiscretizationMethod;
- using LambdaFactory = TensorLambdaFactory<TypeTag, DiffusionMethod>;
-
static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
- // get instance of the interaction volume-local assembler
- using IVTraits = typename InteractionVolume::Traits;
- using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
- IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
-
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
@@ -402,33 +318,8 @@ private:
if (phaseIdx == compIdx)
continue;
- // solve the local system subject to the diffusion tensor (if uses mpfa)
- if (DiffusionMethod == DiscretizationMethods::CCMpfa)
- {
- // 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);
- }
+ // fill data in the handle
+ fillDiffusionHandle(iv, handle, forceUpdateAll, phaseIdx, compIdx);
// fill diffusion caches
for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
@@ -487,35 +378,8 @@ private:
using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
using HeatConductionFiller = typename HeatConductionType::Cache::Filler;
- static constexpr auto HeatConductionMethod = HeatConductionType::myDiscretizationMethod;
- using LambdaFactory = TensorLambdaFactory<TypeTag, HeatConductionMethod>;
-
- // maybe solve the local system subject to fourier coefficient
- if (HeatConductionMethod == DiscretizationMethods::CCMpfa)
- {
- // get instance of the interaction volume-local assembler
- using IVTraits = typename InteractionVolume::Traits;
- using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
- 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);
- }
+ // prepare data in handle
+ fillHeatConductionHandle(iv, handle, forceUpdateAll);
// fill heat conduction caches
for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
@@ -556,6 +420,184 @@ private:
bool forceUpdateAll = false)
{}
+ //! prepares the quantities necessary for advective fluxes in the handle
+ template< class InteractionVolume,
+ class DataHandle,
+ class AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType),
+ typename std::enable_if_t<AdvectionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillAdvectionHandle(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll)
+ {
+ using LambdaFactory = TensorLambdaFactory<TypeTag, DiscretizationMethods::CCMpfa>;
+
+ // get instance of the interaction volume-local assembler
+ using IVTraits = typename InteractionVolume::Traits;
+ using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
+ 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);
+ }
+ }
+
+ //! prepares the quantities necessary for diffusive fluxes in the handle
+ template< class InteractionVolume,
+ class DataHandle,
+ class DiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType),
+ typename std::enable_if_t<DiffusionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillDiffusionHandle(InteractionVolume& iv,
+ DataHandle& handle,
+ bool forceUpdateAll,
+ int phaseIdx, int compIdx)
+ {
+ using LambdaFactory = TensorLambdaFactory<TypeTag, DiscretizationMethods::CCMpfa>;
+
+ // get instance of the interaction volume-local assembler
+ using IVTraits = typename InteractionVolume::Traits;
+ using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
+ IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
+
+ // solve the local system subject to the tensor and update the 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);
+ }
+
+ //! prepares the quantities necessary for conductive fluxes in the handle
+ template< class InteractionVolume,
+ class DataHandle,
+ class HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType),
+ typename std::enable_if_t<HeatConductionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillHeatConductionHandle(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll)
+ {
+ using LambdaFactory = TensorLambdaFactory<TypeTag, DiscretizationMethods::CCMpfa>;
+
+ // get instance of the interaction volume-local assembler
+ using IVTraits = typename InteractionVolume::Traits;
+ using IvLocalAssembler = InteractionVolumeAssembler< IVTraits, InteractionVolume::MpfaMethod >;
+ 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);
+ }
+
+ //! fill handle only when advection uses mpfa
+ template< class InteractionVolume,
+ class DataHandle,
+ class AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType),
+ typename std::enable_if_t<AdvectionType::myDiscretizationMethod != DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillAdvectionHandle(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll) {}
+
+ //! fill handle only when diffusion uses mpfa
+ template< class InteractionVolume,
+ class DataHandle,
+ class DiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType),
+ typename std::enable_if_t<DiffusionType::myDiscretizationMethod != DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillDiffusionHandle(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll, int phaseIdx, int compIdx) {}
+
+ //! fill handle only when heat conduction uses mpfa
+ template< class InteractionVolume,
+ class DataHandle,
+ class HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType),
+ typename std::enable_if_t<HeatConductionType::myDiscretizationMethod != DiscretizationMethods::CCMpfa, int> = 0 >
+ void fillHeatConductionHandle(InteractionVolume& iv, DataHandle& handle, bool forceUpdateAll) {}
+
const Problem* problemPtr_;
const Element* elementPtr_;
const FVElementGeometry* fvGeometryPtr_;