CahnHilliard preconditioner updated

AMDiS/src/parallel/MeshDistributor.cc

@@ -61,6 +61,7 @@
 #include "ProblemInstat.h"
 #include "RefinementManager3d.h"
 #include "Debug.h"
+#include "Timer.h"
 
 namespace AMDiS { namespace Parallel {
 
@@ -719,9 +720,7 @@ namespace AMDiS { namespace Parallel {
     if (mesh->getDim() != 3)
       return;
 
-    MPI::COMM_WORLD.Barrier();
-    double first = MPI::Wtime();
-
+    Timer t;
 
     // === Create set of all edges and all macro element indices in rank's ===
     // === subdomain.                                                      ===
@@ -829,8 +828,7 @@ namespace AMDiS { namespace Parallel {
     }
 
 
-    MPI::COMM_WORLD.Barrier();
-    MSG("Fix 3D mesh refinement needed %.5f seconds\n", MPI::Wtime() - first);
+    MSG("Fix 3D mesh refinement needed %.5f seconds\n", t.elapsed());
   }
 
 
@@ -1070,7 +1068,7 @@ namespace AMDiS { namespace Parallel {
 
     MPI::Intracomm &mpiComm = MPI::COMM_WORLD;
     mpiComm.Barrier();
-    double first = MPI::Wtime();
+    Timer t;
 
     int skip = 0;
     Parameters::get("parallel->debug->skip check mesh change", skip);
@@ -1155,7 +1153,7 @@ namespace AMDiS { namespace Parallel {
     }
 
     mpiComm.Barrier();
-    MSG("Parallel mesh adaption needed %.5f seconds\n", MPI::Wtime() - first);
+    MSG("Parallel mesh adaption needed %.5f seconds\n", t.elapsed());
 
 
     // === Update the DOF numbering and mappings. ===
@@ -1419,8 +1417,7 @@ namespace AMDiS { namespace Parallel {
     if (repartitioning == 0)
       return;
 
-    mpiComm.Barrier();
-    double timePoint = MPI::Wtime();
+    Timer t;
 
 #if (DEBUG != 0)
     ParallelDebug::testDoubleDofs(mesh);
@@ -1473,7 +1470,7 @@ namespace AMDiS { namespace Parallel {
       mpiComm.Barrier();
       repartitioningFailed = repartitioningWaitAfterFail;;
       MSG("Mesh partitioner created empty partition!\n");
-      MSG("Mesh repartitioning needed %.5f seconds\n", MPI::Wtime() - timePoint); 
+      MSG("Mesh repartitioning needed %.5f seconds\n", t.elapsed()); 
       return;
     }
 
@@ -1487,7 +1484,7 @@ namespace AMDiS { namespace Parallel {
       if (!partitioningSucceed || !partitioner->meshChanged()) {
 	mpiComm.Barrier();
 	repartitioningFailed = repartitioningWaitAfterFail;;
-	MSG("Mesh repartitioning needed %.5f seconds\n", MPI::Wtime() - timePoint);   
+	MSG("Mesh repartitioning needed %.5f seconds\n", t.elapsed());   
 	return;
       }
     }
@@ -1706,7 +1703,7 @@ namespace AMDiS { namespace Parallel {
 
 
     mpiComm.Barrier();
-    MSG("Mesh repartitioning needed %.5f seconds\n", MPI::Wtime() - timePoint);   
+    MSG("Mesh repartitioning needed %.5f seconds\n", t.elapsed());   
   }
 
 
@@ -1746,9 +1743,9 @@ namespace AMDiS { namespace Parallel {
   {
     FUNCNAME("MeshDistributor::createBoundaryDofs()");
 
+    Timer *t = NULL;
     if (printTimings)
-      MPI::COMM_WORLD.Barrier();
-    double first = MPI::Wtime();
+      t = new Timer();
 
     // === Create DOF communicator. ===
 
@@ -1805,9 +1802,9 @@ namespace AMDiS { namespace Parallel {
       }
     }
 
-    if (printTimings) {
-      MPI::COMM_WORLD.Barrier();
-      MSG("MESHDIST 01 (createBoundaryDofs) needed %.5f seconds\n", MPI::Wtime() - first);
+    if (printTimings && t) {
+      MSG("MESHDIST 01 (createBoundaryDofs) needed %.5f seconds\n", t->elapsed());
+      delete t;
     }
   }
 
@@ -1828,8 +1825,7 @@ namespace AMDiS { namespace Parallel {
   {
     FUNCNAME("MeshDistributor::updateLocalGlobalNumbering()");
 
-    MPI::COMM_WORLD.Barrier();
-    double first = MPI::Wtime();
+    Timer t;
     
     mesh->dofCompress();
 
@@ -1927,8 +1923,7 @@ namespace AMDiS { namespace Parallel {
 				    debugOutputDir + "mpi-dbg", "dat");
 #endif
 
-    MPI::COMM_WORLD.Barrier();
-    MSG("Rebuild of parallel data structures needed %.5f seconds\n", MPI::Wtime() - first);
+    MSG("Rebuild of parallel data structures needed %.5f seconds\n", t.elapsed());
   }
 
 
@@ -1936,9 +1931,9 @@ namespace AMDiS { namespace Parallel {
   {
     FUNCNAME("MeshDistributor::updateParallelDofMappings()");
 
+    Timer *t = NULL;
     if (printTimings)
-      MPI::COMM_WORLD.Barrier();
-    double first = MPI::Wtime();
+      t = new Timer();
 
     TEST_EXIT(dofMaps.size())("No DOF mapping defined!\n");
 
@@ -1956,9 +1951,9 @@ namespace AMDiS { namespace Parallel {
 #endif
     }
 
-    if (printTimings) {
-      MPI::COMM_WORLD.Barrier();
-      MSG("MESHDIST 02 (update DOF mappings) needed %.5f seconds\n", MPI::Wtime() - first);
+    if (printTimings && t) {
+      MSG("MESHDIST 02 (update DOF mappings) needed %.5f seconds\n", t->elapsed());
+      delete t;
     }
   }
 

AMDiS/src/parallel/PetscSolverCahnHilliard.cc

@@ -25,21 +25,7 @@
 #include "parallel/PetscSolverGlobalMatrix.h"
 
 namespace AMDiS { namespace Parallel {
-  
-  class Multiplier: public AbstractFunction<double, double>
-  {
-  public:
-    Multiplier(double d_)   : AbstractFunction<double, double>(2),d(d_)  {};
-    
-    double operator()(const double& x) const {
-      return d*x;
-    }
-  protected:
-    double d;
-  };
-  
-  
-  
+
   using namespace std;
   
   /// solve Cahn-Hilliard Preconditioner
@@ -84,12 +70,13 @@ namespace AMDiS { namespace Parallel {
   
   PetscSolverCahnHilliard::PetscSolverCahnHilliard(string name, double *epsPtr, double *deltaPtr)
   : PetscSolverGlobalBlockMatrix(name),
-  massMatrixSolver(NULL),
-  laplaceMatrixSolver(NULL),
-  deltaKMatrixSolver(NULL),
-  useOldInitialGuess(false),
-  eps(epsPtr),
-  delta(deltaPtr) { 	
+    massMatrixSolver(NULL),
+    laplaceMatrixSolver(NULL),
+    deltaKMatrixSolver(NULL),
+    useOldInitialGuess(false),
+    eps(epsPtr),
+    delta(deltaPtr) 
+  { 	
     Parameters::get(initFileStr + "->use old initial guess", useOldInitialGuess);
   }
   
@@ -120,15 +107,15 @@ namespace AMDiS { namespace Parallel {
   
   void PetscSolverCahnHilliard::initSolver(KSP &ksp)
   {
-    FUNCNAME("PetscSolverCahnHilliard::initSolver()");
-    
     // Create FGMRES based outer solver
     KSPCreate(meshDistributor->getMpiComm(0), &ksp);
     KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
-    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
-    petsc_helper::setSolver(ksp, "ch_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 1000);
+    if (getInfo() >= 10)
+      KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL);
+    else if (getInfo() >= 20)
+      KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
+    petsc_helper::setSolver(ksp, "ch_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
     KSPSetFromOptions(ksp);
-    
   }
   
   
@@ -137,6 +124,8 @@ namespace AMDiS { namespace Parallel {
     FUNCNAME("PetscSolverCahnHilliard::initPreconditioner()");
     MSG("PetscSolverCahnHilliard::initPreconditioner()\n");
     
+    TEST_EXIT(eps && delta)("eps and/or delta pointers not set!\n");
+    
     //     KSPSetUp(kspInterior);
     
     PCSetType(pc, PCSHELL);
@@ -155,14 +144,13 @@ namespace AMDiS { namespace Parallel {
     DOFMatrix deltaKMatrix(feSpace, feSpace);
     Operator laplaceOp2(feSpace, feSpace);
     
-    if (phase) 
-    {  
+    if (phase) {  
       VecAtQP_ZOT zot(phase, NULL);
       massOp.addTerm(&zot);
       laplaceOp.addTerm(&zot); // M
-      VecAtQP_SOT sot2(phase, new Multiplier((*eps)*sqrt(*delta)));
+      VecAtQP_SOT sot2(phase, NULL, (*eps)*sqrt(*delta));
       laplaceOp.addTerm(&sot2); // eps*sqrt(delta)*K
-      VecAtQP_SOT sot(phase, new Multiplier(-(*delta)));
+      VecAtQP_SOT sot(phase, NULL, -(*delta));
       laplaceOp2.addTerm(&sot); // -delta*K
       massMatrix.assembleOperator(massOp);
       massMatrixSolver = createSubSolver(0, "mass_");
@@ -173,16 +161,13 @@ namespace AMDiS { namespace Parallel {
       laplaceMatrixSolver = createSubSolver(0, "laplace_");
       laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
       
-      
       // === matrix (-delta*K) ===
       deltaKMatrix.assembleOperator(laplaceOp2);
       deltaKMatrixSolver = createSubSolver(0, "laplace2_");
       deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
       
       
-    }       
-    else 
-    { 
+    } else { 
       Simple_ZOT zot;
       massOp.addTerm(&zot);
       laplaceOp.addTerm(&zot); // M
@@ -200,7 +185,6 @@ namespace AMDiS { namespace Parallel {
       laplaceMatrixSolver = createSubSolver(0, "laplace_");
       laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
       
-      
       // === matrix (-delta*K) ===
       deltaKMatrix.assembleOperator(laplaceOp2);
       deltaKMatrixSolver = createSubSolver(0, "laplace2_");
@@ -228,7 +212,7 @@ namespace AMDiS { namespace Parallel {
     //       PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS);
     //     }
     
-    petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
+    petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 2);
     //     petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1);
     //     {
       //             PC pc;

AMDiS/src/parallel/PetscSolverCahnHilliard.h

@@ -63,7 +63,7 @@ namespace AMDiS { namespace Parallel {
     void setChData(double *epsPtr, double *deltaPtr, SystemVector* vec, DOFVector<double> *p=NULL)
     {
       eps = epsPtr;
-      delta = deltaPtr;
+      delta = deltaPtr; // sqrt(tau)
       solution = vec;
       phase = p;
     }