diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
index 01a9f3651fbc47913ad0c8cd4eee1d2423205ea0..c024ebbf6ae14ef534f1ea1cf6c8c1610d576bbd 100644
--- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
@@ -45,7 +45,7 @@ NEW_PROP_TAG(ProblemEnableGravity);
/*!
* \ingroup DarcysLaw
- * \brief Specialization of Darcy's Law for the CCTpfa method.
+ * \brief Specialization of Darcy's Law for the CCMpfa method.
*/
template <class TypeTag>
class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
diff --git a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
index 0d8352e5e0da8f5dc234182539d630d42c4e4947..9069ebeff39e87f9634115bfdd2e6e38b936a617 100644
--- a/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
+++ b/dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh
@@ -18,7 +18,7 @@
*****************************************************************************/
/*!
* \file
- * \brief The global object of flux var caches
+ * \brief The local object of flux var caches
*/
#ifndef DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_FLUXVARSCACHE_HH
#define DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_FLUXVARSCACHE_HH
@@ -32,7 +32,8 @@ namespace Dumux
/*!
* \ingroup ImplicitModel
- * \brief Base class for the stencil local flux variables cache
+ * \brief Base class for the local flux variables cache.
+ * Prepares the cache on all the faces in the stencil.
*/
template<class TypeTag, bool EnableGlobalFluxVariablesCache>
class CCMpfaElementFluxVariablesCache;
@@ -135,7 +136,8 @@ public:
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
- clear();
+ fluxVarsCache_.clear();
+ globalScvfIndices_.clear();
const auto& problem = globalFluxVarsCache().problem_();
const auto& globalFvGeometry = problem.model().globalFvGeometry();
@@ -143,8 +145,10 @@ public:
const auto& assemblyMap = problem.model().localJacobian().assemblyMap();
const auto globalI = problem.elementMapper().index(element);
- // reserve initial guess of memory (won't be enough though - several scvfs per neighbor will be required)
- globalScvfIndices_.reserve(fvGeometry.numScvf() + assemblyMap[globalI].size());
+ // reserve memory
+ auto numNeighborScvfs = 0;
+ for (auto&& facesInNeighbor : assemblyMap[globalI]) numNeighborScvfs += facesInNeighbor.size();
+ globalScvfIndices_.reserve(fvGeometry.numScvf() + numNeighborScvfs);
// first add all the indices inside the element
for (auto&& scvf : scvfs(fvGeometry))
@@ -155,17 +159,11 @@ public:
for (auto fluxVarIdx : assemblyMap[globalI][j])
globalScvfIndices_.push_back(fluxVarIdx);
- // make global indices unique
- std::sort(globalScvfIndices_.begin(), globalScvfIndices_.end());
- globalScvfIndices_.erase(std::unique(globalScvfIndices_.begin(), globalScvfIndices_.end()), globalScvfIndices_.end());
-
// prepare all the caches of the scvfs inside the corresponding interaction volumes using helper class
fluxVarsCache_.resize(globalScvfIndices_.size());
for (auto&& scvf : scvfs(fvGeometry))
- {
if (!(*this)[scvf].isUpdated())
FluxVariablesCacheFiller::fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, *this);
- }
// prepare the caches in the remaining neighbors
unsigned int j = 0;
@@ -215,12 +213,6 @@ public:
private:
const GlobalFluxVariablesCache* globalFluxVarsCachePtr_;
- void clear()
- {
- fluxVarsCache_.clear();
- globalScvfIndices_.clear();
- }
-
// This function updates the transmissibilities after the solution has been deflected during jacobian assembly
void update(const Element& element,
const FVElementGeometry& fvGeometry,
@@ -230,16 +222,14 @@ private:
(*this)[scvf].setUpdateStatus(false);
for (auto&& scvf : scvfs(fvGeometry))
- {
if (!(*this)[scvf].isUpdated())
FluxVariablesCacheFiller::updateFluxVarCache(globalFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf, *this);
- }
}
- // get index of scvf in the local container
+ // get index of an scvf in the local container
int getLocalScvfIdx_(const int scvfIdx) const
{
- auto it = std::lower_bound(globalScvfIndices_.begin(), globalScvfIndices_.end(), scvfIdx);
+ auto it = std::find(globalScvfIndices_.begin(), globalScvfIndices_.end(), scvfIdx);
assert(globalScvfIndices_[std::distance(globalScvfIndices_.begin(), it)] == scvfIdx && "Could not find the flux vars cache for scvfIdx");
return std::distance(globalScvfIndices_.begin(), it);
}
diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
index 3ace085a2de4b641c060d69a7d418da99655e4a6..63319c7e7375a620bbfa813e9a8a8514d71f25e9 100644
--- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
+++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
@@ -18,7 +18,7 @@
*****************************************************************************/
/*!
* \file
- * \brief The local (stencil) volume variables class for cell centered models
+ * \brief The local (stencil) volume variables class for cell centered mpfa models
*/
#ifndef DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_VOLUMEVARIABLES_HH
#define DUMUX_DISCRETIZATION_CCMPFA_ELEMENT_VOLUMEVARIABLES_HH
@@ -31,7 +31,7 @@ namespace Dumux
/*!
* \ingroup ImplicitModel
- * \brief Base class for the volume variables vector
+ * \brief Base class for the local volume variables vector
*/
template<class TypeTag, bool enableGlobalVolVarsCache>
class CCMpfaElementVolumeVariables
@@ -145,25 +145,10 @@ public:
++localIdx;
}
- // if the element is connected to a boundary, additionally prepare BCs
- bool boundary = element.hasBoundaryIntersections();
- if (!boundary)
+ // eventually prepare boundary volume variables
+ auto estimate = boundaryVolVarsEstimate_(element, fvGeometry);
+ if (estimate > 0)
{
- for (auto&& scvf : scvfs(fvGeometry))
- {
- if (globalFvGeometry.scvfTouchesBoundary(scvf))
- {
- boundary = true;
- break;
- }
- }
- }
-
- if (boundary)
- {
- // reserve memory (12 is the case of 3d hexahedron with one face on boundary)
- std::size_t estimate = 12;
- std::vector<IndexType> finishedBoundaries;
volumeVariables_.reserve(numDofs+estimate);
volVarIndices_.reserve(numDofs+estimate);
@@ -190,7 +175,7 @@ public:
// use the inside volume variables for neumann boundaries
else if (!useTpfaBoundary)
{
- volumeVariables_.emplace_back(VolumeVariables(volumeVariables_[0]));
+ volumeVariables_.emplace_back(volumeVariables_[0]);
volVarIndices_.push_back(scvf.outsideScvIdx());
}
}
@@ -208,12 +193,10 @@ public:
for (auto scvfIdx : ivSeed.globalScvfIndices())
{
auto&& ivScvf = fvGeometry.scvf(scvfIdx);
-
// only proceed for scvfs on the boundary and not in the inside element
if (!ivScvf.boundary() || ivScvf.insideScvIdx() == eIdx)
continue;
- // that means we are on a not yet handled boundary scvf
auto insideScvIdx = ivScvf.insideScvIdx();
auto insideElement = globalFvGeometry.element(insideScvIdx);
@@ -232,20 +215,15 @@ public:
// use the inside volume variables for neumann boundaries
else if (!useTpfaBoundary)
{
- volumeVariables_.emplace_back(VolumeVariables((*this)[insideScvIdx]));
+ volumeVariables_.emplace_back((*this)[insideScvIdx]);
volVarIndices_.push_back(ivScvf.outsideScvIdx());
}
}
}
-
- // free unused memory
- volumeVariables_.shrink_to_fit();
- volVarIndices_.shrink_to_fit();
}
}
// Binding of an element, prepares only the volume variables of the element
- // specialization for cc models
void bindElement(const Element& element,
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
@@ -279,7 +257,25 @@ public:
private:
const GlobalVolumeVariables* globalVolVarsPtr_;
- const int getLocalIdx_(const int volVarIdx) const
+ //! checks whether an scvf touches the boundary and returns an estimate of how many
+ //! boundary vol vars will be necessary. In 2d, this is the sum of faces touching
+ //! the boundary, which should be correct. In 3d, we count each face double - probably
+ //! too much for hexahedrons but might be even too little for simplices.
+ int boundaryVolVarsEstimate_(const Element& element,
+ const FVElementGeometry& fvGeometry)
+ {
+ int bVolVarEstimate = 0;
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ bool boundary = scvf.boundary();
+ if (boundary || (!boundary && fvGeometry.globalFvGeometry().scvfTouchesBoundary(scvf)))
+ bVolVarEstimate += dim-1;
+ }
+
+ return bVolVarEstimate;
+ }
+
+ int getLocalIdx_(const int volVarIdx) const
{
auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
assert(it != volVarIndices_.end() && "Could not find the current volume variables for volVarIdx!");
diff --git a/dumux/discretization/cellcentered/mpfa/facetypes.hh b/dumux/discretization/cellcentered/mpfa/facetypes.hh
index 8da25323436a11e6e15353e9f3853f38ef69c04e..87d526e167db11fbf695e0a603483bd62b4e1e74 100644
--- a/dumux/discretization/cellcentered/mpfa/facetypes.hh
+++ b/dumux/discretization/cellcentered/mpfa/facetypes.hh
@@ -18,7 +18,7 @@
*****************************************************************************/
/*!
* \file
- * \brief Face types of the sub control volume faces in the mpfa method.
+ * \brief Face types of the sub control volume faces in cell-centered mpfa methods.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FACETYPES_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_FACETYPES_HH
diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
index 7bf247d3bb64ece91d140d7a288ae479824d7ce1..fb27d52806dec57ffc6cfc30be9e20a92d33bbf1 100644
--- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
@@ -20,7 +20,7 @@
* \file
* \brief This file contains the data which is required to calculate
* molar and mass fluxes of a component in a fluid phase over a face of a finite volume by means
- * of Fick's Law. Specializations are provided for the different discretization methods.
+ * of Fick's Law for cell-centered MPFA models.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FICKS_LAW_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_FICKS_LAW_HH
diff --git a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
index 35f26093fc02f3debc0165922b17951aec406fe4..fcca346bcc23cf14e7aad66705c0f4f5cce965aa 100644
--- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
@@ -19,7 +19,7 @@
/*!
* \file
* \brief This file contains the data which is required to calculate
- * heat conduction fluxes with Fourier's law.
+ * heat conduction fluxes with Fourier's law for cell-centered MPFA models.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FOURIERS_LAW_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_FOURIERS_LAW_HH
diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index 8a6de9b1aaa444d725a8ae964ac28eb491aa3a2c..b7e74e6b686f246f5f5aec62a14690e2f890070a 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -18,7 +18,7 @@
*****************************************************************************/
/*!
* \file
- * \brief Base class for a local finite volume geometry for cell-centered TPFA models
+ * \brief Base class for a local finite volume geometry for cell-centered MPFA models
* This builds up the sub control volumes and sub control volume faces
* for each element in the local scope we are restricting to, e.g. stencil or element.
*/
@@ -34,7 +34,7 @@ namespace Dumux
{
/*!
* \ingroup ImplicitModel
- * \brief Base class for the finite volume geometry vector for cell-centered TPFA models
+ * \brief Base class for the finite volume geometry vector for cell-centered MPFA models
* This builds up the sub control volumes and sub control volume faces
* for each element.
*/
@@ -201,19 +201,18 @@ public:
auto it = std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx);
if (it != scvfIndices_.end())
{
- assert(flipScvfIndices_[std::distance(scvfIndices_.begin(), it)][outsideScvIdx] != -1);
- return neighborScvfs_[ flipScvfIndices_[std::distance(scvfIndices_.begin(), it)][outsideScvIdx] ];
+ const auto localScvfIdx = std::distance(scvfIndices_.begin(), it);
+ return neighborScvfs_[flipScvfIndices_[localScvfIdx][outsideScvIdx]];
}
else
{
const auto neighborScvfIdxLocal = findLocalIndex(scvfIdx, neighborScvfIndices_);
auto&& scvf = neighborScvfs_[neighborScvfIdxLocal];
- assert(neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx] != -1);
if (scvf.outsideScvIdx(outsideScvIdx) == scvIndices_[0])
- return scvfs_[ neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx] ];
+ return scvfs_[neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx]];
else
- return neighborScvfs_[ neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx] ];
+ return neighborScvfs_[neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx]];
}
}
@@ -263,7 +262,7 @@ public:
// reserve memory
const auto numNeighbors = neighborStencil.size();
- const auto estimate = numNeighbors*dim*(4*dim-4); // estimate holds for quadrilaterals in 2D/3D, overestimates else
+ const auto estimate = neighborScvfEstimate_(element, numNeighbors);
neighborScvs_.reserve(numNeighbors);
neighborScvIndices_.reserve(numNeighbors);
neighborScvfIndices_.reserve(estimate);
@@ -274,23 +273,14 @@ public:
unsigned int j = 0;
for (auto globalJ : neighborStencil)
{
- // get data on the neighbor
- auto elementJ = globalFvGeometry().element(globalJ);
- const auto& scvfIndices = assemblyMap[globalI][j];
-
// make neighbor geometries
- makeNeighborGeometries(elementJ, globalJ, scvfIndices);
+ auto elementJ = globalFvGeometry().element(globalJ);
+ makeNeighborGeometries(elementJ, globalJ, assemblyMap[globalI][j]);
// increment counter
j++;
}
- // free unused memory
- neighborScvs_.shrink_to_fit();
- neighborScvIndices_.shrink_to_fit();
- neighborScvfIndices_.shrink_to_fit();
- neighborScvfs_.shrink_to_fit();
-
// set up the scvf flip indices for network grids
if (dim < dimWorld)
makeFlipIndexSet();
@@ -310,6 +300,35 @@ public:
private:
+ //! Estimates the number of neighboring scvfs that have to be prepared
+ //! The estimate holds for inner elements, overestimats on the boundary
+ int neighborScvfEstimate_(const Element& element, int numNeighbors)
+ {
+ auto type = element.geometry().type();
+
+ if (type == Dune::GeometryType(Dune::GeometryType::simplex, dim))
+ {
+ if (dim == 2)
+ return 6 + ( (numNeighbors-3 > 0) ? (numNeighbors-3)*2 : 0 );
+ if (dim == 3)
+ return 12 + ( (numNeighbors-4 > 0) ? (numNeighbors-4)*3 : 0 );
+ }
+ else if (type == Dune::GeometryType(Dune::GeometryType::cube, dim))
+ {
+ if (dim == 2)
+ return 16;
+ if (dim == 3)
+ return 120;
+ }
+ else if (type == Dune::GeometryType(Dune::GeometryType::pyramid, dim))
+ return 16 + ( (numNeighbors-4 > 0) ? (numNeighbors-4)*3 : 0 );
+ else if (type == Dune::GeometryType(Dune::GeometryType::prism, dim))
+ return 18 + ( (numNeighbors-4 > 0) ? (numNeighbors-4)*3 : 0 );
+ else
+ DUNE_THROW(Dune::NotImplemented, "Mpfa for given geometry type.");
+
+ }
+
//! create scvs and scvfs of the bound element
void makeElementGeometries(const Element& element)
{
@@ -503,7 +522,7 @@ private:
const auto vIdxGlobal = scvf.vertexIndex();
const auto insideScvIdx = scvf.insideScvIdx();
- flipScvfIndices_[insideFace].resize(numOutsideScvs, -1);
+ flipScvfIndices_[insideFace].resize(numOutsideScvs);
for (unsigned int nIdx = 0; nIdx < numOutsideScvs; ++nIdx)
{
const auto outsideScvIdx = scvf.outsideScvIdx(nIdx);
@@ -535,7 +554,7 @@ private:
const auto vIdxGlobal = scvf.vertexIndex();
const auto insideScvIdx = scvf.insideScvIdx();
- neighborFlipScvfIndices_[outsideFace].resize(numOutsideScvs, -1);
+ neighborFlipScvfIndices_[outsideFace].resize(numOutsideScvs);
for (unsigned int nIdx = 0; nIdx < numOutsideScvs; ++nIdx)
{
const auto neighborScvIdx = scvf.outsideScvIdx(nIdx);
diff --git a/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh b/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh
index 9c9074500f6a8b90e44ef1822ae23e61022dc016..1201543f6a34508261fdc46c981b46a636dde271 100644
--- a/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh
+++ b/dumux/discretization/cellcentered/mpfa/globalfvgeometrybase.hh
@@ -95,6 +95,10 @@ public:
std::size_t numBoundaryScvf() const
{ return numBoundaryScvf_; }
+ //! The total number of vertices on the boundary
+ std::size_t numBoundaryVertices() const
+ { return numBoundaryVertices_; }
+
// Get an element from a sub control volume contained in it
Element element(const SubControlVolume& scv) const
{ return elementMap_.element(scv.elementIndex()); }
@@ -161,6 +165,7 @@ public:
// Build the SCVs and SCV faces
IndexType scvfIdx = 0;
numBoundaryScvf_ = 0;
+ numBoundaryVertices_ = 0;
for (const auto& element : elements(gridView_))
{
auto eIdx = problem.elementMapper().index(element);
@@ -235,7 +240,10 @@ public:
// store info on which vertices are on the domain boundary
if (boundary)
+ {
boundaryVertices_[vIdxGlobal] = true;
+ numBoundaryVertices_++;
+ }
// is vertex on a branching point?
if (dim < dimWorld && outsideIndices[indexInInside].size() > 1)
@@ -292,9 +300,6 @@ public:
scvfIndicesOfScv_[eIdx] = scvfIndexSet;
}
- // in parallel problems we might have reserved more scvfs than we actually use
- scvfs_.shrink_to_fit();
-
// Make the flip index set for network and surface grids
if (dim < dimWorld)
{
@@ -349,9 +354,7 @@ public:
{ return scvfs_[scvfIdx]; }
const SubControlVolumeFace& flipScvf(IndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
- {
- return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]];
- }
+ { return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]]; }
//! Get the sub control volume face indices of an scv by global index
const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
@@ -378,7 +381,8 @@ private:
std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
std::vector<bool> boundaryVertices_;
std::vector<bool> branchingVertices_;
- IndexType numBoundaryScvf_;
+ std::size_t numBoundaryScvf_;
+ std::size_t numBoundaryVertices_;
// needed for embedded surface and network grids (dim < dimWorld)
std::vector<std::vector<IndexType>> flipScvfIndices_;
@@ -435,6 +439,10 @@ public:
std::size_t numBoundaryScvf() const
{ return numBoundaryScvf_; }
+ //! The total number of vertices on the boundary
+ std::size_t numBoundaryVertices() const
+ { return numBoundaryVertices_; }
+
// Get an element from a sub control volume contained in it
Element element(const SubControlVolume& scv) const
{ return elementMap_.element(scv.elementIndex()); }
@@ -482,6 +490,7 @@ public:
numScvs_ = gridView_.size(0);
numScvf_ = 0;
numBoundaryScvf_ = 0;
+ numBoundaryVertices_ = 0;
elementMap_.resize(numScvs_);
scvfIndicesOfScv_.resize(numScvs_);
neighborVolVarIndices_.resize(numScvs_);
@@ -565,7 +574,10 @@ public:
// store info on which vertices are on the domain boundary
if (boundary)
+ {
boundaryVertices_[vIdxGlobal] = true;
+ numBoundaryVertices_++;
+ }
// is vertex on a branching point?
if (dim < dimWorld && outsideIndices[indexInInside].size() > 1)
@@ -628,9 +640,10 @@ private:
GridView gridView_;
// Information on the global number of geometries
- IndexType numScvs_;
- IndexType numScvf_;
- IndexType numBoundaryScvf_;
+ std::size_t numScvs_;
+ std::size_t numScvf_;
+ std::size_t numBoundaryScvf_;
+ std::size_t numBoundaryVertices_;
// vectors that store the global data
Dumux::ElementMap<GridView> elementMap_;
diff --git a/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh
index 87af94939467c06d328631046c2f2655949dabf4..e49919f688cdfa72b984b2eff8f925af00ee22b5 100644
--- a/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh
+++ b/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh
@@ -33,14 +33,7 @@ namespace Dumux
//! By default we simply inherit from the base class
//! Actual implementations for other methods have to be provided below
template<class TypeTag, MpfaMethods method>
-class CCMpfaGlobalInteractionVolumeSeedsImplementation : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
-{
- using ParentType = CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-
-public:
- CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView gridView) : ParentType(gridView) {}
-};
+class CCMpfaGlobalInteractionVolumeSeedsImplementation;
/*!
* \ingroup Mpfa
@@ -52,6 +45,8 @@ using CCMpfaGlobalInteractionVolumeSeeds = CCMpfaGlobalInteractionVolumeSeedsImp
} // end namespace
// the specializations of this class differing from the default have to be included here
+#include <dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh>
#include <dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh>
+#include <dumux/discretization/cellcentered/mpfa/omethodfps/globalinteractionvolumeseeds.hh>
#endif
diff --git a/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh b/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh
index 402e1a8645e43a5f42257877755949bb24c14adf..cadbd608266eda58b1d1f3acafdcf1d5480d21c0 100644
--- a/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh
+++ b/dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh
@@ -36,8 +36,6 @@ template<class TypeTag>
class CCMpfaGlobalInteractionVolumeSeedsBase
{
using Implementation = typename GET_PROP_TYPE(TypeTag, GlobalInteractionVolumeSeeds);
- // the actual implementation needs to be friend
- friend Implementation;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
@@ -51,21 +49,18 @@ class CCMpfaGlobalInteractionVolumeSeedsBase
using IndexType = typename GridView::IndexSet::IndexType;
public:
- CCMpfaGlobalInteractionVolumeSeedsBase(const GridView gridView) : gridView_(gridView) {}
+ CCMpfaGlobalInteractionVolumeSeedsBase(const GridView& gridView) : gridView_(gridView) {}
// initializes the interaction volumes or the seeds
- void update(const Problem& problem, std::vector<bool>& boundaryVertices)
+ void update(const Problem& p, const std::vector<bool>& boundaryVertices)
{
- problemPtr_ = &problem;
- seeds_.clear();
- boundarySeeds_.clear();
-
- auto numScvf = problem_().model().globalFvGeometry().numScvf();
- scvfIndexMap_.resize(numScvf);
- std::vector<bool> isFaceHandled(numScvf, false);
+ problemPtr_ = &p;
// initialize the seeds according to the mpfa method
- asImp_().initializeSeeds_(boundaryVertices, isFaceHandled);
+ asImp_().initializeSeeds(boundaryVertices,
+ scvfIndexMap_,
+ seeds_,
+ boundarySeeds_);
}
const InteractionVolumeSeed& seed(const SubControlVolumeFace& scvf) const
@@ -74,70 +69,13 @@ public:
const BoundaryInteractionVolumeSeed& boundarySeed(const SubControlVolumeFace& scvf) const
{ return boundarySeeds_[scvfIndexMap_[scvf.index()]]; }
-private:
+ const Problem& problem() const
+ { return *problemPtr_; }
- void initializeSeeds_(std::vector<bool>& boundaryVertices,
- std::vector<bool>& isFaceHandled)
- {
- const auto numScvf = problem_().model().globalFvGeometry().numScvf();
- const auto numBoundaryScvf = problem_().model().globalFvGeometry().numBoundaryScvf();
-
- // reserve memory
- seeds_.reserve(numScvf - numBoundaryScvf);
- boundarySeeds_.reserve(numBoundaryScvf);
-
- IndexType boundarySeedIndex = 0;
- IndexType seedIndex = 0;
- for (const auto& element : elements(gridView_))
- {
- auto fvGeometry = localView(problem_().model().globalFvGeometry());
- fvGeometry.bindElement(element);
- for (const auto& scvf : scvfs(fvGeometry))
- {
- // skip the rest if we already handled this face
- if (isFaceHandled[scvf.index()])
- continue;
-
- if (boundaryVertices[scvf.vertexIndex()])
- {
- // make the boundary interaction volume seed
- boundarySeeds_.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(problem_(), element, fvGeometry, scvf));
-
- // update the index map entries for the global scv faces in the interaction volume
- for (auto scvfIdxGlobal : boundarySeeds_.back().globalScvfIndices())
- {
- scvfIndexMap_[scvfIdxGlobal] = boundarySeedIndex;
- isFaceHandled[scvfIdxGlobal] = true;
- }
-
- // increment counter
- boundarySeedIndex++;
- }
- else
- {
- // make the inner interaction volume seed
- seeds_.emplace_back(Helper::makeInnerInteractionVolumeSeed(problem_(), element, fvGeometry, scvf));
-
- // update the index map entries for the global scv faces in the interaction volume
- for (auto scvfIdxGlobal : seeds_.back().globalScvfIndices())
- {
- scvfIndexMap_[scvfIdxGlobal] = seedIndex;
- isFaceHandled[scvfIdxGlobal] = true;
- }
-
- // increment counter
- seedIndex++;
- }
- }
- }
-
- // shrink vectors to actual size
- seeds_.shrink_to_fit();
- boundarySeeds_.shrink_to_fit();
- }
+ const GridView& gridView() const
+ { return gridView_; }
- const Problem& problem_() const
- { return *problemPtr_; }
+private:
const Implementation& asImp_() const
{ return *static_cast<Implementation*>(this); }
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh
index 86778212839be40182a94e40aa358a1e2a97f14f..263cee21036da10481f5b013e0068ad0c4d9eccd 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh
@@ -33,11 +33,10 @@ namespace Dumux
* \brief Specialization of the class for the mpfa-l method.
*/
template<class TypeTag>
-class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::lMethod> : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
+class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::lMethod>
+ : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
{
using ParentType = CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>;
- // the parent needs to be friend to access the private initialization method
- friend ParentType;
using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
@@ -50,25 +49,37 @@ class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::lMe
using Element = typename GridView::template Codim<0>::Entity;
public:
- CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView gridView) : ParentType(gridView) {}
+ CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView& gridView) : ParentType(gridView) {}
-private:
- void initializeSeeds_(std::vector<bool>& boundaryVertices, std::vector<bool>& isFaceHandled)
+ template<typename SeedVector, typename BoundarySeedVector>
+ void initializeSeeds(const std::vector<bool>& boundaryVertices,
+ std::vector<IndexType>& scvfIndexMap,
+ SeedVector& seeds,
+ BoundarySeedVector& boundarySeeds)
{
- const auto numScvf = this->problem_().model().globalFvGeometry().numScvf();
- const auto numBoundaryScvf = this->problem_().model().globalFvGeometry().numBoundaryScvf();
+ seeds.clear();
+ boundarySeeds.clear();
+ scvfIndexMap.clear();
// reserve memory
- this->seeds_.reserve(numScvf - numBoundaryScvf);
- this->boundarySeeds_.reserve(numBoundaryScvf);
+ const auto numScvf = this->problem().model().globalFvGeometry().numScvf();
+ const auto numBoundaryScvf = this->problem().model().globalFvGeometry().numBoundaryScvf();
+
+ seeds.reserve( std::size_t((numScvf-numBoundaryScvf)/2) );
+ boundarySeeds.reserve(numBoundaryScvf);
+ scvfIndexMap.resize(numScvf);
+
+ // Keep track of which faces have been handled already
+ std::vector<bool> isFaceHandled(numScvf, false);
IndexType boundarySeedIndex = 0;
IndexType seedIndex = 0;
- for (const auto& element : elements(this->gridView_))
+ for (const auto& element : elements(this->gridView()))
{
- auto fvGeometry = localView(this->problem_().model().globalFvGeometry());
+ auto fvGeometry = localView(this->problem().model().globalFvGeometry());
fvGeometry.bindElement(element);
- for (const auto& scvf : scvfs(fvGeometry))
+
+ for (auto&& scvf : scvfs(fvGeometry))
{
// skip the rest if we already handled this face
if (isFaceHandled[scvf.index()])
@@ -77,12 +88,15 @@ private:
if (boundaryVertices[scvf.vertexIndex()])
{
// make the boundary interaction volume seed
- this->boundarySeeds_.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(this->problem_(), element, fvGeometry, scvf));
+ boundarySeeds.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(this->problem(),
+ element,
+ fvGeometry,
+ scvf));
// update the index map entries for the global scv faces in the interaction volume
- for (auto scvfIdxGlobal : this->boundarySeeds_.back().globalScvfIndices())
+ for (auto scvfIdxGlobal : boundarySeeds.back().globalScvfIndices())
{
- this->scvfIndexMap_[scvfIdxGlobal] = boundarySeedIndex;
+ scvfIndexMap[scvfIdxGlobal] = boundarySeedIndex;
isFaceHandled[scvfIdxGlobal] = true;
}
@@ -92,14 +106,17 @@ private:
else
{
// make the inner interaction volume seed only if we are on highest level of all connected elements
- if (isLocalMaxLevel(element, scvf))
+ if (isLocalMaxLevel_(element, scvf))
{
- this->seeds_.emplace_back(Helper::makeInnerInteractionVolumeSeed(this->problem_(), element, fvGeometry, scvf));
+ seeds.emplace_back(Helper::makeInnerInteractionVolumeSeed(this->problem(),
+ element,
+ fvGeometry,
+ scvf));
// update the index map entries for the global scv faces in the interaction volume
- for (auto scvfIdxGlobal : this->seeds_.back().globalScvfIndices())
+ for (auto scvfIdxGlobal : seeds.back().globalScvfIndices())
{
- this->scvfIndexMap_[scvfIdxGlobal] = seedIndex;
+ scvfIndexMap[scvfIdxGlobal] = seedIndex;
isFaceHandled[scvfIdxGlobal] = true;
}
@@ -109,18 +126,13 @@ private:
}
}
}
-
- // shrink vectors to actual size
- this->seeds_.shrink_to_fit();
- this->boundarySeeds_.shrink_to_fit();
}
- bool isLocalMaxLevel(const Element& element, const SubControlVolumeFace& scvf) const
+ bool isLocalMaxLevel_(const Element& element, const SubControlVolumeFace& scvf) const
{
auto inLevel = element.level();
- for (auto outsideIdx : scvf.outsideScvIndices())
- if (this->problem_().model().globalFvGeometry().element(outsideIdx).level() > inLevel)
- return false;
+ if (this->problem().model().globalFvGeometry().element(scvf.outsideScvIdx()).level() > inLevel)
+ return false;
return true;
}
};
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh
new file mode 100644
index 0000000000000000000000000000000000000000..48a59d368120eb7be8f0953de27d8c56d7118d30
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh
@@ -0,0 +1,136 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for the global interaction volumes of the mpfa-o method.
+ */
+#ifndef DUMUX_DISCRETIZATION_MPFA_O_GLOBALINTERACTIONVOLUMESEEDS_HH
+#define DUMUX_DISCRETIZATION_MPFA_O_GLOBALINTERACTIONVOLUMESEEDS_HH
+
+#include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mpfa
+ * \brief Specialization of the class for the mpfa-o method.
+ */
+template<class TypeTag>
+class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::oMethod>
+ : public CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>
+{
+ using ParentType = CCMpfaGlobalInteractionVolumeSeedsBase<TypeTag>;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
+ using InteractionVolumeSeed = typename InteractionVolume::Seed;
+ using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
+ using BoundaryInteractionVolumeSeed = typename BoundaryInteractionVolume::Seed;
+
+ using IndexType = typename GridView::IndexSet::IndexType;
+
+ static const int dim = GridView::dimension;
+
+public:
+ CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView& gridView) : ParentType(gridView) {}
+
+ template<typename SeedVector, typename BoundarySeedVector>
+ void initializeSeeds(const std::vector<bool>& boundaryVertices,
+ std::vector<IndexType>& scvfIndexMap,
+ SeedVector& seeds,
+ BoundarySeedVector& boundarySeeds)
+ {
+ seeds.clear();
+ boundarySeeds.clear();
+ scvfIndexMap.clear();
+
+ // reserve memory
+ const auto numScvf = this->problem().model().globalFvGeometry().numScvf();
+ const auto numBoundaryVertices = this->problem().model().globalFvGeometry().numBoundaryVertices();
+ const auto numInteriorVertices = this->gridView().size(dim) - numBoundaryVertices;
+
+ if (numInteriorVertices > 0)
+ seeds.reserve(numInteriorVertices);
+ boundarySeeds.reserve(numBoundaryVertices);
+ scvfIndexMap.resize(numScvf);
+
+ // Keep track of which faces have been handled already
+ std::vector<bool> isFaceHandled(numScvf, false);
+
+ IndexType boundarySeedIndex = 0;
+ IndexType seedIndex = 0;
+ for (const auto& element : elements(this->gridView()))
+ {
+ auto fvGeometry = localView(this->problem().model().globalFvGeometry());
+ fvGeometry.bindElement(element);
+
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ // skip the rest if we already handled this face
+ if (isFaceHandled[scvf.index()])
+ continue;
+
+ if (boundaryVertices[scvf.vertexIndex()])
+ {
+ // make the boundary interaction volume seed
+ boundarySeeds.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(this->problem(),
+ element,
+ fvGeometry,
+ scvf));
+
+ // update the index map entries for the global scv faces in the interaction volume
+ for (auto scvfIdxGlobal : boundarySeeds.back().globalScvfIndices())
+ {
+ scvfIndexMap[scvfIdxGlobal] = boundarySeedIndex;
+ isFaceHandled[scvfIdxGlobal] = true;
+ }
+
+ // increment counter
+ boundarySeedIndex++;
+ }
+ else
+ {
+ // make the inner interaction volume seed
+ seeds.emplace_back(Helper::makeInnerInteractionVolumeSeed(this->problem(),
+ element,
+ fvGeometry,
+ scvf));
+
+ // update the index map entries for the global scv faces in the interaction volume
+ for (auto scvfIdxGlobal : seeds.back().globalScvfIndices())
+ {
+ scvfIndexMap[scvfIdxGlobal] = seedIndex;
+ isFaceHandled[scvfIdxGlobal] = true;
+ }
+
+ // increment counter
+ seedIndex++;
+ }
+ }
+ }
+ }
+};
+} // end namespace
+
+
+#endif
diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/globalinteractionvolumeseeds.hh b/dumux/discretization/cellcentered/mpfa/omethodfps/globalinteractionvolumeseeds.hh
new file mode 100644
index 0000000000000000000000000000000000000000..120c7acc355c5d0a3b322d27d55ac1b71a013d0d
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/omethodfps/globalinteractionvolumeseeds.hh
@@ -0,0 +1,47 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for the global interaction volumes of the mpfa-o-fps method.
+ */
+#ifndef DUMUX_DISCRETIZATION_MPFA_O_FPS_GLOBALINTERACTIONVOLUMESEEDS_HH
+#define DUMUX_DISCRETIZATION_MPFA_O_FPS_GLOBALINTERACTIONVOLUMESEEDS_HH
+
+#include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mpfa
+ * \brief Specialization of the class for the mpfa-o-fps method.
+ */
+template<class TypeTag>
+class CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::oMethodFps>
+ : public CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::oMethod>
+{
+ using ParentType = CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, MpfaMethods::oMethod>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+public:
+ CCMpfaGlobalInteractionVolumeSeedsImplementation(const GridView& gridView) : ParentType(gridView) {}
+};
+} // end namespace
+
+
+#endif
diff --git a/dumux/porousmediumflow/implicit/fluxvariablescache.hh b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
index e59ad8dbf76c0d36d9d96800264759dc53a889ac..ee9551e7d86b1394d33fd3be33ab76a69d07cd79 100644
--- a/dumux/porousmediumflow/implicit/fluxvariablescache.hh
+++ b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
@@ -155,49 +155,36 @@ private:
Scalar tij_;
};
+///////////////////////////////////////////////////////////////////////////////////////////////////////////////
+//! Classes building up the porous medium flux variables cache for mpfa methods
+//! The cache is dependent on the active physical processes (advection, diffusion, heat conduction)
+//! For each type of process there is a base cache storing the data required to compute the respective fluxes
+//! Specializations of the overall cache are provided for combinations of processes
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
// forward declaration of the base class of the mpfa flux variables cache
template<class TypeTag, bool EnableAdvection, bool EnableMolecularDiffusion, bool EnableEnergyBalance>
-class MpfaPorousMediumFluxVariablesCache
-{};
+class MpfaPorousMediumFluxVariablesCache {};
-// specialization for cell centered mpfa methods
+//! Base class for the advective cache in mpfa methods
template<class TypeTag>
-class PorousMediumFluxVariablesCacheImplementation<TypeTag, DiscretizationMethods::CCMpfa>
- : public MpfaPorousMediumFluxVariablesCache<TypeTag,
- GET_PROP_VALUE(TypeTag, EnableAdvection),
- GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion),
- GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>
-{};
-
-// specialization for the case of pure advection
-template<class TypeTag>
-class MpfaPorousMediumFluxVariablesCache<TypeTag, true, false, false>
+class MpfaAdvectionCache
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
+ using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
// We always use the dynamic types here to be compatible on the boundary
+ using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
using TransmissibilityVector = typename BoundaryInteractionVolume::Vector;
using PositionVector = typename BoundaryInteractionVolume::PositionVector;
public:
- //! the constructor
- MpfaPorousMediumFluxVariablesCache() : isUpdated_(false)
- {
- // We have to initialize the neumann fluxes to zero (for inner interaction volumes)
- phaseNeumannFluxes_.fill(0.0);
- }
+ MpfaAdvectionCache() { phaseNeumannFluxes_.fill(0.0); }
//! update cached objects
template<typename InteractionVolume>
@@ -241,64 +228,33 @@ public:
Scalar advectionNeumannFlux(const unsigned int phaseIdx) const
{ return phaseNeumannFluxes_[phaseIdx]; }
- //! Returns whether or not this cache has been updated
- bool isUpdated() const
- { return isUpdated_; }
-
- //! Sets the update status from outside. Allows an update of the cache specific
- //! to processes that have solution dependent parameters, e.g. only updating
- //! the diffusion transmissibilities leaving the advective ones untouched
- void setUpdateStatus(const bool status)
- {
- isUpdated_ = status;
- }
-
private:
- bool isUpdated_;
+ // Quantities associated with advection
Stencil volVarsStencil_;
PositionVector volVarsPositions_;
TransmissibilityVector tij_;
std::array<Scalar, numPhases> phaseNeumannFluxes_;
};
-// specialization for the case of advection & diffusion
+//! Base class for the diffusive cache in mpfa methods
template<class TypeTag>
-class MpfaPorousMediumFluxVariablesCache<TypeTag, true, true, false>
+class MpfaDiffusionCache
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
+ using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
static const int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
// We always use the dynamic types here to be compatible on the boundary
+ using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
using TransmissibilityVector = typename BoundaryInteractionVolume::Vector;
using PositionVector = typename BoundaryInteractionVolume::PositionVector;
public:
- // the constructor
- MpfaPorousMediumFluxVariablesCache() : isUpdated_(false) {}
-
- // update cached objects for the advective fluxes
- template<typename InteractionVolume>
- void updateAdvection(const SubControlVolumeFace &scvf,
- const InteractionVolume& interactionVolume)
- {
- const auto& localFaceData = interactionVolume.getLocalFaceData(scvf);
- advectionVolVarsStencil_ = interactionVolume.volVarsStencil();
- advectionVolVarsPositions_ = interactionVolume.volVarsPositions();
- advectionTij_ = interactionVolume.getTransmissibilities(localFaceData);
- }
-
// update cached objects for the diffusive fluxes
template<typename InteractionVolume>
void updateDiffusion(const SubControlVolumeFace &scvf,
@@ -311,29 +267,6 @@ public:
diffusionTij_[phaseIdx][compIdx] = interactionVolume.getTransmissibilities(localFaceData);
}
- //! This method is here for compatibility reasons
- //! TODO: How to implement neumann fluxes for !useTpfa when diffusion is active?
- template<typename InteractionVolume>
- void updatePhaseNeumannFlux(const SubControlVolumeFace &scvf,
- const InteractionVolume& interactionVolume,
- const unsigned int phaseIdx) {}
-
- //! Returns the volume variables indices necessary for flux computation
- //! This includes all participating boundary volume variables. Since we
- //! do not allow mixed BC for the mpfa this is the same for all phases.
- const Stencil& advectionVolVarsStencil(const unsigned int phaseIdx) const
- { return advectionVolVarsStencil_; }
-
- //! Returns the position on which the volume variables live. This is
- //! necessary as we need to evaluate gravity also for the boundary volvars
- const PositionVector& advectionVolVarsPositions(const unsigned int phaseIdx) const
- { return advectionVolVarsPositions_; }
-
- //! Returns the transmissibilities associated with the volume variables
- //! All phases flow through the same rock, thus, tij are equal for all phases
- const TransmissibilityVector& advectionTij(const unsigned int phaseIdx) const
- { return advectionTij_; }
-
//! Returns the volume variables indices necessary for diffusive flux
//! computation. This includes all participating boundary volume variables
//! and it can be different for the phases & components.
@@ -347,77 +280,27 @@ public:
const unsigned int compIdx) const
{ return diffusionTij_[phaseIdx][compIdx]; }
- //! This method is needed for compatibility reasons
- //! TODO: How to implement neumann fluxes for !useTpfa when diffusion is active?
- Scalar advectionNeumannFlux(const unsigned int phaseIdx) const
- { return 0.0; }
-
- //! Returns whether or not this cache has been updated
- bool isUpdated() const
- { return isUpdated_; }
-
- //! Sets the update status from outside. Allows an update of the cache specific
- //! to processes that have solution dependent parameters, e.g. only updating
- //! the diffusion transmissibilities leaving the advective ones untouched
- void setUpdateStatus(const bool status)
- {
- isUpdated_ = status;
- }
-
private:
- bool isUpdated_;
- // Quantities associated with advection
- Stencil advectionVolVarsStencil_;
- PositionVector advectionVolVarsPositions_;
- TransmissibilityVector advectionTij_;
-
// Quantities associated with molecular diffusion
std::array< std::array<Stencil, numComponents>, numPhases> diffusionVolVarsStencils_;
std::array< std::array<TransmissibilityVector, numComponents>, numPhases> diffusionTij_;
};
-// specialization for the case of advection & heat conduction
+//! Base class for the heat conduction cache in mpfa methods
template<class TypeTag>
-class MpfaPorousMediumFluxVariablesCache<TypeTag, true, false, true>
+class MpfaHeatConductionCache
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
-
- static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
// We always use the dynamic types here to be compatible on the boundary
+ using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
using TransmissibilityVector = typename BoundaryInteractionVolume::Vector;
- using PositionVector = typename BoundaryInteractionVolume::PositionVector;
-
public:
- // the constructor
- MpfaPorousMediumFluxVariablesCache()
- : isUpdated_(false),
- heatNeumannFlux_(0.0)
- {
- // We have to initialize the neumann fluxes to zero (for inner interaction volumes)
- phaseNeumannFluxes_.fill(0.0);
- }
-
- // update cached objects for the advective fluxes
- template<typename InteractionVolume>
- void updateAdvection(const SubControlVolumeFace &scvf,
- const InteractionVolume& interactionVolume)
- {
- const auto& localFaceData = interactionVolume.getLocalFaceData(scvf);
- advectionVolVarsStencil_ = interactionVolume.volVarsStencil();
- advectionVolVarsPositions_ = interactionVolume.volVarsPositions();
- advectionTij_ = interactionVolume.getTransmissibilities(localFaceData);
- }
+ MpfaHeatConductionCache() : heatNeumannFlux_(0.0) {}
// update cached objects for heat conduction
template<typename InteractionVolume>
@@ -429,42 +312,6 @@ public:
heatConductionTij_ = interactionVolume.getTransmissibilities(localFaceData);
}
- //! update cached neumann boundary flux
- template<typename InteractionVolume>
- void updatePhaseNeumannFlux(const SubControlVolumeFace &scvf,
- const InteractionVolume& interactionVolume,
- const unsigned int phaseIdx)
- {
- const auto& localFaceData = interactionVolume.getLocalFaceData(scvf);
- phaseNeumannFluxes_[phaseIdx] = interactionVolume.getNeumannFlux(localFaceData);
- }
-
- //! update cached neumann boundary flux
- template<typename InteractionVolume>
- void updateHeatNeumannFlux(const SubControlVolumeFace &scvf,
- const InteractionVolume& interactionVolume,
- const unsigned int phaseIdx)
- {
- const auto& localFaceData = interactionVolume.getLocalFaceData(scvf);
- heatNeumannFlux_ = interactionVolume.getNeumannFlux(localFaceData);
- }
-
- //! Returns the volume variables indices necessary for flux computation
- //! This includes all participating boundary volume variables. Since we
- //! do not allow mixed BC for the mpfa this is the same for all phases.
- const Stencil& advectionVolVarsStencil(const unsigned int phaseIdx) const
- { return advectionVolVarsStencil_; }
-
- //! Returns the position on which the volume variables live. This is
- //! necessary as we need to evaluate gravity also for the boundary volvars
- const PositionVector& advectionVolVarsPositions(const unsigned int phaseIdx) const
- { return advectionVolVarsPositions_; }
-
- //! Returns the transmissibilities associated with the volume variables
- //! All phases flow through the same rock, thus, tij are equal for all phases
- const TransmissibilityVector& advectionTij(const unsigned int phaseIdx) const
- { return advectionTij_; }
-
//! Returns the volume variables indices necessary for heat conduction flux
//! computation. This includes all participating boundary volume variables
//! and it can be different for the phases & components.
@@ -476,16 +323,40 @@ public:
const TransmissibilityVector& heatConductionTij() const
{ return heatConductionTij_; }
- //! If the useTpfaBoundary property is set to false, the boundary conditions
- //! are put into the local systems leading to possible contributions on all faces
- Scalar advectionNeumannFlux(const unsigned int phaseIdx) const
- { return phaseNeumannFluxes_[phaseIdx]; }
-
//! If the useTpfaBoundary property is set to false, the boundary conditions
//! are put into the local systems leading to possible contributions on all faces
Scalar heatNeumannFlux() const
{ return heatNeumannFlux_; }
+private:
+ // Quantities associated with heat conduction
+ Stencil heatConductionVolVarsStencil_;
+ TransmissibilityVector heatConductionTij_;
+ Scalar heatNeumannFlux_;
+};
+
+// specialization of the flux variables cache for cell centered mpfa methods
+template<class TypeTag>
+class PorousMediumFluxVariablesCacheImplementation<TypeTag, DiscretizationMethods::CCMpfa>
+ : public MpfaPorousMediumFluxVariablesCache<TypeTag,
+ GET_PROP_VALUE(TypeTag, EnableAdvection),
+ GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion),
+ GET_PROP_VALUE(TypeTag, EnableEnergyBalance)> {};
+
+// specialization for the case of pure advection
+template<class TypeTag>
+class MpfaPorousMediumFluxVariablesCache<TypeTag, true, false, false>
+ : public MpfaAdvectionCache<TypeTag>
+{
+ using AdvectionCache = MpfaAdvectionCache<TypeTag>;
+
+public:
+ //! the constructor
+ MpfaPorousMediumFluxVariablesCache()
+ : AdvectionCache(),
+ isUpdated_(false)
+ {}
+
//! Returns whether or not this cache has been updated
bool isUpdated() const
{ return isUpdated_; }
@@ -500,148 +371,123 @@ public:
private:
bool isUpdated_;
- // Quantities associated with advection
- Stencil advectionVolVarsStencil_;
- PositionVector advectionVolVarsPositions_;
- TransmissibilityVector advectionTij_;
- std::array<Scalar, numPhases> phaseNeumannFluxes_;
-
- // Quantities associated with heat conduction
- Stencil heatConductionVolVarsStencil_;
- TransmissibilityVector heatConductionTij_;
- Scalar heatNeumannFlux_;
};
-// specialization for the case of advection, diffusion & heat conduction
+// specialization for the case of advection & diffusion
template<class TypeTag>
-class MpfaPorousMediumFluxVariablesCache<TypeTag, true, true, true>
+class MpfaPorousMediumFluxVariablesCache<TypeTag, true, true, false>
+ : public MpfaAdvectionCache<TypeTag>,
+ public MpfaDiffusionCache<TypeTag>
{
+ using AdvectionCache = MpfaAdvectionCache<TypeTag>;
+ using DiffusionCache = MpfaDiffusionCache<TypeTag>;
+
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
-
- static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
- static const int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
-
- // We always use the dynamic types here to be compatible on the boundary
- using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
- using TransmissibilityVector = typename BoundaryInteractionVolume::Vector;
- using PositionVector = typename BoundaryInteractionVolume::PositionVector;
public:
// the constructor
MpfaPorousMediumFluxVariablesCache()
- : isUpdated_(false),
- heatNeumannFlux_(0.0)
+ : AdvectionCache(),
+ DiffusionCache(),
+ isUpdated_(false)
{}
- // update cached objects for the advective fluxes
+ //! For compositional problems, neumann fluxes are not associated with a phase anymore
+ //! TODO: How to implement neumann fluxes for !useTpfa
template<typename InteractionVolume>
- void updateAdvection(const SubControlVolumeFace &scvf,
- const InteractionVolume& interactionVolume)
- {
- const auto& localFaceData = interactionVolume.getLocalFaceData(scvf);
- advectionVolVarsStencil_ = interactionVolume.volVarsStencil();
- advectionVolVarsPositions_ = interactionVolume.volVarsPositions();
- advectionTij_ = interactionVolume.getTransmissibilities(localFaceData);
- }
+ void updatePhaseNeumannFlux(const SubControlVolumeFace &scvf,
+ const InteractionVolume& interactionVolume,
+ const unsigned int phaseIdx) {}
- // update cached objects for the diffusive fluxes
- template<typename InteractionVolume>
- void updateDiffusion(const SubControlVolumeFace &scvf,
- const InteractionVolume& interactionVolume,
- const unsigned int phaseIdx,
- const unsigned int compIdx)
- {
- const auto& localFaceData = interactionVolume.getLocalFaceData(scvf);
- diffusionVolVarsStencils_[phaseIdx][compIdx] = interactionVolume.volVarsStencil();
- diffusionTij_[phaseIdx][compIdx] = interactionVolume.getTransmissibilities(localFaceData);
- }
+ //! TODO: How to implement neumann fluxes for !useTpfa
+ Scalar advectionNeumannFlux(const unsigned int phaseIdx) const
+ { return 0.0; }
- // update cached objects for heat conduction
- template<typename InteractionVolume>
- void updateHeatConduction(const SubControlVolumeFace &scvf,
- const InteractionVolume& interactionVolume)
+ //! Returns whether or not this cache has been updated
+ bool isUpdated() const
+ { return isUpdated_; }
+
+ //! Sets the update status from outside. Allows an update of the cache specific
+ //! to processes that have solution dependent parameters, e.g. only updating
+ //! the diffusion transmissibilities leaving the advective ones untouched
+ void setUpdateStatus(const bool status)
{
- const auto& localFaceData = interactionVolume.getLocalFaceData(scvf);
- heatConductionVolVarsStencil_ = interactionVolume.volVarsStencil();
- heatConductionTij_ = interactionVolume.getTransmissibilities(localFaceData);
+ isUpdated_ = status;
}
- //! update cached neumann boundary flux
- template<typename InteractionVolume>
- void updateHeatNeumannFlux(const SubControlVolumeFace &scvf,
- const InteractionVolume& interactionVolume,
- const unsigned int phaseIdx)
+private:
+ bool isUpdated_;
+};
+
+// specialization for the case of advection & heat conduction
+template<class TypeTag>
+class MpfaPorousMediumFluxVariablesCache<TypeTag, true, false, true>
+ : public MpfaAdvectionCache<TypeTag>,
+ public MpfaHeatConductionCache<TypeTag>
+{
+ using AdvectionCache = MpfaAdvectionCache<TypeTag>;
+ using HeatConductionCache = MpfaHeatConductionCache<TypeTag>;
+
+public:
+ // the constructor
+ MpfaPorousMediumFluxVariablesCache()
+ : AdvectionCache(),
+ HeatConductionCache(),
+ isUpdated_(false)
+ {}
+
+ //! Returns whether or not this cache has been updated
+ bool isUpdated() const
+ { return isUpdated_; }
+
+ //! Sets the update status from outside. Allows an update of the cache specific
+ //! to processes that have solution dependent parameters, e.g. only updating
+ //! the diffusion transmissibilities leaving the advective ones untouched
+ void setUpdateStatus(const bool status)
{
- const auto& localFaceData = interactionVolume.getLocalFaceData(scvf);
- heatNeumannFlux_ = interactionVolume.getNeumannFlux(localFaceData);
+ isUpdated_ = status;
}
- //! This method is here for compatibility reasons
+private:
+ bool isUpdated_;
+};
+
+// specialization for the case of advection, diffusion & heat conduction
+template<class TypeTag>
+class MpfaPorousMediumFluxVariablesCache<TypeTag, true, true, true>
+ : public MpfaAdvectionCache<TypeTag>,
+ public MpfaDiffusionCache<TypeTag>,
+ public MpfaHeatConductionCache<TypeTag>
+{
+ using AdvectionCache = MpfaAdvectionCache<TypeTag>;
+ using DiffusionCache = MpfaDiffusionCache<TypeTag>;
+ using HeatConductionCache = MpfaHeatConductionCache<TypeTag>;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
+public:
+ // the constructor
+ MpfaPorousMediumFluxVariablesCache()
+ : AdvectionCache(),
+ DiffusionCache(),
+ HeatConductionCache(),
+ isUpdated_(false)
+ {}
+
//! TODO: How to implement neumann fluxes for !useTpfa when diffusion/heat conduction is active?
template<typename InteractionVolume>
void updatePhaseNeumannFlux(const SubControlVolumeFace &scvf,
const InteractionVolume& interactionVolume,
const unsigned int phaseIdx) {}
- //! Returns the volume variables indices necessary for flux computation
- //! This includes all participating boundary volume variables. Since we
- //! do not allow mixed BC for the mpfa this is the same for all phases.
- const Stencil& advectionVolVarsStencil(const unsigned int phaseIdx) const
- { return advectionVolVarsStencil_; }
-
- //! Returns the position on which the volume variables live. This is
- //! necessary as we need to evaluate gravity also for the boundary volvars
- const PositionVector& advectionVolVarsPositions(const unsigned int phaseIdx) const
- { return advectionVolVarsPositions_; }
-
- //! Returns the transmissibilities associated with the volume variables
- //! All phases flow through the same rock, thus, tij are equal for all phases
- const TransmissibilityVector& advectionTij(const unsigned int phaseIdx) const
- { return advectionTij_; }
-
- //! Returns the volume variables indices necessary for diffusive flux
- //! computation. This includes all participating boundary volume variables
- //! and it can be different for the phases & components.
- const Stencil& diffusionVolVarsStencil(const unsigned int phaseIdx,
- const unsigned int compIdx) const
- { return diffusionVolVarsStencils_[phaseIdx][compIdx]; }
-
- //! Returns the transmissibilities associated with the volume variables
- //! This can be different for the phases & components.
- const TransmissibilityVector& diffusionTij(const unsigned int phaseIdx,
- const unsigned int compIdx) const
- { return diffusionTij_[phaseIdx][compIdx]; }
-
- //! Returns the volume variables indices necessary for heat conduction flux
- //! computation. This includes all participating boundary volume variables
- //! and it can be different for the phases & components.
- const Stencil& heatConductionVolVarsStencil() const
- { return heatConductionVolVarsStencil_; }
-
- //! Returns the transmissibilities associated with the volume variables
- //! This can be different for the phases & components.
- const TransmissibilityVector& heatConductionTij() const
- { return heatConductionTij_; }
//! This method is needed for compatibility reasons
//! TODO: How to implement neumann fluxes for !useTpfa when diffusion is active?
Scalar advectionNeumannFlux(const unsigned int phaseIdx) const
{ return 0.0; }
- //! If the useTpfaBoundary property is set to false, the boundary conditions
- //! are put into the local systems leading to possible contributions on all faces
- Scalar heatNeumannFlux() const
- { return heatNeumannFlux_; }
-
//! Returns whether or not this cache has been updated
bool isUpdated() const
{ return isUpdated_; }
@@ -656,19 +502,6 @@ public:
private:
bool isUpdated_;
- // Quantities associated with advection
- Stencil advectionVolVarsStencil_;
- PositionVector advectionVolVarsPositions_;
- TransmissibilityVector advectionTij_;
-
- // Quantities associated with molecular diffusion
- std::array< std::array<Stencil, numComponents>, numPhases> diffusionVolVarsStencils_;
- std::array< std::array<TransmissibilityVector, numComponents>, numPhases> diffusionTij_;
-
- // Quantities associated with heat conduction
- Stencil heatConductionVolVarsStencil_;
- TransmissibilityVector heatConductionTij_;
- Scalar heatNeumannFlux_;
};
} // end namespace