NSCH preconditioner hinzugefuegt

AMDiS/CMakeLists.txt

@@ -268,6 +268,8 @@ if(ENABLE_PARALLEL_DOMAIN)
 			${SOURCE_DIR}/parallel/PetscSolverGlobalMatrix.cc
 			${SOURCE_DIR}/parallel/PetscSolverGlobalBlockMatrix.cc
 			${SOURCE_DIR}/parallel/PetscSolverNavierStokes.cc
+                        ${SOURCE_DIR}/parallel/PetscSolverCahnHilliard2.cc
+                        ${SOURCE_DIR}/parallel/PetscSolverCahnHilliard.cc                        
 			${SOURCE_DIR}/parallel/PetscSolverSchur.cc)
 	elseif(ENABLE_PARALLEL_DOMAIN STREQUAL "PMTL")
 		set(MTL_INCLUDE_DIR "")

AMDiS/src/parallel/PetscSolverCahnHilliard.cc

+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+#include "PetscSolverCahnHilliard.h"
+#include "parallel/PetscHelper.h"
+#include "parallel/PetscSolverGlobalMatrix.h"
+
+namespace AMDiS {
+
+  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);
+
+  }
+  
+
+  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);
+      }
+
+
+  void PetscSolverCahnHilliard::solvePetscMatrix(SystemVector &vec, 
+						 AdaptInfo *adaptInfo)
+  {
+    FUNCNAME("PetscSolverCahnHilliard::solvePetscMatrix()");
+
+/*    if (useOldInitialGuess) {
+      VecSet(getVecSolInterior(), 0.0);
+      
+      for (int i = 0; i < solution->getSize(); i++)
+      {
+	Vec tmp;
+        VecNestGetSubVec(getVecSolInterior(), i, &tmp);
+	setDofVector(tmp, solution->getDOFVector(i));
+      }
+      
+      vecSolAssembly();
+      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);
+	 
+  }
+
+
+  void PetscSolverCahnHilliard::initPreconditioner(PC pc)
+  {
+    FUNCNAME("PetscSolverCahnHilliard::initPreconditioner()");
+    MSG("PetscSolverCahnHilliard::initPreconditioner()\n");
+
+//     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 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);
+
+    
+    // === Setup solver ===
+    matShellContext.kspMass = massMatrixSolver->getSolver();
+    matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
+    matShellContext.matM = massMatrixSolver->getMatInterior();    
+    matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior();   
+    matShellContext.eps = eps;
+    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.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);
+  }
+
+  
+  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;
+  }
+  
+  void PetscSolverCahnHilliard::exitPreconditioner(PC pc)
+  {
+    FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");
+    
+    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

+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ==  http://www.amdis-fem.org                                              ==
+// ==                                                                        ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file PetscSolverCahnHilliard.h */
+
+#ifndef AMDIS_PETSC_SOLVER_CAHN_HILLIARD_H
+#define AMDIS_PETSC_SOLVER_CAHN_HILLIARD_H
+
+#include "parallel/PetscSolverGlobalBlockMatrix.h"
+
+namespace AMDiS {
+
+  using namespace std;
+
+  struct CahnHilliardData {
+    KSP kspMass, kspLaplace;
+    Mat matM, matMinusDeltaK;
+    double *eps, *delta; 
+    MPI::Intracomm *mpiCommGlobal;
+
+    
+  };
+
+  class PetscSolverCahnHilliard : public PetscSolverGlobalBlockMatrix
+  {
+  public:
+    /// Creator class
+    class Creator : public OEMSolverCreator
+    {
+    public:
+      virtual ~Creator() {}
+
+      /// Returns a new PetscSolver object.
+      OEMSolver* create() 
+      { 
+	return new PetscSolverCahnHilliard(this->name); 
+      }
+    };
+    
+    void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo);
+    
+    PetscSolverCahnHilliard(string name, double *epsPtr = NULL, double *deltaPtr = NULL);
+
+    void setChData(double *epsPtr, double *deltaPtr, SystemVector* vec)
+    {
+      eps = epsPtr;
+      delta = deltaPtr;
+      solution = vec;
+    }
+
+  protected:
+    void initSolver(KSP &ksp);
+
+    void initPreconditioner(PC pc);
+    void exitPreconditioner(PC pc);
+    
+    PetscSolver* createSubSolver(int component, string kspPrefix);
+
+  private:
+    int pressureComponent;
+
+    PetscSolver *massMatrixSolver, *laplaceMatrixSolver, *deltaKMatrixSolver;
+
+    CahnHilliardData matShellContext;
+
+    bool useOldInitialGuess;
+    
+    SystemVector* solution;
+
+    double *eps, *delta;
+  };
+
+}
+
+#endif // AMDIS_PETSC_SOLVER_CAHN_HILLIARD_H

