Bugfix in parallel code.

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;
   };
 

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())

AMDiS/src/parallel/ParallelDomainBase.cc

 #include <algorithm>
 #include <iostream>
 #include <fstream>
-
+#include <boost/lexical_cast.hpp>
 #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<ParallelDomainBase>(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<std::string>(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<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
+      }
+
+      if (probStat->getSolution()->getDOFVector(i)->getBoundaryManager())
+	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat->getSolution()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));
+
+      if (probStat->getRHS()->getDOFVector(i)->getBoundaryManager())
+	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(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()");
@@ -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<DegreeOfFreedom> 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)

AMDiS/src/parallel/ParallelDomainBase.h

@@ -26,21 +26,21 @@
 #include <map>
 #include <set>
 #include <vector>
+#include <mpi.h>
 
+#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,10 +90,8 @@ namespace AMDiS {
     typedef std::vector<MeshStructure> MeshCodeVec;
 
   public:
-    ParallelDomainBase(ProblemIterationInterface *iterationIF,
-		       ProblemTimeInterface *timeIF,
-		       FiniteElemSpace *feSpace,
-		       RefinementManager *refineManager);
+    ParallelDomainBase(ProblemVec *problemStat,
+		       ProblemInstatVec *problemInstat);
 		       
     virtual ~ParallelDomainBase() {}
 
@@ -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;
 

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<ParallelDomainVec>(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<std::string>(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<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
-      }
-
-      if (probVec->getSolution()->getDOFVector(i)->getBoundaryManager())
-	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probVec->getSolution()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));
-
-      if (probVec->getRHS()->getDOFVector(i)->getBoundaryManager())
-	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(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)