From a1f6ce9c68e4d355a432b584dd7358566004b706 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 13 Feb 2012 09:37:00 +0000
Subject: [PATCH] Added some nice to have features for using PETSc solver.

---
 AMDiS/src/parallel/DofComm.cc                 | 16 ++++++----
 AMDiS/src/parallel/DofComm.h                  | 29 ++++++++++---------
 AMDiS/src/parallel/PetscProblemStat.cc        | 20 ++++++++-----
 AMDiS/src/parallel/PetscSolver.cc             | 12 ++++++++
 AMDiS/src/parallel/PetscSolver.h              | 21 +++++++++-----
 AMDiS/src/parallel/PetscSolverGlobalMatrix.cc |  8 ++++-
 6 files changed, 70 insertions(+), 36 deletions(-)

diff --git a/AMDiS/src/parallel/DofComm.cc b/AMDiS/src/parallel/DofComm.cc
index dd812a66..61f8dfa9 100644
--- a/AMDiS/src/parallel/DofComm.cc
+++ b/AMDiS/src/parallel/DofComm.cc
@@ -38,14 +38,18 @@ namespace AMDiS {
   }
 
   
-  void DofComm::Iterator::setNextFeMap()
+  bool DofComm::Iterator::setNextFeMap()
   {
+    FUNCNAME("DofComm::Iterator::setNextFeMap()");
+
     if (dataIter != dofComm.data.end()) {
+      TEST_EXIT_DBG(dataIter->second.size())("Should not happen!\n");
+
       feMapIter = dataIter->second.begin();
       
       if (traverseFeSpace != NULL) {
-	TEST_EXIT_DBG(dataIter->second.count(traverseFeSpace))
-	  ("Should not happen!\n");
+	if ((dataIter->second.count(traverseFeSpace) == 0))
+	  return false;
 	
 	while (feMapIter->first != traverseFeSpace &&
 	       feMapIter != dataIter->second.end())
@@ -56,11 +60,11 @@ namespace AMDiS {
 	  ("Should not happen!\n");
       }
       
-      if (feMapIter != dataIter->second.end())
-	dofIter = feMapIter->second.begin();
-      
+      dofIter = feMapIter->second.begin();      
       dofCounter = 0;
     }
+
+    return true;
   }
 
 }
diff --git a/AMDiS/src/parallel/DofComm.h b/AMDiS/src/parallel/DofComm.h
index 5ef1cc87..58d5dd25 100644
--- a/AMDiS/src/parallel/DofComm.h
+++ b/AMDiS/src/parallel/DofComm.h
@@ -77,7 +77,8 @@ namespace AMDiS {
 
 	dataIter = dofComm.data.begin();
 
-	setNextFeMap();
+	while (setNextFeMap() == false)
+	  ++dataIter;
       }
 
       inline bool end()
