pfc preconditioner improved

extensions/ExtendedProblemStat.h

@@ -66,6 +66,7 @@ public:
   }
   
   
+  //////////////////////////////////////////////////////////////////////////////
   ~ExtendedProblemStat()
   {
     for (int i = 0; i < nComponents; ++i) {
@@ -81,6 +82,7 @@ public:
   }
 
   
+  //////////////////////////////////////////////////////////////////////////////
   void initialize(Flag initFlag,
 		  ProblemStatSeq *adoptProblem = NULL,
 		  Flag adoptFlag = INIT_NOTHING)
@@ -97,7 +99,7 @@ public:
   }
 
   
-  ///
+  //////////////////////////////////////////////////////////////////////////////
   void setExactSolution(DOFVector<double> *dof, int component) 
   {
     TEST_EXIT(exactSolutions[component] != NULL)
@@ -106,7 +108,7 @@ public:
   }
 
   
-  ///
+  //////////////////////////////////////////////////////////////////////////////
   inline DOFVector<double> *getExactSolution(int i = 0) 
   {
     TEST_EXIT(exactSolutions[i] != NULL)
@@ -115,6 +117,7 @@ public:
   }
 
   
+  //////////////////////////////////////////////////////////////////////////////
   void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
 			  bool asmMatrix, bool asmVector)
   {  
@@ -156,6 +159,7 @@ public:
   }
 
   
+  //////////////////////////////////////////////////////////////////////////////
   void solve(AdaptInfo *adaptInfo,
 	      bool createMatrixData = true,
 	      bool storeMatrixData = false)
@@ -167,6 +171,7 @@ public:
 
   
   /// Add arbitrary boundary condition to system
+  //////////////////////////////////////////////////////////////////////////////
   void addBoundaryCondition(BoundaryCondition *bc, int row, int col)
   {
     boundaryConditionSet = true;
@@ -181,6 +186,7 @@ public:
   //============================================================================
   
   
+  //////////////////////////////////////////////////////////////////////////////
   template<typename ValueContainer>
   void addDirichletBC(BoundaryType type, int row, int col,
 		      ValueContainer *values)
@@ -300,6 +306,7 @@ public:
   // ===========================================================================
   
   
+  //////////////////////////////////////////////////////////////////////////////
   void addManualPeriodicBC(int row,
 			    BoundaryType nr, AbstractFunction<bool, WorldVector<double> >* meshIndicator,
 			    AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap)
@@ -308,6 +315,7 @@ public:
   }
 
   
+  //////////////////////////////////////////////////////////////////////////////
   void addManualPeriodicBC(int row,
 			    AbstractFunction<bool, WorldVector<double> >* meshIndicator,
 			    AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap)

extensions/preconditioner/BlockPreconditioner.h

