PetscSolverCahnHilliard geht jetzt auch mit diffuse domain

AMDiS/src/parallel/PetscSolverCahnHilliard.cc

@@ -16,139 +16,189 @@
 #include "parallel/PetscSolverGlobalMatrix.h"
 
 namespace AMDiS {
-
+  
+  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
   PetscErrorCode pcChShell(PC pc, Vec b, Vec x) // solve Px=b
   {FUNCNAME("PCApply()");
   
-   void *ctx;
-    PCShellGetContext(pc, &ctx);
-    CahnHilliardData* data = static_cast<CahnHilliardData*>(ctx);
-    
-    Vec b1, b2, x1, x2;
-    VecNestGetSubVec(b, 0, &b1);
-    VecNestGetSubVec(b, 1, &b2);
-
-    VecNestGetSubVec(x, 0, &x1);
-    VecNestGetSubVec(x, 1, &x2);
-
-    Vec y1, y2;
-    VecDuplicate(b1, &y1);
-    VecDuplicate(b2, &y2);
-    
-//     MatGetDiagonal(data->matM, y2);
-//     VecReciprocal(y2);
-//     VecPointwiseMult(y1, y2, b1);
-    KSPSolve(data->kspMass, b1, y1); 			// M*y1 = b1    
-    MatMultAdd(data->matMinusDeltaK, y1, b2, x1); 	// -> x1 := b2-delta*K*y1
-    
-    KSPSolve(data->kspLaplace, x1, y2); 			// (M+eps*sqrt(delta))*y2 = x1
-    MatMult(data->matM, y2, x1); 			// x1 := M*y2
-    
-    KSPSolve(data->kspLaplace, x1, x2);			// (M+eps*sqrt(delta))*x2 = x1
-    double factor = (*data->eps)/sqrt(*data->delta);
-    VecCopy(x2, x1); 					// x1 := x2
-    VecAXPBYPCZ(x1, 1.0, factor, -factor, y1, y2);	// x1 = 1*y1 + factor*y2 - factor*x1
-    
-    VecDestroy(&y1);
-    VecDestroy(&y2);
-
+  void *ctx;
+  PCShellGetContext(pc, &ctx);
+  CahnHilliardData* data = static_cast<CahnHilliardData*>(ctx);
+  
+  Vec b1, b2, x1, x2;
+  VecNestGetSubVec(b, 0, &b1);
+  VecNestGetSubVec(b, 1, &b2);
+  
+  VecNestGetSubVec(x, 0, &x1);
+  VecNestGetSubVec(x, 1, &x2);
+  
+  Vec y1, y2;
+  VecDuplicate(b1, &y1);
+  VecDuplicate(b2, &y2);
+  
+  //     MatGetDiagonal(data->matM, y2);
+  //     VecReciprocal(y2);
+  //     VecPointwiseMult(y1, y2, b1);
+  KSPSolve(data->kspMass, b1, y1); 			// M*y1 = b1    
+  MatMultAdd(data->matMinusDeltaK, y1, b2, x1); 	// -> x1 := b2-delta*K*y1
+  
+  KSPSolve(data->kspLaplace, x1, y2); 			// (M+eps*sqrt(delta))*y2 = x1
+  MatMult(data->matM, y2, x1); 			// x1 := M*y2
+  
+  KSPSolve(data->kspLaplace, x1, x2);			// (M+eps*sqrt(delta))*x2 = x1
+  double factor = (*data->eps)/sqrt(*data->delta);
+  VecCopy(x2, x1); 					// x1 := x2
+  VecAXPBYPCZ(x1, 1.0, factor, -factor, y1, y2);	// x1 = 1*y1 + factor*y2 - factor*x1
+  
+  VecDestroy(&y1);
+  VecDestroy(&y2);
+  
   }
   
-
+  
   PetscSolverCahnHilliard::PetscSolverCahnHilliard(string name, double *epsPtr, double *deltaPtr)
-    : PetscSolverGlobalBlockMatrix(name),
-      massMatrixSolver(NULL),
-      laplaceMatrixSolver(NULL),
-      deltaKMatrixSolver(NULL),
-      useOldInitialGuess(false),
-      eps(epsPtr),
-      delta(deltaPtr) { 	
-	    Parameters::get(initFileStr + "->use old initial guess", useOldInitialGuess);
-      }
-
-
+  : PetscSolverGlobalBlockMatrix(name),
+  massMatrixSolver(NULL),
+  laplaceMatrixSolver(NULL),
+  deltaKMatrixSolver(NULL),
+  useOldInitialGuess(false),
+  eps(epsPtr),
+  delta(deltaPtr) { 	
+    Parameters::get(initFileStr + "->use old initial guess", useOldInitialGuess);
+  }
+  
+  
   void PetscSolverCahnHilliard::solvePetscMatrix(SystemVector &vec, 
 						 AdaptInfo *adaptInfo)
   {
     FUNCNAME("PetscSolverCahnHilliard::solvePetscMatrix()");
-
-/*    if (useOldInitialGuess) {
-      VecSet(getVecSolInterior(), 0.0);
+    
+  /*  if (useOldInitialGuess) {
+      if (getVecSolInterior())
+      {VecSet(getVecSolInterior(), 0.0);
       
       for (int i = 0; i < solution->getSize(); i++)
       {
 	Vec tmp;
-        VecNestGetSubVec(getVecSolInterior(), i, &tmp);
+	VecNestGetSubVec(getVecSolInterior(), i, &tmp);
 	setDofVector(tmp, solution->getDOFVector(i));
       }
       
       vecSolAssembly();
       KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
     }
-*/ KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
+    KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
+    }*/
     PetscSolverGlobalBlockMatrix::solvePetscMatrix(vec, adaptInfo);
   }
