diff --git a/AMDiS/src/ProblemTimeInterface.h b/AMDiS/src/ProblemTimeInterface.h index 602d6e769083c8a6af7d38a4c86606ad9e4f1fe2..5a68d00a1b5c4a1fb987e0397fbb7cc374e79543 100644 --- a/AMDiS/src/ProblemTimeInterface.h +++ b/AMDiS/src/ProblemTimeInterface.h @@ -44,39 +44,25 @@ namespace AMDiS { public: virtual ~ProblemTimeInterface() {}; - /** \brief - * Executes all needed operations when the simulation time changes. - */ + /// Executes all needed operations when the simulation time changes. virtual void setTime(AdaptInfo *adaptInfo) = 0; - /** \brief - * Called at the beginning of each timestep - */ + /// Called at the beginning of each timestep virtual void initTimestep(AdaptInfo *adaptInfo) = 0; - /** \brief - * Called at the end of each timestep. - */ + /// Called at the end of each timestep. virtual void closeTimestep(AdaptInfo *adaptInfo) = 0; - /** \brief - * Solves the initial problem. - */ + /// Solves the initial problem. virtual void solveInitialProblem(AdaptInfo *adaptInfo) = 0; - /** \brief - * Solves the initial problem. - */ + /// Solves the initial problem. virtual void transferInitialSolution(AdaptInfo *adaptInfo) = 0; - /** \brief - * Function that serializes the problem plus information about the iteration. - */ + /// Function that serializes the problem plus information about the iteration. virtual void serialize(std::ostream &out) = 0; - /** \brief - * Function that deserializes the problem plus information about the iteration. - */ + /// Function that deserializes the problem plus information about the iteration. virtual void deserialize(std::istream &in) = 0; }; diff --git a/AMDiS/src/parallel/InteriorBoundary.h b/AMDiS/src/parallel/InteriorBoundary.h index 3dc8d9572a751a65cb2867e73a95583e67038ec0..aa19f87598e3a4453c3453ec910625de30a0df8d 100644 --- a/AMDiS/src/parallel/InteriorBoundary.h +++ b/AMDiS/src/parallel/InteriorBoundary.h @@ -177,6 +177,9 @@ namespace AMDiS { inline void nextNonempty() { + if (mapIt == bound.boundary.end()) + return; + while (mapIt->second.size() == 0) { ++mapIt; if (mapIt == bound.boundary.end()) diff --git a/AMDiS/src/parallel/ParallelDomainBase.cc b/AMDiS/src/parallel/ParallelDomainBase.cc index 572d6e0241c3208e7fe8828d3b8b997430fb713b..3057c509cf4e42c0b53f7ff82324b99b266f86e2 100644 --- a/AMDiS/src/parallel/ParallelDomainBase.cc +++ b/AMDiS/src/parallel/ParallelDomainBase.cc @@ -1,7 +1,7 @@ #include #include #include - +#include #include "parallel/ParallelDomainBase.h" #include "parallel/StdMpi.h" #include "ParMetisPartitioner.h" @@ -21,12 +21,16 @@ #include "ElementFileWriter.h" #include "VertexVector.h" #include "MeshStructure.h" +#include "ProblemVec.h" +#include "ProblemInstat.h" #include "Debug.h" #include "petscksp.h" namespace AMDiS { + using boost::lexical_cast; + PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *) { if (iter % 10 == 0 && MPI::COMM_WORLD.Get_rank() == 0) @@ -40,22 +44,22 @@ namespace AMDiS { return (*dof1 < *dof2); } - ParallelDomainBase::ParallelDomainBase(ProblemIterationInterface *iIF, - ProblemTimeInterface *tIF, - FiniteElemSpace *fe, - RefinementManager *refinementManager) - : iterationIF(iIF), - timeIF(tIF), - name(iIF->getName()), - feSpace(fe), - mesh(fe->getMesh()), - refineManager(refinementManager), + ParallelDomainBase::ParallelDomainBase(ProblemVec *problemStat, + ProblemInstatVec *problemInstat) + : iterationIF(problemStat), + timeIF(problemInstat), + probStat(problemStat), + name(problemStat->getName()), + feSpace(problemStat->getFESpace(0)), + mesh(feSpace->getMesh()), + refineManager(problemStat->getRefinementManager(0)), + info(problemStat->getInfo()), initialPartitionMesh(true), d_nnz(NULL), o_nnz(NULL), nRankDofs(0), rstart(0), - nComponents(1), + nComponents(problemStat->getNumComponents()), deserialized(false), lastMeshChangeIndex(0) { @@ -70,8 +74,34 @@ namespace AMDiS { mpiSize = MPI::COMM_WORLD.Get_size(); mpiComm = MPI::COMM_WORLD; partitioner = new ParMetisPartitioner(mesh, &mpiComm); + + // Test if all fe spaces are equal. Yet, different fe spaces are not supported for + // domain parallelization. + const FiniteElemSpace *fe = probStat->getFESpace(0); + for (int i = 0; i < nComponents; i++) + TEST_EXIT(fe == probStat->getFESpace(i)) + ("Parallelization does not supported different FE spaces!\n"); + + // Create parallel serialization file writer, if needed. + int writeSerialization = 0; + GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization); + if (writeSerialization) + probStat->getFileWriterList().push_back(new Serializer(this)); + + int readSerialization = 0; + GET_PARAMETER(0, name + "->input->read serialization", "%d", &readSerialization); + if (readSerialization) { + std::string filename = ""; + GET_PARAMETER(0, name + "->input->serialization filename", &filename); + filename += ".p" + lexical_cast(mpiRank); + MSG("Start serialization with %s\n", filename.c_str()); + std::ifstream in(filename.c_str()); + deserialize(in); + in.close(); + } } + void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo) { FUNCNAME("ParallelDomainBase::initParallelization()"); @@ -157,6 +187,46 @@ namespace AMDiS { createPeriodicMap(); } + + + /// === Set DOF rank information to all matrices and vectors. === + + for (int i = 0; i < nComponents; i++) { + for (int j = 0; j < nComponents; j++) + if (probStat->getSystemMatrix(i, j)) + probStat->getSystemMatrix(i, j)->setRankDofs(isRankDof); + + TEST_EXIT_DBG(probStat->getRHS()->getDOFVector(i))("No rhs vector!\n"); + TEST_EXIT_DBG(probStat->getSolution()->getDOFVector(i))("No solution vector!\n"); + + probStat->getRHS()->getDOFVector(i)->setRankDofs(isRankDof); + probStat->getSolution()->getDOFVector(i)->setRankDofs(isRankDof); + } + + // === Remove periodic boundary conditions in sequential problem definition. === + + // Remove periodic boundaries in boundary manager on matrices and vectors. + for (int i = 0; i < nComponents; i++) { + for (int j = 0; j < nComponents; j++) { + DOFMatrix* mat = probStat->getSystemMatrix(i, j); + if (mat && mat->getBoundaryManager()) + removeBoundaryCondition(const_cast(mat->getBoundaryManager()->getBoundaryConditionMap())); + } + + if (probStat->getSolution()->getDOFVector(i)->getBoundaryManager()) + removeBoundaryCondition(const_cast(probStat->getSolution()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap())); + + if (probStat->getRHS()->getDOFVector(i)->getBoundaryManager()) + removeBoundaryCondition(const_cast(probStat->getRHS()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap())); + } + + // Remove periodic boundaries on elements in mesh. + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); + while (elInfo) { + elInfo->getElement()->deleteElementData(PERIODIC); + elInfo = stack.traverseNext(elInfo); + } } @@ -726,6 +796,18 @@ namespace AMDiS { } + void ParallelDomainBase::removeBoundaryCondition(BoundaryIndexMap& boundaryMap) + { + BoundaryIndexMap::iterator it = boundaryMap.begin(); + while (it != boundaryMap.end()) { + if (it->second->isPeriodic()) + boundaryMap.erase(it++); + else + ++it; + } + } + + void ParallelDomainBase::checkMeshChange() { FUNCNAME("ParallelDomainBase::checkMeshChange()"); @@ -929,7 +1011,7 @@ namespace AMDiS { void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data) - { + { int vecSize = data.size(); SerUtil::serialize(out, vecSize); for (int i = 0; i < vecSize; i++) { @@ -1058,6 +1140,15 @@ namespace AMDiS { } + void ParallelDomainBase::solveInitialProblem(AdaptInfo *adaptInfo) + { + if (timeIF) + timeIF->solveInitialProblem(adaptInfo); + + synchVector(*(probStat->getSolution())); + } + + void ParallelDomainBase::createInteriorBoundaryInfo() { FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()"); @@ -1346,7 +1437,8 @@ namespace AMDiS { // Stores all vertices at rank's interior boundaries. std::set rankBoundaryDofs; - // First, traverse the rank owned elements af the interior boundaries. + + // First, traverse the rank owned elements at the interior boundaries. for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) { Element *el = it.getBound().rankObj.el; basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); @@ -2262,6 +2354,29 @@ namespace AMDiS { } + void ParallelDomainBase::solve() + { + FUNCNAME("ParallelDomainBase::solve()"); + +#ifdef _OPENMP + double wtime = omp_get_wtime(); +#endif + clock_t first = clock(); + + fillPetscMatrix(probStat->getSystemMatrix(), probStat->getRHS()); + solvePetscMatrix(*(probStat->getSolution())); + +#ifdef _OPENMP + INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n", + TIME_USED(first, clock()), + omp_get_wtime() - wtime); +#else + INFO(info, 8)("solution of discrete system needed %.5f seconds\n", + TIME_USED(first, clock())); +#endif + } + + void ParallelDomainBase::serialize(std::ostream &out) { SerUtil::serialize(out, elemWeights); @@ -2286,6 +2401,8 @@ namespace AMDiS { SerUtil::serialize(out, rstart); SerUtil::serialize(out, nRankRows); SerUtil::serialize(out, nOverallRows); + + probStat->serialize(out); } @@ -2332,6 +2449,8 @@ namespace AMDiS { SerUtil::deserialize(in, nOverallRows); deserialized = true; + + probStat->deserialize(in); } void ParallelDomainBase::serialize(std::ostream &out, PeriodicDofMap &data) diff --git a/AMDiS/src/parallel/ParallelDomainBase.h b/AMDiS/src/parallel/ParallelDomainBase.h index 326843016a42fbb1c35f672ba4385d065b72553c..a3136f9fb4313977f30cf10caf36604239f71ed6 100644 --- a/AMDiS/src/parallel/ParallelDomainBase.h +++ b/AMDiS/src/parallel/ParallelDomainBase.h @@ -26,21 +26,21 @@ #include #include #include +#include +#include "Global.h" #include "ProblemTimeInterface.h" #include "ProblemIterationInterface.h" #include "FiniteElemSpace.h" #include "AdaptInfo.h" -#include "InteriorBoundary.h" +#include "parallel/InteriorBoundary.h" #include "Serializer.h" +#include "BoundaryManager.h" #include "AMDiS_fwd.h" #include "petsc.h" #include "petscsys.h" #include "petscao.h" -#include "mpi.h" - -#include "Global.h" namespace AMDiS { @@ -90,11 +90,9 @@ namespace AMDiS { typedef std::vector MeshCodeVec; public: - ParallelDomainBase(ProblemIterationInterface *iterationIF, - ProblemTimeInterface *timeIF, - FiniteElemSpace *feSpace, - RefinementManager *refineManager); - + ParallelDomainBase(ProblemVec *problemStat, + ProblemInstatVec *problemInstat); + virtual ~ParallelDomainBase() {} virtual void initParallelization(AdaptInfo *adaptInfo); @@ -143,11 +141,7 @@ namespace AMDiS { timeIF->closeTimestep(adaptInfo); } - virtual void solveInitialProblem(AdaptInfo *adaptInfo) - { - if (timeIF) - timeIF->solveInitialProblem(adaptInfo); - } + virtual void solveInitialProblem(AdaptInfo *adaptInfo); virtual void transferInitialSolution(AdaptInfo *adaptInfo) { @@ -169,7 +163,7 @@ namespace AMDiS { iterationIF->endIteration(adaptInfo); } - virtual void solve() {} + virtual void solve(); virtual int getNumProblems() { @@ -378,6 +372,9 @@ namespace AMDiS { */ void synchVector(SystemVector &vec); + // Removes all periodic boundaries from a given boundary map. + void removeBoundaryCondition(BoundaryIndexMap& boundaryMap); + bool fitElementToMeshCode(MeshStructure &code, Element *el, int ithSide, @@ -489,6 +486,9 @@ namespace AMDiS { /// ProblemTimeInterface *timeIF; + /// + ProblemVec *probStat; + /// The rank of the current process. int mpiRank; diff --git a/AMDiS/src/parallel/ParallelDomainVec.cc b/AMDiS/src/parallel/ParallelDomainVec.cc index 286cbf9e80da94154af1af696c585eaa4ffaf7eb..3ba7af68e180f5f13b4aee496d66800e783c3ba4 100644 --- a/AMDiS/src/parallel/ParallelDomainVec.cc +++ b/AMDiS/src/parallel/ParallelDomainVec.cc @@ -12,132 +12,9 @@ namespace AMDiS { ParallelDomainVec::ParallelDomainVec(ProblemVec *problem, ProblemInstatVec *problemInstat) - : ParallelDomainBase(problem, - problemInstat, - problem->getFESpace(0), - problem->getRefinementManager(0)), - probVec(problem) + : ParallelDomainBase(problem, problemInstat) { FUNCNAME("ParallelDomainVec::ParallelDomainVec()"); - info = probVec->getInfo(); - nComponents = probVec->getNumComponents(); - - // Test if all fe spaces are equal. Yet, different fe spaces are not supported for - // domain parallelization. - const FiniteElemSpace *fe = probVec->getFESpace(0); - for (int i = 0; i < nComponents; i++) - TEST_EXIT(fe == probVec->getFESpace(i)) - ("Parallelization does not supported different FE spaces!\n"); - - // Create parallel serialization file writer, if needed. - int writeSerialization = 0; - GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization); - if (writeSerialization) - problem->getFileWriterList().push_back(new Serializer(this)); - - int readSerialization = 0; - GET_PARAMETER(0, name + "->input->read serialization", "%d", &readSerialization); - if (readSerialization) { - std::string filename = ""; - GET_PARAMETER(0, name + "->input->serialization filename", &filename); - filename += ".p" + lexical_cast(mpiRank); - MSG("Start serialization with %s\n", filename.c_str()); - std::ifstream in(filename.c_str()); - deserialize(in); - in.close(); - } - } - - - void ParallelDomainVec::initParallelization(AdaptInfo *adaptInfo) - { - FUNCNAME("ParallelDomainVec::initParallelization()"); - - ParallelDomainBase::initParallelization(adaptInfo); - - for (int i = 0; i < nComponents; i++) { - for (int j = 0; j < nComponents; j++) - if (probVec->getSystemMatrix(i, j)) - probVec->getSystemMatrix(i, j)->setRankDofs(isRankDof); - - TEST_EXIT_DBG(probVec->getRHS()->getDOFVector(i))("No rhs vector!\n"); - TEST_EXIT_DBG(probVec->getSolution()->getDOFVector(i))("No solution vector!\n"); - - probVec->getRHS()->getDOFVector(i)->setRankDofs(isRankDof); - probVec->getSolution()->getDOFVector(i)->setRankDofs(isRankDof); - } - - // === Remove periodic boundary conditions in sequential problem definition. === - - // Remove periodic boundaries in boundary manager on matrices and vectors. - for (int i = 0; i < nComponents; i++) { - for (int j = 0; j < nComponents; j++) { - DOFMatrix* mat = probVec->getSystemMatrix(i, j); - if (mat && mat->getBoundaryManager()) - removeBoundaryCondition(const_cast(mat->getBoundaryManager()->getBoundaryConditionMap())); - } - - if (probVec->getSolution()->getDOFVector(i)->getBoundaryManager()) - removeBoundaryCondition(const_cast(probVec->getSolution()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap())); - - if (probVec->getRHS()->getDOFVector(i)->getBoundaryManager()) - removeBoundaryCondition(const_cast(probVec->getRHS()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap())); - } - - // Remove periodic boundaries on elements in mesh. - TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); - while (elInfo) { - elInfo->getElement()->deleteElementData(PERIODIC); - elInfo = stack.traverseNext(elInfo); - } } - - void ParallelDomainVec::exitParallelization(AdaptInfo *adaptInfo) - { - ParallelDomainBase::exitParallelization(adaptInfo); - } - - void ParallelDomainVec::solveInitialProblem(AdaptInfo *adaptInfo) - { - ParallelDomainBase::solveInitialProblem(adaptInfo); - - synchVector(*(probVec->getSolution())); - } - - void ParallelDomainVec::solve() - { - FUNCNAME("ParallelDomainVec::solve()"); - -#ifdef _OPENMP - double wtime = omp_get_wtime(); -#endif - clock_t first = clock(); - - fillPetscMatrix(probVec->getSystemMatrix(), probVec->getRHS()); - solvePetscMatrix(*(probVec->getSolution())); - -#ifdef _OPENMP - INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n", - TIME_USED(first, clock()), - omp_get_wtime() - wtime); -#else - INFO(info, 8)("solution of discrete system needed %.5f seconds\n", - TIME_USED(first, clock())); -#endif - } - - - void ParallelDomainVec::removeBoundaryCondition(BoundaryIndexMap& boundaryMap) - { - BoundaryIndexMap::iterator it = boundaryMap.begin(); - while (it != boundaryMap.end()) { - if (it->second->isPeriodic()) - boundaryMap.erase(it++); - else - ++it; - } - } - } diff --git a/AMDiS/src/parallel/ParallelDomainVec.h b/AMDiS/src/parallel/ParallelDomainVec.h index 79c5f7330e9386450f35c6a3e2ffed38a756aad1..c234c07990c1d36a788b57cbb714642b0b43a9e8 100644 --- a/AMDiS/src/parallel/ParallelDomainVec.h +++ b/AMDiS/src/parallel/ParallelDomainVec.h @@ -33,37 +33,6 @@ namespace AMDiS { public: ParallelDomainVec(ProblemVec *problem, ProblemInstatVec *problemInstat); - - void initParallelization(AdaptInfo *adaptInfo); - - void exitParallelization(AdaptInfo *adaptInfo); - - virtual void solveInitialProblem(AdaptInfo *adaptInfo); - - // Writes all data of this object to an output stream. - virtual void serialize(std::ostream &out) - { - probVec->serialize(out); - ParallelDomainBase::serialize(out); - } - - // Reads the object data from an input stream. - virtual void deserialize(std::istream &in) - { - probVec->deserialize(in); - ParallelDomainBase::deserialize(in); - } - - protected: - /// Starts the solution of the linear system using Petsc. - void solve(); - - // Removes all periodic boundaries from a given boundary map. - void removeBoundaryCondition(BoundaryIndexMap& boundaryMap); - - protected: - /// Pointer to the stationary problem. - ProblemVec *probVec; }; }