PetscSolverCahnHilliard can now handle Dirichlet BC for phi

AMDiS/src/parallel/PetscSolverCahnHilliard.cc

@@ -56,7 +56,7 @@ namespace AMDiS { namespace Parallel {
   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
+  KSPSolve(data->kspLaplace2, 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
@@ -144,6 +144,8 @@ namespace AMDiS { namespace Parallel {
     DOFMatrix deltaKMatrix(feSpace, feSpace);
     Operator laplaceOp2(feSpace, feSpace);
     
+    DOFMatrix massMatrix2(feSpace, feSpace);
+    
     if (phase) {  
       VecAtQP_ZOT zot(phase, NULL);
       massOp.addTerm(&zot);
@@ -156,6 +158,10 @@ namespace AMDiS { namespace Parallel {
       massMatrixSolver = createSubSolver(0, "mass_");
       massMatrixSolver->fillPetscMatrix(&massMatrix);
       
+      massMatrix2.assembleOperator(massOp);        
+      massMatrixSolver2 = createSubSolver(0, "laplace_");       
+      massMatrixSolver2->fillPetscMatrix(&massMatrix2);  
+      
       // === matrix (M + eps*sqrt(delta)*K) ===
       laplaceMatrix.assembleOperator(laplaceOp);
       laplaceMatrixSolver = createSubSolver(0, "laplace_");
@@ -180,6 +186,10 @@ namespace AMDiS { namespace Parallel {
       massMatrixSolver = createSubSolver(0, "mass_");
       massMatrixSolver->fillPetscMatrix(&massMatrix);
       
+      massMatrix2.assembleOperator(massOp);        
+      massMatrixSolver2 = createSubSolver(0, "laplace_");       
+      massMatrixSolver2->fillPetscMatrix(&massMatrix2);  
+      
       // === matrix (M + eps*sqrt(delta)*K) ===
       laplaceMatrix.assembleOperator(laplaceOp);
       laplaceMatrixSolver = createSubSolver(0, "laplace_");
@@ -190,8 +200,10 @@ namespace AMDiS { namespace Parallel {
       deltaKMatrixSolver = createSubSolver(0, "laplace2_");
       deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
     }
-    
-    
+        
+    Mat matK0MalMinusEps2;
+    MatNestGetSubMat(mat[0][0], 0, 1, &(matK0MalMinusEps2));
+    MatAXPY(massMatrixSolver2->getMatInterior(), -sqrt(*delta) / (*eps), matK0MalMinusEps2, DIFFERENT_NONZERO_PATTERN);    
     
     
     // === Setup solver ===
@@ -199,6 +211,7 @@ namespace AMDiS { namespace Parallel {
     matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
     matShellContext.matM = massMatrixSolver->getMatInterior();    
     matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior();   
+    matShellContext.kspLaplace2 = massMatrixSolver2->getSolver();  
     matShellContext.eps = eps;
     matShellContext.delta = delta;
     
@@ -213,6 +226,8 @@ namespace AMDiS { namespace Parallel {
     //     }
     
     petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 2);
+    
+    petsc_helper::setSolver(matShellContext.kspLaplace2, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 2);
     //     petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1);
     //     {
       //             PC pc;
@@ -250,15 +265,19 @@ namespace AMDiS { namespace Parallel {
     massMatrixSolver->destroyMatrixData();
     laplaceMatrixSolver->destroyMatrixData();
     deltaKMatrixSolver->destroyMatrixData();
+    massMatrixSolver2->destroyMatrixData();
     
     massMatrixSolver->destroyVectorData();
     laplaceMatrixSolver->destroyVectorData();
     deltaKMatrixSolver->destroyVectorData();
-    
+    massMatrixSolver2->destroyVectorData();
     
     delete massMatrixSolver;
     massMatrixSolver = NULL;
     
+    delete massMatrixSolver2;    
+    massMatrixSolver2 = NULL;
+    
     delete laplaceMatrixSolver;
     laplaceMatrixSolver = NULL;
     

AMDiS/src/parallel/PetscSolverCahnHilliard.h

@@ -32,7 +32,7 @@ namespace AMDiS { namespace Parallel {
   using namespace std;
 
   struct CahnHilliardData {
-    KSP kspMass, kspLaplace;
+    KSP kspMass, kspLaplace, kspLaplace2;
     Mat matM, matMinusDeltaK;
     double *eps, *delta; 
     MPI::Intracomm *mpiCommGlobal;
@@ -79,7 +79,7 @@ namespace AMDiS { namespace Parallel {
   private:
     int pressureComponent;
 
-    PetscSolver *massMatrixSolver, *laplaceMatrixSolver, *deltaKMatrixSolver;
+    PetscSolver *massMatrixSolver, *massMatrixSolver2, *laplaceMatrixSolver, *deltaKMatrixSolver;
 
     CahnHilliardData matShellContext;