@@ -87,9 +88,9 @@ namespace AMDiS {
       
       inline void nextRank()
       {
-	++dataIter;
-
-	setNextFeMap();
+	do {
+	  ++dataIter;
+	} while (setNextFeMap() == false);
       }
 
       inline void nextFeSpace()
@@ -101,9 +102,9 @@ namespace AMDiS {
       {
 	++feMapIter;
 	if (feMapIter == dataIter->second.end()) {
-	  ++dataIter;
-
-	  setNextFeMap();
+	  do {
+	    ++dataIter;
+	  } while (setNextFeMap() == false);
 	} else {
 	  dofIter = feMapIter->second.begin();
 	  dofCounter = 0;
@@ -122,15 +123,17 @@ namespace AMDiS {
 	    ++feMapIter;
 	}
 
-	TEST_EXIT_DBG(feMapIter != dataIter->second.end())
-	  ("Should not happen!\n");
-
-	dofIter = feMapIter->second.begin();
-	dofCounter = 0;
+	if (feMapIter != dataIter->second.end()) {
+	  dofIter = feMapIter->second.begin();
+	  dofCounter = 0;
+	}
       }
 
       inline bool endDofIter()
       {
+	if (feMapIter == dataIter->second.end())
+	  return true;
+
 	return (dofIter == feMapIter->second.end());
       }
       
@@ -171,7 +174,7 @@ namespace AMDiS {
       }
 
     protected:
-      void setNextFeMap();
+      bool setNextFeMap();
 
     protected:
       DofComm &dofComm;
diff --git a/AMDiS/src/parallel/PetscProblemStat.cc b/AMDiS/src/parallel/PetscProblemStat.cc
index 16f4df52..a2afb7d5 100644
--- a/AMDiS/src/parallel/PetscProblemStat.cc
+++ b/AMDiS/src/parallel/PetscProblemStat.cc
@@ -31,26 +31,30 @@ namespace AMDiS {
   {
     FUNCNAME("PetscProblemStat::PetscProblemStat()");
 
-    string name("");
-    Parameters::get("parallel->solver", name);
+    string tmp("");
+    Parameters::get("parallel->solver", tmp);
     
-    if (name == "petsc-schur") {
+    if (tmp == "petsc-schur") {
       petscSolver = new PetscSolverSchur();
-    } else if (name == "petsc-feti") {
+    } else if (tmp == "petsc-feti") {
       petscSolver = new PetscSolverFeti();
-    } else if (name == "petsc-block") {
+    } else if (tmp == "petsc-block") {
       petscSolver = new PetscSolverGlobalBlockMatrix();
-    } else if (name == "petsc" || name == "") {
+    } else if (tmp == "petsc" || tmp == "") {
       petscSolver = new PetscSolverGlobalMatrix();
-    } else if (name == "bddcml") {
+    } else if (tmp == "bddcml") {
 #ifdef HAVE_BDDC_ML
       petscSolver = new BddcMlSolver();
 #else
       ERROR_EXIT("AMDiS was compiled without BDDC-ML support!\n");
 #endif
     } else {
-      ERROR_EXIT("No parallel solver %s available!\n", name.c_str());
+      ERROR_EXIT("No parallel solver %s available!\n", tmp.c_str());
     }
+
+    tmp = "";
+    Parameters::get(nameStr + "->solver->petsc prefix", tmp);
+    petscSolver->setKspPrefix(tmp);    
   }
 
 
diff --git a/AMDiS/src/parallel/PetscSolver.cc b/AMDiS/src/parallel/PetscSolver.cc
index 6c5b2109..ad2d7cc9 100644
--- a/AMDiS/src/parallel/PetscSolver.cc
+++ b/AMDiS/src/parallel/PetscSolver.cc
@@ -18,6 +18,18 @@ namespace AMDiS {
 
   using namespace std;
 
+  PetscSolver::PetscSolver()
+    : meshDistributor(NULL),
+      mpiRank(-1),
+      kspPrefix("")
+  {
+    string kspStr = "";
+    Parameters::get("parallel->solver->petsc->ksp", kspStr);
+    if (kspStr != "")
+      PetscOptionsInsertString(kspStr.c_str());
+  }
+
+
   void PetscSolver::printSolutionInfo(AdaptInfo *adaptInfo,
 				      bool iterationCounter, 
 				      bool residual)
diff --git a/AMDiS/src/parallel/PetscSolver.h b/AMDiS/src/parallel/PetscSolver.h
index cb02f068..3c2e3e1d 100644
--- a/AMDiS/src/parallel/PetscSolver.h
+++ b/AMDiS/src/parallel/PetscSolver.h
@@ -46,10 +46,7 @@ namespace AMDiS {
   class PetscSolver
   {
   public:
-    PetscSolver()
-      : meshDistributor(NULL),
-	mpiRank(-1)
-    {}
+    PetscSolver();
 
     virtual ~PetscSolver() {}
 
@@ -96,6 +93,11 @@ namespace AMDiS {
       return pc; 
     }
 
+    void setKspPrefix(std::string s)
+    {
+      kspPrefix = s;
+    }
+
   protected:
     void printSolutionInfo(AdaptInfo* adaptInfo,
 			   bool iterationCounter = true,
@@ -135,15 +137,18 @@ namespace AMDiS {
     /// Petsc's matrix structure.
     Mat petscMatrix;
 
-    /** \brief
-     * PETSc's vector structures for the rhs vector, the solution vector and a
-     * temporary vector for calculating the final residuum.
-     */
+    /// PETSc's vector structures for the rhs vector, the solution vector and a
+    /// temporary vector for calculating the final residuum.
     Vec petscRhsVec, petscSolVec, petscTmpVec;
 
+    /// PETSc solver object
     KSP solver;
 
+    /// PETSc preconditioner object
     PC pc;
+
+    /// KSP database prefix
+    string kspPrefix;
   };
 
 
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 8470ad97..cd21ab43 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -10,6 +10,7 @@
 // See also license.opensource.txt in the distribution.
 
 
+#include "AMDiS.h"
 #include "parallel/PetscSolverGlobalMatrix.h"
 #include "parallel/StdMpi.h"
 #include "parallel/MpiHelper.h"
@@ -113,6 +114,7 @@ namespace AMDiS {
     KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); 
     KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
     KSPSetType(solver, KSPBCGS);
+    KSPSetOptionsPrefix(solver, kspPrefix.c_str());
     KSPSetFromOptions(solver);
     PCSetFromOptions(pc);
 
@@ -171,7 +173,11 @@ namespace AMDiS {
     }
 
     // PETSc.
-    KSPSolve(solver, petscRhsVec, petscSolVec);
+    PetscErrorCode solverError = KSPSolve(solver, petscRhsVec, petscSolVec);
+    if (solverError != 0) {
+      AMDiS::finalize();
+      exit(-1);
+    }
 
     // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
     PetscScalar *vecPointer;
-- 
GitLab