From be06847b30df4dce8dc92ec58fa635fe86936372 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 9 Jun 2009 16:20:40 +0000
Subject: [PATCH] Work in progress on parallel domain decomposition.
---
AMDiS/src/DataCollector.cc | 106 ++++++++++-------------------
AMDiS/src/DataCollector.h | 14 ++--
AMDiS/src/ParallelDomainProblem.cc | 78 +++++++++++++++++----
AMDiS/src/ParallelDomainProblem.h | 4 +-
AMDiS/src/VtkWriter.cc | 3 +-
AMDiS/src/VtkWriter.h | 26 +++----
AMDiS/src/VtkWriter.hh | 11 ++-
7 files changed, 128 insertions(+), 114 deletions(-)
diff --git a/AMDiS/src/DataCollector.cc b/AMDiS/src/DataCollector.cc
index b533bc57..30feb017 100644
--- a/AMDiS/src/DataCollector.cc
+++ b/AMDiS/src/DataCollector.cc
@@ -6,6 +6,7 @@
#include "Projection.h"
namespace AMDiS {
+
DataCollector::DataCollector(const FiniteElemSpace *feSpace,
DOFVector<double> *values,
int level,
@@ -61,20 +62,17 @@ namespace AMDiS {
void DataCollector::fillAllData()
{
- if (!elementDataCollected_) {
+ if (!elementDataCollected_)
startCollectingElementData();
- }
- if (!periodicDataCollected_) {
+ if (!periodicDataCollected_)
startCollectingPeriodicData();
- }
- if (!valueDataCollected_) {
+ if (!valueDataCollected_)
startCollectingValueData();
- }
}
- int DataCollector::startCollectingElementData()
+ void DataCollector::startCollectingElementData()
{
FUNCNAME("DataCollector::startCollectingElementData()");
@@ -99,26 +97,24 @@ namespace AMDiS {
// Traverse elements to create element information
elInfo = stack.traverseFirst(mesh_, level_, flag);
+
while (elInfo) {
- if (!writeElem_ || writeElem_(elInfo)) {
+ if (!writeElem_ || writeElem_(elInfo))
addElementData(elInfo);
- }
+
elInfo = stack.traverseNext(elInfo);
}
elementDataCollected_ = true;
-
- return(0);
}
- int DataCollector::startCollectingValueData()
+ void DataCollector::startCollectingValueData()
{
FUNCNAME("DataCollector::startCollectingValueData()");
DOFVector<int>::Iterator intPointIt(interpPointInd_, USED_DOFS);
- for (intPointIt.reset(); !intPointIt.end(); ++intPointIt) {
+ for (intPointIt.reset(); !intPointIt.end(); ++intPointIt)
(*intPointIt) = -1;
- }
interpPoints_.clear();
@@ -130,13 +126,11 @@ namespace AMDiS {
// Traverse elements to add value information and to mark
// interpolation points.
- ElInfo *elInfo = stack.traverseFirst(mesh_,
- level_,
+ ElInfo *elInfo = stack.traverseFirst(mesh_, level_,
traverseFlag_ | Mesh::FILL_COORDS);
while (elInfo) {
- if (!writeElem_ || writeElem_(elInfo)) {
+ if (!writeElem_ || writeElem_(elInfo))
addValueData(elInfo);
- }
elInfo = stack.traverseNext(elInfo);
}
@@ -146,11 +140,9 @@ namespace AMDiS {
// Remove all interpolation marks and, instead, set to each
// interpolation point its continous index starting from 0.
int i = 0;
- for (intPointIt.reset(); !intPointIt.end(); ++intPointIt) {
- if (*intPointIt == -3) {
+ for (intPointIt.reset(); !intPointIt.end(); ++intPointIt)
+ if (*intPointIt == -3)
*intPointIt = i++;
- }
- }
// Traverse elements to create interpolation values.
elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_ | Mesh::FILL_COORDS);
@@ -163,11 +155,9 @@ namespace AMDiS {
delete [] localDOFs_;
valueDataCollected_ = true;
-
- return 0;
}
- int DataCollector::startCollectingPeriodicData()
+ void DataCollector::startCollectingPeriodicData()
{
FUNCNAME("DataCollector::startCollectingPeriodicData()");
@@ -210,11 +200,9 @@ namespace AMDiS {
}
periodicDataCollected_ = true;
-
- return(0);
}
- int DataCollector::addElementData(ElInfo* elInfo)
+ void DataCollector::addElementData(ElInfo* elInfo)
{
FUNCNAME("DataCollector::addElementData()");
@@ -280,19 +268,15 @@ namespace AMDiS {
outputIndices_[elInfo->getNeighbour(i)->getIndex()] :
-1;
}
-
- if (dim_ == 3) {
+ if (dim_ == 3)
elementInfo.type = (dynamic_cast<ElInfo3d*>(elInfo)->getType());
- }
// remember element info
elements_.push_back(elementInfo);
-
- return(0);
}
- int DataCollector::addValueData(ElInfo *elInfo)
+ void DataCollector::addValueData(ElInfo *elInfo)
{
FUNCNAME("DataCollector::addValueData()");
@@ -345,28 +329,24 @@ namespace AMDiS {
nInterpPoints_++;
}
}
- }
-
- return(0);
+ }
}
- int DataCollector::addInterpData(ElInfo *elInfo)
+ void DataCollector::addInterpData(ElInfo *elInfo)
{
FUNCNAME("DataCollector::addInterpData()");
basisFcts_->getLocalIndices(elInfo->getElement(), localAdmin_, localDOFs_);
std::vector<int> elemInterpPoints(0);
- for (int i = mesh_->getGeo(VERTEX); i < nBasisFcts_; i++) {
+ for (int i = mesh_->getGeo(VERTEX); i < nBasisFcts_; i++)
elemInterpPoints.push_back((*interpPointInd_)[localDOFs_[i]]);
- }
- interpPoints_.push_back(elemInterpPoints);
-
- return(0);
+ interpPoints_.push_back(elemInterpPoints);
}
- int DataCollector::addPeriodicData(ElInfo *elInfo) {
+ void DataCollector::addPeriodicData(ElInfo *elInfo)
+ {
FUNCNAME("DataCollector::addPeriodicData");
LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
@@ -423,69 +403,60 @@ namespace AMDiS {
}
}
}
-
- return(0);
}
std::list<ElementInfo>* DataCollector::getElementInfos()
{
- if (!elementDataCollected_) {
+ if (!elementDataCollected_)
startCollectingElementData();
- }
return &elements_;
}
DOFVector< std::list<VertexInfo> >* DataCollector::getVertexInfos()
{
- if (!elementDataCollected_) {
+ if (!elementDataCollected_)
startCollectingElementData();
- }
return vertexInfos_;
}
std::list<PeriodicInfo>* DataCollector::getPeriodicInfos()
{
- if (!periodicDataCollected_) {
+ if (!periodicDataCollected_)
startCollectingPeriodicData();
- }
return &periodicInfos_;
}
int DataCollector::getNumberVertices()
{
- if (!elementDataCollected_) {
+ if (!elementDataCollected_)
startCollectingElementData();
- }
return nVertices_;
}
int DataCollector::getNumberElements()
{
- if (!elementDataCollected_) {
+ if (!elementDataCollected_)
startCollectingElementData();
- }
return nElements_;
}
int DataCollector::getNumberInterpPoints()
{
- if (!valueDataCollected_) {
+ if (!valueDataCollected_)
startCollectingValueData();
- }
return nInterpPoints_;
}
int DataCollector::getNumberConnections()
{
- if (!periodicDataCollected_) {
+ if (!periodicDataCollected_)
startCollectingPeriodicData();
- }
return nConnections_;
}
@@ -502,45 +473,40 @@ namespace AMDiS {
DOFVector<double>* DataCollector::getValues()
{
- if (!valueDataCollected_) {
+ if (!valueDataCollected_)
startCollectingValueData();
- }
return values_;
}
DOFVector< std::list<WorldVector<double> > >* DataCollector::getDofCoords()
{
- if (!valueDataCollected_) {
+ if (!valueDataCollected_)
startCollectingValueData();
- }
return dofCoords_;
}
DOFVector<int>* DataCollector::getInterpPointInd()
{
- if (!valueDataCollected_) {
+ if (!valueDataCollected_)
startCollectingValueData();
- }
return interpPointInd_;
}
DOFVector< std::list<WorldVector<double> > >* DataCollector::getInterpPointCoords()
{
- if (!valueDataCollected_) {
+ if (!valueDataCollected_)
startCollectingValueData();
- }
return interpPointCoords_;
}
std::vector< std::vector<int> >* DataCollector::getInterpPoints()
{
- if (!valueDataCollected_) {
+ if (!valueDataCollected_)
startCollectingValueData();
- }
return &interpPoints_;
}
diff --git a/AMDiS/src/DataCollector.h b/AMDiS/src/DataCollector.h
index d46af551..57234325 100644
--- a/AMDiS/src/DataCollector.h
+++ b/AMDiS/src/DataCollector.h
@@ -101,25 +101,25 @@ namespace AMDiS {
protected:
/// Start collecting element and vertex data of the problem.
- int startCollectingElementData();
+ void startCollectingElementData();
/// Start collecting value data of the problem.
- int startCollectingValueData();
+ void startCollectingValueData();
/// Start collecting periodic data of the problem.
- int startCollectingPeriodicData();
+ void startCollectingPeriodicData();
/// Adds information about one element and its vertices.
- int addElementData(ElInfo* elInfo);
+ void addElementData(ElInfo* elInfo);
/// Adds value information of one element.
- int addValueData(ElInfo *elInfo);
+ void addValueData(ElInfo *elInfo);
/// Adds information about interpolation points of vertices.
- int addInterpData(ElInfo *elInfo);
+ void addInterpData(ElInfo *elInfo);
/// Adds value information of one element.
- int addPeriodicData(ElInfo *elInfo);
+ void addPeriodicData(ElInfo *elInfo);
/// Vector with vertex values
DOFVector<double> *values_;
diff --git a/AMDiS/src/ParallelDomainProblem.cc b/AMDiS/src/ParallelDomainProblem.cc
index 7cf96050..126ba828 100644
--- a/AMDiS/src/ParallelDomainProblem.cc
+++ b/AMDiS/src/ParallelDomainProblem.cc
@@ -10,6 +10,9 @@
#include "PartitionElementData.h"
#include "DOFMatrix.h"
#include "DOFVector.h"
+#include "VtkWriter.h"
+
+#include "petscksp.h"
namespace AMDiS {
@@ -138,9 +141,8 @@ namespace AMDiS {
PartitionElementData *partitionData =
dynamic_cast<PartitionElementData*>
((*it)->getElement()->getElementData(PARTITION_ED));
- if (partitionData->getPartitionStatus() != IN) {
+ if (partitionData->getPartitionStatus() != IN)
macrosToRemove.push_back(*it);
- }
}
mesh->removeMacroElements(macrosToRemove);
@@ -151,17 +153,15 @@ namespace AMDiS {
int *lOrder = (int*)(malloc(sizeof(int) * rankDofs.size()));
for (std::vector<const DegreeOfFreedom*>::iterator it = rankDofs.begin();
- it != rankDofs.end(); ++it) {
+ it != rankDofs.end(); ++it)
gOrder[nRankDOFs++] = (*it)[0];
- }
int rstart = 0;
MPI_Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
rstart -= nRankDOFs;
- for (int i = 0; i < nRankDOFs; i++) {
+ for (int i = 0; i < nRankDOFs; i++)
lOrder[i] = rstart + i;
- }
AOCreateBasic(PETSC_COMM_WORLD, nRankDOFs, gOrder, lOrder, &applicationOrdering);
@@ -275,6 +275,10 @@ namespace AMDiS {
ierr = VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
ierr = VecSetSizes(petscRhsVec, rankDofs.size(), partitionDofs.size());
ierr = VecSetType(petscRhsVec, VECMPI);
+
+ ierr = VecCreate(PETSC_COMM_WORLD, &petscSolVec);
+ ierr = VecSetSizes(petscSolVec, rankDofs.size(), partitionDofs.size());
+ ierr = VecSetType(petscSolVec, VECMPI);
}
void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo)
@@ -285,14 +289,27 @@ namespace AMDiS {
void ParallelDomainProblemBase::fillPetscMatrix(DOFMatrix *mat,
DOFVector<double> *vec)
{
- /* DOFMatrix::Iterator rowIt(mat, USED_DOFS);
- for (rowIt.reset(); !rowIt.end(); ++rowIt) {
- for (int i = 0; i < static_cast<int>((*rowIt).size()); i++) {
- if ((*rowIt)[i].col >= 0) {
- MatSetValues(petscMatrix, 1, &i, 1, &((*rowIt)[i].col), &((*rowIt)[i].entry), ADD_VALUES);
+ using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
+ namespace traits= mtl::traits;
+ typedef DOFMatrix::base_matrix_type Matrix;
+
+ traits::row<Matrix>::type row(mat->getBaseMatrix());
+ traits::col<Matrix>::type col(mat->getBaseMatrix());
+ traits::const_value<Matrix>::type value(mat->getBaseMatrix());
+
+ typedef traits::range_generator<major, Matrix>::type cursor_type;
+ typedef traits::range_generator<nz, cursor_type>::type icursor_type;
+
+ std::cout.precision(10);
+ for (cursor_type cursor = begin<major>(mat->getBaseMatrix()), cend = end<major>(mat->getBaseMatrix()); cursor != cend; ++cursor)
+ for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
+ if (value(*icursor) != 0.0) {
+ int r = row(*icursor);
+ int c = col(*icursor);
+ double v = value(*icursor);
+ MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
}
- }
- }
+
MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);
@@ -303,7 +320,38 @@ namespace AMDiS {
double value = *dofIt;
VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
- }*/
+ }
+ }
+
+ void ParallelDomainProblemBase::solvePetscMatrix(DOFVector<double> *vec)
+ {
+ KSP ksp;
+ PC pc;
+
+ KSPCreate(PETSC_COMM_WORLD, &ksp);
+ KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
+ KSPGetPC(ksp, &pc);
+ PCSetType(pc, PCJACOBI);
+ KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
+ KSPSetType(ksp, KSPBCGS);
+ KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, 0);
+ KSPSolve(ksp, petscRhsVec, petscSolVec);
+
+ PetscScalar *vecPointer;
+ VecGetArray(petscSolVec, &vecPointer);
+
+ DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
+ int counter = 0;
+ for (dofIt.reset(); !dofIt.end(); ++dofIt)
+ *vec = vecPointer[counter++];
+
+ VecRestoreArray(petscSolVec, &vecPointer);
+
+ std::stringstream oss;
+ oss << "test" << MPI::COMM_WORLD.Get_rank() << ".vtu";
+ VtkWriter::writeFile(vec, oss.str());
+
+ std::cout << "USED SIZE = " << vec->getUsedSize() << std::endl;
}
double ParallelDomainProblemBase::setElemWeights(AdaptInfo *adaptInfo)
@@ -370,6 +418,8 @@ namespace AMDiS {
fillPetscMatrix(probScal->getSystemMatrix(), probScal->getRHS());
+ solvePetscMatrix(probScal->getSolution());
+
// if (toDo.isSet(SOLVE))
// iterationIF->getProblem()->solve(adaptInfo, false);
diff --git a/AMDiS/src/ParallelDomainProblem.h b/AMDiS/src/ParallelDomainProblem.h
index 3da0ec43..c5ad751d 100644
--- a/AMDiS/src/ParallelDomainProblem.h
+++ b/AMDiS/src/ParallelDomainProblem.h
@@ -105,6 +105,8 @@ namespace AMDiS {
void fillPetscMatrix(DOFMatrix *mat, DOFVector<double> *vec);
+ void solvePetscMatrix(DOFVector<double> *vec);
+
virtual ProblemStatBase *getProblem(int number = 0) {}
virtual void serialize(std::ostream&) {}
@@ -168,7 +170,7 @@ namespace AMDiS {
Vec petscRhsVec;
- Vec petscSolutionVec;
+ Vec petscSolVec;
/// Number of DOFs in the rank mesh.
int nRankDOFs;
diff --git a/AMDiS/src/VtkWriter.cc b/AMDiS/src/VtkWriter.cc
index c995213e..151ef4d4 100644
--- a/AMDiS/src/VtkWriter.cc
+++ b/AMDiS/src/VtkWriter.cc
@@ -75,8 +75,7 @@ namespace AMDiS {
return 0;
}
- void VtkWriter::writeFile(DOFVector<double> *values,
- const char *filename)
+ void VtkWriter::writeFile(DOFVector<double> *values, const std::string filename)
{
DataCollector *dc = new DataCollector(values->getFESpace(), values);
diff --git a/AMDiS/src/VtkWriter.h b/AMDiS/src/VtkWriter.h
index 65d71fca..645ce6d9 100644
--- a/AMDiS/src/VtkWriter.h
+++ b/AMDiS/src/VtkWriter.h
@@ -51,11 +51,11 @@ namespace AMDiS {
int writeFile(const std::string name);
/// May be used to simply write ParaView files.
- static void writeFile(DOFVector<double> *values,
- const char *filename);
+ static void writeFile(DOFVector<double> *values, const std::string filename);
/// Set a compressing method for file output.
- void setCompression(FileCompression c) {
+ void setCompression(FileCompression c)
+ {
compress = c;
}
@@ -95,25 +95,25 @@ namespace AMDiS {
#ifdef HAVE_BOOST
/// Writes a world coordinate to a given file.
inline void writeCoord(boost::iostreams::filtering_ostream &file,
- WorldVector<double> coord) {
- for (int i = 0; i < Global::getGeo(WORLD); i++) {
+ WorldVector<double> coord)
+ {
+ for (int i = 0; i < Global::getGeo(WORLD); i++)
file << " " << std::scientific << coord[i];
- }
- for (int i = Global::getGeo(WORLD); i < 3; i++) {
+ for (int i = Global::getGeo(WORLD); i < 3; i++)
file << " 0.0";
- }
+
file << "\n";
}
#endif
/// Writes a world coordinate to a given file.
- inline void writeCoord(std::ofstream &file, WorldVector<double> coord) {
- for (int i = 0; i < Global::getGeo(WORLD); i++) {
+ inline void writeCoord(std::ofstream &file, WorldVector<double> coord)
+ {
+ for (int i = 0; i < Global::getGeo(WORLD); i++)
file << " " << std::scientific << coord[i];
- }
- for (int i = Global::getGeo(WORLD); i < 3; i++) {
+ for (int i = Global::getGeo(WORLD); i < 3; i++)
file << " 0.0";
- }
+
file << "\n";
}
diff --git a/AMDiS/src/VtkWriter.hh b/AMDiS/src/VtkWriter.hh
index 42f39a44..856285ac 100644
--- a/AMDiS/src/VtkWriter.hh
+++ b/AMDiS/src/VtkWriter.hh
@@ -97,7 +97,7 @@ namespace AMDiS {
for (it2 = it->begin(); it2 != it->end(); ++it2) {
it2->outputIndex = counter++;
writeCoord(file, it2->coords);
- }
+ }
}
// For the second dim case, write also the interpolation points.
@@ -107,9 +107,8 @@ namespace AMDiS {
for (pointIt.reset(); !pointIt.end(); ++pointIt) {
std::list< WorldVector<double> >::iterator it2;
- for (it2 = pointIt->begin(); it2 != pointIt->end(); ++it2) {
+ for (it2 = pointIt->begin(); it2 != pointIt->end(); ++it2)
writeCoord(file, *it2);
- }
}
}
}
@@ -147,9 +146,8 @@ namespace AMDiS {
++intPointIt, ++valueIt, ++coordIt) {
if (*intPointIt == -2) {
- for (int i = 0; i < static_cast<int>(coordIt->size()); i++) {
+ for (int i = 0; i < static_cast<int>(coordIt->size()); i++)
file << " " << (fabs(*valueIt) < 1e-40 ? 0.0 : *valueIt) << "\n";
- }
}
}
@@ -163,9 +161,8 @@ namespace AMDiS {
++intPointIt, ++valueIt, ++interpCoordIt) {
if (*intPointIt >= 0) {
- for (int i = 0; i < static_cast<int>(interpCoordIt->size()); i++) {
+ for (int i = 0; i < static_cast<int>(interpCoordIt->size()); i++)
file << " " << (fabs(*valueIt) < 1e-40 ? 0.0 : *valueIt) << "\n";
- }
}
}
}
--
GitLab