From ca47f527bf82cd4b7fb769e3f8463db2d7533610 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= <dennis.glaeser@iws.uni-stuttgart.de>
Date: Thu, 15 Nov 2018 16:57:47 +0100
Subject: [PATCH] [mpfa][ivdatahandle] get rid of resize method
Resizing of the matrices place in the local assembler directly wherr
needed. This has a small negative impact on performance but saves a lot of code.
---
.../mpfa/fluxvariablescachefiller.hh | 2 -
.../mpfa/interactionvolumedatahandle.hh | 122 +-----------------
.../mpfa/omethod/localassembler.hh | 80 +++++++++---
3 files changed, 61 insertions(+), 143 deletions(-)
diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
index dfac3b1cb1..0a59d20489 100644
--- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
+++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
@@ -125,7 +125,6 @@ public:
// prepare the corresponding data handle
fluxVarsCacheContainer.secondaryDataHandles().emplace_back();
secondaryIvDataHandle_ = &fluxVarsCacheContainer.secondaryDataHandles().back();
- secondaryIvDataHandle_->resize(*secondaryIv_);
// fill the caches for all the scvfs in the interaction volume
fillCachesInInteractionVolume_(fluxVarsCacheContainer, *secondaryIv_, *secondaryIvDataHandle_, ivIndexInContainer, true);
@@ -156,7 +155,6 @@ public:
// prepare the corresponding data handle
fluxVarsCacheContainer.primaryDataHandles().emplace_back();
primaryIvDataHandle_ = &fluxVarsCacheContainer.primaryDataHandles().back();
- primaryIvDataHandle_->resize(*primaryIv_);
// fill the caches for all the scvfs in the interaction volume
fillCachesInInteractionVolume_(fluxVarsCacheContainer, *primaryIv_, *primaryIvDataHandle_, ivIndexInContainer, true);
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
index e768aa65ac..863bc739d8 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
@@ -29,18 +29,12 @@
#include <dune/common/dynvector.hh>
#include <dumux/common/parameters.hh>
-#include <dumux/common/matrixvectorhelper.hh>
namespace Dumux
{
//! Empty data handle class
-class EmptyDataHandle
-{
-public:
- template< class InteractionVolume >
- void resize(const InteractionVolume& iv) {}
-};
+class EmptyDataHandle {};
//! Data handle for quantities related to advection
template<class MatVecTraits, class PhysicsTraits, bool EnableAdvection>
@@ -58,47 +52,6 @@ class AdvectionDataHandle
static constexpr int numPhases = PhysicsTraits::numPhases;
public:
-
- //! prepare data handle for subsequent fill
- template< class IV >
- void resize(const IV& iv)
- {
- using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
-
- // resize matrices
- resizeMatrix(T_, iv.numFaces(), iv.numKnowns());
- resizeMatrix(A_, iv.numUnknowns(), iv.numUnknowns());
- resizeMatrix(CA_, iv.numFaces(), iv.numUnknowns());
-
- // lambdas to resize a vector
- auto resizeToNumFaces = [&iv] (auto& v) { resizeVector(v, iv.numFaces()); };
- auto resizeToNumKnowns = [&iv] (auto& v) { resizeVector(v, iv.numKnowns()); };
-
- // resize pressure/gravity vector
- std::for_each(p_.begin(), p_.end(), resizeToNumKnowns);
- std::for_each(g_.begin(), g_.end(), resizeToNumFaces);
-
- // on surface grids, resize additional containers
- static constexpr int dim = IV::Traits::GridView::dimension;
- static constexpr int dimWorld = IV::Traits::GridView::dimensionworld;
- if (dim < dimWorld)
- {
- outsideT_.resize(iv.numFaces());
- std::for_each(outsideG_.begin(), outsideG_.end(), resizeToNumFaces);
-
- // resize the entries for each face
- for (LocalIndexType i = 0; i < iv.numFaces(); ++i)
- {
- const auto& localScvf = iv.localScvf(i);
- const auto numNeighbors = localScvf.neighboringLocalScvIndices().size()-1;
- outsideT_[i].resize(numNeighbors);
-
- std::for_each(outsideT_[i].begin(), outsideT_[i].end(), resizeToNumKnowns);
- std::for_each(outsideG_.begin(), outsideG_.end(), [&] (auto& v) { resizeVector(v[i], numNeighbors); });
- }
- }
- }
-
//! Access to the iv-wide pressure of one phase
const CellVector& pressures(unsigned int pIdx) const { return p_[pIdx]; }
CellVector& pressures(unsigned int pIdx) { return p_[pIdx]; }
@@ -111,7 +64,6 @@ public:
const std::array< FaceVector, numPhases >& gravity() const { return g_; }
std::array< FaceVector, numPhases >& gravity() { return g_; }
-
//! Projection matrix for Neumann/gravity contribution computation
const CMatrix& advectionCA() const { return CA_; }
CMatrix& advectionCA() { return CA_; }
@@ -162,37 +114,6 @@ public:
void setPhaseIndex(unsigned int phaseIdx) { phaseIdx_ = phaseIdx; }
void setComponentIndex(unsigned int compIdx) { compIdx_ = compIdx; }
- //! prepare data handle for subsequent fill
- template< class InteractionVolume >
- void resize(const InteractionVolume& iv)
- {
- for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
- {
- for (unsigned int cIdx = 0; cIdx < numComponents; ++cIdx)
- {
- // resize transmissibility matrix & mole fraction vector
- resizeMatrix(T_[pIdx][cIdx], iv.numFaces(), iv.numKnowns());
- resizeVector(xj_[pIdx][cIdx], iv.numKnowns());
-
- // resize outsideTij on surface grids
- using GridView = typename InteractionVolume::Traits::GridView;
- if (int(GridView::dimension) < int(GridView::dimensionworld))
- {
- outsideT_[pIdx][cIdx].resize(iv.numFaces());
-
- using LocalIndexType = typename InteractionVolume::Traits::IndexSet::LocalIndexType;
- for (LocalIndexType i = 0; i < iv.numFaces(); ++i)
- {
- const auto numNeighbors = iv.localScvf(i).neighboringLocalScvIndices().size() - 1;
- outsideT_[pIdx][cIdx][i].resize(numNeighbors);
- for (LocalIndexType j = 0; j < numNeighbors; ++j)
- resizeVector(outsideT_[pIdx][cIdx][i][j], iv.numKnowns());
- }
- }
- }
- }
- }
-
//! Access to the iv-wide mole fractions of a component in one phase
const CellVector& moleFractions(unsigned int pIdx, unsigned int compIdx) const { return xj_[pIdx][compIdx]; }
CellVector& moleFractions(unsigned int pIdx, unsigned int compIdx) { return xj_[pIdx][compIdx]; }
@@ -222,31 +143,6 @@ class HeatConductionDataHandle
using CellVector = typename MatVecTraits::CellVector;
public:
- //! prepare data handle for subsequent fill
- template< class InteractionVolume >
- void resize(const InteractionVolume& iv)
- {
- //! resize transmissibility matrix & temperature vector
- resizeMatrix(T_, iv.numFaces(), iv.numKnowns());
- resizeVector(Tj_, iv.numKnowns());
-
- //! resize outsideTij on surface grids
- using GridView = typename InteractionVolume::Traits::GridView;
- if (int(GridView::dimension) < int(GridView::dimensionworld))
- {
- outsideT_.resize(iv.numFaces());
-
- using LocalIndexType = typename InteractionVolume::Traits::IndexSet::LocalIndexType;
- for (LocalIndexType i = 0; i < iv.numFaces(); ++i)
- {
- const auto numNeighbors = iv.localScvf(i).neighboringLocalScvIndices().size() - 1;
- outsideT_[i].resize(numNeighbors);
- for (LocalIndexType j = 0; j < numNeighbors; ++j)
- resizeVector(outsideT_[i][j], iv.numKnowns());
- }
- }
- }
-
//! Access to the iv-wide temperatures
const CellVector& temperatures() const { return Tj_; }
CellVector& temperatures() { return Tj_; }
@@ -286,21 +182,7 @@ class InteractionVolumeDataHandle
: public AdvectionDataHandle<MVT, PT, PT::enableAdvection>
, public DiffusionDataHandle<MVT, PT, PT::enableMolecularDiffusion>
, public HeatConductionDataHandle<MVT, PT, PT::enableHeatConduction>
-{
- using AdvectionHandle = AdvectionDataHandle<MVT, PT, PT::enableAdvection>;
- using DiffusionHandle = DiffusionDataHandle<MVT, PT, PT::enableMolecularDiffusion>;
- using HeatConductionHandle = HeatConductionDataHandle<MVT, PT, PT::enableHeatConduction>;
-
-public:
- //! prepare data handles for subsequent fills
- template< class InteractionVolume >
- void resize(const InteractionVolume& iv)
- {
- AdvectionHandle::resize(iv);
- DiffusionHandle::resize(iv);
- HeatConductionHandle::resize(iv);
- }
-};
+{};
} // end namespace Dumux
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
index 82aa7b8273..151f6f49a7 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
@@ -111,6 +111,23 @@ public:
iv.B().leftmultiply(iv.A());
T += multiplyMatrices(iv.C(), iv.B());
+ // bring outsideT vector to the right size
+ outsideTij.resize(iv.numFaces());
+ using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
+ for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
+ {
+ // gravitational acceleration on this face
+ const auto& curLocalScvf = iv.localScvf(faceIdx);
+ const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex());
+ const auto numOutsideFaces = curGlobalScvf.boundary() ? 0 : curGlobalScvf.numOutsideScvs();
+
+ // resize each face entry to the right number of outside faces
+ outsideTij[faceIdx].resize(numOutsideFaces);
+ std::for_each(outsideTij[faceIdx].begin(),
+ outsideTij[faceIdx].end(),
+ [&iv](auto& v) { resizeVector(v, iv.numKnowns()); });
+ }
+
// compute outside transmissibilities
for (const auto& localFaceData : iv.localFaceData())
{
@@ -240,6 +257,23 @@ public:
A = iv.A();
CA = iv.C().rightmultiply(A);
+ // bring outsideT vector to the right size
+ outsideTij.resize(iv.numFaces());
+ using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
+ for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
+ {
+ // gravitational acceleration on this face
+ const auto& curLocalScvf = iv.localScvf(faceIdx);
+ const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex());
+ const auto numOutsideFaces = curGlobalScvf.boundary() ? 0 : curGlobalScvf.numOutsideScvs();
+
+ // resize each face entry to the right number of outside faces
+ outsideTij[faceIdx].resize(numOutsideFaces);
+ std::for_each(outsideTij[faceIdx].begin(),
+ outsideTij[faceIdx].end(),
+ [&iv](auto& v) { resizeVector(v, iv.numKnowns()); });
+ }
+
// compute outside transmissibilities
for (const auto& localFaceData : iv.localFaceData())
{
@@ -298,6 +332,7 @@ public:
void assemble(typename IV::Traits::MatVecTraits::CellVector& u, const IV& iv, const GetU& getU)
{
// put the cell pressures first
+ resizeVector(u, iv.numKnowns());
using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
for (LocalIndexType i = 0; i < iv.numScvs(); ++i)
u[i] = getU( iv.localScv(i).gridScvIndex() );
@@ -337,8 +372,10 @@ public:
//! - first index adresses the respective phases
//! - second index adresses the face within the interaction volume
- // make sure g vector and CA matrix have the correct sizes already
- assert(std::all_of(g.begin(), g.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); }));
+ // resize the g vector
+ std::for_each(g.begin(), g.end(), [&iv] (auto& v) { resizeVector(v, iv.numFaces()); });
+
+ // make sure CA matrix has the correct size already
assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns());
//! For each face, we...
@@ -469,10 +506,13 @@ public:
//! - second index adresses the face within the interaction volume
//! - third index adresses the i-th "outside" face of the current face
- // we require the CA matrix and the gravity containers to have the correct size already
+ // resize the g vector
+ auto resizeToNumFaces = [&iv] (auto& v) { resizeVector(v, iv.numFaces()); };
+ std::for_each(g.begin(), g.end(), resizeToNumFaces);
+ std::for_each(outsideG.begin(), outsideG.end(), resizeToNumFaces);
+
+ // we require the CA matrix to have the correct size already
assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns());
- assert(std::all_of(g.begin(), g.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); }));
- assert(std::all_of(outsideG.begin(), outsideG.end(), [&iv](const auto& v) { return v.size() == iv.numFaces(); }));
//! For each face, we...
//! - arithmetically average the phase densities
@@ -490,7 +530,6 @@ public:
resizeVector(sum_alphas[pIdx], iv.numUnknowns());
sum_alphas[pIdx] = 0.0;
g[pIdx] = 0.0;
- std::for_each(outsideG[pIdx].begin(), outsideG[pIdx].end(), [] (auto& v) { v = 0.0; });
}
for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
@@ -513,6 +552,13 @@ public:
std::vector< Scalar > alpha_outside(numOutsideFaces);
std::vector< Scalar > rho(numPhases);
+ // resize & reset entries in outside gravity vector
+ for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
+ {
+ outsideG[pIdx][faceIdx].resize(numOutsideFaces);
+ outsideG[pIdx][faceIdx] = 0.0;
+ }
+
if (!curLocalScvf.isDirichlet())
{
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
@@ -641,15 +687,12 @@ private:
static constexpr int dim = IV::Traits::GridView::dimension;
static constexpr int dimWorld = IV::Traits::GridView::dimensionworld;
- // Matrix D is assumed to have the right size already
- assert(D.rows() == iv.numFaces() && D.cols() == iv.numKnowns());
-
// if only Dirichlet faces are present in the iv,
// the matrices A, B & C are undefined and D = T
if (iv.numUnknowns() == 0)
{
- // reset matrix beforehand
- D = 0.0;
+ // resize & reset D matrix
+ resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0;
// Loop over all the faces, in this case these are all dirichlet boundaries
for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
@@ -680,16 +723,11 @@ private:
}
else
{
- // we require the matrices A,B,C to have the correct size already
- assert(A.rows() == iv.numUnknowns() && A.cols() == iv.numUnknowns());
- assert(B.rows() == iv.numUnknowns() && B.cols() == iv.numKnowns());
- assert(C.rows() == iv.numFaces() && C.cols() == iv.numUnknowns());
-
- // reset matrices
- A = 0.0;
- B = 0.0;
- C = 0.0;
- D = 0.0;
+ // resize & reset matrices
+ resizeMatrix(A, iv.numUnknowns(), iv.numUnknowns()); A = 0.0;
+ resizeMatrix(B, iv.numUnknowns(), iv.numKnowns()); B = 0.0;
+ resizeMatrix(C, iv.numFaces(), iv.numUnknowns()); C = 0.0;
+ resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0;
auto& wijk = iv.omegas();
for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
--
GitLab