-      
+  
   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);
-         KSPSetFromOptions(ksp);
-	 
+    KSPSetFromOptions(ksp);
+    
   }
-
-
+  
+  
   void PetscSolverCahnHilliard::initPreconditioner(PC pc)
   {
     FUNCNAME("PetscSolverCahnHilliard::initPreconditioner()");
     MSG("PetscSolverCahnHilliard::initPreconditioner()\n");
-
-//     KSPSetUp(kspInterior);
-
+    
+    //     KSPSetUp(kspInterior);
+    
     PCSetType(pc, PCSHELL);
     PCShellSetApply(pc, pcChShell);
     PCShellSetContext(pc, &matShellContext);
-
+    
     const FiniteElemSpace *feSpace = componentSpaces[0];
     
     // === matrix M ===	
-    DOFMatrix massMatrix(feSpace, feSpace);
-    Operator massOp(feSpace, feSpace);
-    Simple_ZOT zot;
-    massOp.addTerm(&zot);
-    massMatrix.assembleOperator(massOp);
-    massMatrixSolver = createSubSolver(0, "mass_");
-    massMatrixSolver->fillPetscMatrix(&massMatrix);
-    
-    
-    // === matrix (M + eps*sqrt(delta)*K) ===
     DOFMatrix laplaceMatrix(feSpace, feSpace);
     Operator laplaceOp(feSpace, feSpace);
-    laplaceOp.addTerm(&zot); // M
-    Simple_SOT sot2((*eps)*sqrt(*delta));
-    laplaceOp.addTerm(&sot2); // eps*sqrt(delta)*K
-    laplaceMatrix.assembleOperator(laplaceOp);
-    laplaceMatrixSolver = createSubSolver(0, "laplace_");
-    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
-
     
-    // === matrix (-delta*K) ===
+    DOFMatrix massMatrix(feSpace, feSpace);
+    Operator massOp(feSpace, feSpace);
+    
     DOFMatrix deltaKMatrix(feSpace, feSpace);
     Operator laplaceOp2(feSpace, feSpace);
-    Simple_SOT sot(-(*delta));
-    laplaceOp2.addTerm(&sot); // -delta*K
-    deltaKMatrix.assembleOperator(laplaceOp2);
-    deltaKMatrixSolver = createSubSolver(0, "laplace2_");
-    deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
-
+    
+    if (phase) 
+    {  
+      VecAtQP_ZOT zot(phase, NULL);
+      massOp.addTerm(&zot);
+      laplaceOp.addTerm(&zot); // M
+      VecAtQP_SOT sot2(phase, new Multiplier((*eps)*sqrt(*delta)));
+      laplaceOp.addTerm(&sot2); // eps*sqrt(delta)*K
+      VecAtQP_SOT sot(phase, new Multiplier(-(*delta)));
+      laplaceOp2.addTerm(&sot); // -delta*K
+      massMatrix.assembleOperator(massOp);
+      massMatrixSolver = createSubSolver(0, "mass_");
+      massMatrixSolver->fillPetscMatrix(&massMatrix);
+      
+      // === matrix (M + eps*sqrt(delta)*K) ===
+      laplaceMatrix.assembleOperator(laplaceOp);
+      laplaceMatrixSolver = createSubSolver(0, "laplace_");
+      laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
+      
+      
+      // === matrix (-delta*K) ===
+      deltaKMatrix.assembleOperator(laplaceOp2);
+      deltaKMatrixSolver = createSubSolver(0, "laplace2_");
+      deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
+      
+      
+    }       
+    else 
+    { 
+      Simple_ZOT zot;
+      massOp.addTerm(&zot);
+      laplaceOp.addTerm(&zot); // M
+      Simple_SOT sot2((*eps)*sqrt(*delta));
+      laplaceOp.addTerm(&sot2); // eps*sqrt(delta)*K
+      Simple_SOT sot(-(*delta));    
+      laplaceOp2.addTerm(&sot); // -delta*K
+      
+      massMatrix.assembleOperator(massOp);
+      massMatrixSolver = createSubSolver(0, "mass_");
+      massMatrixSolver->fillPetscMatrix(&massMatrix);
+      
+      // === matrix (M + eps*sqrt(delta)*K) ===
+      laplaceMatrix.assembleOperator(laplaceOp);
+      laplaceMatrixSolver = createSubSolver(0, "laplace_");
+      laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
+      
+      
+      // === matrix (-delta*K) ===
+      deltaKMatrix.assembleOperator(laplaceOp2);
+      deltaKMatrixSolver = createSubSolver(0, "laplace2_");
+      deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
+    }
+    
+    
+    
     
     // === Setup solver ===
     matShellContext.kspMass = massMatrixSolver->getSolver();
@@ -159,43 +209,43 @@ namespace AMDiS {
     matShellContext.delta = delta;
     
     matShellContext.mpiCommGlobal= &(meshDistributor->getMpiComm(0));
-
+    
     petsc_helper::setSolver(matShellContext.kspMass, "", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
-//     petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1);
-//     {
-//             PC pc;
-//       KSPGetPC(matShellContext.kspMass, &pc);
-//       PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS);
-//     }
-
+    //     petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1);
+    //     {
+      //             PC pc;
+    //       KSPGetPC(matShellContext.kspMass, &pc);
+    //       PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS);
+    //     }
+    
     petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
-//     petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1);
-//     {
-//             PC pc;
-//       KSPGetPC(matShellContext.kspLaplace, &pc);
-//       PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS);
-//     }
-     PCSetFromOptions(pc);
+    //     petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1);
+    //     {
+      //             PC pc;
+    //       KSPGetPC(matShellContext.kspLaplace, &pc);
+    //       PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS);
+    //     }
+    PCSetFromOptions(pc);
   }
-
+  
   
   PetscSolver* PetscSolverCahnHilliard::createSubSolver(int component,
 							string kspPrefix)
   {
     FUNCNAME("PetscSolverCahnHilliard::createSubSolver()");
-
+    
     vector<const FiniteElemSpace*> fe;
     fe.push_back(componentSpaces[component]);
-
+    
     PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
     subSolver->setKspPrefix(kspPrefix.c_str());
     subSolver->setMeshDistributor(meshDistributor, 0);
     subSolver->init(fe, fe);
-
+    
     ParallelDofMapping &subDofMap = subSolver->getDofMap();
     subDofMap[0] = dofMap[component];
     subDofMap.update();
-
+    
     return subSolver;
   }
   
@@ -206,21 +256,21 @@ namespace AMDiS {
     massMatrixSolver->destroyMatrixData();
     laplaceMatrixSolver->destroyMatrixData();
     deltaKMatrixSolver->destroyMatrixData();
-
+    
     massMatrixSolver->destroyVectorData();
     laplaceMatrixSolver->destroyVectorData();
     deltaKMatrixSolver->destroyVectorData();
-
+    
     
     delete massMatrixSolver;
     massMatrixSolver = NULL;
-
+    
     delete laplaceMatrixSolver;
     laplaceMatrixSolver = NULL;
-
+    
     delete deltaKMatrixSolver;
     deltaKMatrixSolver = NULL;
   }
-
-
+  
+  
 }

AMDiS/src/parallel/PetscSolverCahnHilliard.h

@@ -58,11 +58,12 @@ namespace AMDiS {
     
     PetscSolverCahnHilliard(string name, double *epsPtr = NULL, double *deltaPtr = NULL);
 
-    void setChData(double *epsPtr, double *deltaPtr, SystemVector* vec)
+    void setChData(double *epsPtr, double *deltaPtr, SystemVector* vec, DOFVector<double> *p=NULL)
     {
       eps = epsPtr;
       delta = deltaPtr;
       solution = vec;
+      phase = p;
     }
 
   protected:
@@ -83,6 +84,8 @@ namespace AMDiS {
     bool useOldInitialGuess;
     
     SystemVector* solution;
+    
+    DOFVector<double> *phase;
 
     double *eps, *delta;
   };