Skip to content
Snippets Groups Projects
Commit d2fc21e3 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[mpfa][localassembly] put gravity assembly in base class 

This should be identical for all mpfa schemes
parent baa881eb
No related branches found
No related tags found
2 merge requests!1337WIP Fix/dirichlet caching v2,!1298Restructure mpfa flux caches
Hide whitespace changes
Inline Side-by-side
dumux/discretization/cellcentered/mpfa/localassembler.hh
+ 140
− 1
View file @ d2fc21e3
......@@ -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:
......
dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
+ 0
− 153
View file @ d2fc21e3
......@@ -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
......
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