Fix a billion of very small bugs, most related to performance issues.

AMDiS/src/CreatorMap.cc

@@ -92,7 +92,7 @@ namespace AMDiS {
     #ifdef HAVE_BDDC_ML
       addCreator("bddcml", new BddcMlSolver::Creator);
     #endif
-    addCreator("petsc-stokes", new PetscSolverNavierStokes::Creator);
+    addCreator("petsc-navierstokes", new PetscSolverNavierStokes::Creator);
 #endif
   }
 

AMDiS/src/Element.cc

@@ -653,12 +653,8 @@ namespace AMDiS {
   {
     FUNCNAME("Element::getAllDofs()");
 
-    int sb = dofs.size();
-
     getNodeDofs(feSpace, bound, dofs, baseDofPtr);
 
-    //    MSG("    -> getNodeDofs %d\n", dofs.size() - sb);
-
     if (dofGeoIndex != NULL) {
       // In the case dofGeoIndex should be filled, set all node DOFs to be 
       // vertex DOFs.
@@ -667,14 +663,8 @@ namespace AMDiS {
 	(*dofGeoIndex)[i] = VERTEX;
     }
 
-    if (feSpace->getBasisFcts()->getDegree() > 1) {
-
-      sb = dofs.size();
-
+    if (feSpace->getBasisFcts()->getDegree() > 1)
       getHigherOrderDofs(feSpace, bound, dofs, baseDofPtr, dofGeoIndex);
-      
-      //      MSG("     -> getHODofs %d\n", dofs.size() - sb);
-    }
 
     if (dofGeoIndex) {
       TEST_EXIT_DBG(dofs.size() == dofGeoIndex->size())

AMDiS/src/io/ArhReader.cc

@@ -108,7 +108,7 @@ namespace AMDiS {
     string arhPrefix = "";
     map<int, int> elInRank;
     map<int, int> elCodeSize;
-    readMetaData(filename, elInRank, elCodeSize, arhPrefix);
+    int nProc = readMetaData(filename, elInRank, elCodeSize, arhPrefix);
     
 
     // === Check which arh files must be read by current rank. ===
@@ -133,8 +133,12 @@ namespace AMDiS {
 
     for (std::set<int>::iterator it = readArhFiles.begin();
 	 it != readArhFiles.end(); ++it) {
-      string arhFilename = 
-	directory.native() + "/" + arhPrefix + "-p" + boost::lexical_cast<string>(*it) + "-.arh";
+      string arhFilename = directory.native() + "/" + arhPrefix;
+      if (nProc == 1)
+	arhFilename += ".arh";
+      else
+	arhFilename += "-p" + boost::lexical_cast<string>(*it) + "-.arh";
+
       MSG("ARH file read from: %s\n", arhFilename.c_str());
       readFile(arhFilename, mesh, vecs);
     }

AMDiS/src/parallel/MeshDistributor.cc

@@ -81,9 +81,11 @@ namespace AMDiS {
       writeSerializationFile(false),
       repartitioningAllowed(false),
       repartitionIthChange(20),
+      repartitioningWaitAfterFail(20),
       nMeshChangesAfterLastRepartitioning(0),
       repartitioningCounter(0),
       repartitioningFailed(0),
+
       debugOutputDir(""),
       lastMeshChangeIndex(0),
       createBoundaryDofFlag(0),
@@ -101,6 +103,7 @@ namespace AMDiS {
     Parameters::get(name + "->repartitioning", repartitioningAllowed);
     Parameters::get(name + "->debug output dir", debugOutputDir);
     Parameters::get(name + "->repartition ith change", repartitionIthChange);
+    Parameters::get(name + "->repartition wait after fail", repartitioningWaitAfterFail);
     Parameters::get(name + "->mesh adaptivity", meshAdaptivity);
     
     string partStr = "parmetis";
@@ -366,8 +369,9 @@ namespace AMDiS {
       // and now partition the mesh    
       bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
       TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
-      partitioner->createPartitionMap(partitionMap);
     }
+
+    partitioner->createPartitionMap(partitionMap);
   }
 
 
@@ -1372,7 +1376,7 @@ namespace AMDiS {
       partitioner->partition(elemWeights, ADAPTIVE_REPART);
     if (!partitioningSucceed) {
       MPI::COMM_WORLD.Barrier();
-      repartitioningFailed = 20;
+      repartitioningFailed = repartitioningWaitAfterFail;;
       MSG("Mesh partitioner created empty partition!\n");
       MSG("Mesh repartitioning needed %.5f seconds\n", MPI::Wtime() - timePoint); 
       return;
@@ -1395,7 +1399,7 @@ namespace AMDiS {
       
       if (!partitioningSucceed || !partitioner->meshChanged()) {
 	MPI::COMM_WORLD.Barrier();
-	repartitioningFailed = 20;
+	repartitioningFailed = repartitioningWaitAfterFail;;
 	MSG("Mesh repartitioning needed %.5f seconds\n", MPI::Wtime() - timePoint);   
 	return;
       }
@@ -1562,7 +1566,9 @@ namespace AMDiS {
     // update the number of elements, vertices, etc. of the mesh.
     mesh->removeMacroElements(deleteMacroElements, feSpaces);
 
+
     // === Remove double DOFs. ===
+
     MeshManipulation meshManipulation(mesh);
     meshManipulation.deleteDoubleDofs(feSpaces, newMacroEl, elObjDb);
 

AMDiS/src/parallel/MeshDistributor.h

@@ -536,6 +536,9 @@ namespace AMDiS {
     /// repartitionings.
     int repartitionIthChange;
 
+    ///
+    int repartitioningWaitAfterFail;
+
     /// Counts the number of mesh changes after the last mesh repartitioning 
     /// was done.
     int nMeshChangesAfterLastRepartitioning;

AMDiS/src/parallel/MeshManipulation.cc

@@ -58,7 +58,7 @@ namespace AMDiS {
     for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
 	 it != newMacroEl.end(); ++it)
       macroIndexMap[(*it)->getIndex()] = *it;
-
+   
     // === Traverse mesh and put all "old" macro element to macrosProcessed  ===
     // === that stores all macro elements which are really connected to the  ===
     // === overall mesh structure.                                           ===
@@ -83,7 +83,6 @@ namespace AMDiS {
     feSpace->getMesh()->getDofIndexCoords(coords);
 #endif
 
-
     // === Traverse all new elements and connect them to the overall mesh     ===
     // === structure by deleting all double DOFs on macro element's vertices, ===
     // === edges and faces.                                                   ===
@@ -105,8 +104,8 @@ namespace AMDiS {
 	  if (elIt->elIndex == (*it)->getIndex())
 	    continue;
 
-	  if (macrosProcessed.count(elIt->elIndex) == 1) {
-	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex) == 1)
+	  if (macrosProcessed.count(elIt->elIndex)) {
+	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))
 	      ("Should not happen!\n");
 
 	    Element *el0 = (*it)->getElement();
@@ -129,12 +128,12 @@ namespace AMDiS {
 	vector<ElementObjectData> &edgeEl = 
 	  elObjDb.getElementsEdge((*it)->getIndex(), i);
 
-	for (vector<ElementObjectData>::iterator elIt = edgeEl.begin();
+      	for (vector<ElementObjectData>::iterator elIt = edgeEl.begin();
 	     elIt != edgeEl.end(); ++elIt) {
 	  if (elIt->elIndex == (*it)->getIndex())
 	    continue;
 
-	  if (macrosProcessed.count(elIt->elIndex) == 1) {
+	  if (macrosProcessed.count(elIt->elIndex)) {
 	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))
 	      ("Should not happen!\n");
 	    
@@ -150,9 +149,9 @@ namespace AMDiS {
 	    vector<GeoIndex> dofGeoIndex0, dofGeoIndex1;
 	    el0->getAllDofs(feSpace, b0, dofs0, true, &dofGeoIndex0);
 	    el1->getAllDofs(feSpace, b1, dofs1, true, &dofGeoIndex1);
-
+	    
 #if (DEBUG != 0)
-	    if (feSpaces.size() == 1)
+	    if (feSpaces.size())
 	      debug::testDofsByCoords(coords, dofs0, dofs1);
 	    else 
 	      TEST_EXIT_DBG(dofs0.size() == dofs1.size())
@@ -184,7 +183,7 @@ namespace AMDiS {
 	  if (elIt->elIndex == (*it)->getIndex())
 	    continue;
 
-	  if (macrosProcessed.count(elIt->elIndex) == 1) {
+	  if (macrosProcessed.count(elIt->elIndex)) {
 	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))
 	      ("Should not happen!\n");
 
@@ -202,7 +201,7 @@ namespace AMDiS {
 	    el1->getAllDofs(feSpace, b1, dofs1, true, &dofGeoIndex1);
 
 #if (DEBUG != 0)
-	    if (feSpaces.size() == 1)
+	    if (feSpaces.size())
 	      debug::testDofsByCoords(coords, dofs0, dofs1);
 	    else 
 	      TEST_EXIT_DBG(dofs0.size() == dofs1.size())
@@ -355,7 +354,8 @@ namespace AMDiS {
     // === To the mesh structure code.                                          ===
 
     TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
+    ElInfo *elInfo = 
+      stack.traverseFirstOneMacro(mesh, boundEl.elIndex, -1, traverseFlag);
 
     if (s0 != -1) {
       while (elInfo && elInfo->getElement() != child0)

AMDiS/src/parallel/ParallelDebug.cc

@@ -833,8 +833,7 @@ namespace AMDiS {
 
     map<int, double> vec;    
     TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, 
-					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
+    ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, Mesh::CALL_LEAF_EL);
     
     while (elInfo) {		  
       int index = elInfo->getElement()->getIndex();

AMDiS/src/parallel/PetscHelper.cc

@@ -250,7 +250,7 @@ namespace AMDiS {
 
 
     void setSolver(KSP ksp, 
-		   string kspPrefix,
+		   const char* kspPrefix,
 		   KSPType kspType, 
 		   PCType pcType, 
 		   PetscReal rtol,
@@ -259,8 +259,9 @@ namespace AMDiS {
     {
       KSPSetType(ksp, kspType);
       KSPSetTolerances(ksp, rtol, atol, PETSC_DEFAULT, maxIt);
-      if (kspPrefix != "")
-	KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
+      //if (kspPrefix != "")
+      if (strcmp(kspPrefix, "") != 0)
+	KSPSetOptionsPrefix(ksp, kspPrefix);
       KSPSetFromOptions(ksp);
 
       PC pc;

AMDiS/src/parallel/PetscHelper.h

@@ -92,7 +92,7 @@ namespace AMDiS {
     void matNestConvert(Mat matNest, Mat &mat);
 
     void setSolver(KSP ksp, 
-		   string kspPrefix,
+		   const char* kspPrefix,
 		   KSPType kspType, 
 		   PCType pcType, 
 		   PetscReal rtol = PETSC_DEFAULT,

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -501,20 +501,20 @@ namespace AMDiS {
 	   isNames[i].c_str());
       }
 
-      createFieldSplit(pc, isNames[i], blockComponents);
+      createFieldSplit(pc, isNames[i].c_str(), blockComponents);
     }
   }
 
 
   void PetscSolverGlobalMatrix::createFieldSplit(PC pc, 
-						 string splitName, 
+						 const char* splitName, 
 						 vector<int> &components)
   {
     FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()");
 
     IS is;
     interiorMap->createIndexSet(is, components[0], components.size());
-    PCFieldSplitSetIS(pc, splitName.c_str(), is);
+    PCFieldSplitSetIS(pc, splitName, is);
     ISDestroy(&is);
   }
 
@@ -823,7 +823,7 @@ namespace AMDiS {
     fe.push_back(componentSpaces[component]);
 
     PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
-    subSolver->setKspPrefix(kspPrefix.c_str());
+    subSolver->setKspPrefix(kspPrefix);
     subSolver->setMeshDistributor(meshDistributor,
 				  mpiCommGlobal,
 				  mpiCommLocal);

AMDiS/src/parallel/PetscSolverGlobalMatrix.h

@@ -88,10 +88,10 @@ namespace AMDiS {
      * \param[in] components  System component numbers of the field split. At
      *                        the moment only continuous splits are allowed.
      */
-    void createFieldSplit(PC pc, string splitName, vector<int> &components);
+    void createFieldSplit(PC pc, const char* splitName, vector<int> &components);
 
     /// Wrapper to create field split from only one component.
-    void createFieldSplit(PC pc, string splitName, int component)
+    void createFieldSplit(PC pc, const char* splitName, int component)
     {
       vector<int> components;
       components.push_back(component);

AMDiS/src/parallel/PetscSolverNavierStokes.cc

@@ -46,6 +46,7 @@ namespace AMDiS {
     : PetscSolverGlobalMatrix(name),
       pressureComponent(-1),
       pressureNullSpace(true),
+      useOldInitialGuess(false),
       massMatrixSolver(NULL),
       laplaceMatrixSolver(NULL),
       nu(NULL),
@@ -61,6 +62,28 @@ namespace AMDiS {
     Parameters::get(initFileStr + "->navierstokes->pressure null space",
 		    pressureNullSpace);
     TEST_EXIT(pressureNullSpace)("This is not yet tested, may be wrong!\n");
+
+    Parameters::get(initFileStr + "->navierstokes->use old initial guess", 
+		    useOldInitialGuess);
+  }
+
+
+  void PetscSolverNavierStokes::solvePetscMatrix(SystemVector &vec, 
+						 AdaptInfo *adaptInfo)
+  {
+    FUNCNAME("PetscSolverNavierStokes::solvePetscMatrix()");
+
+    if (useOldInitialGuess) {
+      VecSet(getVecSolInterior(), 0.0);
+      
+      for (int i = 0; i < solution->getSize(); i++)
+	setDofVector(getVecSolInterior(), solution->getDOFVector(i), i, true);
+      
+      vecSolAssembly();
+      KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
+    }
+
+    PetscSolverGlobalMatrix::solvePetscMatrix(vec, adaptInfo);
   }
 
 
@@ -151,7 +174,7 @@ namespace AMDiS {
     else
       laplaceOp.addTerm(new VecAtQP_SOT(phase, &idFct));
     laplaceMatrix.assembleOperator(laplaceOp);
-    laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_");
+    laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_"));
     laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
 
 

AMDiS/src/parallel/PetscSolverNavierStokes.h

@@ -120,8 +120,11 @@ namespace AMDiS {
 	return new PetscSolverNavierStokes(this->name); 
       }
     };
+
     PetscSolverNavierStokes(string name);
 
+    void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo);
+
     void setStokesData(double *nuPtr, double *invTauPtr, SystemVector *vec)
     {
       nu = nuPtr;
@@ -144,6 +147,9 @@ namespace AMDiS {
 
     bool pressureNullSpace;
 
+    /// If true, old solution is used for initial guess in solver phase.
+    bool useOldInitialGuess;
+
     PetscSolver *massMatrixSolver, *laplaceMatrixSolver, *conDifMatrixSolver;
 
     NavierStokesSchurData matShellContext;