Commit 26dfeee8 authored by Dennis Gläser's avatar Dennis Gläser

[mpfa] reintroduce stencils in flux var caches

They are needed for analytic assembly
parent 00535bdf
......@@ -93,6 +93,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
using GridIndexType = typename GridView::IndexSet::IndexType;
// In the current implementation of the flux variables cache we cannot
// make a disctinction between dynamic (mpfa-o) and static (mpfa-l)
......@@ -106,24 +107,31 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
// of the fluxVarsCacheContainer
using Vector = Dune::DynamicVector< Scalar >;
using Matrix = Dune::DynamicMatrix< Scalar >;
using Stencil = std::vector< GridIndexType >;
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
using PrimaryStencil = typename PrimaryInteractionVolume::Traits::Stencil;
static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value,
"The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" );
static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value,
"The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" );
static_assert( std::is_convertible<PrimaryStencil*, Stencil*>::value,
"The stencil type used in primary interaction volumes is not convertible to std::vector<GridIndexType>!" );
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
using SecondaryStencil = typename SecondaryInteractionVolume::Traits::Stencil;
static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value,
"The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" );
static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value,
"The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" );
static_assert( std::is_convertible<SecondaryStencil*, Stencil*>::value,
"The stencil type used in secondary interaction volumes is not convertible to std::vector<GridIndexType>!" );
public:
// export the filler type
......@@ -147,6 +155,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const DataHandle& dataHandle,
const SubControlVolumeFace &scvf)
{
stencil_ = &iv.stencil();
switchFluxSign_ = localFaceData.isOutside();
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
......@@ -188,6 +197,9 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
}
}
//! The stencil corresponding to the transmissibilities
const Stencil& advectionStencil() const { return *stencil_; }
//! Coefficients for the cell (& Dirichlet) unknowns in flux expressions
const Vector& advectionTij() const { return *Tij_; }
......@@ -205,6 +217,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
private:
bool switchFluxSign_;
const Stencil* stencil_; //!< The stencil, i.e. the grid indices j
const Vector* Tij_; //!< The transmissibilities such that f = Tij*pj
std::array< Scalar, numPhases > g_; //!< Gravitational flux contribution on this face
std::array<const Vector*, numPhases> pj_; //!< The interaction-volume wide phase pressures pj
......
......@@ -100,26 +100,34 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
// is implemented. One idea to overcome the performance drop could be only
// storing the iv-local index here and obtain tij always from the datahandle
// of the fluxVarsCacheContainer
using GridIndexType = typename GridView::IndexSet::IndexType;
using Vector = Dune::DynamicVector< Scalar >;
using Matrix = Dune::DynamicMatrix< Scalar >;
using Stencil = std::vector< GridIndexType >;
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
using PrimaryStencil = typename PrimaryInteractionVolume::Traits::Stencil;
static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value,
"The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" );
static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value,
"The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" );
static_assert( std::is_convertible<PrimaryStencil*, Stencil*>::value,
"The stencil type used in primary interaction volumes is not convertible to std::vector<GridIndexType>!" );
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
using SecondaryStencil = typename SecondaryInteractionVolume::Traits::Stencil;
static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value,
"The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" );
static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value,
"The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" );
static_assert( std::is_convertible<SecondaryStencil*, Stencil*>::value,
"The stencil type used in secondary interaction volumes is not convertible to std::vector<GridIndexType>!" );
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
......@@ -148,6 +156,7 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SubControlVolumeFace &scvf,
unsigned int phaseIdx, unsigned int compIdx)
{
stencil_[phaseIdx][compIdx] = &iv.stencil();
switchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutside();
// store pointer to the mole fraction vector of this iv
......@@ -176,10 +185,15 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const Vector& moleFractions(unsigned int phaseIdx, unsigned int compIdx) const
{ return *xj_[phaseIdx][compIdx]; }
//! The stencils corresponding to the transmissibilities
const Stencil& diffusionStencil(unsigned int phaseIdx, unsigned int compIdx) const
{ return *stencil_[phaseIdx][compIdx]; }
private:
std::array< std::array<bool, numComponents>, numPhases > switchFluxSign_;
std::array< std::array<const Vector*, numComponents>, numPhases > Tij_; //!< The transmissibilities such that f = Tij*xj
std::array< std::array<const Vector*, numComponents>, numPhases > xj_; //!< The interaction-volume wide mole fractions xj
std::array< std::array<const Stencil*, numComponents>, numPhases > stencil_; //!< The stencils, i.e. the grid indices j
std::array< std::array<const Vector*, numComponents>, numPhases > Tij_; //!< The transmissibilities such that f = Tij*xj
std::array< std::array<const Vector*, numComponents>, numPhases > xj_; //!< The interaction-volume wide mole fractions xj
};
public:
......
......@@ -95,26 +95,34 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
// is implemented. One idea to overcome the performance drop could be only
// storing the iv-local index here and obtain tij always from the datahandle
// of the fluxVarsCacheContainer
using GridIndexType = typename GridView::IndexSet::IndexType;
using Vector = Dune::DynamicVector< Scalar >;
using Matrix = Dune::DynamicMatrix< Scalar >;
using Stencil = std::vector< GridIndexType >;
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
using PrimaryStencil = typename PrimaryInteractionVolume::Traits::Stencil;
static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value,
"The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" );
static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value,
"The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" );
static_assert( std::is_convertible<PrimaryStencil*, Stencil*>::value,
"The stencil type used in primary interaction volumes is not convertible to std::vector<GridIndexType>!" );
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
using SecondaryStencil = typename SecondaryInteractionVolume::Traits::Stencil;
static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value,
"The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" );
static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value,
"The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" );
static_assert( std::is_convertible<SecondaryStencil*, Stencil*>::value,
"The stencil type used in secondary interaction volumes is not convertible to std::vector<GridIndexType>!" );
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
......@@ -141,6 +149,7 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const DataHandle& dataHandle,
const SubControlVolumeFace &scvf)
{
stencil_ = &iv.stencil();
switchFluxSign_ = localFaceData.isOutside();
// store pointer to the temperature vector of this iv
......@@ -157,6 +166,9 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
//! Coefficients for the cell (& Dirichlet) unknowns in flux expressions
const Vector& heatConductionTij() const { return *Tij_; }
//! The stencil corresponding to the transmissibilities
const Stencil& heatConductionStencil() const { return *stencil_; }
//! 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.
......@@ -168,8 +180,9 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
private:
bool switchFluxSign_;
const Vector* Tij_; //!< The transmissibilities such that f = Tij*Tj
const Vector* Tj_; //!< The interaction-volume wide temperature Tj
const Stencil* stencil_; //!< The stencil, i.e. the grid indices j
const Vector* Tij_; //!< The transmissibilities such that f = Tij*Tj
const Vector* Tj_; //!< The interaction-volume wide temperature Tj
};
public:
......
......@@ -57,6 +57,8 @@ namespace Dumux
* using Matrix = ...;
* //! export the type used for iv-local vectors
* using Vector = ...;
* //! export the type used for the iv-stencils
* using Stencil = ...;
* //! export the data handle type for this iv
* using DataHandle = InteractionVolumeDataHandle< TypeTag, InteractionVolumeType >;
* \endcode
......@@ -112,6 +114,9 @@ public:
//! returns a reference to the container with the local face data
const typename Traits::LocalFaceDataContainer& localFaceData() const { asImp().localFaceData(); }
//! returns the cell-stencil of this interaction volume
const typename Traits::Stencil& stencil() const { return asImp().stencil(); }
//! returns the local scvf entity corresponding to a given iv-local scvf idx
const typename Traits::LocalScvfType& localScvf(typename Traits::LocalIndexType ivLocalScvfIdx) const
{ return asImp().localScvf(ivLocalScvfIdx); }
......
......@@ -83,6 +83,8 @@ public:
using Matrix = Dune::DynamicMatrix< Scalar >;
//! export the type used for iv-local vectors
using Vector = Dune::DynamicVector< Scalar >;
//! export the type used for the iv-stencils
using Stencil = std::vector< GridIndexType >;
//! export the data handle type for this iv
using DataHandle = InteractionVolumeDataHandle< TypeTag, InteractionVolumeType >;
};
......@@ -116,6 +118,7 @@ class CCMpfaInteractionVolumeImplementation< TypeTag, MpfaMethods::oMethod >
using IndexSet = typename TraitsType::IndexSet;
using GridIndexType = typename TraitsType::GridIndexType;
using LocalIndexType = typename TraitsType::LocalIndexType;
using Stencil = typename TraitsType::Stencil;
//! Data attached to scvf touching Dirichlet boundaries.
//! For the default o-scheme, we only store the corresponding vol vars index.
......@@ -144,6 +147,10 @@ public:
const Problem& problem,
const FVElementGeometry& fvGeometry)
{
// for the o-scheme, the stencil is equal to the scv
// index set of the dual grid's nodal index set
stencil_ = &indexSet.nodalIndexSet().globalScvIndices();
// number of interaction-volume-local (= node-local for o-scheme) scvs/scvf
numFaces_ = indexSet.numFaces();
const auto numLocalScvs = indexSet.numScvs();
......@@ -258,6 +265,9 @@ public:
//! returns the number of scvfs embedded in this interaction volume
std::size_t numScvfs() const { return scvfs_.size(); }
//! returns the cell-stencil of this interaction volume
const Stencil& stencil() const { return *stencil_; }
//! returns the grid element corresponding to a given iv-local scv idx
const Element& element(const LocalIndexType ivLocalScvIdx) const { return elements_[ivLocalScvIdx]; }
......@@ -313,6 +323,9 @@ public:
}
private:
// pointer to cell stencil (in iv index set)
const Stencil* stencil_;
// Variables defining the local scope
std::vector<Element> elements_;
std::vector<LocalScvType> scvs_;
......
......@@ -124,7 +124,7 @@ public:
"1p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!");
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
const auto& volVarIndices = fluxVarsCache.advectionVolVarsStencil();
const auto& volVarIndices = fluxVarsCache.advectionStencil();
const auto& tij = fluxVarsCache.advectionTij();
// we know the "upwind factor" is constant, get inner one here and compute derivatives
......@@ -211,7 +211,7 @@ public:
const SubControlVolumeFace& scvf) const
{
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
const auto& volVarIndices = fluxVarsCache.advectionVolVarsStencil();
const auto& volVarIndices = fluxVarsCache.advectionStencil();
const auto& tij = fluxVarsCache.advectionTij();
// we know the "upwind factor" is constant, get inner one here and compute derivatives
......
Markdown is supported
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