Commit e8cfacab authored by Thomas Witkowski's avatar Thomas Witkowski

First work in repartitioning in parallel computations. A lot of clean up in code structure.

parent f6559482
......@@ -75,6 +75,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/DOFIndexed.cc
${SOURCE_DIR}/CreatorMap.cc
${SOURCE_DIR}/ProblemInterpolScal.cc
${SOURCE_DIR}/ProblemInterpolVec.cc
${SOURCE_DIR}/MacroInfo.cc
${SOURCE_DIR}/MacroReader.cc
${SOURCE_DIR}/ValueReader.cc
${SOURCE_DIR}/Projection.cc
......@@ -165,10 +166,6 @@ if(ENABLE_PARMETIS)
CMAKE_FORCE_CXX_COMPILER(mpiCC "MPI C++ compiler")
include_directories( ${LIB_DIR}/ParMetis-3.1)
SET(PARALLEL_AMDIS_SRC ${SOURCE_DIR}/ConditionalEstimator.cc ${SOURCE_DIR}/ConditionalMarker.cc
${SOURCE_DIR}/ParallelProblem.cc ${SOURCE_DIR}/ParMetisPartitioner.cc
${SOURCE_DIR}/PollutionError.cc)
SET(COMPILEFLAGS "${COMPILEFLAGS} -DHAVE_PARALLEL_AMDIS=1")
INSTALL(FILES ${LIB_DIR}/ParMetis-3.1/parmetis.h
${LIB_DIR}/ParMetis-3.1/libparmetis.a
${LIB_DIR}/ParMetis-3.1/libmetis.a
......@@ -179,6 +176,7 @@ endif(ENABLE_PARMETIS)
if(ENABLE_PARALLEL_DOMAIN)
include_directories(${PETSC_DIR}/include ${PETSC_DIR}/${PETSC_ARCH}/include)
SET(PARALLEL_DOMAIN_AMDIS_SRC
${SOURCE_DIR}/parallel/ParMetisPartitioner.cc
${SOURCE_DIR}/parallel/MeshDistributor.cc
${SOURCE_DIR}/parallel/StdMpi.cc
${SOURCE_DIR}/parallel/ParallelDebug.cc
......@@ -247,7 +245,7 @@ SET(COMPOSITE_FEM_SRC ${COMPOSITE_SOURCE_DIR}/CFE_Integration.cc ${COMPOSITE_SOU
include_directories(${MTL_DIR})
include_directories(${SOURCE_DIR})
add_library(amdis SHARED ${AMDIS_SRC} ${PARALLEL_AMDIS_SRC} ${PARALLEL_DOMAIN_AMDIS_SRC})
add_library(amdis SHARED ${AMDIS_SRC} ${PARALLEL_DOMAIN_AMDIS_SRC})
add_library(compositeFEM SHARED ${COMPOSITE_FEM_SRC})
target_link_libraries(compositeFEM amdis)
......
......@@ -10,25 +10,13 @@ AMDIS_INCLUDES = -I$(SOURCE_DIR)
libamdis_la_CXXFLAGS =
if USE_PARALLEL_AMDIS
PARALLEL_AMDIS_SOURCES = \
$(PARALLEL_DIR)/ConditionalEstimator.h $(PARALLEL_DIR)/ConditionalEstimator.cc \
$(PARALLEL_DIR)/ConditionalMarker.h \ $(PARALLEL_DIR)/ConditionalMarker.cc \
$(PARALLEL_DIR)/ParallelError.h $(PARALLEL_DIR)/ParallelError.hh \
$(PARALLEL_DIR)/ParallelProblem.h $(PARALLEL_DIR)/ParallelProblem.cc \
$(PARALLEL_DIR)/ParMetisPartitioner.h $(PARALLEL_DIR)/ParMetisPartitioner.cc \
$(PARALLEL_DIR)/PartitionElementData.h \
$(PARALLEL_DIR)/PollutionError.h $(PARALLEL_DIR)/PollutionError.cc
PARALLEL_INCLUDES = -I$(MPI_DIR)/include -I$(PARMETIS_DIR)
libamdis_la_CXXFLAGS += -DHAVE_PARALLEL_AMDIS=1
else
PARALLEL_AMDIS_SOURCES =
PARALLEL_INCLUDES =
endif
PARALLEL_AMDIS_SOURCES =
PARALLEL_INCLUDES =
if USE_PARALLEL_DOMAIN_AMDIS
PARALLEL_AMDIS_SOURCES += \
$(SOURCE_DIR)/parallel/StdMpi.h $(SOURCE_DIR)/parallel/StdMpi.cc \
$(SOURCE_DIR)/parallel/ParMetisPartitioner.h $(SOURCE_DIR)/parallel/ParMetisPartitioner.cc \
$(SOURCE_DIR)/parallel/MeshDistributor.h $(SOURCE_DIR)/parallel/MeshDistributor.cc \
$(SOURCE_DIR)/parallel/ParallelDebug.h $(SOURCE_DIR)/parallel/ParallelDebug.cc \
$(SOURCE_DIR)/parallel/ParallelProblemStatBase.h \
......@@ -98,6 +86,7 @@ $(SOURCE_DIR)/ProblemInterpolVec.h $(SOURCE_DIR)/ProblemInterpolVec.cc \
$(SOURCE_DIR)/Serializable.h \
$(SOURCE_DIR)/BallProject.h \
$(SOURCE_DIR)/CylinderProject.h \
$(SOURCE_DIR)/MacroInfo.h $(SOURCE_DIR)/MacroInfo.cc \
$(SOURCE_DIR)/MacroReader.h $(SOURCE_DIR)/MacroReader.cc \
$(SOURCE_DIR)/ValueReader.h $(SOURCE_DIR)/ValueReader.cc \
$(SOURCE_DIR)/Projection.h $(SOURCE_DIR)/Projection.cc \
......
This diff is collapsed.
......@@ -75,7 +75,7 @@ namespace AMDiS {
}
if (mesh->getNumberOfDOFs(CENTER) && !mesh->queryCoarseDOFs()) {
parent->setDOF(mesh->getNode(CENTER), mesh->getDOF(CENTER));
parent->setDOF(mesh->getNode(CENTER), mesh->getDof(CENTER));
}
/*--------------------------------------------------------------------------*/
......
......@@ -84,7 +84,7 @@ namespace AMDiS {
// get new dof on el at the midpoint of the coarsening edge
if (!el->getDOF(node + 2)) {
el->setDOF(node + 2, mesh->getDOF(EDGE));
el->setDOF(node + 2, mesh->getDof(EDGE));
if (neigh)
neigh->setDOF(node + 2, const_cast<int*>(el->getDOF(node + 2)));
}
......
......@@ -324,7 +324,7 @@ namespace AMDiS {
/****************************************************************************/
node = mesh->getNode(EDGE);
if (!(dof = const_cast<int*>( el->getDOF(node))))
dof = mesh->getDOF(EDGE);
dof = mesh->getDof(EDGE);
} else {
dof = NULL;
}
......
#include "ConditionalEstimator.h"
#include "ElInfo.h"
#include "PartitionElementData.h"
#include "Element.h"
#include "Traverse.h"
#include "mpi.h"
namespace AMDiS {
ConditionalEstimator::ConditionalEstimator(Estimator* decorated)
: decoratedEstimator_(decorated),
elementCount_(0),
row_(decorated ? decorated->getRow() : 0)
{
FUNCNAME("ConditionalEstimator::ConditionalEstimator(Estimator* decorated)");
if (decorated!=NULL) {
name = decorated->getName();
int swap;
int retVal =
Parameters::getGlobalParameter(0, name + "->estimate outer", "%d", &swap);
if (retVal > 0)
estimateOut = (bool)swap;
else
estimateOut = false;
}
}
double ConditionalEstimator::estimate(double ts)
{
if (decoratedEstimator_) {
double partition_sum = 0.0;
elementCount_ = 0;
decoratedEstimator_->init(ts);
mesh = decoratedEstimator_->getMesh();
traverseFlag = decoratedEstimator_->getTraverseFlag();
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
while (elInfo) {
PartitionElementData *elData = dynamic_cast<PartitionElementData*>
(elInfo->getElement()->getElementData(PARTITION_ED));
TEST_EXIT_DBG(elInfo->getElement()->isLeaf())("element not leaf\n");
TEST_EXIT_DBG(elData)("no partition data on leaf element %d (rank %d)\n",
elInfo->getElement()->getIndex(),
MPI::COMM_WORLD.Get_rank());
PartitionStatus status = elData->getPartitionStatus();
if (status == IN || status == OVERLAP || (estimateOut && status==OUT))
decoratedEstimator_->estimateElement(elInfo);
else
elInfo->getElement()->setMark(0);
elInfo = stack.traverseNext(elInfo);
}
elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
while (elInfo) {
PartitionElementData *elData = dynamic_cast<PartitionElementData*>
(elInfo->getElement()->getElementData(PARTITION_ED));
PartitionStatus status = elData->getPartitionStatus();
if (status == IN || status==OVERLAP || (estimateOut && status==OUT)) {
elementCount_++;
partition_sum += elInfo->getElement()->getEstimation(row_);
}
elInfo = stack.traverseNext(elInfo);
}
// !!! test !!!
double total_sum = 0.0;
MPI::COMM_WORLD.Allreduce(&partition_sum,
&total_sum,
1,
MPI_DOUBLE,
MPI_SUM);
decoratedEstimator_->setErrorSum(total_sum);
double total_max = 0.0;
est_max = decoratedEstimator_->getErrorMax();
MPI::COMM_WORLD.Allreduce(&est_max,
&total_max,
1,
MPI_DOUBLE,
MPI_MAX);
decoratedEstimator_->setErrorMax(total_max);
decoratedEstimator_->exit();
est_sum = decoratedEstimator_->getErrorSum();
est_max = decoratedEstimator_->getErrorMax();
// !!! test !!!
#if 0
decoratedEstimator_->exit();
est_sum = sqrt(partition_sum); //decoratedEstimator_->getErrorSum();
est_t_sum = decoratedEstimator_->getTimeEst();
est_max = decoratedEstimator_->getErrorMax();
#endif
// MSG("rank %d , estimate %e (total %e) elements %d (total %d)\n",
// MPI::COMM_WORLD.Get_rank(), est_sum, total_sum,
// elementCount_,
// mesh->getNumberOfLeaves());
return est_sum;
} else {
return 0.0;
}
}
}
#include "ConditionalMarker.h"
namespace AMDiS
{
ConditionalMarker::ConditionalMarker(const std::string name_,
int row,
Marker *decoratedMarker,
int globalCoarseGridLevel,
int localCoarseGridLevel)
: Marker(name_, row),
decoratedMarker_(decoratedMarker),
globalCoarseGridLevel_(globalCoarseGridLevel),
localCoarseGridLevel_(localCoarseGridLevel)
{
if (decoratedMarker != NULL)
name = decoratedMarker->getName();
else
std::cout << "no decorated Marker, mark outer path is" << name << "->mark outer" << std::endl;
int swap;
int retVal = Parameters::getGlobalParameter(0, name + "->mark outer", "%d", &swap);
if (retVal>0)
markOuter = (bool)swap;
else
markOuter = false;
MSG("markOuter:" + markOuter);
}
void ConditionalMarker::initMarking(AdaptInfo *adaptInfo, Mesh *mesh)
{
if (decoratedMarker_)
decoratedMarker_->initMarking(adaptInfo, mesh);
}
void ConditionalMarker::finishMarking(AdaptInfo *adaptInfo)
{
if (decoratedMarker_) {
int tmp = decoratedMarker_->getElMarkRefine();
MPI::COMM_WORLD.Allreduce(&tmp,
&elMarkRefine,
1,
MPI_INT,
MPI_SUM);
tmp = decoratedMarker_->getElMarkCoarsen();
MPI::COMM_WORLD.Allreduce(&tmp,
&elMarkCoarsen,
1,
MPI_INT,
MPI_SUM);
decoratedMarker_->finishMarking(adaptInfo);
}
}
void ConditionalMarker::markElement(AdaptInfo *adaptInfo, ElInfo *elInfo)
{
FUNCNAME("ConditionalMarker::markElement()");
if (decoratedMarker_) {
PartitionElementData *elData =
dynamic_cast<PartitionElementData*>(elInfo->getElement()->
getElementData(PARTITION_ED));
TEST_EXIT_DBG(elData)("no partition data\n");
decoratedMarker_->markElement(adaptInfo, elInfo);
if (!markOuter && elData->getPartitionStatus() == OUT) {
Element *element = elInfo->getElement();
if (element->getMark() > 0)
element->setMark(0);
}
int minLevel =
elData->getPartitionStatus() != OUT ?
localCoarseGridLevel_ :
globalCoarseGridLevel_;
if (elData->getLevel() + elInfo->getElement()->getMark() < minLevel)
elInfo->getElement()->setMark(-(elData->getLevel() - minLevel));
}
}
}
// ============================================================================
// == ==
// == AMDiS - Adaptive multidimensional simulations ==
// == ==
// ============================================================================
// == ==
// == crystal growth group ==
// == ==
// == Stiftung caesar ==
// == Ludwig-Erhard-Allee 2 ==
// == 53175 Bonn ==
// == germany ==
// == ==
// ============================================================================
// == ==
// == http://www.caesar.de/cg/AMDiS ==
// == ==
// ============================================================================
/** \file ConditionalMarker.h */
#ifndef AMDIS_CONDITIONALMARKER_H
#define AMDIS_CONDITIONALMARKER_H
#include "mpi.h"
#include "Marker.h"
#include "Marker.h"
#include "PartitionElementData.h"
#include "Traverse.h"
#include "ElInfo.h"
namespace AMDiS {
/**
* \ingroup Adaption
*
* \brief
*
*/
class ConditionalMarker : public Marker
{
public:
/// Constructor.
ConditionalMarker(const std::string name,
int row,
Marker *decoratedMarker,
int globalCoarseGridLevel,
int localCoarseGridLevel);
/// Implementation of MarkScal::initMarking().
virtual void initMarking(AdaptInfo *adaptInfo, Mesh *mesh);
virtual void finishMarking(AdaptInfo *adaptInfo);
void markElement(AdaptInfo *adaptInfo, ElInfo *elInfo);
protected:
Marker *decoratedMarker_;
int globalCoarseGridLevel_;
int localCoarseGridLevel_;
bool markOuter;
};
}
#endif
......@@ -32,7 +32,7 @@ namespace AMDiS {
public:
static FixVec<T,d1>& convertVec(FixVec<T,d2>& rhs, Mesh* mesh)
{
TEST_EXIT(mesh->getGeo(d1)==mesh->getGeo(d2))("Incompatible dimensions.\n");
TEST_EXIT(mesh->getGeo(d1) == mesh->getGeo(d2))("Incompatible dimensions!\n");
return reinterpret_cast<FixVec<T,d1>&>(rhs);
}
};
......
#include "MacroInfo.h"
#include "Mesh.h"
#include "MacroReader.h"
#include "FixVecConvert.h"
#include "SurfaceRegion_ED.h"
#include "ElementRegion_ED.h"
namespace AMDiS {
void MacroInfo::fill(Mesh *pmesh, int nElements, int nVertices)
{
FUNCNAME("MacroInfo::fill()");
TEST_EXIT(pmesh)("no mesh\n");
int dim = pmesh->getDim();
mesh = pmesh;
mesh->setNumberOfElements(nElements);
mesh->setNumberOfLeaves(nElements);
mesh->setNumberOfVertices(nVertices);
for (int i = 0; i < nElements; i++) {
MacroElement *newMacro = new MacroElement(mesh->getDim());
mel.push_back(newMacro);
mesh->addMacroElement(mel[i]);
}
dof = new DegreeOfFreedom*[nVertices];
coords = new WorldVector<double>[nVertices];
mel_vertex = new int*[nElements];
for (int i = 0; i < nElements; i++)
mel_vertex[i] = new int[mesh->getGeo(VERTEX)];
for (int i = 0; i < nVertices; i++)
dof[i] = mesh->getDof(VERTEX);
for (int i = 0; i < nElements; i++) {
mel[i]->element = mesh->createNewElement();
mel[i]->index = i;
if (dim == 3)
mel[i]->elType = 0;
}
neigh_set = false;
bound_set = false;
}
void MacroInfo::clear()
{
for (int i = 0; i < mesh->getNumberOfMacros(); i++)
delete [] mel_vertex[i];
delete [] mel_vertex;
delete [] coords;
coords = NULL;
delete [] dof;
dof = NULL;
mesh = NULL;
neigh_set = false;
}
/****************************************************************************/
/* read_indices() reads dim+1 indices from file into id[0-dim], */
/* returns true if dim+1 inputs arguments could be read successfully by */
/* fscanf(), else false */
/****************************************************************************/
int MacroInfo::read_indices(FILE *file, DimVec<int> &id)
{
int dim = mesh->getDim();
for (int i = 0; i <= dim; i++)
if (fscanf(file, "%d", &id[i]) != 1)
return(false);
return(true);
}
#define N_KEYS 14
#define N_MIN_KEYS 7
static const char *keys[N_KEYS] = {
"DIM", // 0
"DIM_OF_WORLD", // 1
"number of vertices", // 2
"number of elements", // 3
"vertex coordinates", // 4
"element vertices", // 5
"element boundaries", // 6
"element neighbours", // 7
"element type", // 8
"projections", // 9
"element region", // 10
"surface region", // 11
"mesh name", // 12
"time" // 13
};
static int get_key_no(const char *key)
{
for (int i = 0; i < N_KEYS; i++)
if (!strcmp(keys[i], key))
return i;
return -1;
}
#include <ctype.h>
static const char *read_key(const char *line)
{
static char key[100];
char *k = key;
while (isspace(*line))
line++;
while ((*k++ = *line++) != ':');
*--k = '\0';
return const_cast<const char *>(key);
}
void MacroInfo::readAMDiSMacro(std::string filename, Mesh* mesh)
{
FUNCNAME("MacroInfo::readAMDiSMacro()");
int dim, dow;
int nElements, nVertices;
int j, k;
double dbl;
char line[256];
int line_no, n_keys, sort_key[N_KEYS], nv_key, ne_key;
int key_def[N_KEYS] = {0,0,0,0,0,0,0,0,0,0,0,0};
const char *key;
DimVec<int> *ind = NULL;
TEST_EXIT(filename != "")("No filename specified!\n");
TEST_EXIT(filename.length() < static_cast<unsigned int>(127))
("can only handle filenames up to 127 characters\n");
FILE *file = fopen(filename.c_str(), "r");
TEST_EXIT(file)("cannot open file %s\n", filename.c_str());
// === Looking for all keys in the macro file. ===
line_no = n_keys = 0;
while (fgets(line, 255, file)) {
line_no++;
if (!strchr(line, ':'))
continue;
key = read_key(line);
int i_key = get_key_no(key);
TEST_EXIT(i_key >= 0)
("macro file %s must not contain key %s on line %d\n",
filename.c_str(), key, line_no);
TEST_EXIT(!key_def[i_key])("key %s defined second time on line %d in file %s\n");
sort_key[n_keys++] = i_key;
key_def[i_key] = true;
}
fclose(file);
// === Test, if there is data for every key and if all is defined in ===
// === right order. ===
for (int i_key = 0; i_key < N_MIN_KEYS; i_key++) {
for (j = 0; j < n_keys; j++)
if (sort_key[j] == i_key)
break;
TEST_EXIT(j < n_keys)("You do not have specified data for %s in %s\n",
keys[i_key], filename.c_str());
for (j = 0; j < n_keys; j++)
if (sort_key[j] == 2) break;
nv_key = j;
for (j = 0; j < n_keys; j++)
if (sort_key[j] == 3) break;
ne_key = j;
switch (i_key) {
case 0:
case 1:
TEST_EXIT(sort_key[i_key] < 2)
("You have to specify DIM or mesh->getGeo(WORLD) before all other data\n");
break;
case 4:
TEST_EXIT(nv_key < i_key)
("Before reading data for %s, you have to specify the %s in file\n",
keys[4], keys[2], filename.c_str());
break;
case 5:
TEST_EXIT(nv_key < i_key && ne_key < i_key)
("Before reading data for %s, you have to specify the %s and %s in file %s\n",
keys[5], keys[3], keys[2], filename.c_str());
case 6:
case 7:
case 8:
TEST_EXIT(ne_key < i_key)
("Before reading data for %s, you have to specify the %s in file %s\n",
keys[i_key], keys[3], filename.c_str());
}
}
for (int i_key = 0; i_key < N_KEYS; i_key++)
key_def[i_key] = false;
// === And now, reading data. ===
file = fopen(filename.c_str(), "r");
TEST_EXIT(file)("cannot open file %s\n", filename.c_str());
int result;
for (int i_key = 0; i_key < n_keys; i_key++) {
switch (sort_key[i_key]) {
case 0:
// line "DIM"
result = fscanf(file, "%*s %d", &dim);
TEST_EXIT(result == 1)("cannot read DIM correctly in file %s\n", filename.c_str());
ind = new DimVec<int>(dim, NO_INIT);
key_def[0] = true;
break;
case 1:
// line "DIM_OF_WORLD"
result = fscanf(file, "%*s %d", &dow);
TEST_EXIT(result == 1)
("cannot read Global::getGeo(WORLD) correctly in file %s\n", filename.c_str());