@@ -62,9 +62,27 @@ namespace AMDiS {
     template<typename SolverType, typename RunnerType>
     void createSubSolver(std::string param, SolverType*& solver, RunnerType*& runner, std::string solverType = "0")
     {
-      // === create solver ===
-      string initFileStr = param + "->solver";
+      // definition of standard-backends
+#if defined HAVE_PARALLEL_PETSC
+      string backend("p_petsc");
+#elif defined HAVE_PARALLEL_MTL
+      string backend("p_mtl");
+#elif defined HAVE_PETSC
+      string backend("petsc");
+#else
+      string backend("mtl");
+#endif
+    
+      // === read backend-name ===
+      string initFileStr = param + "->solver";    
+      Parameters::get(initFileStr + "->backend", backend);
+
+      // === read solver-name ===
       Parameters::get(initFileStr, solverType);
+      
+      if (backend != "0" && backend != "no" && backend != "")
+	solverType = backend + "_" + solverType;
+    
       LinearSolverCreator *solverCreator = 
 	dynamic_cast<LinearSolverCreator*>(CreatorMap<LinearSolver>::getCreator(solverType, initFileStr));
       TEST_EXIT(solverCreator)
@@ -72,7 +90,6 @@ namespace AMDiS {
       solverCreator->setName(initFileStr);
       solver = dynamic_cast<SolverType*>(solverCreator->create());
       assert(solver != NULL);
-      solver->initParameters();
       
       runner = dynamic_cast<RunnerType*>(solver->getRunner());
       assert(runner != NULL);

extensions/preconditioner/MTLPreconPfc.h

@@ -73,7 +73,7 @@ struct MTLPreconPfc : BlockPreconditioner
     tau(NULL),
     solverM(NULL), solverM2(NULL), solverMpL(NULL), runnerM(NULL), runnerM2(NULL), runnerMpL(NULL),
     Pid(NULL), Pdiag(NULL),
-    maxIterS(30), restartS(30)
+    maxIterS(30), restartS(100)
   {
     createSubSolver(name + "->subsolver M", solverM, runnerM, "cg");
     createSubSolver(name + "->subsolver M2", solverM2, runnerM2, "cg");
@@ -120,7 +120,6 @@ struct MTLPreconPfc : BlockPreconditioner
   void solveM2(VectorX& x, const VectorB& b) const
   {
     runnerM2->solve(getM(), x, b);
-//     umfpack_solve(getM(), x, b);
   }
   
   template <typename VectorX, typename VectorB>
@@ -134,7 +133,7 @@ struct MTLPreconPfc : BlockPreconditioner
   { 
     itl::basic_iteration<double> iter(b, maxIterS, 0, 0);
     x = 0.0;
-    itl::fgmres(S, x, b, *Pid, *Pdiag, iter, 30);	// S*x = b
+    itl::fgmres_householder(S, x, b, *Pdiag, iter, restartS);	// S*x = b
   }
   
   const MTLMatrix& getM() const { return A->getSubMatrix(2,2); }

extensions/preconditioner/MTLPreconPfc.hh

@@ -16,8 +16,7 @@
  ******************************************************************************/
 
 void MTLPreconPfc::init(const BlockMatrix& A_, const MTLMatrix& fullMatrix_)
-{ FUNCNAME("MTLPreconPfc::init()");
-
+{
   assert(tau != NULL);
   super::init(A_, fullMatrix_);
   
@@ -33,30 +32,13 @@ void MTLPreconPfc::init(const BlockMatrix& A_, const MTLMatrix& fullMatrix_)
   y1.change_dim(num_rows(getM()));
   tmp.change_dim(num_rows(getM()));  
   
+  // init sub-solvers
   Pid = new itl::pc::identity<MTLTypes::MTLMatrix, double>(getM());
-  Pdiag = new itl::pc::diagonal<MTLTypes::MTLMatrix, double>(getM());
+  Pdiag = new itl::pc::diagonal<MTLTypes::MTLMatrix, double>(MpL);
   
   runnerM->init(A_, getM());
   runnerM2->init(A_, getM());
   runnerMpL->init(A_, MpL);
-  
-#ifndef NDEBUG
-  mtl::io::matrix_market_ostream out("A.mtx");
-  out << fullMatrix_;
-  out.close();
-  
-  mtl::io::matrix_market_ostream out2("M.mtx");
-  out2 << getM();
-  out2.close();
-  
-  mtl::io::matrix_market_ostream out3("L.mtx");
-  out3 << getL();
-  out3.close();
-  
-  mtl::io::matrix_market_ostream out4("MpL.mtx");
-  out4 << MpL;
-  out4.close();
-#endif
 }
 
 
@@ -106,8 +88,6 @@ void PfcSchurMatrix<MTLPreconPfc>::mult(const VectorIn& b, VectorOut& x, Assign)
   y1 = super->getL() * b;
   super->solveM(y2, y1);
   y3 = super->getL() * y2;
-//   super->solveM(y2, y3);
-//   y3 = super->getL() * y2;
 
   y2 = super->getM() * b;
   y2-= (2.0*delta) * y1;

extensions/preconditioner/PetscPreconPfc.cc

@@ -33,8 +33,8 @@ namespace AMDiS {
     MatMult(data->matK, b, y1);		// y1 := K*b
     KSPSolve(data->kspM, y1, y2);	// M*y2 = y1
     MatMult(data->matK, y2, y3);	// y3 := K*y2
-    KSPSolve(data->kspM, y3, y2);	// M*y2 = y3
-    MatMult(data->matK, y2, y3);	// y3 := K*y2
+//     KSPSolve(data->kspM, y3, y2);	// M*y2 = y3
+//     MatMult(data->matK, y2, y3);	// y3 := K*y2
     
     VecScale(y3, data->delta);		// y3 = delta*y3
     MatMultAdd(data->matM, b, y3, x);	// x = M*b + y3
@@ -83,7 +83,7 @@ namespace AMDiS {
 
     MatMult(data->matK, x2, tmp); 			// tmp := K*x2
     VecAYPX(tmp, -1.0, b3);				// tmp := b3 - tmp
-    KSPSolve(data->kspM, tmp, x3); 			// M*x3 = tmp    
+    KSPSolve(data->kspM2, tmp, x3); 			// M*x3 = tmp    
 
     VecDestroy(&y1);
     VecDestroy(&y2);
@@ -93,11 +93,11 @@ namespace AMDiS {
 
   void PetscPreconPfc::init(const BlockMatrix& A, const Mat& matrix)
   {
-    FUNCNAME("PetscPreconPfc::initSolver()");
-
+    assert(tau != NULL);
     dynamic_cast<PetscSolver<PetscPreconPfc>*>(&oem)->setNestedVectors(true);
     PetscRunner::init(A, matrix);
     
+    // init shell preconditioner
     PCSetType(getPc(), PCSHELL);
     PCShellSetApply(getPc(), pcPfcShell);
     PCShellSetContext(getPc(), &data);
@@ -108,16 +108,7 @@ namespace AMDiS {
     MatDuplicate(data.matM, MAT_COPY_VALUES, &MpK);
     MatAXPY(MpK, sqrt(*tau), data.matK, SAME_NONZERO_PATTERN);
     
-    createSubSolver(data.kspM, data.matM, "mass_");
-    createSubSolver(data.kspMpK, MpK, "MpK_");
-    
-    // === Setup preconditioner data ===    
-    data.delta = sqrt(*tau);
-    data.tau = (*tau);
-
-    setSolver(data.kspM, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 5);
-    setSolver(data.kspMpK, "MpK_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 5);
-    
+    // schur-komplete-Matrix
     PetscInt numRows, numCols;
     MatGetSize(data.matM, &numRows, &numCols);
     MatCreateShell(PETSC_COMM_SELF, numRows, numCols, 
@@ -129,6 +120,19 @@ namespace AMDiS {
     KSPSetOperators(data.kspSchur, matSchur, MpK, SAME_NONZERO_PATTERN);
     setSolver(data.kspSchur, "schur_", KSPFGMRES, PCJACOBI, 1e-6, 1e-8, 4);
     KSPSetFromOptions(data.kspSchur);
+    
+    // init sub-solvers
+    createSubSolver(data.kspM, data.matM, "mass_");
+    createSubSolver(data.kspM2, data.matM, "mass2_");
+    createSubSolver(data.kspMpK, MpK, "MpK_");
+    
+    data.delta = sqrt(*tau);
+    data.tau = (*tau);
+
+    setSolver(data.kspM, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 5);
+    setSolver(data.kspM2, "mass2_", KSPCG, PCJACOBI, 0.0, 1e-14, 15);
+    setSolver(data.kspMpK, "MpK_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 5);
+    
   }
 
 
@@ -140,6 +144,7 @@ namespace AMDiS {
     MatDestroy(&matSchur);
     
     KSPDestroy(&data.kspM);
+    KSPDestroy(&data.kspM2);
     KSPDestroy(&data.kspMpK);
     KSPDestroy(&data.kspSchur);
     

extensions/preconditioner/PetscPreconPfc.h

@@ -25,7 +25,7 @@ namespace AMDiS {
   using namespace std;
 
   struct PfcData {
-    KSP kspM, kspMpK, kspSchur;
+    KSP kspM, kspM2, kspMpK, kspSchur;
     Mat matM, matK;
     double delta, tau;
   };

extensions/preconditioner/PhaseFieldCrystal_.cc

@@ -51,8 +51,7 @@ void PhaseFieldCrystal_::fillOperators()
   opMTemp->addZeroOrderTerm(new Simple_ZOT());
   // -2*rho_old^3
   Operator *opMPowExpl = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
-  //opMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1),new Pow3Functor(-2.0)));
-  opMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1),new AMDiS::Pow<3>(1.0, 3*degree)));
+  opMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1),new AMDiS::Pow<3>(-2.0, 3*degree)));
   // -3*rho_old^2
   Operator *opMPowImpl = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
   opMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1),new AMDiS::Pow<2>(-3.0, 2*degree)));
@@ -68,11 +67,10 @@ void PhaseFieldCrystal_::fillOperators()
   // mu-2*nu-laplace(nu)-(1+r)*rho = (rho_old^3) + ExtPot - eps^2/(rho_old+0.9)
   // ----------------------------------------------------------------------
   prob->addMatrixOperator(opMTemp, 0, 1, getTempParameter(), getTempParameter()); // -phi*(1+r)*rho
-  //prob->addMatrixOperator(opMPowImpl, 0, 1); // -3*rho*rho_old^2
+  prob->addMatrixOperator(opMPowImpl, 0, 1); // -3*rho*rho_old^2
   prob->addMatrixOperator(opL, 0, 1, &two, &two); // -2*phi*laplace(rho) * psi
   prob->addMatrixOperator(opMnew, 0, 0); // phi*mu * psi
   prob->addMatrixOperator(opL, 0, 2); // phi*grad(nu) * grad(psi)
-  // prob.addMatrixOperator(opMnew,2,2,&minus2);
   // . . . vectorOperators . . . . . . . . . . . . . . .
   prob->addVectorOperator(opMPowExpl, 0); // -2*phi^old*rho_old^3
   

extensions/preconditioner/PhaseFieldCrystal_.h

@@ -18,7 +18,7 @@ public: // typedefs
   typedef BaseProblem<ExtendedProblemStat> super;
 
 public:
-
+ 
   PhaseFieldCrystal_(const std::string &name_, bool createProblem = true);
   ~PhaseFieldCrystal_() {}