Skip to content
Snippets Groups Projects
Commit d921fcb7 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[mpfa][fickslaw] store pointer to handle 

instead of storing pointers to the transmissibilities and pressures
etc, we now only store a pointer to the data handle and reconstruct
the fluxes in the FicksLaw implementation. This dramatically reduces
the number of lines needed for the cache implementation.
parent d8a1a2c9
No related branches found
No related tags found
2 merge requests!1337WIP Fix/dirichlet caching v2,!1298Restructure mpfa flux caches
Hide whitespace changes
Inline Side-by-side
dumux/discretization/cellcentered/mpfa/fickslaw.hh
+ 67
− 116
View file @ d921fcb7
......@@ -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,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment