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

[mpfa-l] port to new flux var cache concept

parent 962b60ac
No related branches found
No related tags found
2 merge requests!617[WIP] Next,!371port mpfa l method to new flux var cache concept
Hide whitespace changes
Inline Side-by-side
dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
+ 7
− 8
View file @ b1e657e7
......@@ -98,7 +98,6 @@ public:
const Matrix& M1,
const Matrix& M2)
{
Scalar eps = 1e-10;
Scalar t11, t12, t21, t22;
t11 = M1[I1.contiFaceLocalIdx][0];
t21 = M2[I2.contiFaceLocalIdx][0];
......@@ -116,7 +115,7 @@ public:
Scalar s1 = std::abs(t11-t12);
Scalar s2 = std::abs(t22-t21);
if (s1 < s2 + eps*s1)
if (s1 < s2)
return 0;
else
return 1;
......@@ -132,16 +131,16 @@ private:
const SubControlVolumeFace& scvf)
{
// make the first scv seed, we know this element will NOT be on the lowest local level
auto scvfVector = Implementation::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry);
auto localScvfIdx = Implementation::getLocalFaceIndex(scvf, scvfVector);
const auto scvfVector = Implementation::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry);
const auto localScvfIdx = Implementation::getLocalFaceIndex(scvf, scvfVector);
scvSeeds.emplace_back( GlobalIndexSet({scvfVector[0]->index(), scvfVector[1]->index()}),
scvf.insideScvIdx(),
localScvfIdx );
// get the surrounding elements and "outside" data
LocalIndexType otherScvfIdx = 1-localScvfIdx;
auto e2 = problem.model().globalFvGeometry().element(scvf.outsideScvIdx());
auto e3 = problem.model().globalFvGeometry().element(scvfVector[otherScvfIdx]->outsideScvIdx());
const auto e2 = problem.model().globalFvGeometry().element(scvf.outsideScvIdx());
const auto e3 = problem.model().globalFvGeometry().element(scvfVector[otherScvfIdx]->outsideScvIdx());
auto e2Geometry = localView(problem.model().globalFvGeometry());
auto e3Geometry = localView(problem.model().globalFvGeometry());
......@@ -166,10 +165,10 @@ private:
e3Scvfs[localScvfIdx]->index());
// Outer seed for outside of e2, we know the local scvf index here will be localScvfIdx in 2d
auto e4 = problem.model().globalFvGeometry().element(e2Scvfs[localScvfIdx]->outsideScvIdx());
const auto e4 = problem.model().globalFvGeometry().element(e2Scvfs[localScvfIdx]->outsideScvIdx());
auto e4Geometry = localView(problem.model().globalFvGeometry());
e4Geometry.bindElement(e4);
auto e4Scvfs = Implementation::getCommonAndNextScvFace(*e2Scvfs[localScvfIdx], e4Geometry, /*clockwise?*/localScvfIdx == 1);
const auto e4Scvfs = Implementation::getCommonAndNextScvFace(*e2Scvfs[localScvfIdx], e4Geometry, /*clockwise?*/localScvfIdx == 1);
outerScvSeeds.emplace_back(e4Scvfs[otherScvfIdx]->insideScvIdx(),
e4Scvfs[otherScvfIdx]->index());
......
dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh
+ 12
− 12
View file @ b1e657e7
......@@ -58,22 +58,22 @@ public:
std::vector<Scalar> detX;
std::vector<Element> elements;
std::array<GlobalIndexType, 2> globalScvfs;
std::array<GlobalIndexType, 2> globalScvfIndices;
// Constructor signature for dim == 2
template<class ScvSeedType, class OuterScvSeedType>
template<class ScvSeed, class OuterScvSeed>
InteractionRegion(const Problem& problem,
const FVElementGeometry& fvGeometry,
const ScvSeedType& scvSeed,
const OuterScvSeedType& outerSeed1,
const OuterScvSeedType& outerSeed2,
const ScvSeed& scvSeed,
const OuterScvSeed& outerSeed1,
const OuterScvSeed& outerSeed2,
const Element& element1,
const Element& element2,
const Element& element3)
{
contiFaceLocalIdx = scvSeed.contiFaceLocalIdx();
globalScvfs[0] = scvSeed.globalScvfIndices()[contiFaceLocalIdx];
globalScvfs[1] = contiFaceLocalIdx == 0 ? outerSeed1.globalScvfIndex() : outerSeed2.globalScvfIndex();
globalScvfIndices[0] = scvSeed.globalScvfIndices()[contiFaceLocalIdx];
globalScvfIndices[1] = contiFaceLocalIdx == 0 ? outerSeed1.globalScvfIndex() : outerSeed2.globalScvfIndex();
elements.resize(3);
elements[0] = element1;
......@@ -81,11 +81,11 @@ public:
elements[2] = element3;
// The participating sub control entities
auto&& scv1 = fvGeometry.scv(scvSeed.globalIndex());
auto&& scv2 = fvGeometry.scv(outerSeed1.globalIndex());
auto&& scv3 = fvGeometry.scv(outerSeed2.globalIndex());
auto&& scvf1 = fvGeometry.scvf(scvSeed.globalScvfIndices()[0]);
auto&& scvf2 = fvGeometry.scvf(scvSeed.globalScvfIndices()[1]);
const auto& scv1 = fvGeometry.scv(scvSeed.globalIndex());
const auto& scv2 = fvGeometry.scv(outerSeed1.globalIndex());
const auto& scv3 = fvGeometry.scv(outerSeed2.globalIndex());
const auto& scvf1 = fvGeometry.scvf(scvSeed.globalScvfIndices()[0]);
const auto& scvf2 = fvGeometry.scvf(scvSeed.globalScvfIndices()[1]);
// The necessary coordinates and normals
GlobalPosition v = scvf1.vertexCorner();
......
dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
+ 69
− 39
View file @ b1e657e7
......@@ -80,6 +80,7 @@ class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::lMethod> : pub
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using InteriorBoundaryData = typename GET_PROP_TYPE(TypeTag, InteriorBoundaryData);
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
......@@ -93,6 +94,7 @@ class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::lMethod> : pub
public:
using Matrix = typename Traits::Matrix;
using typename ParentType::GlobalLocalFaceDataPair;
using typename ParentType::LocalIndexType;
using typename ParentType::LocalIndexSet;
using typename ParentType::LocalFaceData;
......@@ -108,7 +110,8 @@ public:
const Problem& problem,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
: problemPtr_(&problem),
: seedPtr_(&seed),
problemPtr_(&problem),
fvGeometryPtr_(&fvGeometry),
elemVolVarsPtr_(&elemVolVars),
regionUnique_(seed.isUnique()),
......@@ -125,7 +128,9 @@ public:
{
const auto& ir = interactionRegions_[0];
T_ = assembleMatrix_(getTensor, ir);
postSolve_(ir);
interactionRegionId_ = 0;
updateGlobalLocalFaceData();
systemSolved_ = true;
}
else
{
......@@ -133,22 +138,25 @@ public:
M[0] = assembleMatrix_(getTensor, interactionRegions_[0]);
M[1] = assembleMatrix_(getTensor, interactionRegions_[1]);
auto id = MpfaHelper::selectionCriterion(interactionRegions_[0], interactionRegions_[1], M[0], M[1]);
postSolve_(interactionRegions_[id]);
T_ = M[id];
interactionRegionId_ = MpfaHelper::selectionCriterion(interactionRegions_[0], interactionRegions_[1], M[0], M[1]);
updateGlobalLocalFaceData();
systemSolved_ = true;
T_ = M[interactionRegionId_];
}
}
//! returns the local index pair. This holds the scvf's local index and a boolean whether or not the flux has to be inverted
LocalFaceData getLocalFaceData(const SubControlVolumeFace& scvf) const
const LocalFaceData getLocalFaceData(const SubControlVolumeFace& scvf) const
{
assert(systemSolved_ && "Scvf Indices not set yet. You have to call solveLocalSystem() beforehand.");
assert((scvf.index() == globalScvfIndices_[0] || scvf.index() == globalScvfIndices_[1]) && "The provided scvf is not the flux face of the interaction volume.");
assert((scvf.index() == globalLocalScvfPairedData_[0].first->index()
|| scvf.index() == globalLocalScvfPairedData_[1].first->index())
&& "The provided scvf is not the flux face of the interaction volume.");
if (globalScvfIndices_[0] == scvf.index())
return LocalFaceData(contiFaceLocalIdx_, /*dummy*/0, false);
if (globalLocalScvfPairedData_[0].first->index() == scvf.index())
return globalLocalScvfPairedData_[0].second;
else
return LocalFaceData(contiFaceLocalIdx_, /*dummy*/0, true);
return globalLocalScvfPairedData_[1].second;
}
//! returns the transmissibilities corresponding to the bound scvf
......@@ -164,21 +172,18 @@ public:
const GlobalIndexSet& volVarsStencil() const
{
assert(systemSolved_ && "volVarsStencil not set yet. You have to call solveLocalSystem() beforehand.");
return volVarsStencil_;
assert(systemSolved_ && "stencil not set yet. You have to call solveLocalSystem() beforehand.");
return interactionRegion().scvIndices;
}
const PositionVector& volVarsPositions() const
{
assert(systemSolved_ && "globalScvfs not set yet. You have to call solveLocalSystem() beforehand.");
return volVarsPositions_;
assert(systemSolved_ && "stencil not set yet. You have to call solveLocalSystem() beforehand.");
return interactionRegion().scvCenters;
}
const GlobalIndexSet& globalScvfs() const
{
assert(systemSolved_ && "globalScvfs not set yet. You have to call solveLocalSystem() beforehand.");
return globalScvfIndices_;
}
const std::vector<GlobalLocalFaceDataPair>& globalLocalScvfPairedData() const
{ return globalLocalScvfPairedData_; }
//! Boundaries will be treated by a different mpfa method (e.g. o method). Thus, on
//! faces in l-method interaction volumes there will never be a Neumann flux contribution.
......@@ -192,6 +197,13 @@ public:
const Matrix& matrix() const
{ return T_; }
const InteractionRegion& interactionRegion() const
{ return interactionRegions_[interactionRegionId_]; }
// The l-method can't treat interior boundaries
const std::vector<InteriorBoundaryData> interiorBoundaryData() const
{ return std::vector<InteriorBoundaryData>(); }
private:
//! Assembles and solves the local equation system
//! Specialization for dim = 2
......@@ -210,9 +222,9 @@ private:
const auto scv3 = fvGeometry_().scv(ir.scvIndices[2]);
// Get diffusion tensors in the three scvs
const auto T1 = getTensor(e1, elemVolVars_()[scv1], scv1);
const auto T2 = getTensor(e2, elemVolVars_()[scv2], scv2);
const auto T3 = getTensor(e3, elemVolVars_()[scv3], scv3);
const auto T1 = getTensor(problem_(), e1, elemVolVars_()[scv1], fvGeometry_(), scv1);
const auto T2 = getTensor(problem_(), e2, elemVolVars_()[scv2], fvGeometry_(), scv2);
const auto T3 = getTensor(problem_(), e3, elemVolVars_()[scv3], fvGeometry_(), scv3);
// required omega factors
Scalar w111 = calculateOmega_(ir.normal[0], ir.nu[0], ir.detX[0], T1);
......@@ -248,7 +260,7 @@ private:
B[1][1] = 0.0;
B[1][2] = -w235 - w236;
// T = CA^-1B + D
// T = C*A^-1*B + D
A.invert();
auto T = (A.leftmultiply(C)).rightmultiplyany(B);
T[0][0] += w111 + w112;
......@@ -274,8 +286,10 @@ private:
auto e1 = globalFvGeometry.element(scvSeed1.globalIndex());
auto e2 = globalFvGeometry.element(scvSeed2.globalIndex());
auto e3 = globalFvGeometry.element(scvSeed3.globalIndex());
interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, scvSeed2, scvSeed3, e1, e2, e3);
if (scvSeed1.contiFaceLocalIdx() == 0)
interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, scvSeed2, scvSeed3, e1, e2, e3);
else
interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed1, scvSeed3, scvSeed2, e1, e3, e2);
}
else
{
......@@ -301,17 +315,36 @@ private:
interactionRegions_.emplace_back(problem_(), fvGeometry, scvSeed2, OuterScvSeedType(scvSeed1), outerScvSeed2, e2, e1, e4);
}
}
// as an initial guess we set the first interaction region as the one used
// this is updated after solution of the local system
interactionRegionId_ = 0;
updateGlobalLocalFaceData();
}
void postSolve_(const InteractionRegion& ir)
//! set the global&local face data according to the current choice for the interaction region
void updateGlobalLocalFaceData()
{
globalScvfIndices_.resize(2);
globalScvfIndices_[0] = ir.globalScvfs[0];
globalScvfIndices_[1] = ir.globalScvfs[1];
volVarsStencil_ = ir.scvIndices;
volVarsPositions_ = ir.scvCenters;
contiFaceLocalIdx_ = ir.contiFaceLocalIdx;
systemSolved_ = true;
const auto numPairs = globalLocalScvfPairedData_.size();
// if no data has been stored yet, initialize
if (numPairs == 0)
{
globalLocalScvfPairedData_.clear();
globalLocalScvfPairedData_.reserve(2);
globalLocalScvfPairedData_.emplace_back(&fvGeometry_().scvf(seed_().globalScvfIndices()[0]),
LocalFaceData(interactionRegion().contiFaceLocalIdx, /*dummy*/0, false));
globalLocalScvfPairedData_.emplace_back(&fvGeometry_().scvf(seed_().globalScvfIndices()[1]),
LocalFaceData(interactionRegion().contiFaceLocalIdx, /*dummy*/0, true));
}
// if order was swapped, change bools indicating which face is the "outside" version of the local face
else if (globalLocalScvfPairedData_[1].first->index() == interactionRegion().globalScvfIndices[0])
{
globalLocalScvfPairedData_[0].second.localScvfIndex = interactionRegion().contiFaceLocalIdx;
globalLocalScvfPairedData_[1].second.localScvfIndex = interactionRegion().contiFaceLocalIdx;
globalLocalScvfPairedData_[0].second.isOutside = true;
globalLocalScvfPairedData_[1].second.isOutside = false;
}
}
Scalar calculateOmega_(const GlobalPosition& normal,
......@@ -347,7 +380,7 @@ private:
}
const Seed& seed_() const
{ return seedPtr_; }
{ return *seedPtr_; }
const Problem& problem_() const
{ return *problemPtr_; }
......@@ -366,12 +399,9 @@ private:
bool regionUnique_;
bool systemSolved_;
LocalIndexType contiFaceLocalIdx_;
GlobalIndexSet globalScvfIndices_;
GlobalIndexSet volVarsStencil_;
PositionVector volVarsPositions_;
LocalIndexType interactionRegionId_;
std::vector<InteractionRegion> interactionRegions_;
std::vector<GlobalLocalFaceDataPair> globalLocalScvfPairedData_;
Matrix T_;
};
......
dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh
+ 2
− 2
View file @ b1e657e7
......@@ -61,8 +61,8 @@ public:
Scalar q,
bool boundary)
: ParentType(helper,
std::move(corners),
std::move(unitOuterNormal),
std::forward<std::vector<GlobalPosition>>(corners),
std::forward<GlobalPosition>(unitOuterNormal),
vertexIndex,
localIndex,
scvfIndex,
......
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