Skip to content

GitLab

Commit aded0540 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Revision of parallel domain decomposition code.

parent b4439d35
Hide whitespace changes
Inline Side-by-side
Showing with 338 additions and 247 deletions
AMDiS/libtool
View file @ aded0540
......@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
# Libtool was configured on host deimos103:
# Libtool was configured on host p2q071:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -82,13 +82,13 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
LTCC="gcc"
LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
# A language-specific compiler.
CC="gcc"
CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# Is the compiler the GNU C compiler?
with_gcc=yes
......@@ -171,7 +171,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
link_static_flag="-static"
link_static_flag=""
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
......@@ -6760,7 +6760,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
# Libtool was configured on host deimos103:
# Libtool was configured on host p2q071:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -6798,13 +6798,13 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
LTCC="gcc"
LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
# A language-specific compiler.
CC="g++"
CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpiCC"
# Is the compiler the GNU C compiler?
with_gcc=yes
......@@ -6887,7 +6887,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
link_static_flag="-static"
link_static_flag=""
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
......@@ -6954,11 +6954,11 @@ predeps=""
# Dependencies to place after the objects being linked to create a
# shared library.
postdeps="-lstdc++ -lm -lgcc_s -lc -lgcc_s"
postdeps="-lmpi_cxx -lmpi -lopen-rte -lopen-pal -libverbs -lrt -lnuma -ldl -lnsl -lutil -ldl -lstdc++ -lm -lgcc_s -lpthread -lc -lgcc_s"
# The library search path used internally by the compiler when linking
# a shared library.
compiler_lib_search_path="-L/usr/lib64/gcc/x86_64-suse-linux/4.1.2 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64 -L/lib/../lib64 -L/usr/lib/../lib64 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../.."
compiler_lib_search_path="-L/usr/lib64 -L/licsoft/libraries/openmpi/1.2.6/64bit/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64 -L/lib/../lib64 -L/usr/lib/../lib64 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../.."
# Method to check whether dependent libraries are shared objects.
deplibs_check_method="pass_all"
......@@ -7065,7 +7065,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
# Libtool was configured on host deimos103:
# Libtool was configured on host p2q071:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -7103,7 +7103,7 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
LTCC="gcc"
LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
......
AMDiS/src/Element.cc
View file @ aded0540
......@@ -566,71 +566,4 @@ namespace AMDiS {
return result;
}
bool fitElementToMeshCode(RefinementManager *refineManager,
MeshStructure &code,
Element *el,
int ithSide,
int elType)
{
FUNCNAME("fitElementToMeshCode()");
if (code.empty())
return false;
int s1 = el->getSideOfChild(0, ithSide, elType);
int s2 = el->getSideOfChild(1, ithSide, elType);
TEST_EXIT_DBG(s1 != -1 || s2 != -1)("This should not happen!\n");
if (s1 != -1 && s2 != -1)
return fitElementToMeshCode2(refineManager, code, el, ithSide, elType);
if (el->isLeaf()) {
if (code.getNumElements() == 1 && code.isLeafElement())
return false;
el->setMark(1);
refineManager->refineElement(el->getMesh(), el);
}
if (s1 != -1)
return fitElementToMeshCode2(refineManager, code,
el->getFirstChild(), s1, el->getChildType(elType));
else
return fitElementToMeshCode2(refineManager, code,
el->getSecondChild(), s2, el->getChildType(elType));
}
bool fitElementToMeshCode2(RefinementManager *refineManager, MeshStructure &code,
Element *el, int ithSide, int elType)
{
if (code.isLeafElement())
return false;
bool value = false;
if (el->isLeaf()) {
el->setMark(1);
refineManager->refineElement(el->getMesh(), el);
value = true;
}
int s1 = el->getSideOfChild(0, ithSide, elType);
int s2 = el->getSideOfChild(1, ithSide, elType);
if (s1 != -1) {
code.nextElement();
value |= fitElementToMeshCode2(refineManager, code, el->getFirstChild(),
s1, el->getChildType(elType));
}
if (s2 != -1) {
code.nextElement();
value |= fitElementToMeshCode2(refineManager, code, el->getSecondChild(),
s2, el->getChildType(elType));
}
return value;
}
}
AMDiS/src/Element.h
View file @ aded0540
......@@ -574,18 +574,6 @@ namespace AMDiS {
friend class Mesh;
};
bool fitElementToMeshCode(RefinementManager *refineManager,
MeshStructure &code,
Element *el,
int ithSide,
int elType);
bool fitElementToMeshCode2(RefinementManager *refineManager,
MeshStructure &code,
Element *el,
int ithSide,
int elType);
}
#endif // AMDIS_ELEMENT_H
......
AMDiS/src/InteriorBoundary.cc
View file @ aded0540
......@@ -76,12 +76,14 @@ namespace AMDiS {
SerUtil::serialize(out, bound.rankObj.subObj);
SerUtil::serialize(out, bound.rankObj.ithObj);
SerUtil::serialize(out, bound.rankObj.reverseMode);
serializeExcludeList(out, bound.rankObj.excludedSubstructures);
SerUtil::serialize(out, bound.neighObj.elIndex);
SerUtil::serialize(out, bound.neighObj.elType);
SerUtil::serialize(out, bound.neighObj.subObj);
SerUtil::serialize(out, bound.neighObj.ithObj);
SerUtil::serialize(out, bound.neighObj.reverseMode);
serializeExcludeList(out, bound.neighObj.excludedSubstructures);
}
}
}
......@@ -107,16 +109,48 @@ namespace AMDiS {
SerUtil::deserialize(in, bound.rankObj.subObj);
SerUtil::deserialize(in, bound.rankObj.ithObj);
SerUtil::deserialize(in, bound.rankObj.reverseMode);
deserializeExcludeList(in, bound.rankObj.excludedSubstructures);
SerUtil::deserialize(in, bound.neighObj.elIndex);
SerUtil::deserialize(in, bound.neighObj.elType);
SerUtil::deserialize(in, bound.neighObj.subObj);
SerUtil::deserialize(in, bound.neighObj.ithObj);
SerUtil::deserialize(in, bound.neighObj.reverseMode);
deserializeExcludeList(in, bound.neighObj.excludedSubstructures);
bound.rankObj.el = elIndexMap[bound.rankObj.elIndex];
bound.neighObj.el = NULL;
}
}
}
void InteriorBoundary::serializeExcludeList(std::ostream &out, ExcludeList &list)
{
int size = list.size();
SerUtil::serialize(out, size);
for (int i = 0; i < size; i++) {
SerUtil::serialize(out, list[i].first);
SerUtil::serialize(out, list[i].second);
}
}
void InteriorBoundary::deserializeExcludeList(std::istream &in, ExcludeList &list)
{
int size = 0;
SerUtil::deserialize(in, size);
list.resize(0);
list.reserve(size);
for (int i = 0; i < size; i++) {
GeoIndex a;
int b;
SerUtil::deserialize(in, a);
SerUtil::deserialize(in, b);
list.push_back(std::make_pair(a, b));
}
}
}
AMDiS/src/InteriorBoundary.h
View file @ aded0540
......@@ -31,12 +31,14 @@
namespace AMDiS {
typedef std::vector<std::pair<GeoIndex, int> > ExcludeList;
/// Defines the geometrical objects that forms the boundary;
struct BoundaryObject {
BoundaryObject()
: reverseMode(false),
notIncludedSubStructures(0)
excludedSubstructures(0)
{}
BoundaryObject(Element *e, int eType, GeoIndex sObj, int iObj)
......@@ -46,7 +48,7 @@ namespace AMDiS {
subObj(sObj),
ithObj(iObj),
reverseMode(false),
notIncludedSubStructures(0)
excludedSubstructures(0)
{}
void setReverseMode(BoundaryObject &otherBound, FiniteElemSpace *feSpace);
......@@ -83,7 +85,15 @@ namespace AMDiS {
bool reverseMode;
std::vector<std::pair<GeoIndex, int> > notIncludedSubStructures;
/** \brief
* In many situations it may be necessary to exclude some parts of the element to
* be part of the boundary. In 3d, when a face is part of the boundary, an edge or
* an vertex may be exludeded. In 2d only vertices may be exluded to be part of
* an edge boundary. This list contains pairs of exludeded structures. The first
* component of every pair denotes if it is a vertex or an edge, and the second
* component denotes the local index of the structure.
*/
ExcludeList excludedSubstructures;
};
/** \brief
......@@ -117,6 +127,11 @@ namespace AMDiS {
/// Reads the state of an interior boundary from a file.
void deserialize(std::istream &in, std::map<int, Element*> &elIndexMap);
protected:
void serializeExcludeList(std::ostream &out, ExcludeList &list);
void deserializeExcludeList(std::istream &in, ExcludeList &list);
public:
RankToBoundMap boundary;
};
......
AMDiS/src/ParMetisPartitioner.cc
View file @ aded0540
......@@ -287,14 +287,12 @@ namespace AMDiS {
switch(mode) {
case INITIAL:
// TODO: Removed weight because it does not work correctly for very large
// macro meshes.
ParMETIS_V3_PartKway(parMetisMesh->getElementDist(),
parMetisGraph.getXAdj(),
parMetisGraph.getAdjncy(),
wgts,
NULL,
NULL,
NULL,
&wgtflag,
&numflag,
&ncon,
&nparts,
......
AMDiS/src/ParallelDomainBase.cc
View file @ aded0540
#include <algorithm>
#include <iostream>
#include <fstream>
#include "ParallelDomainBase.h"
#include "ParMetisPartitioner.h"
......@@ -49,6 +51,8 @@ namespace AMDiS {
mesh(fe->getMesh()),
refineManager(refinementManager),
initialPartitionMesh(true),
d_nnz(NULL),
o_nnz(NULL),
nRankDofs(0),
rstart(0),
nComponents(1),
......@@ -364,38 +368,25 @@ namespace AMDiS {
}
void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
void ParallelDomainBase::createPetscNnzStructure(Matrix<DOFMatrix*> *mat)
{
FUNCNAME("ParallelDomainBase::fillPetscMatrix()");
clock_t first = clock();
// === Create PETSc vector (rhs, solution and a temporary vector). ===
FUNCNAME("ParallelDomainBase::createPetscNnzStructure()");
VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
VecSetType(petscRhsVec, VECMPI);
VecCreate(PETSC_COMM_WORLD, &petscSolVec);
VecSetSizes(petscSolVec, nRankRows, nOverallRows);
VecSetType(petscSolVec, VECMPI);
TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n");
TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n");
VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
VecSetType(petscTmpVec, VECMPI);
d_nnz = new int[nRankRows];
o_nnz = new int[nRankRows];
for (int i = 0; i < nRankRows; i++) {
d_nnz[i] = 0;
o_nnz[i] = 0;
}
using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
namespace traits = mtl::traits;
typedef DOFMatrix::base_matrix_type Matrix;
typedef std::vector<std::pair<int, int> > MatrixNnzEntry;
int d_nnz[nRankRows];
int o_nnz[nRankRows];
for (int i = 0; i < nRankRows; i++) {
d_nnz[i] = 0;
o_nnz[i] = 0;
}
// Stores to each rank a list of nnz entries (i.e. pairs of row and column index)
// that this rank will send to. This nnz entries will be assembled on this rank,
// but because the row DOFs are not DOFs of this rank they will be send to the
......@@ -531,12 +522,39 @@ namespace AMDiS {
}
}
}
}
void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
{
FUNCNAME("ParallelDomainBase::fillPetscMatrix()");
clock_t first = clock();
// === Create PETSc vector (rhs, solution and a temporary vector). ===
VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
VecSetType(petscRhsVec, VECMPI);
VecCreate(PETSC_COMM_WORLD, &petscSolVec);
VecSetSizes(petscSolVec, nRankRows, nOverallRows);
VecSetType(petscSolVec, VECMPI);
VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
VecSetType(petscTmpVec, VECMPI);
if (!d_nnz)
createPetscNnzStructure(mat);
// === Create PETSc matrix with the computed nnz data structure. ===
MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
0, d_nnz, 0, o_nnz, &petscMatrix);
INFO(info, 8)("Fill petsc matrix 1 needed %.5f seconds\n", TIME_USED(first, clock()));
#if (DEBUG != 0)
int a, b;
MatGetOwnershipRange(petscMatrix, &a, &b);
......@@ -551,6 +569,8 @@ namespace AMDiS {
if ((*mat)[i][j])
setDofMatrix((*mat)[i][j], nComponents, i, j);
INFO(info, 8)("Fill petsc matrix 2 needed %.5f seconds\n", TIME_USED(first, clock()));
MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);
......@@ -727,10 +747,20 @@ namespace AMDiS {
do {
// To check the interior boundaries, the ownership of the boundaries is not
// important. Therefore, we add all boundaries to one boundary container.
RankToBoundMap allBound = myIntBoundary.boundary;
RankToBoundMap allBound;
for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin();
it != myIntBoundary.boundary.end(); ++it)
for (std::vector<AtomicBoundary>::iterator bit = it->second.begin();
bit != it->second.end(); ++bit)
if (bit->rankObj.subObj == FACE)
allBound[it->first].push_back(*bit);
for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
it != otherIntBoundary.boundary.end(); ++it)
allBound[it->first] = it->second;
for (std::vector<AtomicBoundary>::iterator bit = it->second.begin();
bit != it->second.end(); ++bit)
if (bit->rankObj.subObj == FACE)
allBound[it->first].push_back(*bit);
// === Check the boundaries and adapt mesh if necessary. ===
......@@ -741,6 +771,7 @@ namespace AMDiS {
int sendValue = static_cast<int>(!meshChanged);
recvAllValues = 0;
mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
MSG("LOOP\n");
} while (recvAllValues != 0);
......@@ -754,6 +785,19 @@ namespace AMDiS {
// === Because the mesh has been changed, update the DOF numbering and mappings. ===
updateLocalGlobalNumbering();
// === If there is a non zero pattern computed for Petsc matrix, delete it. So ===
// === it will be recomputed after next assembling. ===
if (d_nnz) {
delete [] d_nnz;
d_nnz = NULL;
}
if (o_nnz) {
delete [] o_nnz;
o_nnz = NULL;
}
}
......@@ -761,6 +805,8 @@ namespace AMDiS {
{
FUNCNAME("ParallelDomainBase::checkAndAdaptBoundary()");
MSG("CHECK_AND_ADAPT 1!\n");
// === Create mesh structure codes for all ranks boundary elements. ===
std::map<int, MeshCodeVec> sendCodes;
......@@ -774,14 +820,21 @@ namespace AMDiS {
}
}
MSG("CHECK_AND_ADAPT 2!\n");
StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
stdMpi.send(sendCodes);
stdMpi.recv(allBound);
stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
MSG("CHECK_AND_ADAPT 3!\n");
clock_t first = clock();
// === Compare received mesh structure codes. ===
bool meshFitTogether = true;
int superCounter = 0;
for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
......@@ -790,18 +843,27 @@ namespace AMDiS {
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end(); ++boundIt) {
MeshStructure elCode;
elCode.init(boundIt->rankObj);
if (elCode.getCode() != recvCodes[i].getCode()) {
TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
bool b = fitElementToMeshCode(refineManager,
recvCodes[i],
superCounter++;
clock_t fff = clock();
int refCounter = 0;
bool b = fitElementToMeshCode(recvCodes[i],
boundIt->rankObj.el,
boundIt->rankObj.ithObj,
boundIt->rankObj.elType);
boundIt->rankObj.elType, refCounter);
MSG("SIZE OF ELINFO3d = %d\n", sizeof(ElInfo3d));
MSG("Code length %d with %d refs needed %.5f seconds\n",
recvCodes[i].getNumElements(),
refCounter,
TIME_USED(fff, clock()));
if (b)
meshFitTogether = false;
}
......@@ -810,9 +872,84 @@ namespace AMDiS {
}
}
MSG("time for %d needed %.5f seconds\n", superCounter, TIME_USED(first, clock()));
MSG("CHECK_AND_ADAPT 4!\n");
return meshFitTogether;
}
bool ParallelDomainBase::fitElementToMeshCode(MeshStructure &code,
Element *el,
int ithSide,
int elType,
int &refCounter)
{
FUNCNAME("fitElementToMeshCode()");
if (code.empty())
return false;
int s1 = el->getSideOfChild(0, ithSide, elType);
int s2 = el->getSideOfChild(1, ithSide, elType);
TEST_EXIT_DBG(s1 != -1 || s2 != -1)("This should not happen!\n");
if (s1 != -1 && s2 != -1)
return fitElementToMeshCode2(code, el, ithSide, elType, refCounter);
if (el->isLeaf()) {
if (code.getNumElements() == 1 && code.isLeafElement())
return false;
el->setMark(1);
refineManager->refineElement(el->getMesh(), el);
refCounter++;
}
if (s1 != -1)
return fitElementToMeshCode2(code,
el->getFirstChild(), s1, el->getChildType(elType), refCounter);
else
return fitElementToMeshCode2(code,
el->getSecondChild(), s2, el->getChildType(elType), refCounter);
}
bool ParallelDomainBase::fitElementToMeshCode2(MeshStructure &code,
Element *el,
int ithSide,
int elType, int &refCounter)
{
if (code.isLeafElement())
return false;
bool value = false;
if (el->isLeaf()) {
el->setMark(1);
refineManager->refineElement(el->getMesh(), el);
value = true;
refCounter++;
}
int s1 = el->getSideOfChild(0, ithSide, elType);
int s2 = el->getSideOfChild(1, ithSide, elType);
if (s1 != -1) {
code.nextElement();
value |= fitElementToMeshCode2(code, el->getFirstChild(),