From a6e8512bca9409833256074ee120abda0f6fd02e Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Tue, 10 May 2011 16:20:51 +0000
Subject: [PATCH] make test_impes work in parallel (kind of)
there seem to be issues with the accuracy of the linear solver,
i.e. the residual reduction needs to be quite large...
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@5788 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
.../decoupled/2p/diffusion/fv/fvpressure2p.hh | 9 -
dumux/linear/borderindex.hh | 5 +
dumux/linear/boxbicgstabilu0solver.hh | 5 +-
dumux/linear/domesticoverlapfrombcrsmatrix.hh | 13 +-
dumux/linear/elementborderlistfromgrid.hh | 39 ++-
dumux/linear/foreignoverlapfrombcrsmatrix.hh | 227 +++++++-----------
dumux/linear/globalindices.hh | 117 +++++----
dumux/linear/impetbicgstabilu0solver.hh | 17 +-
dumux/linear/overlappingbcrsmatrix.hh | 106 +++++++-
dumux/linear/overlappingblockvector.hh | 11 +-
dumux/linear/seqsolverbackend.hh | 2 +-
dumux/linear/vertexborderlistfromgrid.hh | 12 +-
test/decoupled/2p/test_impes.cc | 20 +-
13 files changed, 358 insertions(+), 225 deletions(-)
diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
index b77cc1e32b..a96a958580 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
@@ -348,15 +348,6 @@ void FVPressure2P<TypeTag>::assemble(bool first)
// cell index
int globalIdxI = problem_.variables().index(*eIt);
-#if HAVE_MPI
- if (eIt->partitionType() != Dune::InteriorEntity)
- {
- A_[globalIdxI][globalIdxI]=1.0;
- f_[globalIdxI]=0;
- continue;
- }
-#endif
-
// get global coordinate of cell center
const GlobalPosition& globalPos = eIt->geometry().center();
diff --git a/dumux/linear/borderindex.hh b/dumux/linear/borderindex.hh
index 02bf61860a..8646a130da 100644
--- a/dumux/linear/borderindex.hh
+++ b/dumux/linear/borderindex.hh
@@ -41,6 +41,11 @@ struct BorderIndex
//! Distance to the process border for the peer (in hops)
int borderDistance;
+
+ //! True if and only if the entity which corrosponds to the index
+ //! is in the interior of more than one process (i.e. the entity
+ //! is completely on the border)
+ bool isShared;
};
};
diff --git a/dumux/linear/boxbicgstabilu0solver.hh b/dumux/linear/boxbicgstabilu0solver.hh
index 4e5acab45e..13c17cbd68 100644
--- a/dumux/linear/boxbicgstabilu0solver.hh
+++ b/dumux/linear/boxbicgstabilu0solver.hh
@@ -66,7 +66,7 @@ class BoxBiCGStabILU0Solver
typedef Dune::BiCGSTABSolver<OverlappingVector> Solver;
public:
- BoxBiCGStabILU0Solver(const Problem &problem, int overlapSize=1)
+ BoxBiCGStabILU0Solver(const Problem &problem, int overlapSize=3)
: problem_(problem)
, overlapSize_(overlapSize)
{
@@ -152,7 +152,8 @@ private:
// create the overlapping Jacobian matrix
overlapMatrix_ = new OverlappingMatrix (M,
- borderListCreator.borderList(),
+ borderListCreator.foreignBorderList(),
+ borderListCreator.domesticBorderList(),
overlapSize_);
// create the overlapping vectors for the residual and the
diff --git a/dumux/linear/domesticoverlapfrombcrsmatrix.hh b/dumux/linear/domesticoverlapfrombcrsmatrix.hh
index 1cd64366b0..5859f4e17e 100644
--- a/dumux/linear/domesticoverlapfrombcrsmatrix.hh
+++ b/dumux/linear/domesticoverlapfrombcrsmatrix.hh
@@ -86,8 +86,9 @@ public:
*/
DomesticOverlapFromBCRSMatrix(const BCRSMatrix &A,
const BorderList &borderList,
+ const BorderList &domesticBorderList,
int overlapSize)
- : foreignOverlap_(A, borderList, overlapSize)
+ : foreignOverlap_(A, borderList, domesticBorderList, overlapSize)
, globalIndices_(foreignOverlap_)
{
myRank_ = 0;
@@ -257,6 +258,16 @@ public:
return foreignOverlap_.masterOf(domesticIdx) == myRank_;
};
+ /*!
+ * \brief Return true iff a given index is shared by more than one process
+ */
+ bool isShared(int domesticIdx) const
+ {
+ if (!isLocal(domesticIdx))
+ return false;
+ return foreignOverlap_.isShared(domesticIdx);
+ };
+
/*!
* \brief Return true iff a given rank is the master of a given
* domestic index.
diff --git a/dumux/linear/elementborderlistfromgrid.hh b/dumux/linear/elementborderlistfromgrid.hh
index 67b1fd3a01..045cd4cd19 100644
--- a/dumux/linear/elementborderlistfromgrid.hh
+++ b/dumux/linear/elementborderlistfromgrid.hh
@@ -68,9 +68,14 @@ public:
const ElementMapper &map)
: gv_(gv), map_(map)
{
+ forwardComm_ = true;
+ gv.communicate(*this,
+ Dune::InteriorBorder_All_Interface,
+ Dune::ForwardCommunication);
+
+ forwardComm_ = false;
gv.communicate(*this,
Dune::InteriorBorder_All_Interface,
-// Dune::Overlap_All_Interface,
Dune::BackwardCommunication);
};
@@ -100,19 +105,39 @@ public:
bIdx.localIdx = map_.map(e);
buff.read(bIdx.peerRank);
buff.read(bIdx.peerIdx);
- bIdx.borderDistance = 1;
+ bIdx.borderDistance = 0;
+ // this index is not shared between at least two processes
+ // because an element is always in the interior of exactly one
+ // process
+ bIdx.isShared = false;
- borderList_.push_back(bIdx);
+ if (forwardComm_) {
+ // add the index to the list of indices which we intersect
+ // with remote processes.
+ domesticBorderList_.push_back(bIdx);
+ }
+ else {
+ // add the index to the list of remote processes
+ // intersecting with us.
+ foreignBorderList_.push_back(bIdx);
+ }
};
- // Access to the initial seed list.
- const BorderList &borderList() const
- { return borderList_; }
+ // Access to the initial border list.
+ const BorderList &foreignBorderList() const
+ { return foreignBorderList_; }
+
+ // Access to the indices which are on the border but which are in
+ // the interior of a remote process
+ const BorderList &domesticBorderList() const
+ { return domesticBorderList_; }
private:
const GridView gv_;
const ElementMapper &map_;
- BorderList borderList_;
+ bool forwardComm_;
+ BorderList foreignBorderList_;
+ BorderList domesticBorderList_;
};
} // namespace Dumux
diff --git a/dumux/linear/foreignoverlapfrombcrsmatrix.hh b/dumux/linear/foreignoverlapfrombcrsmatrix.hh
index f07728eef9..fa5414a74f 100644
--- a/dumux/linear/foreignoverlapfrombcrsmatrix.hh
+++ b/dumux/linear/foreignoverlapfrombcrsmatrix.hh
@@ -85,9 +85,11 @@ public:
* an initial list of border indices.
*/
ForeignOverlapFromBCRSMatrix(const BCRSMatrix &A,
- const BorderList &borderList,
+ const BorderList &foreignBorderList,
+ const BorderList &domesticBorderList,
int overlapSize)
- : borderList_(borderList)
+ : foreignBorderList_(foreignBorderList)
+ , domesticBorderList_(domesticBorderList)
{
overlapSize_ = overlapSize;
@@ -101,7 +103,7 @@ public:
// calculate the border list. From this, create an initial
// seed list of indices which are in the overlap.
SeedList initialSeedList;
- borderListToSeedList_(initialSeedList, borderList);
+ borderListToSeedList_(initialSeedList, foreignBorderList);
// find the set of processes which have an overlap with the
// local processes. (i.e. the set of processes which we will
@@ -112,89 +114,24 @@ public:
// i.e. find the distance of each row from the seed set.
foreignOverlapByIndex_.resize(A.N());
extendForeignOverlap_(A, initialSeedList, overlapSize);
+
+ updateMasterRanks_(foreignBorderList, domesticBorderList);
// group foreign overlap by peer process rank
groupForeignOverlapByRank_();
}
- /*!
- * \brief Return the border list.
- *
- * The border list is the list of (local index, peer index, peer
- * rank) triples for all indices on a process border.
- */
- const BorderList& borderList() const
- { return borderList_; };
-
/*!
* \brief Returns true iff a local index is a border index.
*/
bool isBorder(int localIdx) const
{ return borderIndices_.count(localIdx) > 0; };
- /*!
- * \brief Returns true iff a local index is a border index with a given peer.
- */
- bool isBorderFor(int peerRank, int localIdx) const
- {
- typedef std::map<ProcessRank, BorderDistance> BorderDistMap;
- const BorderDistMap &borderDist = foreignOverlapByIndex_[localIdx];
- BorderDistMap::const_iterator bdIt = borderDist.find(peerRank);
-
- if (bdIt == borderDist.end())
- return false; // this index is not seen by the peer
-
- BorderDistance bDist = bdIt->second;
- if (bDist == 0)
- // the index is on the border to the peer
- return true;
-
- // at this point, the local index is in the interior of the
- // foreign overlap for the peer, so it is a remote index
- return false;
- };
-
- /*!
- * \brief Returns true iff a local index is a front index for a given peer.
- */
- bool isFrontFor(int peerRank, int localIdx) const
- {
- typedef std::map<ProcessRank, BorderDistance> BorderDistMap;
- const BorderDistMap &borderDist = foreignOverlapByIndex_[localIdx];
- BorderDistMap::const_iterator bdIt = borderDist.find(peerRank);
-
- if (bdIt == borderDist.end())
- return false; // this index is not seen by the peer
-
- BorderDistance bDist = bdIt->second;
- if (bDist == overlapSize_)
- // the index is on the front of the peer
- return true;
- return false;
- };
-
/*!
* \brief Return the rank of the master process for a local index.
*/
int masterOf(int localIdx) const
- {
- if (!isBorder(localIdx))
- return myRank_; // interior index
-
- // if the local index is a border index, loop over all ranks
- // for which this index is also a border index. the lowest
- // rank wins!
- typedef typename std::map<ProcessRank, BorderDistance>::const_iterator iterator;
- iterator it = foreignOverlapByIndex_[localIdx].begin();
- iterator endIt = foreignOverlapByIndex_[localIdx].end();
- LocalIndex masterIdx = myRank_;
- for (; it != endIt; ++it) {
- if (it->second == 0)
- masterIdx = std::min(masterIdx, it->first);
- }
-
- return masterIdx;
- };
+ { return masterRank_[localIdx]; };
/*!
* \brief Return true if the current rank is the "master" of an
@@ -207,70 +144,21 @@ public:
* rank.
*/
bool iAmMasterOf(int localIdx) const
- {
- if (!isBorder(localIdx))
- return true; // interior index
-
- // if the local index is a border index, loop over all ranks
- // for which this index is also a border index. the lowest
- // rank wins!
- typedef typename std::map<ProcessRank, BorderDistance>::const_iterator iterator;
- iterator it = foreignOverlapByIndex_[localIdx].begin();
- iterator endIt = foreignOverlapByIndex_[localIdx].end();
- LocalIndex masterIdx = myRank_;
- for (; it != endIt; ++it) {
- if (it->first < myRank_ && it->second == 0)
- return false;
- }
-
- return masterIdx == myRank_;
- };
+ { return masterRank_[localIdx] == myRank_; };
/*!
- * \brief Return true if a given local index is a remote index for
- * a peer.
+ * \brief Returns the list of indices which intersect the process
+ * border and are in the interior of the local process.
*/
- bool isRemoteIndexFor(ProcessRank peerRank, Index localIdx) const
- {
- typedef std::map<ProcessRank, BorderDistance> BorderDistMap;
- const BorderDistMap &borderDist = foreignOverlapByIndex_[localIdx];
- BorderDistMap::const_iterator bdIt = borderDist.find(peerRank);
-
- if (bdIt == borderDist.end())
- return false; // this index is not seen by the peer
-
- BorderDistance bDist = bdIt->second;
- if (bDist == 0)
- // the index is on the border to the peer, so for the peer
- // it is a local index.
- return false;
-
- // at this point, the local index is in the interior of the
- // foreign overlap for the peer, so it is a remote index
- return true;
- };
+ const BorderList &foreignBorderList() const
+ { return foreignBorderList_; }
/*!
- * \brief Return true if a given local index is also a local index
- * for a peer.
- */
- bool isLocalIndexFor(ProcessRank peerRank, Index localIdx) const
- {
- if (!isLocal(localIdx))
- // our own remote indices do not count!
- return false;
-
- typedef std::map<ProcessRank, BorderDistance> BorderDistMap;
- const BorderDistMap &borderDist = foreignOverlapByIndex_[localIdx];
- BorderDistMap::const_iterator bdIt = borderDist.find(peerRank);
- if (bdIt == borderDist.end())
- return false; // this index is not seen by the peer
-
- // the index is also local for the peer if it an index on the
- // border.
- BorderDistance bDist = bdIt->second;
- return bDist == 0;
- };
+ * \brief Returns the list of indices which intersect the process
+ * border and are in the interior of some remote process.
+ */
+ const BorderList &domesticBorderList() const
+ { return domesticBorderList_; }
/*!
* \brief Return true if a given local index is a domestic index
@@ -329,7 +217,15 @@ public:
* \brief Returns true iff a domestic index is local
*/
bool isLocal(int domesticIdx) const
- { return domesticIdx < numLocal(); };
+ {
+ return domesticIdx < numLocal();
+ };
+
+ /*!
+ * \brief Returns true iff a local index is shared with an other process.
+ */
+ bool isShared(int domesticIdx) const
+ { return isLocal(domesticIdx) && isShared_[domesticIdx]; }
/*!
* \brief Return the number of peer ranks for which a given local
@@ -382,11 +278,60 @@ protected:
initialSeedList.push_back(IndexRankDist(it->localIdx,
it->peerRank,
it->borderDistance));
- if (it->borderDistance == 0)
- borderIndices_.insert(it->localIdx);
+ borderIndices_.insert(it->localIdx);
};
};
+ // given a list of border indices and provided that
+ // borderListToSeedList_() was already called, calculate the
+ // master process of each local index.
+ void updateMasterRanks_(const BorderList &foreignBorderList,
+ const BorderList &domesticBorderList)
+ {
+ // determine the minimum rank for all indices
+ masterRank_.resize(numLocal_, myRank_);
+ for (int localIdx = 0; localIdx < masterRank_.size(); ++localIdx) {
+ int masterRank = myRank_;
+ if (isBorder(localIdx)) {
+ // if the local index is a border index, loop over all ranks
+ // for which this index is also a border index. the lowest
+ // rank wins!
+ typedef typename std::map<ProcessRank, BorderDistance>::const_iterator iterator;
+ iterator it = foreignOverlapByIndex_[localIdx].begin();
+ iterator endIt = foreignOverlapByIndex_[localIdx].end();
+ for (; it != endIt; ++it) {
+ if (it->second == 0) {
+ // if the border distance is zero, the rank with the minimum
+ masterRank = std::min(masterRank, it->first);
+ }
+ }
+ }
+ masterRank_[localIdx] = masterRank;
+ }
+
+ // overwrite the master rank of the non-shared border indices
+ isShared_.resize(numLocal_, false);
+ BorderList::const_iterator it = foreignBorderList.begin();
+ BorderList::const_iterator endIt = foreignBorderList.end();
+ for (; it != endIt; ++it) {
+ if (!it->isShared)
+ masterRank_[it->localIdx] = myRank_;
+ else
+ isShared_[it->localIdx] = true;
+ };
+
+ // overwrite the master rank of the non-shared border on the
+ // domestic overlap
+ it = domesticBorderList_.begin();
+ endIt = domesticBorderList_.end();
+ for (; it != endIt; ++it) {
+ if (!it->isShared) {
+ masterRank_[it->localIdx] = it->peerRank;
+ }
+ };
+ }
+
+
// extend the foreign overlaps by one level. this uses a greedy
// algorithm.
void extendForeignOverlap_(const BCRSMatrix &A,
@@ -502,11 +447,21 @@ protected:
PeerSet peerSet_;
// the list of indices on the border
- const BorderList &borderList_;
+ const BorderList &foreignBorderList_;
+ const BorderList &domesticBorderList_;
+
+ // an array which contains the rank of the master process for each
+ // index
+ std::vector<ProcessRank> masterRank_;
+
+ // an array which stores whether an index is also in the interior
+ // of some other process
+ std::vector<bool> isShared_;
- // set of all local indices which are on the border
+ // set of all local indices which are on the border of some remote
+ // process
std::set<LocalIndex> borderIndices_;
-
+
// stores the set of process ranks which are in the overlap for a
// given row index "owned" by the current rank. The second value
// store the distance from the nearest process border.
diff --git a/dumux/linear/globalindices.hh b/dumux/linear/globalindices.hh
index a5acd134dc..960b61d6dd 100644
--- a/dumux/linear/globalindices.hh
+++ b/dumux/linear/globalindices.hh
@@ -205,11 +205,8 @@ public:
*/
void print() const
{
- int myRank = 0;
-#if HAVE_MPI
- MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
-#endif // HAVE_MPI
- std::cout << "(domestic index, global index, domestic->global->domestic) list for rank " << myRank << "\n";
+ std::cout << "(domestic index, global index, domestic->global->domestic) list for rank " <<
+ myRank_ << "\n";
for (int domIdx = 0; domIdx < domesticToGlobal_.size(); ++ domIdx) {
std::cout << "(" << domIdx
@@ -272,61 +269,87 @@ protected:
0, // tag
MPI_COMM_WORLD); // communicator
};
+
+ typename PeerSet::const_iterator peerIt;
+ typename PeerSet::const_iterator peerEndIt = peerSet_().end();
+ // receive the border indices from the lower ranks
+ peerIt = peerSet_().begin();
+ for (; peerIt != peerEndIt; ++peerIt) {
+ if (*peerIt < myRank_)
+ receiveBorderFrom_(*peerIt);
+ }
- // retrieve the global indices for which we are not master
- // from the processes with lower rank
- PeerSet::const_iterator peerIt = peerSet_().begin();
- PeerSet::const_iterator peerEndIt = peerSet_().end();
+ // send the border indices to the higher ranks
+ peerIt = peerSet_().begin();
for (; peerIt != peerEndIt; ++peerIt) {
- int peerRank = *peerIt;
- if (peerRank > myRank_)
- continue; // ignore processes with higher rank
-
- // receive (local index on myRank, global index) pairs and
- // update maps
- BorderList::const_iterator borderIt = borderList_().begin();
- BorderList::const_iterator borderEndIt = borderList_().end();
- for (; borderIt != borderEndIt; ++borderIt) {
- int borderPeer = borderIt->peerRank;
- if (borderPeer != peerRank)
- continue;
- if (borderIt->borderDistance == 0)
- receiveBorderIndex(peerRank);
- }
+ if (*peerIt > myRank_)
+ sendBorderTo_(*peerIt);
}
- // send the global indices for which we are master
- // to the processes with higher rank
+ // receive the border indices from the higher ranks
peerIt = peerSet_().begin();
- peerEndIt = peerSet_().end();
for (; peerIt != peerEndIt; ++peerIt) {
- int peerRank = *peerIt;
- if (peerRank < myRank_)
- continue; // ignore processes with lower rank
-
- // send (local index on myRank, global index) pairs to the
- // peers
- BorderList::const_iterator borderIt = borderList_().begin();
- BorderList::const_iterator borderEndIt = borderList_().end();
- for (; borderIt != borderEndIt; ++borderIt) {
- int borderPeer = borderIt->peerRank;
- if (borderPeer != peerRank)
- continue;
-
- int localIdx = borderIt->localIdx;
- int peerIdx = borderIt->peerIdx;
- if (borderIt->borderDistance == 0)
- sendBorderIndex(peerRank, localIdx, peerIdx);
- }
+ if (*peerIt > myRank_)
+ receiveBorderFrom_(*peerIt);
+ }
+
+ // send the border indices to the lower ranks
+ peerIt = peerSet_().begin();
+ for (; peerIt != peerEndIt; ++peerIt) {
+ if (*peerIt < myRank_)
+ sendBorderTo_(*peerIt);
}
#endif // HAVE_MPI
}
+ void sendBorderTo_(ProcessRank peerRank)
+ {
+#if HAVE_MPI
+ // send (local index on myRank, global index) pairs to the
+ // peers
+ BorderList::const_iterator borderIt = foreignBorderList_().begin();
+ BorderList::const_iterator borderEndIt = foreignBorderList_().end();
+ for (; borderIt != borderEndIt; ++borderIt) {
+ int borderPeer = borderIt->peerRank;
+ if (borderPeer != peerRank)
+ continue;
+
+ int localIdx = borderIt->localIdx;
+ int peerIdx = borderIt->peerIdx;
+ if (foreignOverlap_.iAmMasterOf(borderIt->localIdx)) {
+ sendBorderIndex(borderPeer, localIdx, peerIdx);
+ }
+ }
+#endif // HAVE_MPI
+ }
+
+ void receiveBorderFrom_(ProcessRank peerRank)
+ {
+#if HAVE_MPI
+ // retrieve the global indices for which we are not master
+ // from the processes with lower rank
+ BorderList::const_iterator borderIt = domesticBorderList_().begin();
+ BorderList::const_iterator borderEndIt = domesticBorderList_().end();
+ for (; borderIt != borderEndIt; ++borderIt) {
+ int borderPeer = borderIt->peerRank;
+ if (borderPeer != peerRank)
+ continue;
+
+ if (foreignOverlap_.masterOf(borderIt->localIdx) == borderPeer) {
+ receiveBorderIndex(borderPeer);
+ }
+ }
+#endif // HAVE_MPI
+ };
+
const PeerSet &peerSet_() const
{ return foreignOverlap_.peerSet(); }
- const BorderList &borderList_() const
- { return foreignOverlap_.borderList(); }
+ const BorderList &foreignBorderList_() const
+ { return foreignOverlap_.foreignBorderList(); }
+
+ const BorderList &domesticBorderList_() const
+ { return foreignOverlap_.domesticBorderList(); }
int myRank_;
diff --git a/dumux/linear/impetbicgstabilu0solver.hh b/dumux/linear/impetbicgstabilu0solver.hh
index d371610872..0320347368 100644
--- a/dumux/linear/impetbicgstabilu0solver.hh
+++ b/dumux/linear/impetbicgstabilu0solver.hh
@@ -115,9 +115,9 @@ public:
const Vector &b)
{
int verbosityLevel = GET_PROP_VALUE(TypeTag, PTAG(LSVerbosity));
- static const int maxIter = GET_PROP_VALUE(TypeTag, PTAG(LSMaxIterations));
- static const double residReduction = GET_PROP_VALUE(TypeTag, PTAG(LSResidualReduction));
- static const double relaxation = GET_PROP_VALUE(TypeTag, PTAG(PreconditionerRelaxation));
+ static const int maxIter = GET_PROP_VALUE(TypeTag, PTAG(LSMaxIterations));
+ static const double residReduction = GET_PROP_VALUE(TypeTag, PTAG(LSResidualReduction));
+ static const double relaxation = GET_PROP_VALUE(TypeTag, PTAG(PreconditionerRelaxation));
if (!overlapMatrix_) {
// make sure that the overlapping matrix and block vectors
@@ -128,10 +128,16 @@ public:
// copy the values of the non-overlapping linear system of
// equations to the overlapping one. On ther border, we add up
// the values of all processes (using the assignAdd() methods)
- overlapMatrix_->assignAdd(M);
+ overlapMatrix_->assignCopy(M);
overlapb_->assignAdd(b);
(*overlapx_) = 0.0;
+ /*
+ overlapMatrix_->print();
+ overlapb_->print();
+ exit(1);
+ */
+
// create sequential and overlapping preconditioners
//SeqPreconditioner seqPreCond(*overlapMatrix_, 1, 1.0);
SeqPreconditioner seqPreCond(*overlapMatrix_, relaxation);
@@ -168,7 +174,8 @@ private:
// create the overlapping Jacobian matrix
overlapMatrix_ = new OverlappingMatrix (M,
- borderListCreator.borderList(),
+ borderListCreator.foreignBorderList(),
+ borderListCreator.domesticBorderList(),
overlapSize_);
// create the overlapping vectors for the residual and the
diff --git a/dumux/linear/overlappingbcrsmatrix.hh b/dumux/linear/overlappingbcrsmatrix.hh
index 4e35c90713..5c90d020ab 100644
--- a/dumux/linear/overlappingbcrsmatrix.hh
+++ b/dumux/linear/overlappingbcrsmatrix.hh
@@ -77,10 +77,11 @@ public:
};
OverlappingBCRSMatrix(const BCRSMatrix &M,
- const BorderList &borderList,
+ const BorderList &foreignBorderList,
+ const BorderList &domesticBorderList,
int overlapSize)
{
- overlap_ = std::shared_ptr<Overlap>(new Overlap(M, borderList, overlapSize));
+ overlap_ = std::shared_ptr<Overlap>(new Overlap(M, foreignBorderList, domesticBorderList, overlapSize));
myRank_ = 0;
#if HAVE_MPI
MPI_Comm_rank(MPI_COMM_WORLD, &myRank_);
@@ -150,7 +151,40 @@ public:
}
// communicate and add the contents of overlapping rows
- sync_();
+ syncAdd_();
+ }
+
+ /*!
+ * \brief Assign and syncronize the overlapping matrix from a
+ * non-overlapping one.
+ *
+ * The non-master entries are copied from the master
+ */
+ void assignCopy(const BCRSMatrix &M)
+ {
+ // first, set everything to 0
+ BCRSMatrix::operator=(0.0);
+
+ // assign the local rows
+ for (int rowIdx = 0; rowIdx < M.N(); ++rowIdx) {
+ ConstColIterator colIt = M[rowIdx].begin();
+ ConstColIterator colEndIt = M[rowIdx].end();
+ ColIterator myColIt = (*this)[rowIdx].begin();
+ for (; colIt != colEndIt; ++colIt) {
+ while (true) {
+ if (myColIt.index() == colIt.index())
+ break;
+
+ ++ myColIt;
+ }
+ assert(myColIt.index() == colIt.index());
+
+ (*myColIt) = *colIt;
+ }
+ }
+
+ // communicate and add the contents of overlapping rows
+ syncCopy_();
}
void print() const
@@ -319,8 +353,9 @@ private:
ColIt colEndIt = M[rowIdx].end();
int j = 0;
for (; colIt != colEndIt; ++colIt) {
- if (overlap_->isDomesticIndexFor(peerRank, colIt.index()))
+ if (overlap_->isDomesticIndexFor(peerRank, colIt.index())) {
++ j;
+ }
}
(*rowSizesSendBuff_[peerRank])[i] = j;
@@ -409,7 +444,7 @@ private:
};
// communicates and adds up the contents of overlapping rows
- void sync_()
+ void syncAdd_()
{
// first, send all entries to the peers
const PeerSet &peerSet = overlap_->foreignOverlap().peerSet();
@@ -438,6 +473,37 @@ private:
}
}
+ // communicates and copies the contents of overlapping rows from
+ // the master
+ void syncCopy_()
+ {
+ // first, send all entries to the peers
+ const PeerSet &peerSet = overlap_->foreignOverlap().peerSet();
+ typename PeerSet::const_iterator peerIt = peerSet.begin();
+ typename PeerSet::const_iterator peerEndIt = peerSet.end();
+ for (; peerIt != peerEndIt; ++peerIt) {
+ int peerRank = *peerIt;
+
+ sendEntries_(peerRank);
+ }
+
+ // then, receive entries from the peers
+ peerIt = peerSet.begin();
+ for (; peerIt != peerEndIt; ++peerIt) {
+ int peerRank = *peerIt;
+
+ receiveCopyEntries_(peerRank);
+ }
+
+ // finally, make sure that everything which we send was
+ // received by the peers
+ peerIt = peerSet.begin();
+ for (; peerIt != peerEndIt; ++peerIt) {
+ int peerRank = *peerIt;
+ entryValuesSendBuff_[peerRank]->wait();
+ }
+ }
+
void sendEntries_(int peerRank)
{
#if HAVE_MPI
@@ -494,6 +560,36 @@ private:
#endif // HAVE_MPI
}
+ void receiveCopyEntries_(int peerRank)
+ {
+#if HAVE_MPI
+ MpiBuffer<block_type> &mpiRecvBuff = *entryValuesRecvBuff_[peerRank];
+
+ MpiBuffer<int> &mpiRowIndicesRecvBuff = *rowIndicesRecvBuff_[peerRank];
+ MpiBuffer<int> &mpiRowSizesRecvBuff = *rowSizesRecvBuff_[peerRank];
+ MpiBuffer<int> &mpiColIndicesRecvBuff = *entryIndicesRecvBuff_[peerRank];
+
+ mpiRecvBuff.receive(peerRank);
+
+ // retrieve the values from the receive buffer
+ int k = 0;
+ for (int i = 0; i < mpiRowIndicesRecvBuff.size(); ++i) {
+ Index domRowIdx = mpiRowIndicesRecvBuff[i];
+ for (int j = 0; j < mpiRowSizesRecvBuff[i]; ++j) {
+ Index domColIdx = mpiColIndicesRecvBuff[k];
+
+ if (!overlap_->iAmMasterOf(domRowIdx) ||
+ !overlap_->iAmMasterOf(domColIdx))
+ {
+ (*this)[domRowIdx][domColIdx] = mpiRecvBuff[k];
+ }
+
+ ++ k;
+ }
+ }
+#endif // HAVE_MPI
+ }
+
void globalToDomesticBuff_(MpiBuffer<Index> &idxBuff)
{
for (int i = 0; i < idxBuff.size(); ++i) {
diff --git a/dumux/linear/overlappingblockvector.hh b/dumux/linear/overlappingblockvector.hh
index 888b7866a7..0111d4714d 100644
--- a/dumux/linear/overlappingblockvector.hh
+++ b/dumux/linear/overlappingblockvector.hh
@@ -82,12 +82,17 @@ public:
*/
void assignAdd(const BlockVector &nbv)
{
- Valgrind::SetDefined(nbv[7]);
-
// assign the local rows
int numLocal = overlap_->numLocal();
for (int rowIdx = 0; rowIdx < numLocal; ++rowIdx) {
- (*this)[rowIdx] = nbv[rowIdx];
+ if (overlap_->iAmMasterOf(rowIdx) ||
+ overlap_->isShared(rowIdx))
+ {
+ (*this)[rowIdx] = nbv[rowIdx];
+ }
+ else
+ (*this)[rowIdx] = 0;
+
}
// set the remote indices to 0
diff --git a/dumux/linear/seqsolverbackend.hh b/dumux/linear/seqsolverbackend.hh
index c47f1cb3a2..e8fb037851 100644
--- a/dumux/linear/seqsolverbackend.hh
+++ b/dumux/linear/seqsolverbackend.hh
@@ -59,7 +59,7 @@ SET_PROP_DEFAULT(LSVerbosity)
SET_PROP_DEFAULT(LSResidualReduction)
{public:
- static constexpr double value = 1e-6;
+ static constexpr double value = 1e-13;
};
SET_PROP_DEFAULT(LSMaxIterations)
diff --git a/dumux/linear/vertexborderlistfromgrid.hh b/dumux/linear/vertexborderlistfromgrid.hh
index c8c0e71275..0716e88ee1 100644
--- a/dumux/linear/vertexborderlistfromgrid.hh
+++ b/dumux/linear/vertexborderlistfromgrid.hh
@@ -98,12 +98,20 @@ public:
buff.read(bIdx.peerRank);
buff.read(bIdx.peerIdx);
bIdx.borderDistance = 0;
+ // vertices on the border are always in the interior of more
+ // than one process which means that they are shared.
+ bIdx.isShared = true;
borderList_.push_back(bIdx);
};
- // Access to the initial seed list.
- const BorderList &borderList() const
+ // Access to the foreign border list.
+ const BorderList &foreignBorderList() const
+ { return borderList_; }
+
+ // Access to the domestic border list (same as foreign border list
+ // because all vertices are shared entities)
+ const BorderList &domesticBorderList() const
{ return borderList_; }
private:
diff --git a/test/decoupled/2p/test_impes.cc b/test/decoupled/2p/test_impes.cc
index 5e2e910f73..d2ca40caa9 100644
--- a/test/decoupled/2p/test_impes.cc
+++ b/test/decoupled/2p/test_impes.cc
@@ -91,17 +91,23 @@ int main(int argc, char** argv)
////////////////////////////////////////////////////////////
// create the grid
////////////////////////////////////////////////////////////
- Dune::FieldVector<int,dim> N(6); N[0] = 30;
- Dune::FieldVector<double ,dim> L(0);
- Dune::FieldVector<double,dim> H(60); H[0] = 300;
+ Dune::FieldVector<int,dim> numCells;
+ numCells[0] = 30;
+ numCells[1] = 6;
+ Dune::FieldVector<double,dim> lowerLeft;
+ lowerLeft[0] = 0;
+ lowerLeft[1] = 0;
+ Dune::FieldVector<double,dim> upperRight;
+ upperRight[0] = 300;
+ upperRight[1] = 60;
Dune::FieldVector<bool,dim> periodic(false);
- Grid grid(Dune::MPIHelper::getCommunicator(), H, N, periodic, /*overlap=*/1);
-
+
+ Grid grid(Dune::MPIHelper::getCommunicator(), upperRight, numCells, periodic, /*overlap=*/1);
+ grid.loadBalance();
////////////////////////////////////////////////////////////
// instantiate and run the concrete problem
////////////////////////////////////////////////////////////
-
- Problem problem(grid.leafView(), L, H);
+ Problem problem(grid.leafView(), lowerLeft, upperRight);
// load restart file if necessarry
if (restart)
--
GitLab