From b1e657e7a8317dd8d02f4b2ff1a1079a1b049c23 Mon Sep 17 00:00:00 2001
From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de>
Date: Fri, 3 Feb 2017 18:47:36 +0100
Subject: [PATCH] [mpfa-l] port to new flux var cache concept
---
.../cellcentered/mpfa/lmethod/helper.hh | 15 ++-
.../mpfa/lmethod/interactionregions.hh | 24 ++--
.../mpfa/lmethod/interactionvolume.hh | 108 +++++++++++-------
.../mpfa/lmethod/subcontrolvolumeface.hh | 4 +-
4 files changed, 90 insertions(+), 61 deletions(-)
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh b/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
index f5d88a6897..1ca199a5d2 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
@@ -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());
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh
index fc4c039ea1..17380bc1a9 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionregions.hh
@@ -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();
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
index 51e7929705..e56cb825c5 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
@@ -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_;
};
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh
index 13e90f8edc..9468f8c2fb 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/subcontrolvolumeface.hh
@@ -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,
--
GitLab