From 3484b9be9280aee3c4574d80b80d80b0a23b09e5 Mon Sep 17 00:00:00 2001
From: Sebastian Aland <sebastian.aland@yahoo.com>
Date: Tue, 12 Feb 2013 12:51:57 +0000
Subject: [PATCH] PetscSolverCahnHilliard2 korrigiert
---
.../src/parallel/PetscSolverCahnHilliard2.cc | 209 +++++++++---------
AMDiS/src/parallel/PetscSolverCahnHilliard2.h | 46 +++-
2 files changed, 142 insertions(+), 113 deletions(-)
diff --git a/AMDiS/src/parallel/PetscSolverCahnHilliard2.cc b/AMDiS/src/parallel/PetscSolverCahnHilliard2.cc
index e7a996f7..ae421783 100644
--- a/AMDiS/src/parallel/PetscSolverCahnHilliard2.cc
+++ b/AMDiS/src/parallel/PetscSolverCahnHilliard2.cc
@@ -13,14 +13,15 @@
#include "PetscSolverCahnHilliard2.h"
#include "parallel/PetscHelper.h"
-#include "parallel/PetscSolverGlobalMatrix.h"
+#include "TransformDOF.h"
namespace AMDiS {
+
using namespace std;
-
- /// solve Schurcomplement Pre-Preconditioner
- PetscErrorCode pcChShell2b(PC pc, Vec b, Vec x) // solve Px=b
+
+
+ PetscErrorCode pcChShell2b(PC pc, Vec b, Vec x) // solve Px=b
{FUNCNAME("PCApply()");
void *ctx;
PCShellGetContext(pc, &ctx);
@@ -31,41 +32,26 @@ namespace AMDiS {
VecDuplicate(b, &y2);
KSPSolve(data->kspMplusK, b, y1);
- MatMult(data->matM, y1, y2);
+ MatMult(data->matMass, y1, y2);
KSPSolve(data->kspMplusK, y2, x);
}
/// solve Cahn-Hilliard Preconditioner
- PetscErrorCode pcChShell2(PC pc, Vec b, Vec x) // solve Px=b
+ PetscErrorCode pcChSchurShell(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
+ MatGetDiagonal(data->matMass, x);
+ VecReciprocal(x);
+ MatDiagonalScale(K, x, 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
+ MatAYPX(S, sqr(*data->eps)/(*data->delta), data->matMass, DIFFERENT_NONZERO_PATTERN); // S = eps^2/delta*S + M
/// create new solver for S
KSP kspS;
@@ -79,35 +65,46 @@ namespace AMDiS {
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);
+ KSPSolve(kspS, b, x); // S*x2 = x1
+
MatDestroy(&S);
MatDestroy(&K);
KSPDestroy(&kspS);
}
+
- PetscSolverCahnHilliard2::PetscSolverCahnHilliard2(string name, double *epsPtr, double *deltaPtr)
- : PetscSolverGlobalBlockMatrix(name),
+ PetscSolverCahnHilliard2::PetscSolverCahnHilliard2(string name)
+ : PetscSolverGlobalMatrix(name),
+ useOldInitialGuess(false),
+ laplaceSolutionMode(0),
massMatrixSolver(NULL),
- laplaceMatrixSolver(NULL),
+ laplaceMatrixSolver(NULL),
deltaKMatrixSolver(NULL),
- eps(epsPtr),
- delta(deltaPtr) { }
+ phase(NULL)
+ {
+ Parameters::get(initFileStr + "->cahnhilliard->use old initial guess",
+ useOldInitialGuess);
+ }
+
+
+ void PetscSolverCahnHilliard2::solvePetscMatrix(SystemVector &vec,
+ AdaptInfo *adaptInfo)
+ {
+ FUNCNAME("PetscSolverCahnHilliard2::solvePetscMatrix()");
+
+ if (useOldInitialGuess) {
+ VecSet(getVecSolInterior(), 0.0);
+
+ for (int i = 0; i < solution->getSize(); i++)
+ setDofVector(getVecSolInterior(), solution->getDOFVector(i), i, true);
+
+ vecSolAssembly();
+ KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
+ }
+
+ PetscSolverGlobalMatrix::solvePetscMatrix(vec, adaptInfo);
+ }
void PetscSolverCahnHilliard2::initSolver(KSP &ksp)
@@ -115,27 +112,63 @@ namespace AMDiS {
FUNCNAME("PetscSolverCahnHilliard2::initSolver()");
// Create FGMRES based outer solver
- KSPCreate(meshDistributor->getMpiComm(0), &ksp);
+
+ MSG("CREATE POS 1: %p\n", &ksp);
+ KSPCreate(domainComm, &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);
+ petsc_helper::setSolver(ksp, "ch_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 10000);
+
+
}
void PetscSolverCahnHilliard2::initPreconditioner(PC pc)
{
FUNCNAME("PetscSolverCahnHilliard2::initPreconditioner()");
- MSG("PetscSolverCahnHilliard2::initPreconditioner()\n");
-// KSPSetUp(kspInterior);
+ MPI::COMM_WORLD.Barrier();
+ double wtime = MPI::Wtime();
+
+ int dim = componentSpaces[0]->getMesh()->getDim();
+
+ vector<int> chPotentialComponent;
+ chPotentialComponent.push_back(0);
+ vector<int> chSchurComponent;
+ chSchurComponent.push_back(1);
- PCSetType(pc, PCSHELL);
- PCShellSetApply(pc, pcChShell2);
- PCShellSetContext(pc, &matShellContext);
+ PCSetType(pc, PCFIELDSPLIT);
+ PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
+ PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
- const FiniteElemSpace *feSpace = componentSpaces[0];
-
- // === matrix M ===
+ createFieldSplit(pc, "ch_potential", chPotentialComponent);
+ createFieldSplit(pc, "ch_schur", chSchurComponent);
+ PCSetFromOptions(pc);
+
+ KSPSetUp(kspInterior);
+
+ KSP *subKspCH;
+ int nSubKspCH;
+ PCFieldSplitGetSubKSP(pc, &nSubKspCH, &subKspCH);
+
+
+ TEST_EXIT(nSubKspCH == 2)
+ ("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");
+
+ KSP kspChPotential = subKspCH[0];
+ KSP kspChSchur = subKspCH[1];
+ PetscFree(subKspCH);
+
+ KSPSetType(kspChSchur, KSPPREONLY);
+ PC pcSub;
+ KSPGetPC(kspChSchur, &pcSub);
+ PCSetType(pcSub, PCSHELL);
+ PCShellSetApply(pcSub, pcChSchurShell);
+ PCShellSetContext(pcSub, &matShellContext);
+
+ const FiniteElemSpace *feSpace = componentSpaces[0];
+
+ // === Mass matrix solver ===
DOFMatrix massMatrix(feSpace, feSpace);
Operator massOp(feSpace, feSpace);
Simple_ZOT zot;
@@ -144,8 +177,7 @@ namespace AMDiS {
massMatrixSolver = createSubSolver(0, "mass_");
massMatrixSolver->fillPetscMatrix(&massMatrix);
-
- // === matrix (M + eps*sqrt(delta)*K) ===
+ // === Laplace matrix solver ===
DOFMatrix laplaceMatrix(feSpace, feSpace);
Operator laplaceOp(feSpace, feSpace);
laplaceOp.addTerm(&zot); // M
@@ -166,80 +198,47 @@ namespace AMDiS {
// === 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.kspMass = massMatrixSolver->getSolver();
+ matShellContext.matMass = massMatrixSolver->getMatInterior();
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);
+ petsc_helper::setSolver(kspChPotential, "chPotential", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
+
+ MSG("Setup of Cahn-Hilliard preconditioner needed %.5f seconds\n",
+ MPI::Wtime() - wtime);
}
-
- 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()");
-
+ FUNCNAME("PetscSolverCahnHilliard2::exitPreconditioner()");
+
massMatrixSolver->destroyMatrixData();
+ massMatrixSolver->destroyVectorData();
laplaceMatrixSolver->destroyMatrixData();
+ laplaceMatrixSolver->destroyVectorData();
deltaKMatrixSolver->destroyMatrixData();
+ deltaKMatrixSolver->destroyVectorData();
+
massMatrixSolver->destroyVectorData();
laplaceMatrixSolver->destroyVectorData();
deltaKMatrixSolver->destroyVectorData();
-
delete massMatrixSolver;
massMatrixSolver = NULL;
delete laplaceMatrixSolver;
laplaceMatrixSolver = NULL;
-
+
delete deltaKMatrixSolver;
deltaKMatrixSolver = NULL;
- }
-
+ }
}
diff --git a/AMDiS/src/parallel/PetscSolverCahnHilliard2.h b/AMDiS/src/parallel/PetscSolverCahnHilliard2.h
index 31bf7570..1c4b924c 100644
--- a/AMDiS/src/parallel/PetscSolverCahnHilliard2.h
+++ b/AMDiS/src/parallel/PetscSolverCahnHilliard2.h
@@ -23,7 +23,7 @@
#ifndef AMDIS_PETSC_SOLVER_CAHN_HILLIARD2_H
#define AMDIS_PETSC_SOLVER_CAHN_HILLIARD2_H
-#include "parallel/PetscSolverGlobalBlockMatrix.h"
+#include "parallel/PetscSolverGlobalMatrix.h"
namespace AMDiS {
@@ -31,12 +31,12 @@ namespace AMDiS {
struct CahnHilliardData2 {
KSP kspMass, kspMinusDeltaK, kspMplusK;
- Mat matM, matMinusDeltaK;
+ Mat matMass, matMinusDeltaK;
double *eps, *delta;
MPI::Intracomm *mpiCommGlobal;
};
- class PetscSolverCahnHilliard2 : public PetscSolverGlobalBlockMatrix
+ class PetscSolverCahnHilliard2 : public PetscSolverGlobalMatrix
{
public:
/// Creator class
@@ -51,33 +51,63 @@ namespace AMDiS {
return new PetscSolverCahnHilliard2(this->name);
}
};
-
- PetscSolverCahnHilliard2(string name, double *epsPtr = NULL, double *deltaPtr = NULL);
+ PetscSolverCahnHilliard2(string name);
+
+ void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo);
+
void setChData(double *epsPtr, double *deltaPtr)
{
eps = epsPtr;
delta = deltaPtr;
}
+
+
+ void setPhase(DOFVector<double> *d, double eP3=0)
+ {
+ phase = d;
+ epsPhase3 = eP3;
+ }
protected:
void initSolver(KSP &ksp);
void initPreconditioner(PC pc);
+
void exitPreconditioner(PC pc);
-
- PetscSolver* createSubSolver(int component, string kspPrefix);
private:
int pressureComponent;
+ bool pressureNullSpace;
+
+ /// If true, old solution is used for initial guess in solver phase.
+ bool useOldInitialGuess;
+
+ /// 0: approximate solve 1: direct solver
+ int velocitySolutionMode;
+
+ /// 0: approximate solve 1: direct solver
+ int massSolutionMode;
+
+ /// 0: approximate solve 1: direct solver
+ int laplaceSolutionMode;
+
PetscSolver *massMatrixSolver, *laplaceMatrixSolver, *deltaKMatrixSolver;
CahnHilliardData2 matShellContext;
double *eps, *delta;
+
+
+ double epsPhase3;
+
+ SystemVector* solution;
+
+ DOFVector<double>* phase;
+
};
}
-#endif // AMDIS_PETSC_SOLVER_CAHN_HILLIARD2_H
+#endif
--
GitLab