From 077c544eb3189516441a7615f1010c018ed36008 Mon Sep 17 00:00:00 2001
From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de>
Date: Tue, 30 Jan 2018 15:22:12 +0100
Subject: [PATCH] [mpfa][fluxcachefiller] outsource data handle preparation
This introduces another set of functions specialized via enable_if,
but this is necessary as otherwise compilation fails for the tracer model,
because it is tried to access pressures via the vol vars, which the volvars
of the tracer model do not provide. By ensuring that this is only called
for processes using mpfa, we are on the safe side.
---
.../mpfa/fluxvariablescachefiller.hh | 326 ++++++++++--------
1 file changed, 184 insertions(+), 142 deletions(-)
diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
index 5a90cc053a..b0ced7e512 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_;
--
GitLab