From d2fc21e3d7812be0c28ba5e8b5ed776e5e4d3c6e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= <dennis.glaeser@iws.uni-stuttgart.de>
Date: Mon, 19 Nov 2018 12:35:12 +0100
Subject: [PATCH] [mpfa][localassembly] put gravity assembly in base class
This should be identical for all mpfa schemes
---
.../cellcentered/mpfa/localassembler.hh | 141 +++++++++++++++-
.../mpfa/omethod/localassembler.hh | 153 ------------------
2 files changed, 140 insertions(+), 154 deletions(-)
diff --git a/dumux/discretization/cellcentered/mpfa/localassembler.hh b/dumux/discretization/cellcentered/mpfa/localassembler.hh
index 2dc51bb868..c3247c6b44 100644
--- a/dumux/discretization/cellcentered/mpfa/localassembler.hh
+++ b/dumux/discretization/cellcentered/mpfa/localassembler.hh
@@ -156,7 +156,146 @@ class InteractionVolumeAssemblerBase
template< class DataHandle, class IV, class GetRho >
void assembleGravity(DataHandle& handle, const IV& iv, const GetRho& getRho)
{
- DUNE_THROW(Dune::NotImplemented, "Implementation does not provide a assembleGravity() function");
+ using GridView = typename IV::Traits::GridView;
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+ static constexpr bool isSurfaceGrid = dim < dimWorld;
+
+ // resize the gravity vectors
+ auto& g = handle.g();
+ auto& outsideG = handle.gOutside();
+ resizeVector_(g, iv.numFaces());
+ if (isSurfaceGrid)
+ resizeVector_(outsideG, iv.numFaces());
+
+ // we require the CA matrix to have the correct size already
+ assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns());
+
+ //! For each face, we...
+ //! - arithmetically average the phase densities
+ //! - compute the term \f$ \alpha := \mathbf{A} \rho \ \mathbf{n}^T \mathbf{K} \mathbf{g} \f$ in each neighboring cell
+ //! - compute \f$ \alpha^* = \sum{\alpha_{outside, i}} - \alpha_{inside} \f$
+ using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type;
+ using FaceVector = typename IV::Traits::MatVecTraits::FaceVector;
+ using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
+
+ FaceVector sum_alphas;
+ resizeVector_(sum_alphas, iv.numUnknowns());
+ sum_alphas = 0.0;
+ std::fill(g.begin(), g.end(), 0.0);
+ for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
+ {
+ // gravitational acceleration on this face
+ const auto& curLocalScvf = iv.localScvf(faceIdx);
+ const auto& curGlobalScvf = fvGeometry().scvf(curLocalScvf.gridScvfIndex());
+ const auto& gravity = problem().gravityAtPos(curGlobalScvf.ipGlobal());
+
+ // get permeability tensor in "positive" sub volume
+ const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
+ const auto& posGlobalScv = fvGeometry().scv(iv.localScv(neighborScvIndices[0]).gridScvIndex());
+ const auto& posVolVars = elemVolVars()[posGlobalScv];
+ const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(),
+ posVolVars.permeability(),
+ gravity);
+
+ const auto numOutsideFaces = !curGlobalScvf.boundary() ? curGlobalScvf.numOutsideScvs() : 0;
+ using OutsideAlphaStorage = std::conditional_t< isSurfaceGrid,
+ std::vector<Scalar>,
+ Dune::ReservedVector<Scalar, 1> >;
+ OutsideAlphaStorage alpha_outside; alpha_outside.resize(numOutsideFaces);
+ std::fill(alpha_outside.begin(), alpha_outside.end(), 0.0);
+ Scalar rho;
+
+ if (isSurfaceGrid)
+ {
+ resizeVector_(outsideG[faceIdx], numOutsideFaces);
+ std::fill(outsideG[faceIdx].begin(), outsideG[faceIdx].end(), 0.0);
+ }
+
+ if (!curLocalScvf.isDirichlet())
+ {
+ rho = getRho(posVolVars);
+
+ // arithmetically average density on inside faces
+ const auto localDofIdx = curLocalScvf.localDofIndex();
+ if (!curGlobalScvf.boundary())
+ {
+ for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside)
+ {
+ // obtain outside tensor
+ const auto negLocalScvIdx = neighborScvIndices[idxInOutside+1];
+ const auto& negGlobalScv = fvGeometry().scv(iv.localScv(negLocalScvIdx).gridScvIndex());
+ const auto& negVolVars = elemVolVars()[negGlobalScv];
+ const auto& flipScvf = !isSurfaceGrid ? curGlobalScvf
+ : fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside);
+
+ alpha_outside[idxInOutside] = negVolVars.extrusionFactor()*vtmv(flipScvf.unitOuterNormal(),
+ negVolVars.permeability(),
+ gravity);
+ if (isSurfaceGrid)
+ alpha_outside[idxInOutside] *= -1.0;
+
+ rho += getRho(negVolVars);
+ sum_alphas[localDofIdx] += alpha_outside[idxInOutside];
+ }
+ }
+
+ rho /= numOutsideFaces + 1;
+ sum_alphas[localDofIdx] -= alpha_inside;
+ sum_alphas[localDofIdx] *= rho*curGlobalScvf.area();
+ }
+ // use density resulting from Dirichlet BCs
+ else
+ rho = getRho(elemVolVars()[curGlobalScvf.outsideScvIdx()]);
+
+ // add "inside" & "outside" alphas to gravity containers
+ g[faceIdx] += alpha_inside*rho*curGlobalScvf.area();
+
+ if (isSurfaceGrid)
+ {
+ unsigned int i = 0;
+ for (const auto& alpha : alpha_outside)
+ outsideG[faceIdx][i++] += alpha*rho*curGlobalScvf.area();
+ }
+ }
+
+ // g += CA*sum_alphas
+ // outsideG = wikj*A^-1*sum_alphas + outsideG
+ handle.CA().umv(sum_alphas, g);
+
+ if (isSurfaceGrid)
+ {
+ using FaceVector = typename IV::Traits::MatVecTraits::FaceVector;
+ FaceVector AG;
+ resizeVector_(AG, iv.numUnknowns());
+ handle.A().mv(sum_alphas, AG);
+
+ // compute gravitational accelerations
+ for (const auto& localFaceData : iv.localFaceData())
+ {
+ // continue only for "outside" faces
+ if (!localFaceData.isOutsideFace())
+ continue;
+
+ const auto localScvIdx = localFaceData.ivLocalInsideScvIndex();
+ const auto localScvfIdx = localFaceData.ivLocalScvfIndex();
+ const auto idxInOutside = localFaceData.scvfLocalOutsideScvfIndex();
+ const auto& posLocalScv = iv.localScv(localScvIdx);
+ const auto& wijk = iv.omegas()[localScvfIdx][idxInOutside+1];
+
+ // make sure the given outside gravity container has the right size
+ assert(outsideG[localScvfIdx].size() == iv.localScvf(localScvfIdx).neighboringLocalScvIndices().size()-1);
+
+ // add contributions from all local directions
+ for (LocalIndexType localDir = 0; localDir < dim; localDir++)
+ {
+ // the scvf corresponding to this local direction in the scv
+ const auto& curLocalScvf = iv.localScvf(posLocalScv.localScvfIndex(localDir));
+ if (!curLocalScvf.isDirichlet())
+ outsideG[localScvfIdx][idxInOutside] -= wijk[localDir]*AG[curLocalScvf.localDofIndex()];
+ }
+ }
+ }
}
protected:
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
index 10d1fa3435..4828fc0cff 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
@@ -165,159 +165,6 @@ public:
u[i++] = getU( this->elemVolVars()[data.volVarIndex()] );
}
- /*!
- * \brief Assembles the gravitational flux contributions on the scvfs within an mpfa-o
- * interaction volume.
- *
- * \param handle The data handle in which the vector is stored
- * \param iv The mpfa-o interaction volume
- * \param getRho Lambda to obtain the density from volume variables
- */
- template< class DataHandle, class IV, class GetRho >
- void assembleGravity(DataHandle& handle, const IV& iv, const GetRho& getRho)
- {
- using GridView = typename IV::Traits::GridView;
- static constexpr int dim = GridView::dimension;
- static constexpr int dimWorld = GridView::dimensionworld;
- static constexpr bool isSurfaceGrid = dim < dimWorld;
-
- // resize the gravity vectors
- auto& g = handle.g();
- auto& outsideG = handle.gOutside();
- this->resizeVector_(g, iv.numFaces());
- if (isSurfaceGrid)
- this->resizeVector_(outsideG, iv.numFaces());
-
- // we require the CA matrix to have the correct size already
- assert(CA.rows() == iv.numFaces() && CA.cols() == iv.numUnknowns());
-
- //! For each face, we...
- //! - arithmetically average the phase densities
- //! - compute the term \f$ \alpha := \mathbf{A} \rho \ \mathbf{n}^T \mathbf{K} \mathbf{g} \f$ in each neighboring cell
- //! - compute \f$ \alpha^* = \sum{\alpha_{outside, i}} - \alpha_{inside} \f$
- using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type;
- using FaceVector = typename IV::Traits::MatVecTraits::FaceVector;
- using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
-
- FaceVector sum_alphas;
- this->resizeVector_(sum_alphas, iv.numUnknowns());
- sum_alphas = 0.0;
- std::fill(g.begin(), g.end(), 0.0);
- 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& gravity = this->problem().gravityAtPos(curGlobalScvf.ipGlobal());
-
- // get permeability tensor in "positive" sub volume
- const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
- const auto& posGlobalScv = this->fvGeometry().scv(iv.localScv(neighborScvIndices[0]).gridScvIndex());
- const auto& posVolVars = this->elemVolVars()[posGlobalScv];
- const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(),
- posVolVars.permeability(),
- gravity);
-
- const auto numOutsideFaces = !curGlobalScvf.boundary() ? curGlobalScvf.numOutsideScvs() : 0;
- using OutsideAlphaStorage = std::conditional_t< isSurfaceGrid,
- std::vector<Scalar>,
- Dune::ReservedVector<Scalar, 1> >;
- OutsideAlphaStorage alpha_outside; alpha_outside.resize(numOutsideFaces);
- std::fill(alpha_outside.begin(), alpha_outside.end(), 0.0);
- Scalar rho;
-
- if (isSurfaceGrid)
- {
- this->resizeVector_(outsideG[faceIdx], numOutsideFaces);
- std::fill(outsideG[faceIdx].begin(), outsideG[faceIdx].end(), 0.0);
- }
-
- if (!curLocalScvf.isDirichlet())
- {
- rho = getRho(posVolVars);
-
- // arithmetically average density on inside faces
- const auto localDofIdx = curLocalScvf.localDofIndex();
- if (!curGlobalScvf.boundary())
- {
- for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside)
- {
- // obtain outside tensor
- const auto negLocalScvIdx = neighborScvIndices[idxInOutside+1];
- const auto& negGlobalScv = this->fvGeometry().scv(iv.localScv(negLocalScvIdx).gridScvIndex());
- const auto& negVolVars = this->elemVolVars()[negGlobalScv];
- const auto& flipScvf = !isSurfaceGrid ? curGlobalScvf
- : this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside);
-
- alpha_outside[idxInOutside] = negVolVars.extrusionFactor()*vtmv(flipScvf.unitOuterNormal(),
- negVolVars.permeability(),
- gravity);
- if (isSurfaceGrid)
- alpha_outside[idxInOutside] *= -1.0;
-
- rho += getRho(negVolVars);
- sum_alphas[localDofIdx] += alpha_outside[idxInOutside];
- }
- }
-
- rho /= numOutsideFaces + 1;
- sum_alphas[localDofIdx] -= alpha_inside;
- sum_alphas[localDofIdx] *= rho*curGlobalScvf.area();
- }
- // use density resulting from Dirichlet BCs
- else
- rho = getRho(this->elemVolVars()[curGlobalScvf.outsideScvIdx()]);
-
- // add "inside" & "outside" alphas to gravity containers
- g[faceIdx] += alpha_inside*rho*curGlobalScvf.area();
-
- if (isSurfaceGrid)
- {
- unsigned int i = 0;
- for (const auto& alpha : alpha_outside)
- outsideG[faceIdx][i++] += alpha*rho*curGlobalScvf.area();
- }
- }
-
- // g += CA*sum_alphas
- // outsideG = wikj*A^-1*sum_alphas + outsideG
- handle.CA().umv(sum_alphas, g);
-
- if (isSurfaceGrid)
- {
- using FaceVector = typename IV::Traits::MatVecTraits::FaceVector;
- FaceVector AG;
- this->resizeVector_(AG, iv.numUnknowns());
- handle.A().mv(sum_alphas, AG);
-
- // compute gravitational accelerations
- for (const auto& localFaceData : iv.localFaceData())
- {
- // continue only for "outside" faces
- if (!localFaceData.isOutsideFace())
- continue;
-
- const auto localScvIdx = localFaceData.ivLocalInsideScvIndex();
- const auto localScvfIdx = localFaceData.ivLocalScvfIndex();
- const auto idxInOutside = localFaceData.scvfLocalOutsideScvfIndex();
- const auto& posLocalScv = iv.localScv(localScvIdx);
- const auto& wijk = iv.omegas()[localScvfIdx][idxInOutside+1];
-
- // make sure the given outside gravity container has the right size
- assert(outsideG[localScvfIdx].size() == iv.localScvf(localScvfIdx).neighboringLocalScvIndices().size()-1);
-
- // add contributions from all local directions
- for (LocalIndexType localDir = 0; localDir < dim; localDir++)
- {
- // the scvf corresponding to this local direction in the scv
- const auto& curLocalScvf = iv.localScvf(posLocalScv.localScvfIndex(localDir));
- if (!curLocalScvf.isDirichlet())
- outsideG[localScvfIdx][idxInOutside] -= wijk[localDir]*AG[curLocalScvf.localDofIndex()];
- }
- }
- }
- }
-
private:
/*!
* \brief Assemble the matrices involved in the flux expressions
--
GitLab