Commit d8a1a2c9 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[mpfa][darcyslaw] 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 DarcysLaw implementation. This dramatically reduces
the number of lines needed for the cache implementation.
parent 0b4b7410
......@@ -50,6 +50,9 @@ class DarcysLawImplementation<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 SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
......@@ -74,44 +77,34 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethod::ccmpfa>
{
// get interaction volume related data from the filler class & upate the cache
if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
scvfFluxVarsCache.updateAdvection(problem,
fluxVarsCacheFiller.secondaryInteractionVolume(),
scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.secondaryInteractionVolume(),
fluxVarsCacheFiller.secondaryIvLocalFaceData(),
fluxVarsCacheFiller.secondaryIvDataHandle(),
scvf);
fluxVarsCacheFiller.secondaryIvDataHandle());
else
scvfFluxVarsCache.updateAdvection(problem,
fluxVarsCacheFiller.primaryInteractionVolume(),
scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.primaryInteractionVolume(),
fluxVarsCacheFiller.primaryIvLocalFaceData(),
fluxVarsCacheFiller.primaryIvDataHandle(),
scvf);
fluxVarsCacheFiller.primaryIvDataHandle());
}
};
//! The cache used in conjunction with the mpfa Darcy's Law
class MpfaDarcysLawCache
{
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
static constexpr int numPhases = GetPropType<TypeTag, Properties::ModelTraits>::numPhases();
using DualGridNodalIndexSet = GetPropType<TypeTag, Properties::DualGridNodalIndexSet>;
using Stencil = typename DualGridNodalIndexSet::NodalGridStencilType;
using MpfaHelper = typename FVGridGeometry::MpfaHelper;
static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs();
static constexpr bool considerSecondaryIVs = FVGridGeometry::MpfaHelper::considerSecondaryIVs();
using PrimaryDataHandle = typename ElementFluxVariablesCache::PrimaryIvDataHandle::AdvectionHandle;
using SecondaryDataHandle = typename ElementFluxVariablesCache::SecondaryIvDataHandle::AdvectionHandle;
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;
//! 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; }
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 primary data handles)
void setHandlePointer_(const PrimaryDataHandle& dataHandle)
{ primaryHandlePtr_ = &dataHandle; }
public:
// export the filler type
......@@ -121,168 +114,39 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethod::ccmpfa>
* \brief Update cached objects (transmissibilities and gravity).
* This is used for updates with primary interaction volumes.
*
* \param problem The problem
* \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<class Problem>
void updateAdvection(const Problem& problem,
const PrimaryInteractionVolume& iv,
const PrimaryIvLocalFaceData& localFaceData,
const PrimaryIvDataHandle& dataHandle,
const SubControlVolumeFace &scvf)
template<class IV, class LocalFaceData, class DataHandle>
void updateAdvection(const IV& iv,
const LocalFaceData& localFaceData,
const DataHandle& dataHandle)
{
const auto& handle = dataHandle.advectionHandle();
switchFluxSign_ = localFaceData.isOutsideFace();
stencil_ = &iv.stencil();
static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
// standard grids
if (dim == dimWorld)
{
primaryTij_ = &handle.T()[ivLocalIdx];
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
handle.setPhaseIndex(phaseIdx);
primaryPj_[phaseIdx] = &handle.uj();
if (enableGravity) g_[phaseIdx] = handle.g()[ivLocalIdx];
}
}
// surface grids
else
{
if (!localFaceData.isOutsideFace())
{
primaryTij_ = &handle.T()[ivLocalIdx];
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
handle.setPhaseIndex(phaseIdx);
primaryPj_[phaseIdx] = &handle.uj();
if (enableGravity) g_[phaseIdx] = handle.g()[ivLocalIdx];
}
}
else
{
const auto idxInOutsideFaces = localFaceData.scvfLocalOutsideScvfIndex();
primaryTij_ = &handle.tijOutside()[ivLocalIdx][idxInOutsideFaces];
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
handle.setPhaseIndex(phaseIdx);
primaryPj_[phaseIdx] = &handle.uj();
if (enableGravity) g_[phaseIdx] = handle.gOutside()[ivLocalIdx][idxInOutsideFaces];
}
}
}
}
/*!
* \brief Update cached objects (transmissibilities and gravity).
* This is used for updates with secondary interaction volumes.
*
* \param problem The problem
* \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<class Problem, bool doSecondary = considerSecondaryIVs, std::enable_if_t<doSecondary, int > = 0 >
void updateAdvection(const Problem& problem,
const SecondaryInteractionVolume& iv,
const SecondaryIvLocalFaceData& localFaceData,
const SecondaryIvDataHandle& dataHandle,
const SubControlVolumeFace &scvf)
{
const auto& handle = dataHandle.advectionHandle();
switchFluxSign_ = localFaceData.isOutsideFace();
stencil_ = &iv.stencil();
static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
// standard grids
if (dim == dimWorld)
{
secondaryTij_ = &handle.T()[ivLocalIdx];
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
handle.setPhaseIndex(phaseIdx);
secondaryPj_[phaseIdx] = &handle.uj();
if (enableGravity) g_[phaseIdx] = handle.g()[ivLocalIdx];
}
}
// surface grids
else
{
if (!localFaceData.isOutsideFace())
{
secondaryTij_ = &handle.T()[ivLocalIdx];
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
handle.setPhaseIndex(phaseIdx);
secondaryPj_[phaseIdx] = &handle.uj();
if (enableGravity) g_[phaseIdx] = handle.g()[ivLocalIdx];
}
}
else
{
const auto idxInOutsideFaces = localFaceData.scvfLocalOutsideScvfIndex();
secondaryTij_ = &handle.tijOutside()[ivLocalIdx][idxInOutsideFaces];
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
handle.setPhaseIndex(phaseIdx);
secondaryPj_[phaseIdx] = &handle.uj();
if (enableGravity) g_[phaseIdx] = handle.gOutside()[ivLocalIdx][idxInOutsideFaces];
}
}
}
setHandlePointer_(dataHandle.advectionHandle());
}
//! The stencil corresponding to the transmissibilities (primary type)
const Stencil& advectionStencil() const { return *stencil_; }
//! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (primary type)
const PrimaryIvTij& advectionTijPrimaryIv() const { return *primaryTij_; }
//! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (secondary type)
const SecondaryIvTij& advectionTijSecondaryIv() const { return *secondaryTij_; }
//! The cell (& maybe Dirichlet) pressures within this interaction volume (primary type)
const PrimaryIvCellVector& pressuresPrimaryIv(unsigned int phaseIdx) const { return *primaryPj_[phaseIdx]; }
//! The cell (& maybe Dirichlet) pressures within this interaction volume (secondary type)
const SecondaryIvCellVector& pressuresSecondaryIv(unsigned int phaseIdx) const { return *secondaryPj_[phaseIdx]; }
//! The corresponding data handles
const PrimaryDataHandle& advectionPrimaryDataHandle() const { return *primaryHandlePtr_; }
const SecondaryDataHandle& advectionSecondaryDataHandle() const { return *secondaryHandlePtr_; }
//! The gravitational acceleration for a phase on this scvf
Scalar gravity(unsigned int phaseIdx) const { return g_[phaseIdx]; }
//! 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 advectionSwitchFluxSign() const { return switchFluxSign_; }
private:
bool switchFluxSign_;
//! pointers to the corresponding iv-data handles
const PrimaryDataHandle* primaryHandlePtr_;
const SecondaryDataHandle* secondaryHandlePtr_;
//! The stencil, i.e. the grid indices j
const Stencil* stencil_;
//! The transmissibilities such that f = Tij*pj
const PrimaryIvTij* primaryTij_;
const SecondaryIvTij* secondaryTij_;
//! The interaction-volume wide phase pressures pj
std::array<const PrimaryIvCellVector*, numPhases> primaryPj_;
std::array<const SecondaryIvCellVector*, numPhases> secondaryPj_;
//! Gravitational flux contribution on this face
std::array< Scalar, numPhases > g_;
};
public:
......@@ -303,28 +167,41 @@ public:
{
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
// compute t_ij*p_j
Scalar scvfFlux;
// forward to the private function taking the iv data handle
if (fluxVarsCache.usesSecondaryIv())
{
const auto& tij = fluxVarsCache.advectionTijSecondaryIv();
const auto& pj = fluxVarsCache.pressuresSecondaryIv(phaseIdx);
scvfFlux = tij*pj;
}
return flux_(problem, fluxVarsCache, fluxVarsCache.advectionSecondaryDataHandle(), phaseIdx);
else
{
const auto& tij = fluxVarsCache.advectionTijPrimaryIv();
const auto& pj = fluxVarsCache.pressuresPrimaryIv(phaseIdx);
scvfFlux = tij*pj;
}
return flux_(problem, fluxVarsCache, fluxVarsCache.advectionPrimaryDataHandle(), phaseIdx);
}
private:
template< class Problem, class FluxVarsCache, class DataHandle >
static Scalar flux_(const Problem& problem,
const FluxVarsCache& cache,
const DataHandle& dataHandle,
int phaseIdx)
{
dataHandle.setPhaseIndex(phaseIdx);
const bool switchSign = cache.advectionSwitchFluxSign();
const auto localFaceIdx = cache.ivLocalFaceIndex();
const auto idxInOutside = cache.indexInOutsideFaces();
const auto& pj = dataHandle.uj();
const auto& tij = dim == dimWorld ? dataHandle.T()[localFaceIdx]
: (!switchSign ? dataHandle.T()[localFaceIdx]
: dataHandle.tijOutside()[localFaceIdx][idxInOutside]);
Scalar scvfFlux = tij*pj;
// maybe add gravitational acceleration
static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
if (enableGravity)
scvfFlux += fluxVarsCache.gravity(phaseIdx);
scvfFlux += dim == dimWorld ? dataHandle.g()[localFaceIdx]
: (!switchSign ? dataHandle.g()[localFaceIdx]
: dataHandle.gOutside()[localFaceIdx][idxInOutside]);
// switch the sign if necessary
if (fluxVarsCache.advectionSwitchFluxSign())
if (cache.advectionSwitchFluxSign())
scvfFlux *= -1.0;
return scvfFlux;
......
......@@ -231,6 +231,9 @@ private:
ivFluxVarCaches[i]->setIvIndexInContainer(ivIndexInContainer);
ivFluxVarCaches[i]->setUpdateStatus(true);
ivFluxVarCaches[i]->setSecondaryIvUsage(isSecondary);
ivFluxVarCaches[i]->setIvLocalFaceIndex(d.ivLocalScvfIndex());
if (dim < dimWorld)
ivFluxVarCaches[i]->setIndexInOutsideFaces(d.scvfLocalOutsideScvfIndex());
i++;
}
......
......@@ -220,10 +220,12 @@ class HeatConductionDataHandle<MatVecTraits, PhysicsTraits, false> : public Empt
template<class MVT, class PT>
class InteractionVolumeDataHandle
{
public:
//! export the underlying process-specific handle types
using AdvectionHandle = AdvectionDataHandle<MVT, PT, PT::enableAdvection>;
using DiffusionHandle = DiffusionDataHandle<MVT, PT, PT::enableMolecularDiffusion>;
using HeatConductionHandle = HeatConductionDataHandle<MVT, PT, PT::enableHeatConduction>;
public:
//! return references to the handle containing data related to advection
const AdvectionHandle& advectionHandle() const { return advectionHandle_; }
......
......@@ -126,41 +126,20 @@ public:
/ curElemVolVars[scvf.insideScvIdx()].viscosity();
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
const auto& stencil = fluxVarsCache.advectionStencil();
if (fluxVarsCache.usesSecondaryIv())
{
const auto& tij = fluxVarsCache.advectionTijSecondaryIv();
// We assume same the tij are order as the stencil up to stencil.size()
// any contribution of Dirichlet BCs is assumed to be placed afterwards
assert(stencil.size() <= tij.size());
const auto localFaceIdx = fluxVarsCache.ivLocalFaceIndex();
const auto usesSecondary = fluxVarsCache.usesSecondaryIv();
const auto switchSign = fluxVarsCache.advectionSwitchFluxSign();
// add partial derivatives to the respective given matrices
for (unsigned int i = 0; i < stencil.size();++i)
{
if (fluxVarsCache.advectionSwitchFluxSign())
derivativeMatrices[stencil[i]][conti0EqIdx][pressureIdx] -= tij[i]*up;
else
derivativeMatrices[stencil[i]][conti0EqIdx][pressureIdx] += tij[i]*up;
}
}
else
{
const auto& tij = fluxVarsCache.advectionTijPrimaryIv();
// We assume same the tij are order as the stencil up to stencil.size()
// any contribution of Dirichlet BCs is assumed to be placed afterwards
assert(stencil.size() <= tij.size());
// add partial derivatives to the respective given matrices
for (unsigned int i = 0; i < stencil.size();++i)
{
if (fluxVarsCache.advectionSwitchFluxSign())
derivativeMatrices[stencil[i]][conti0EqIdx][pressureIdx] -= tij[i]*up;
else
derivativeMatrices[stencil[i]][conti0EqIdx][pressureIdx] += tij[i]*up;
}
}
const auto& stencil = fluxVarsCache.advectionStencil();
const auto& tij = usesSecondary ? fluxVarsCache.advectionSecondaryDataHandle().T()[localFaceIdx]
: fluxVarsCache.advectionPrimaryDataHandle().T()[localFaceIdx];
// We assume same the tij are order as the stencil up to stencil.size()
// any contribution of Dirichlet BCs is assumed to be placed afterwards
assert(stencil.size() <= tij.size());
for (unsigned int i = 0; i < stencil.size();++i)
derivativeMatrices[stencil[i]][conti0EqIdx][pressureIdx] += switchSign ? -tij[i]*up
: tij[i]*up;
}
//! flux derivatives for the box scheme
......
......@@ -103,27 +103,32 @@ public:
template< bool doSecondary = considerSecondary, std::enable_if_t<doSecondary, int> = 0>
bool usesSecondaryIv() const { return usesSecondaryIv_; }
//! Returns the index of the iv (this scvf is embedded in) in its container
GridIndexType ivIndexInContainer() const { return ivIndexInContainer_; }
//! Returns interaction volume-local face index
unsigned int ivLocalFaceIndex() const { return ivLocalFaceIdx_; }
//! Returns index of the face among "outside" faces of iv-local "positive" face
unsigned int indexInOutsideFaces() const { return idxInOutsideFaces_; }
//! Sets the update status. When set to true, consecutive updates will be skipped
void setUpdateStatus(bool status) { isUpdated_ = status; }
//! Sets if this cache is associated witha secondary iv
//! Sets if this cache is associated with a secondary iv
void setSecondaryIvUsage(bool status) { usesSecondaryIv_ = status; }
//! Sets the index of the iv (this scvf is embedded in) in its container
void setIvIndexInContainer(GridIndexType ivIndex) { ivIndexInContainer_ = ivIndex; }
//! Returns the index of the iv (this scvf is embedded in) in its container
GridIndexType ivIndexInContainer() const { return ivIndexInContainer_; }
//! Sets the iv-local face index
void setIvLocalFaceIndex(unsigned int idx) { ivLocalFaceIdx_ = idx; }
//! Sets the index of the face among the "positive" face's outside scvfs
void setIndexInOutsideFaces(unsigned int idx) { idxInOutsideFaces_ = idx; }
private:
//! indicates if cache has been fully updated
bool isUpdated_ = false;
//! indicates if cache is associated with secondary interaction volume
bool usesSecondaryIv_ = false;
bool isUpdated_ = false; //!< returns true if cache has been fully updated
bool usesSecondaryIv_ = false; //!< returns true if scvf is embedded in secondary interaction volume
//! the index of the iv (this scvf is embedded in) in its container
GridIndexType ivIndexInContainer_;
GridIndexType ivIndexInContainer_; //!< index of the iv (this scvf is embedded in) in its container
unsigned int ivLocalFaceIdx_; //!< the interaction volume-local face index of this scvf
unsigned int idxInOutsideFaces_; //!< index of scvf among outside scvfs of iv-local "positive" face (only surface grids)
};
} // end namespace Dumux
......
Supports Markdown
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