diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
index 7d13e5c514271edd32753af26be2437caaf6e5cc..3d5e759d4d95babc40e4bd7a5de687110319cde0 100644
--- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
@@ -46,6 +46,9 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::ccmpfa>
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
@@ -79,12 +82,12 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::ccmpfa>
scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.secondaryInteractionVolume(),
fluxVarsCacheFiller.secondaryIvLocalFaceData(),
fluxVarsCacheFiller.secondaryIvDataHandle(),
- scvf, phaseIdx, compIdx);
+ phaseIdx, compIdx);
else
scvfFluxVarsCache.updateDiffusion(fluxVarsCacheFiller.primaryInteractionVolume(),
fluxVarsCacheFiller.primaryIvLocalFaceData(),
fluxVarsCacheFiller.primaryIvDataHandle(),
- scvf, phaseIdx, compIdx);
+ phaseIdx, compIdx);
}
};
@@ -94,24 +97,19 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::ccmpfa>
using DualGridNodalIndexSet = GetPropType<TypeTag, Properties::DualGridNodalIndexSet>;
using Stencil = typename DualGridNodalIndexSet::NodalGridStencilType;
- using MpfaHelper = typename FVGridGeometry::MpfaHelper;
- static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs();
-
- using PrimaryInteractionVolume = GetPropType<TypeTag, Properties::PrimaryInteractionVolume>;
- using PrimaryIvLocalFaceData = typename PrimaryInteractionVolume::Traits::LocalFaceData;
- using PrimaryIvDataHandle = typename ElementFluxVariablesCache::PrimaryIvDataHandle;
- using PrimaryIvCellVector = typename PrimaryInteractionVolume::Traits::MatVecTraits::CellVector;
- using PrimaryIvTij = typename PrimaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type;
+ static constexpr int numPhases = GetPropType<TypeTag, Properties::ModelTraits>::numPhases();
+ static constexpr bool considerSecondaryIVs = FVGridGeometry::MpfaHelper::considerSecondaryIVs();
+ using PrimaryDataHandle = typename ElementFluxVariablesCache::PrimaryIvDataHandle::DiffusionHandle;
+ using SecondaryDataHandle = typename ElementFluxVariablesCache::SecondaryIvDataHandle::DiffusionHandle;
- using SecondaryInteractionVolume = GetPropType<TypeTag, Properties::SecondaryInteractionVolume>;
- using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
- using SecondaryIvDataHandle = typename ElementFluxVariablesCache::SecondaryIvDataHandle;
- using SecondaryIvCellVector = typename SecondaryInteractionVolume::Traits::MatVecTraits::CellVector;
- using SecondaryIvTij = typename SecondaryInteractionVolume::Traits::MatVecTraits::TMatrix::row_type;
+ //! sets the pointer to the data handle (overload for secondary data handles)
+ template< bool doSecondary = considerSecondaryIVs, std::enable_if_t<doSecondary, int> = 0 >
+ void setHandlePointer_(const SecondaryDataHandle& dataHandle)
+ { secondaryHandlePtr_ = &dataHandle; }
- static constexpr int dim = GridView::dimension;
- static constexpr int dimWorld = GridView::dimensionworld;
- static constexpr int numPhases = GetPropType<TypeTag, Properties::ModelTraits>::numPhases();
+ //! sets the pointer to the data handle (overload for primary data handles)
+ void setHandlePointer_(const PrimaryDataHandle& dataHandle)
+ { primaryHandlePtr_ = &dataHandle; }
public:
// export filler type
@@ -124,104 +122,39 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::ccmpfa>
* \param iv The interaction volume this scvf is embedded in
* \param localFaceData iv-local info on this scvf
* \param dataHandle Transmissibility matrix & gravity data of this iv
- * \param scvf The sub-control volume face
*/
- void updateDiffusion(const PrimaryInteractionVolume& iv,
- const PrimaryIvLocalFaceData& localFaceData,
- const PrimaryIvDataHandle& dataHandle,
- const SubControlVolumeFace &scvf,
+ template<class IV, class LocalFaceData, class DataHandle>
+ void updateDiffusion(const IV& iv,
+ const LocalFaceData& localFaceData,
+ const DataHandle& dataHandle,
unsigned int phaseIdx, unsigned int compIdx)
{
- const auto& handle = dataHandle.diffusionHandle();
-
stencil_[phaseIdx][compIdx] = &iv.stencil();
switchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutsideFace();
-
- // store pointer to the mole fraction vector of this iv
- primaryXj_[phaseIdx][compIdx] = &handle.uj();
-
- const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
- if (dim == dimWorld)
- primaryTij_[phaseIdx][compIdx] = &handle.T()[ivLocalIdx];
- else
- primaryTij_[phaseIdx][compIdx] = localFaceData.isOutsideFace()
- ? &handle.tijOutside()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()]
- : &handle.T()[ivLocalIdx];
+ setHandlePointer_(dataHandle.diffusionHandle());
}
- /*!
- * \brief Update cached objects (transmissibilities).
- * This is used for updates with secondary interaction volumes.
- *
- * \param iv The interaction volume this scvf is embedded in
- * \param localFaceData iv-local info on this scvf
- * \param dataHandle Transmissibility matrix & gravity data of this iv
- * \param scvf The sub-control volume face
- */
- template< bool doSecondary = considerSecondaryIVs, std::enable_if_t<doSecondary, int > = 0 >
- void updateDiffusion(const SecondaryInteractionVolume& iv,
- const SecondaryIvLocalFaceData& localFaceData,
- const SecondaryIvDataHandle& dataHandle,
- const SubControlVolumeFace &scvf,
- unsigned int phaseIdx, unsigned int compIdx)
- {
- const auto& handle = dataHandle.diffusionHandle();
-
- stencil_[phaseIdx][compIdx] = &iv.stencil();
- switchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutsideFace();
-
- // store pointer to the mole fraction vector of this iv
- secondaryXj_[phaseIdx][compIdx] = &handle.uj();
+ //! The stencils corresponding to the transmissibilities
+ const Stencil& diffusionStencil(unsigned int phaseIdx, unsigned int compIdx) const
+ { return *stencil_[phaseIdx][compIdx]; }
- const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
- if (dim == dimWorld)
- secondaryTij_[phaseIdx][compIdx] = &handle.T()[ivLocalIdx];
- else
- secondaryTij_[phaseIdx][compIdx] = localFaceData.isOutsideFace()
- ? &handle.tijOutside()[ivLocalIdx][localFaceData.scvfLocalOutsideScvfIndex()]
- : &handle.T()[ivLocalIdx];
- }
+ //! The corresponding data handles
+ const PrimaryDataHandle& diffusionPrimaryDataHandle() const { return *primaryHandlePtr_; }
+ const SecondaryDataHandle& diffusionSecondaryDataHandle() const { return *secondaryHandlePtr_; }
- //! In the interaction volume-local system of eq we have one unknown per face.
- //! On scvfs on this face, but in "outside" (neighbor) elements of it, we have
- //! to take the negative value of the fluxes due to the flipped normal vector.
- //! This function returns whether or not this scvf is an "outside" face in the iv.
+ //! Returns whether or not this scvf is an "outside" face in the scope of the iv.
bool diffusionSwitchFluxSign(unsigned int phaseIdx, unsigned int compIdx) const
{ return switchFluxSign_[phaseIdx][compIdx]; }
- //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (primary type)
- const PrimaryIvTij& diffusionTijPrimaryIv(unsigned int phaseIdx, unsigned int compIdx) const
- { return *primaryTij_[phaseIdx][compIdx]; }
-
- //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (secondary type)
- const SecondaryIvTij& diffusionTijSecondaryIv(unsigned int phaseIdx, unsigned int compIdx) const
- { return *secondaryTij_[phaseIdx][compIdx]; }
-
- //! The cell (& maybe Dirichlet) mole fractions within this interaction volume (primary type)
- const PrimaryIvCellVector& moleFractionsPrimaryIv(unsigned int phaseIdx, unsigned int compIdx) const
- { return *primaryXj_[phaseIdx][compIdx]; }
-
- //! The cell (& maybe Dirichlet) mole fractions within this interaction volume (secondary type)
- const SecondaryIvCellVector& moleFractionsSecondaryIv(unsigned int phaseIdx, unsigned int compIdx) const
- { return *secondaryXj_[phaseIdx][compIdx]; }
-
- //! The stencils corresponding to the transmissibilities
- const Stencil& diffusionStencil(unsigned int phaseIdx, unsigned int compIdx) const
- { return *stencil_[phaseIdx][compIdx]; }
private:
+ //! phase-/component- specific data
std::array< std::array<bool, numComponents>, numPhases > switchFluxSign_;
-
- //! The stencils, i.e. the grid indices j
std::array< std::array<const Stencil*, numComponents>, numPhases > stencil_;
- //! The transmissibilities such that f = Tij*xj
- std::array< std::array<const PrimaryIvCellVector*, numComponents>, numPhases > primaryTij_;
- std::array< std::array<const SecondaryIvCellVector*, numComponents>, numPhases > secondaryTij_;
-
- //! The interaction-volume wide mole fractions xj
- std::array< std::array<const PrimaryIvCellVector*, numComponents>, numPhases > primaryXj_;
- std::array< std::array<const SecondaryIvCellVector*, numComponents>, numPhases > secondaryXj_;
+ //! pointers to the corresponding iv-data handles
+ const PrimaryDataHandle* primaryHandlePtr_;
+ const SecondaryDataHandle* secondaryHandlePtr_;
};
public:
@@ -259,25 +192,17 @@ public:
// calculate the density at the interface
const auto rho = interpolateDensity(elemVolVars, scvf, phaseIdx);
- // calculate Tij*xj
- Scalar flux;
+ // compute the flux
if (fluxVarsCache.usesSecondaryIv())
- {
- const auto& tij = fluxVarsCache.diffusionTijSecondaryIv(phaseIdx, compIdx);
- const auto& xj = fluxVarsCache.moleFractionsSecondaryIv(phaseIdx, compIdx);
- flux = tij*xj;
- }
+ componentFlux[compIdx] = rho*effFactor*computeVolumeFlux(problem,
+ fluxVarsCache,
+ fluxVarsCache.diffusionSecondaryDataHandle(),
+ phaseIdx, compIdx);
else
- {
- const auto& tij = fluxVarsCache.diffusionTijPrimaryIv(phaseIdx, compIdx);
- const auto& xj = fluxVarsCache.moleFractionsPrimaryIv(phaseIdx, compIdx);
- flux = tij*xj;
- }
-
- if (fluxVarsCache.diffusionSwitchFluxSign(phaseIdx, compIdx))
- flux *= -1.0;
-
- componentFlux[compIdx] = flux*rho*effFactor;
+ componentFlux[compIdx] = rho*effFactor*computeVolumeFlux(problem,
+ fluxVarsCache,
+ fluxVarsCache.diffusionPrimaryDataHandle(),
+ phaseIdx, compIdx);
}
// accumulate the phase component flux
@@ -289,6 +214,32 @@ public:
}
private:
+ template< class Problem, class FluxVarsCache, class DataHandle >
+ static Scalar computeVolumeFlux(const Problem& problem,
+ const FluxVarsCache& cache,
+ const DataHandle& dataHandle,
+ int phaseIdx, int compIdx)
+ {
+ dataHandle.setPhaseIndex(phaseIdx);
+ dataHandle.setComponentIndex(compIdx);
+
+ const bool switchSign = cache.diffusionSwitchFluxSign(phaseIdx, compIdx);
+
+ const auto localFaceIdx = cache.ivLocalFaceIndex();
+ const auto idxInOutside = cache.indexInOutsideFaces();
+ const auto& xj = dataHandle.uj();
+ const auto& tij = dim == dimWorld ? dataHandle.T()[localFaceIdx]
+ : (!switchSign ? dataHandle.T()[localFaceIdx]
+ : dataHandle.tijOutside()[localFaceIdx][idxInOutside]);
+ Scalar scvfFlux = tij*xj;
+
+ // switch the sign if necessary
+ if (cache.advectionSwitchFluxSign())
+ scvfFlux *= -1.0;
+
+ return scvfFlux;
+ }
+
//! compute the density at branching facets for network grids as arithmetic mean
static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,