AMDiS/src/parallel/PetscSolverCahnHilliard2.cc

+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+#include "PetscSolverCahnHilliard2.h"
+#include "parallel/PetscHelper.h"
+#include "parallel/PetscSolverGlobalMatrix.h"
+
+namespace AMDiS {
+
+  using namespace std;
+  
+  /// solve Schurcomplement Pre-Preconditioner
+  PetscErrorCode pcChShell2b(PC pc, Vec b, Vec x) // solve Px=b
+  {FUNCNAME("PCApply()");
+    void *ctx;
+    PCShellGetContext(pc, &ctx);
+    CahnHilliardData2* data = static_cast<CahnHilliardData2*>(ctx);
+    
+    Vec y1, y2;
+    VecDuplicate(b, &y1);
+    VecDuplicate(b, &y2);
+    
+    KSPSolve(data->kspMplusK, b, y1);
+    MatMult(data->matM, y1, y2); 
+    KSPSolve(data->kspMplusK, y2, x);
+  }
+
+  /// solve Cahn-Hilliard Preconditioner
+  PetscErrorCode pcChShell2(PC pc, Vec b, Vec x) // solve Px=b
+  {FUNCNAME("PCApply()");
+    void *ctx;
+    PCShellGetContext(pc, &ctx);
+    CahnHilliardData2* data = static_cast<CahnHilliardData2*>(ctx);
+    
+    Vec b1, b2, x1, x2;
+    VecNestGetSubVec(b, 0, &b1);
+    VecNestGetSubVec(b, 1, &b2);
+
+    VecNestGetSubVec(x, 0, &x1);
+    VecNestGetSubVec(x, 1, &x2);
+
+    Vec y1;
+    VecDuplicate(b1, &y1);
+    
+//     VecPointwiseMult(y1, y2, b1);
+    KSPSolve(data->kspMass, b1, y1); 			// M*y1 = b1    
+    
+    /// create S = M + eps^2*delta*K*D^(-1)*K where D=diag(M)
+    Mat K, S;
+    MatDuplicate(data->matMinusDeltaK, MAT_COPY_VALUES, &K);
+    
+    MatGetDiagonal(data->matM, x1);
+    VecReciprocal(x1);
+    MatDiagonalScale(K, x1, PETSC_NULL); 						// K' := -delta*D^(-1)*K
+    MatMatMult(data->matMinusDeltaK, K, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &S); 		// S := -delta*K*K'
+    MatAYPX(S, sqr(*data->eps)/(*data->delta), data->matM, DIFFERENT_NONZERO_PATTERN); 	// S = eps^2/delta*S + M
+    
+    MatMultAdd(data->matMinusDeltaK, y1, b2, x1); 					// x1 := b2-delta*K*y1
+    
+    /// create new solver for S
+    KSP kspS;
+    KSPCreate(*data->mpiCommGlobal, &kspS);
+    KSPSetOperators(kspS, S, S, SAME_NONZERO_PATTERN);    
+    petsc_helper::setSolver(kspS, "S_", KSPFGMRES, PCSHELL, 1e-6, 1e-8, 1);
+    {
+      PC pc;
+      KSPGetPC(kspS, &pc);
+      PCShellSetApply(pc, pcChShell2b);
+      PCShellSetContext(pc, data);
+    }
+    
+    KSPSolve(kspS, x1, x2); 				// S*x2 = x1
+    MatMult(data->matM, x2, y1); 			// y1 := M*x2
+    VecAXPY(y1, -1.0, b2);				// y1 := y1 - b2
+    
+    // Project out constant null space
+    {
+      int vecSize;
+      VecGetSize(y1, &vecSize);
+      PetscScalar vecSum;
+      VecSum(y1, &vecSum);
+      vecSum = vecSum / static_cast<PetscScalar>(-1.0 * vecSize);
+      VecShift(y1, vecSum);
+    }
+    KSPSolve(data->kspMinusDeltaK, y1, x1);		// -delta*K*x1 = y1
+    
+    VecDestroy(&y1);
+    MatDestroy(&S);
+    MatDestroy(&K);
+    KSPDestroy(&kspS);
+  }
+  
+
+  PetscSolverCahnHilliard2::PetscSolverCahnHilliard2(string name, double *epsPtr, double *deltaPtr)
+    : PetscSolverGlobalBlockMatrix(name),
+      massMatrixSolver(NULL),
+      laplaceMatrixSolver(NULL),
+      deltaKMatrixSolver(NULL),
+      eps(epsPtr),
+      delta(deltaPtr) {  }
+
+
+  void PetscSolverCahnHilliard2::initSolver(KSP &ksp)
+  {
+    FUNCNAME("PetscSolverCahnHilliard2::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);
+  }
+
+
+  void PetscSolverCahnHilliard2::initPreconditioner(PC pc)
+  {
+    FUNCNAME("PetscSolverCahnHilliard2::initPreconditioner()");
+    MSG("PetscSolverCahnHilliard2::initPreconditioner()\n");
+
+//     KSPSetUp(kspInterior);
+
+    PCSetType(pc, PCSHELL);
+    PCShellSetApply(pc, pcChShell2);
+    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, "MpK_");
+    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
+    
+    // === matrix (-delta*K) ===
+    DOFMatrix deltaKMatrix(feSpace, feSpace);
+    Operator laplaceOp2(feSpace, feSpace);
+    Simple_SOT sot(-(*delta));
+    laplaceOp2.addTerm(&sot); // -delta*K
+    deltaKMatrix.assembleOperator(laplaceOp2);
+    deltaKMatrixSolver = createSubSolver(0, "laplace_");
+    deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
+
+    
+    // === Setup solver ===
+    matShellContext.kspMass = massMatrixSolver->getSolver();
+    matShellContext.kspMinusDeltaK = deltaKMatrixSolver->getSolver();
+    matShellContext.kspMplusK = laplaceMatrixSolver->getSolver();
+    matShellContext.matM = massMatrixSolver->getMatInterior();    
+    matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior();   
+    matShellContext.eps = eps;
+    matShellContext.delta = delta;
+    matShellContext.mpiCommGlobal = &(meshDistributor->getMpiComm(0));
+
+    petsc_helper::setSolver(matShellContext.kspMass, "mass_", 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.kspMinusDeltaK, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
+    petsc_helper::setSolver(matShellContext.kspMplusK, "MpK_", 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);
+  }
+
+  
+  PetscSolver* PetscSolverCahnHilliard2::createSubSolver(int component,
+							string kspPrefix)
+  {
+    FUNCNAME("PetscSolverCahnHilliard2::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;
+  }
+  
+  
+  
+  
+  void PetscSolverCahnHilliard2::exitPreconditioner(PC pc)
+  {
+    FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");
+    
+    massMatrixSolver->destroyMatrixData();
+    laplaceMatrixSolver->destroyMatrixData();
+    deltaKMatrixSolver->destroyMatrixData();
+
+    massMatrixSolver->destroyVectorData();
+    laplaceMatrixSolver->destroyVectorData();
+    deltaKMatrixSolver->destroyVectorData();
+
+    
+    delete massMatrixSolver;
+    massMatrixSolver = NULL;
+
+    delete laplaceMatrixSolver;
+    laplaceMatrixSolver = NULL;
+
+    delete deltaKMatrixSolver;
+    deltaKMatrixSolver = NULL;
+  }
+