From 7f81c57043323c1efef4a244b5681e8121b7ddd7 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Thu, 25 Oct 2012 08:39:57 +0000
Subject: [PATCH] Fixed some issues for removing periodic boundary conditions.
---
AMDiS/src/ProblemStat.cc | 2 +
AMDiS/src/io/MacroReader.cc | 4 +-
AMDiS/src/parallel/ElementObjectDatabase.cc | 33 +--
AMDiS/src/parallel/ElementObjectDatabase.h | 24 +-
AMDiS/src/parallel/InteriorBoundary.cc | 249 +++++++++++++-----
AMDiS/src/parallel/InteriorBoundary.h | 20 ++
AMDiS/src/parallel/MeshDistributor.cc | 4 +-
AMDiS/src/parallel/MeshDistributor.h | 4 +-
.../src/parallel/ParallelCoarseSpaceMatVec.h | 7 +
AMDiS/src/parallel/ParallelProblemStatBase.h | 4 +-
AMDiS/src/parallel/PetscProblemStat.cc | 9 +-
AMDiS/src/parallel/PetscSolverFeti.cc | 180 +++++--------
AMDiS/src/parallel/PetscSolverFeti.h | 10 +-
AMDiS/src/parallel/PetscSolverFetiMonitor.cc | 2 +-
.../src/parallel/PetscSolverFetiOperators.cc | 16 +-
AMDiS/src/parallel/PetscSolverFetiStructs.h | 2 -
AMDiS/src/parallel/PetscSolverGlobalMatrix.cc | 2 +-
17 files changed, 323 insertions(+), 249 deletions(-)
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index 4c6d4471..e37ed74d 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/AMDiS/src/ProblemStat.cc
@@ -505,6 +505,7 @@ namespace AMDiS {
{
FUNCNAME("ProblemStat::createSolver()");
+#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
// === create solver ===
string solverType("0");
string initFileStr = name + "->solver";
@@ -515,6 +516,7 @@ namespace AMDiS {
solverCreator->setName(initFileStr);
solver = solverCreator->create();
solver->initParameters();
+#endif
}
diff --git a/AMDiS/src/io/MacroReader.cc b/AMDiS/src/io/MacroReader.cc
index 6e6e4b7a..0f079bf2 100644
--- a/AMDiS/src/io/MacroReader.cc
+++ b/AMDiS/src/io/MacroReader.cc
@@ -48,9 +48,7 @@ namespace AMDiS {
DegreeOfFreedom **dof = macroInfo->dof;
// === read periodic data =================================
- if (periodicFile != "") {
- WARNING("Periodic boundaries may lead to errors in small meshes if element neighbours are not set!\n");
-
+ if (periodicFile != "") {
FILE *file = fopen(periodicFile.c_str(), "r");
TEST_EXIT(file)("can't open file %s\n", periodicFile.c_str());
diff --git a/AMDiS/src/parallel/ElementObjectDatabase.cc b/AMDiS/src/parallel/ElementObjectDatabase.cc
index 9ec7f81f..38b147d2 100644
--- a/AMDiS/src/parallel/ElementObjectDatabase.cc
+++ b/AMDiS/src/parallel/ElementObjectDatabase.cc
@@ -105,16 +105,15 @@ namespace AMDiS {
periodicEdgeAssoc[edge0].insert(edge1);
// Add both vertices of the edge to be periodic.
- periodicVertices[make_pair(edge0.first, edge1.first)] = boundaryType;
- periodicVertices[make_pair(edge0.second, edge1.second)] = boundaryType;
+ DegreeOfFreedom mappedDof0 =
+ mesh->getPeriodicAssociations(boundaryType)[edge0.first];
+ DegreeOfFreedom mappedDof1 =
+ mesh->getPeriodicAssociations(boundaryType)[edge0.second];
+ periodicVertices[make_pair(edge0.first, mappedDof0)] = boundaryType;
+ periodicVertices[make_pair(edge0.second, mappedDof1)] = boundaryType;
+
periodicDofAssoc[edge0.first].insert(boundaryType);
- periodicDofAssoc[edge0.second].insert(boundaryType);
-
- TEST_EXIT_DBG(edge0.first ==
- mesh->getPeriodicAssociations(boundaryType)[edge1.first] &&
- edge0.second ==
- mesh->getPeriodicAssociations(boundaryType)[edge1.second])
- ("Should not happen!\n");
+ periodicDofAssoc[edge0.second].insert(boundaryType);
}
}
break;
@@ -134,12 +133,16 @@ namespace AMDiS {
periodicFaces[make_pair(face0, face1)] = elInfo->getBoundary(i);
/// Add all three vertices of the face to be periodic.
- periodicVertices[make_pair(face0.get<0>(), face1.get<0>())] =
- boundaryType;
- periodicVertices[make_pair(face0.get<1>(), face1.get<1>())] =
- boundaryType;
- periodicVertices[make_pair(face0.get<2>(), face1.get<2>())] =
- boundaryType;
+ DegreeOfFreedom mappedDof0 =
+ mesh->getPeriodicAssociations(boundaryType)[face0.get<0>()];
+ DegreeOfFreedom mappedDof1 =
+ mesh->getPeriodicAssociations(boundaryType)[face0.get<1>()];
+ DegreeOfFreedom mappedDof2 =
+ mesh->getPeriodicAssociations(boundaryType)[face0.get<2>()];
+
+ periodicVertices[make_pair(face0.get<0>(), mappedDof0)] = boundaryType;
+ periodicVertices[make_pair(face0.get<1>(), mappedDof1)] = boundaryType;
+ periodicVertices[make_pair(face0.get<2>(), mappedDof2)] = boundaryType;
periodicDofAssoc[face0.get<0>()].insert(boundaryType);
periodicDofAssoc[face0.get<1>()].insert(boundaryType);
diff --git a/AMDiS/src/parallel/ElementObjectDatabase.h b/AMDiS/src/parallel/ElementObjectDatabase.h
index 69fb5224..f781c429 100644
--- a/AMDiS/src/parallel/ElementObjectDatabase.h
+++ b/AMDiS/src/parallel/ElementObjectDatabase.h
@@ -309,24 +309,32 @@ namespace AMDiS {
return faceElements[face];
}
-
-
- /// Returns a map that maps to each rank all macro elements in this rank that
- /// have a given vertex DOF in common.
+ /// Returns, for a given vertex, a map that maps from rank numbers to
+ /// element data objects, which identify on the rank the element which
+ /// contains this vertex. If more than one element in one subdomain contains
+ /// the vertex, the element with the highest element index is given. If the
+ /// vertex is not contained in a rank's subdomain, it will not be considered
+ /// in this mapping.
map<int, ElementObjectData>& getElementsInRank(DegreeOfFreedom vertex)
{
return vertexInRank[vertex];
}
- /// Returns a map that maps to each rank all macro elements in this rank that
- /// have a given edge in common.
+ /// Returns, for a given edge, a map that maps from rank numbers to
+ /// element data objects, which identify on the rank the element which
+ /// contains this edge. If more than one element in one subdomain contains
+ /// the edge, the element with the highest element index is given. If the
+ /// edge is not contained in a rank's subdomain, it will not be considered
+ /// in this mapping.
map<int, ElementObjectData>& getElementsInRank(DofEdge edge)
{
return edgeInRank[edge];
}
- /// Returns a map that maps to each rank all macro elements in this rank that
- /// have a given face in common.
+ /// Returns, for a given face, a map that maps from rank numbers to
+ /// element data objects, which identify on the rank the element which
+ /// contains this face. If the face is not contained in a rank's subdomain,
+ /// it will not be considered in this mapping.
map<int, ElementObjectData>& getElementsInRank(DofFace face)
{
return faceInRank[face];
diff --git a/AMDiS/src/parallel/InteriorBoundary.cc b/AMDiS/src/parallel/InteriorBoundary.cc
index feb0f5cd..24243fb0 100644
--- a/AMDiS/src/parallel/InteriorBoundary.cc
+++ b/AMDiS/src/parallel/InteriorBoundary.cc
@@ -110,7 +110,7 @@ namespace AMDiS {
bound.neighObj.ithObj = it2->second.ithObject;
bound.type = INTERIOR;
-
+
AtomicBoundary& b = getNewOwn(it2->first);
b = bound;
if (geoIndex == EDGE)
@@ -181,6 +181,7 @@ namespace AMDiS {
bound.neighObj.subObj = VERTEX;
bound.neighObj.ithObj = perDofEl1.ithObject;
+
if (removePeriodicBoundary) {
bound.type = INTERIOR;
@@ -191,47 +192,56 @@ namespace AMDiS {
int owner = std::max(elObjDb.getOwner(it->first.first, level),
elObjDb.getOwner(it->first.second, level));
+
+ // Consider the following configuration of four elements partitioned
+ // to four ranks:
+ //
+ // |--------|--------|(*)
+ // | /| /|
+ // | 3 / | 0 / |
+ // | / | / |
+ // | / | / |
+ // | / | / |
+ // | / 2 | / 1 |
+ // | / | / |
+ // |/-------|/-------|
+ //
+ // We are in now interested in vertex (*). For the first, rank 1 owns
+ // the interior boundary with rank 0 on this vertex. When we remove
+ // periodic boundaries, which are applied on the left and right outer
+ // boundary, rank 3 becomes the owner of this vertex. Rank 0 and rank 1
+ // will have an interior boundary with rank 3 on this vertex, but rank
+ // 0 and rank 1 have no more interior boundaries on vertex (*). Thus we
+ // have to search and remove for this scenario.
+
+ if (owner != globalMpiRank) {
+ removeOwn(bound.rankObj);
+ removeOther(bound.rankObj);
+ }
+
+
+ // And than we can add the boundary to either own or other part of the
+ // interior database.
+
if (owner == globalMpiRank) {
- for (map<int, ElementObjectData>::iterator it2 = objData.begin();
- it2 != objData.end(); ++it2) {
- if (it2->first == globalMpiRank)
- continue;
-
- if (!allRanks && levelRanks.count(it2->first) == 0)
- continue;
-
- AtomicBoundary& b = getNewOwn(it2->first);
- b = bound;
- b.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
- b.neighObj.elIndex = it2->second.elIndex;
- b.neighObj.elType = elObjDb.getElementType(it2->second.elIndex);
- b.neighObj.ithObj = it2->second.ithObject;
- }
+ AtomicBoundary& b = getNewOwn(otherElementRank);
+ b = bound;
} else {
+ ElementObjectData& ownerBoundEl = objData[owner];
+ bound.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
+ bound.neighObj.elIndex = ownerBoundEl.elIndex;
+ bound.neighObj.elType = -1;
+ bound.neighObj.ithObj = ownerBoundEl.ithObject;
+
AtomicBoundary& b = getNewOther(owner);
b = bound;
-
- ElementObjectData& ownerBoundEl = objData[owner];
- b.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
- b.neighObj.elIndex = ownerBoundEl.elIndex;
- b.neighObj.elType = -1;
- b.neighObj.ithObj = ownerBoundEl.ithObject;
}
+
} else {
bound.type = it->second;
- if (MeshDistributor::sebastianMode) {
- if (bound.rankObj.elIndex == 3 &&
- bound.rankObj.ithObj == 1 ||
- bound.rankObj.elIndex == 78 &&
- bound.rankObj.ithObj == 2) {
- AtomicBoundary& b = getNewPeriodic(otherElementRank);
- b = bound;
- }
- } else {
- AtomicBoundary& b = getNewPeriodic(otherElementRank);
- b = bound;
- }
+ AtomicBoundary& b = getNewPeriodic(otherElementRank);
+ b = bound;
}
}
}
@@ -276,24 +286,15 @@ namespace AMDiS {
elObjDb.getOwner(it->first.second, level));
ElementObjectData& rankBoundEl = objData[globalMpiRank];
+ // See comments in the same part of code for VERTEX
+ if (owner != globalMpiRank) {
+ removeOwn(bound.rankObj);
+ removeOther(bound.rankObj);
+ }
+
if (owner == globalMpiRank) {
- for (map<int, ElementObjectData>::iterator it2 = objData.begin();
- it2 != objData.end(); ++it2) {
- if (it2->first == globalMpiRank)
- continue;
-
- if (!allRanks && levelRanks.count(it2->first) == 0)
- continue;
-
- AtomicBoundary& b = getNewOwn(it2->first);
- b = bound;
- b.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
- b.neighObj.elIndex = it2->second.elIndex;
- b.neighObj.elType = elObjDb.getElementType(it2->second.elIndex);
- b.neighObj.ithObj = it2->second.ithObject;
- b.neighObj.reverseMode =
- elObjDb.getEdgeReverseMode(rankBoundEl, it2->second);
- }
+ AtomicBoundary& b = getNewOwn(otherElementRank);
+ b = bound;
} else {
AtomicBoundary& b = getNewOther(owner);
b = bound;
@@ -368,23 +369,8 @@ namespace AMDiS {
ElementObjectData& rankBoundEl = objData[globalMpiRank];
if (owner == globalMpiRank) {
- for (map<int, ElementObjectData>::iterator it2 = objData.begin();
- it2 != objData.end(); ++it2) {
- if (it2->first == globalMpiRank)
- continue;
-
- if (!allRanks && levelRanks.count(it2->first) == 0)
- continue;
-
- AtomicBoundary& b = getNewOwn(it2->first);
- b = bound;
- b.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
- b.neighObj.elIndex = it2->second.elIndex;
- b.neighObj.elType = elObjDb.getElementType(it2->second.elIndex);
- b.neighObj.ithObj = it2->second.ithObject;
- b.neighObj.reverseMode =
- elObjDb.getFaceReverseMode(rankBoundEl, it2->second);
- }
+ AtomicBoundary& b = getNewOwn(otherElementRank);
+ b = bound;
} else {
AtomicBoundary& b = getNewOther(owner);
b = bound;
@@ -413,6 +399,7 @@ namespace AMDiS {
}
}
+
// === Once we have this information, we must care about the order of the ===
// === atomic bounds in the three boundary handling object. Eventually ===
// === all the boundaries have to be in the same order on both ranks that ===
@@ -454,19 +441,25 @@ namespace AMDiS {
// If the expected object is not at place, search for it.
- BoundaryObject &recvedBound =
+ BoundaryObject &receivedBound =
stdMpi.getRecvData()[rank][j].rankObj;
- if ((rankIt->second)[j].neighObj != recvedBound) {
+ if ((rankIt->second)[j].neighObj != receivedBound) {
unsigned int k = j + 1;
- for (; k < rankIt->second.size(); k++)
- if ((rankIt->second)[k].neighObj == recvedBound)
+ for (; k < rankIt->second.size(); k++) {
+ if ((rankIt->second)[k].neighObj == receivedBound)
break;
+ }
// The element must always be found, because the list is just in
// another order.
- TEST_EXIT_DBG(k < rankIt->second.size())("Should never happen!\n");
+ TEST_EXIT_DBG(k < rankIt->second.size())
+ ("Cannot find element %d/%d/%d received from rank %d!\n",
+ receivedBound.elIndex,
+ receivedBound.subObj,
+ receivedBound.ithObj,
+ rank);
// Swap the current with the found element.
AtomicBoundary tmpBound = (rankIt->second)[k];
@@ -531,6 +524,11 @@ namespace AMDiS {
}
}
} // periodicBoundary.boundary.size() > 0
+
+
+ // === Run verification procedure. ===
+
+ verifyBoundary();
}
@@ -644,8 +642,37 @@ namespace AMDiS {
}
+ void InteriorBoundary::verifyBoundary()
+ {
+ FUNCNAME("InteriorBoundary::verifyBoundary()");
+
+ // === Check if no other boundery rank object is also included in own ===
+ // === boundary part, which would make no sence. ===
+
+ std::set<BoundaryObject> allOwnBounds;
+ for (map<int, vector<AtomicBoundary> >::iterator it = own.begin();
+ it != own.end(); ++it)
+ for (vector<AtomicBoundary>::iterator bIt = it->second.begin();
+ bIt != it->second.end(); ++bIt)
+ allOwnBounds.insert(bIt->rankObj);
+
+ for (map<int, vector<AtomicBoundary> >::iterator it = other.begin();
+ it != other.end(); ++it)
+ for (vector<AtomicBoundary>::iterator bIt = it->second.begin();
+ bIt != it->second.end(); ++bIt)
+ if (allOwnBounds.count(bIt->rankObj)) {
+ ERROR_EXIT("Boundary %d/%d/%d is included in both, own and other interior boundary part!\n",
+ bIt->rankObj.elIndex,
+ bIt->rankObj.subObj,
+ bIt->rankObj.ithObj);
+ }
+ }
+
+
AtomicBoundary& InteriorBoundary::getNewOwn(int rank)
{
+ FUNCNAME("InteriorBoundary::getNewOwn()");
+
int size = own[rank].size();
own[rank].resize(size + 1);
return own[rank][size];
@@ -654,6 +681,8 @@ namespace AMDiS {
AtomicBoundary& InteriorBoundary::getNewOther(int rank)
{
+ FUNCNAME("InteriorBoundary::getNewOther()");
+
int size = other[rank].size();
other[rank].resize(size + 1);
return other[rank][size];
@@ -662,12 +691,86 @@ namespace AMDiS {
AtomicBoundary& InteriorBoundary::getNewPeriodic(int rank)
{
+ FUNCNAME("InteriorBoundary::getNewPeriodic()");
+
int size = periodic[rank].size();
periodic[rank].resize(size + 1);
return periodic[rank][size];
}
+ bool InteriorBoundary::checkOther(AtomicBoundary& bound, int rank)
+ {
+ FUNCNAME("InteriorBoundary::checkOther()");
+
+ int globalMpiRank = MPI::COMM_WORLD.Get_rank();
+ TEST_EXIT(rank > globalMpiRank)
+ ("Wrong MPI ranks: %d %d\n", rank, globalMpiRank);
+
+ bool found = false;
+
+ if (other.count(rank)) {
+ vector<AtomicBoundary>& otherBounds = other[rank];
+ for (int i = 0; i < static_cast<int>(otherBounds.size()); i++) {
+ if (otherBounds[i].rankObj == bound.rankObj) {
+ if (otherBounds[i].neighObj != bound.neighObj) {
+ ERROR_EXIT("If this really can happen, I have to thing about this!\n");
+ } else {
+ found = true;
+ break;
+ }
+ }
+ }
+ }
+
+ return found;
+ }
+
+
+ bool InteriorBoundary::removeOwn(BoundaryObject& bound)
+ {
+ FUNCNAME("InteriorBoundary::removeOwn()");
+
+ bool removed = false;
+
+ for (map<int, vector<AtomicBoundary> >::iterator it = own.begin();
+ it != own.end(); ++it) {
+ for (vector<AtomicBoundary>::iterator bIt = it->second.begin();
+ bIt != it->second.end(); ++bIt) {
+ if (bIt->rankObj == bound) {
+ it->second.erase(bIt);
+ removed = true;
+ break;
+ }
+ }
+ }
+
+ return removed;
+ }
+
+
+ bool InteriorBoundary::removeOther(BoundaryObject& bound)
+ {
+ FUNCNAME("InteriorBoundary::removeOther()");
+
+ bool removed = false;
+
+ for (map<int, vector<AtomicBoundary> >::iterator it = other.begin();
+ it != other.end(); ++it) {
+ for (vector<AtomicBoundary>::iterator bIt = it->second.begin();
+ bIt != it->second.end(); ++bIt) {
+ if (bIt->rankObj == bound) {
+ it->second.erase(bIt);
+ removed = true;
+ break;
+ }
+ }
+ }
+
+ return removed;
+ }
+
+
void InteriorBoundary::serializeExcludeList(ostream &out,
ExcludeList &list)
{
diff --git a/AMDiS/src/parallel/InteriorBoundary.h b/AMDiS/src/parallel/InteriorBoundary.h
index 2660253f..594c295e 100644
--- a/AMDiS/src/parallel/InteriorBoundary.h
+++ b/AMDiS/src/parallel/InteriorBoundary.h
@@ -86,12 +86,32 @@ namespace AMDiS {
void deserialize(istream &in, Mesh *mesh);
private:
+ /// In this function, we put some verification procedures to check for
+ /// consistency if the created interior boundary data. This function is
+ /// enabled even for optimized mode to check all possible meshes. Thus,
+ /// everything herein should be fast. If the code is stable, we can think
+ /// to remove this function.
+ void verifyBoundary();
+
AtomicBoundary& getNewOwn(int rank);
AtomicBoundary& getNewOther(int rank);
AtomicBoundary& getNewPeriodic(int rank);
+ /// Checks whether the given boundary is already a other boundary with
+ /// given rank number. Returns true, if the bound is already in the other
+ /// boundary database.
+ bool checkOther(AtomicBoundary& bound, int rank);
+
+ /// Removes the given boundary object from all owned boundaries. Returns
+ /// true, if at least one object has been removed, otherwise false.
+ bool removeOwn(BoundaryObject& bound);
+
+ /// Removes the given boundary object from all other boundaries. Returns
+ /// true, if at least one object has been removed, otherwise false.
+ bool removeOther(BoundaryObject& bound);
+
void serialize(ostream &out, RankToBoundMap& boundary);
void deserialize(istream &in,
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 444cb894..c84d4e23 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -68,8 +68,6 @@ namespace AMDiS {
return (*dof1 < *dof2);
}
- bool MeshDistributor::sebastianMode = false;
-
MeshDistributor::MeshDistributor()
: problemStat(0),
initialized(false),
@@ -1696,7 +1694,7 @@ namespace AMDiS {
dofMaps[i]->update();
#if (DEBUG != 0)
- dofMaps[i]->printInfo();
+ // dofMaps[i]->printInfo();
#endif
}
diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h
index db4b23f5..e8a5a791 100644
--- a/AMDiS/src/parallel/MeshDistributor.h
+++ b/AMDiS/src/parallel/MeshDistributor.h
@@ -290,7 +290,7 @@ namespace AMDiS {
void setBoundaryDofRequirement(Flag flag)
{
- createBoundaryDofFlag = flag;
+ createBoundaryDofFlag |= flag;
}
BoundaryDofInfo& getBoundaryDofInfo(const FiniteElemSpace *feSpace,
@@ -584,8 +584,6 @@ namespace AMDiS {
vector<ParallelDofMapping*> dofMaps;
public:
- static bool sebastianMode;
-
/// The boundary DOFs are sorted by subobject entities, i.e., first all
/// face DOFs, edge DOFs and to the last vertex DOFs will be set to
/// communication structure vectors, \ref sendDofs and \ref recvDofs.
diff --git a/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h b/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h
index 3ceb8e73..b4911552 100644
--- a/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h
+++ b/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h
@@ -47,6 +47,7 @@ namespace AMDiS {
*/
class ParallelCoarseSpaceMatVec {
public:
+ /// Constructor
ParallelCoarseSpaceMatVec();
/// Set parallel DOF mapping for the interior DOFs.
@@ -54,6 +55,12 @@ namespace AMDiS {
{
interiorMap = interiorDofs;
}
+
+ /// Returns the parallel DOF mapping for the interior DOFs.
+ ParallelDofMapping* getDofMapping()
+ {
+ return interiorMap;
+ }
/** \brief
* Sets the coarse space for all or a specific component.
diff --git a/AMDiS/src/parallel/ParallelProblemStatBase.h b/AMDiS/src/parallel/ParallelProblemStatBase.h
index 7320196a..6eeb5c67 100644
--- a/AMDiS/src/parallel/ParallelProblemStatBase.h
+++ b/AMDiS/src/parallel/ParallelProblemStatBase.h
@@ -40,8 +40,8 @@ namespace AMDiS {
virtual ~ParallelProblemStatBase() {}
void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
- bool assembleMatrix,
- bool assembleVector);
+ bool assembleMatrix = true,
+ bool assembleVector = true);
void initialize(Flag initFlag,
ProblemStatSeq *adoptProblem = NULL,
diff --git a/AMDiS/src/parallel/PetscProblemStat.cc b/AMDiS/src/parallel/PetscProblemStat.cc
index b9c57b2e..4aa37d37 100644
--- a/AMDiS/src/parallel/PetscProblemStat.cc
+++ b/AMDiS/src/parallel/PetscProblemStat.cc
@@ -33,15 +33,16 @@ namespace AMDiS {
FUNCNAME("PetscProblemStat::PetscProblemStat()");
string tmp("");
- Parameters::get("parallel->solver", tmp);
-
+ string initFileStr = nameStr + "->solver";
+ Parameters::get(initFileStr, tmp);
+
if (tmp == "petsc-schur") {
petscSolver = new PetscSolverSchur();
} else if (tmp == "petsc-feti") {
- petscSolver = new PetscSolverFeti();
+ petscSolver = new PetscSolverFeti(initFileStr);
} else if (tmp == "petsc-block") {
petscSolver = new PetscSolverGlobalBlockMatrix();
- } else if (tmp == "petsc" || tmp == "") {
+ } else if (tmp == "petsc") {
petscSolver = new PetscSolverGlobalMatrix();
} else if (tmp == "bddcml") {
#ifdef HAVE_BDDC_ML
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 72402ac3..e91404cc 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -28,8 +28,9 @@ namespace AMDiS {
using namespace std;
- PetscSolverFeti::PetscSolverFeti()
+ PetscSolverFeti::PetscSolverFeti(string name)
: PetscSolver(),
+ initFileStr(name),
primalDofMap(COMPONENT_WISE),
dualDofMap(COMPONENT_WISE),
interfaceDofMap(COMPONENT_WISE),
@@ -51,8 +52,8 @@ namespace AMDiS {
FUNCNAME("PetscSolverFeti::PetscSolverFeti()");
string preconditionerName = "";
- Parameters::get("parallel->solver->precon", preconditionerName);
- if (preconditionerName == "" || preconditionerName == "none") {
+ Parameters::get(initFileStr + "->left precon", preconditionerName);
+ if (preconditionerName == "" || preconditionerName == "no") {
MSG("Create FETI-DP solver with no preconditioner!\n");
fetiPreconditioner = FETI_NONE;
} else if (preconditionerName == "dirichlet") {
@@ -66,20 +67,34 @@ namespace AMDiS {
preconditionerName.c_str());
}
- Parameters::get("parallel->feti->schur primal solver", schurPrimalSolver);
+ preconditionerName = "";
+ Parameters::get(initFileStr + "->right precon", preconditionerName);
+ if (preconditionerName != "" && preconditionerName != "no") {
+ ERROR_EXIT("FETI-DP does not support right preconditioning! (parameter \"%s->right precon\" has value \"%s\")\n",
+ initFileStr.c_str(), preconditionerName.c_str());
+ }
+
+ Parameters::get(initFileStr + "->feti->schur primal solver", schurPrimalSolver);
TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1)
("Wrong solver \"%d\"for the Schur primal complement!\n",
schurPrimalSolver);
+ Parameters::get(initFileStr + "->feti->stokes mode", stokesMode);
+ if (stokesMode) {
+ Parameters::get(initFileStr + "->feti->pressure component", pressureComponent);
+ TEST_EXIT(pressureComponent >= 0)
+ ("FETI-DP in Stokes mode, no pressure component defined!\n");
+ }
+
+ Parameters::get(initFileStr + "->feti->augmented lagrange", augmentedLagrange);
+
+ Parameters::get(initFileStr + "->feti->symmetric", isSymmetric);
+
Parameters::get("parallel->multi level test", multiLevelTest);
if (multiLevelTest)
meshLevel = 1;
Parameters::get("parallel->print timings", printTimings);
-
- Parameters::get("parallel->feti->augmented lagrange", augmentedLagrange);
-
- Parameters::get("parallel->feti->symmetric", isSymmetric);
}
@@ -92,14 +107,6 @@ namespace AMDiS {
MeshLevelData& levelData = meshDistributor->getMeshLevelData();
- stokesMode = false;
- Parameters::get("parallel->feti->stokes mode", stokesMode);
- if (stokesMode) {
- Parameters::get("parallel->feti->pressure component", pressureComponent);
- TEST_EXIT(pressureComponent >= 0)
- ("FETI-DP in Stokes mode, no pressure component defined!\n");
- }
-
if (subdomain == NULL) {
subdomain = new PetscSolverGlobalMatrix();
subdomain->setSymmetric(isSymmetric);
@@ -1149,8 +1156,19 @@ namespace AMDiS {
KSPSetType(ksp_feti, KSPGMRES);
KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
KSPSetFromOptions(ksp_feti);
- if (stokesMode)
- KSPMonitorSet(ksp_feti, KSPMonitorFetiStokes, &fetiKspData, PETSC_NULL);
+
+
+ // === Set KSP monitor. ===
+
+ bool monitor = false;
+ Parameters::get(initFileStr + "->feti->monitor", monitor);
+ if (monitor) {
+ if (stokesMode)
+ KSPMonitorSet(ksp_feti, KSPMonitorFetiStokes, &fetiKspData, PETSC_NULL);
+ else
+ KSPMonitorSet(ksp_feti, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
+ }
+
// === Create null space objects. ===
@@ -1257,64 +1275,53 @@ namespace AMDiS {
}
if (stokesMode) {
- // === Create H2 vec ===
+ // === Create mass matrix solver ===
const FiniteElemSpace *pressureFeSpace =
componentSpaces[pressureComponent];
- DOFVector<double> tmpvec(pressureFeSpace, "tmpvec");
- createH2vec(tmpvec);
- interfaceDofMap.createVec(fetiInterfaceLumpedPreconData.h2vec);
-
- DofMap& m = interfaceDofMap[pressureComponent].getMap();
- for (DofMap::iterator it = m.begin(); it != m.end(); ++it) {
- if (dofMap[pressureFeSpace].isRankDof(it->first)) {
- int index = interfaceDofMap.getMatIndex(pressureComponent, it->first);
- VecSetValue(fetiInterfaceLumpedPreconData.h2vec,
- index, tmpvec[it->first], INSERT_VALUES);
- }
- }
-
- VecAssemblyBegin(fetiInterfaceLumpedPreconData.h2vec);
- VecAssemblyEnd(fetiInterfaceLumpedPreconData.h2vec);
-
- // === Create mass matrix solver ===
-
- if (!massMatrixSolver) {
- ParallelDofMapping *massMapping = new ParallelDofMapping(COMPONENT_WISE);
+ // Create parallel DOF mapping in pressure space.
+ ParallelDofMapping *massMapping = NULL;
+ if (massMatrixSolver) {
+ massMapping = massMatrixSolver->getDofMapping();
+ } else {
+ massMapping =
+ new ParallelDofMapping(COMPONENT_WISE);
massMapping->init(meshDistributor->getMeshLevelData(),
pressureFeSpace, pressureFeSpace);
massMapping->setDofComm(meshDistributor->getDofComm());
massMapping->setMpiComm(meshDistributor->getMeshLevelData().getMpiComm(0), 0);
massMapping->setComputeMatIndex(true);
- (*massMapping)[0] = interfaceDofMap[pressureComponent];
- massMapping->update();
-
- DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
- Operator op(pressureFeSpace, pressureFeSpace);
- Simple_ZOT zot;
- op.addTerm(&zot);
- massMatrix.assembleOperator(op);
+ }
+ (*massMapping)[0] = interfaceDofMap[pressureComponent];
+ massMapping->update();
+
+ DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
+ Operator op(pressureFeSpace, pressureFeSpace);
+ Simple_ZOT zot;
+ op.addTerm(&zot);
+ massMatrix.assembleOperator(op);
+ if (!massMatrixSolver) {
massMatrixSolver = new PetscSolverGlobalMatrix;
massMatrixSolver->setKspPrefix("mass_");
massMatrixSolver->setMeshDistributor(meshDistributor,
mpiCommGlobal,
mpiCommLocal);
massMatrixSolver->setDofMapping(massMapping);
- massMatrixSolver->fillPetscMatrix(&massMatrix);
-
- int r, c;
- MatGetSize(massMatrixSolver->getMatInterior(), &r, &c);
- MatInfo info;
- MatGetInfo(massMatrixSolver->getMatInterior(), MAT_GLOBAL_SUM, &info);
- MSG("MASS MAT INFO: size = %d x %d nnz = %d\n",
- r, c,
- static_cast<int>(info.nz_used));
-
- fetiInterfaceLumpedPreconData.ksp_mass =
- massMatrixSolver->getSolver();
}
+
+ massMatrixSolver->fillPetscMatrix(&massMatrix);
+
+ int r, c;
+ MatGetSize(massMatrixSolver->getMatInterior(), &r, &c);
+ MatInfo info;
+ MatGetInfo(massMatrixSolver->getMatInterior(), MAT_GLOBAL_SUM, &info);
+ MSG("MASS MAT INFO: size = %d x %d nnz = %d\n",
+ r, c, static_cast<int>(info.nz_used));
+
+ fetiInterfaceLumpedPreconData.ksp_mass = massMatrixSolver->getSolver();
+
// === Create tmp vectors ===
@@ -2374,59 +2381,4 @@ namespace AMDiS {
MeshDistributor::globalMeshDistributor->synchVector(vec);
}
-
- void PetscSolverFeti::createH2vec(DOFVector<double> &vec)
- {
- FUNCNAME("PetscSolverFeti::createH2vec()");
-
- vec.set(0.0);
-
- Mesh *mesh = vec.getFeSpace()->getMesh();
- TEST_EXIT(mesh->getDim() == 2)("This works only in 2D!\n");
- int n0 = vec.getFeSpace()->getAdmin()->getNumberOfPreDofs(VERTEX);
-
- TraverseStack stack;
- ElInfo *elInfo =
- stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
- while (elInfo) {
- Element *el = elInfo->getElement();
- std::vector<double> length(3, 0.0);
- double sumLength = 0.0;
-
- for (int i = 0; i < 3; i++) {
- int idx0 = el->getVertexOfEdge(i, 0);
- int idx1 = el->getVertexOfEdge(i, 1);
- DegreeOfFreedom dof0 = el->getDof(idx0, n0);
- DegreeOfFreedom dof1 = el->getDof(idx1, n0);
-
- WorldVector<double> c0 = elInfo->getCoord(idx0);
- WorldVector<double> c1 = elInfo->getCoord(idx1);
- c0 -= c1;
- length[i] = norm(c0);
- sumLength += length[i];
- }
-
- sumLength *= 0.5;
- double area =
- sqrt(sumLength * (sumLength - length[0]) * (sumLength - length[1]) * (sumLength - length[2]));
-
- for (int i = 0; i < 3; i++) {
- DegreeOfFreedom d = el->getDof(i, n0);
- vec[d] += area;
- }
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- meshDistributor->synchAddVector(vec);
-
- double scalePower = 2.0;
- Parameters::get("lp->scale power", scalePower);
- DOFIterator<double> dofIt(&vec, USED_DOFS);
- for (dofIt.reset(); !dofIt.end(); ++dofIt) {
- DegreeOfFreedom d = dofIt.getDOFIndex();
- vec[d] = 1.0 / pow(vec[d], scalePower);
- }
- }
-
}
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index e10d725b..e4398497 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -41,7 +41,8 @@ namespace AMDiS {
class PetscSolverFeti : public PetscSolver
{
public:
- PetscSolverFeti();
+ /// Constructor of FETI-DP solver class.
+ PetscSolverFeti(string name);
/// Assemble the sequentially created matrices to the global matrices
/// required by the FETI-DP method.
@@ -228,11 +229,10 @@ namespace AMDiS {
return false;
}
- /// Traverse the mesh and create a DOF vector where each DOF is in the
- /// order of: 1/h^2
- void createH2vec(DOFVector<double> &vec);
-
protected:
+ /// Prefix string for parameters in init file.
+ string initFileStr;
+
/// Mapping from primal DOF indices to a global index of primals.
ParallelDofMapping primalDofMap;
diff --git a/AMDiS/src/parallel/PetscSolverFetiMonitor.cc b/AMDiS/src/parallel/PetscSolverFetiMonitor.cc
index 866d5efb..07ee4d7f 100644
--- a/AMDiS/src/parallel/PetscSolverFetiMonitor.cc
+++ b/AMDiS/src/parallel/PetscSolverFetiMonitor.cc
@@ -35,7 +35,7 @@ namespace AMDiS {
VecDestroy(&v);
VecDestroy(&w);
- PetscPrintf(PETSC_COMM_WORLD, "%3D KSP residual norm %1.12e [%1.12e %1.12e] and preconditioned norm [%1.12e]\n",
+ PetscPrintf(PETSC_COMM_WORLD, "%3D KSP residual norm %1.12e [ %1.12e %1.12e ] and preconditioned norm [%1.12e]\n",
n, norm, norm0, norm1, rnorm);
return 0;
diff --git a/AMDiS/src/parallel/PetscSolverFetiOperators.cc b/AMDiS/src/parallel/PetscSolverFetiOperators.cc
index cd9c94c3..2b66bd67 100644
--- a/AMDiS/src/parallel/PetscSolverFetiOperators.cc
+++ b/AMDiS/src/parallel/PetscSolverFetiOperators.cc
@@ -466,21 +466,7 @@ namespace AMDiS {
VecNestGetSubVec(yvec, 0, &y_interface);
VecNestGetSubVec(yvec, 1, &y_lagrange);
- bool useMassMatrix = false;
- Parameters::get("lp->mass matrix", useMassMatrix);
- if (useMassMatrix) {
- KSPSolve(data->ksp_mass, x_interface, y_interface);
- } else {
- VecCopy(x_interface, y_interface);
- double scaleFactor = 1.0;
- Parameters::get("lp->scal", scaleFactor);
- bool mult = false;
- Parameters::get("lp->mult", mult);
- if (mult)
- VecPointwiseMult(y_interface, x_interface, data->h2vec);
- if (scaleFactor != 0.0)
- VecScale(y_interface, scaleFactor);
- }
+ KSPSolve(data->ksp_mass, x_interface, y_interface);
MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b0);
diff --git a/AMDiS/src/parallel/PetscSolverFetiStructs.h b/AMDiS/src/parallel/PetscSolverFetiStructs.h
index 8fa97340..2ed5cedb 100644
--- a/AMDiS/src/parallel/PetscSolverFetiStructs.h
+++ b/AMDiS/src/parallel/PetscSolverFetiStructs.h
@@ -144,8 +144,6 @@ namespace AMDiS {
PetscSolver* subSolver;
- Vec h2vec;
-
Vec tmp_primal;
KSP ksp_mass;
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index b3c0208b..a0ebc076 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -487,8 +487,8 @@ namespace AMDiS {
if (bIt->second && bIt->second->isDirichlet()) {
dirichletRow = true;
if ((dynamic_cast<DirichletBC*>(bIt->second))->applyBoundaryCondition()) {
- TEST_EXIT_DBG(dirichletMainCol == -1)("Should not happen!\n");
dirichletMainCol = j;
+ break;
}
}
}
--
GitLab