Commit 3484b9be authored by Sebastian Aland's avatar Sebastian Aland
Browse files

PetscSolverCahnHilliard2 korrigiert

parent 0e7a58de
...@@ -13,14 +13,15 @@ ...@@ -13,14 +13,15 @@
#include "PetscSolverCahnHilliard2.h" #include "PetscSolverCahnHilliard2.h"
#include "parallel/PetscHelper.h" #include "parallel/PetscHelper.h"
#include "parallel/PetscSolverGlobalMatrix.h" #include "TransformDOF.h"
namespace AMDiS { namespace AMDiS {
using namespace std; 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()"); {FUNCNAME("PCApply()");
void *ctx; void *ctx;
PCShellGetContext(pc, &ctx); PCShellGetContext(pc, &ctx);
...@@ -31,41 +32,26 @@ namespace AMDiS { ...@@ -31,41 +32,26 @@ namespace AMDiS {
VecDuplicate(b, &y2); VecDuplicate(b, &y2);
KSPSolve(data->kspMplusK, b, y1); KSPSolve(data->kspMplusK, b, y1);
MatMult(data->matM, y1, y2); MatMult(data->matMass, y1, y2);
KSPSolve(data->kspMplusK, y2, x); KSPSolve(data->kspMplusK, y2, x);
} }
/// solve Cahn-Hilliard Preconditioner /// 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()"); {FUNCNAME("PCApply()");
void *ctx; void *ctx;
PCShellGetContext(pc, &ctx); PCShellGetContext(pc, &ctx);
CahnHilliardData2* data = static_cast<CahnHilliardData2*>(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) /// create S = M + eps^2*delta*K*D^(-1)*K where D=diag(M)
Mat K, S; Mat K, S;
MatDuplicate(data->matMinusDeltaK, MAT_COPY_VALUES, &K); MatDuplicate(data->matMinusDeltaK, MAT_COPY_VALUES, &K);
MatGetDiagonal(data->matM, x1); MatGetDiagonal(data->matMass, x);
VecReciprocal(x1); VecReciprocal(x);
MatDiagonalScale(K, x1, PETSC_NULL); // K' := -delta*D^(-1)*K MatDiagonalScale(K, x, PETSC_NULL); // K' := -delta*D^(-1)*K
MatMatMult(data->matMinusDeltaK, K, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &S); // S := -delta*K*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 MatAYPX(S, sqr(*data->eps)/(*data->delta), data->matMass, 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 /// create new solver for S
KSP kspS; KSP kspS;
...@@ -79,35 +65,46 @@ namespace AMDiS { ...@@ -79,35 +65,46 @@ namespace AMDiS {
PCShellSetContext(pc, data); PCShellSetContext(pc, data);
} }
KSPSolve(kspS, x1, x2); // S*x2 = x1 KSPSolve(kspS, b, x); // 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(&S);
MatDestroy(&K); MatDestroy(&K);
KSPDestroy(&kspS); KSPDestroy(&kspS);
} }
PetscSolverCahnHilliard2::PetscSolverCahnHilliard2(string name, double *epsPtr, double *deltaPtr) PetscSolverCahnHilliard2::PetscSolverCahnHilliard2(string name)
: PetscSolverGlobalBlockMatrix(name), : PetscSolverGlobalMatrix(name),
useOldInitialGuess(false),
laplaceSolutionMode(0),
massMatrixSolver(NULL), massMatrixSolver(NULL),
laplaceMatrixSolver(NULL), laplaceMatrixSolver(NULL),
deltaKMatrixSolver(NULL), deltaKMatrixSolver(NULL),
eps(epsPtr), phase(NULL)
delta(deltaPtr) { } {
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) void PetscSolverCahnHilliard2::initSolver(KSP &ksp)
...@@ -115,27 +112,63 @@ namespace AMDiS { ...@@ -115,27 +112,63 @@ namespace AMDiS {
FUNCNAME("PetscSolverCahnHilliard2::initSolver()"); FUNCNAME("PetscSolverCahnHilliard2::initSolver()");
// Create FGMRES based outer solver // 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); KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL); 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) void PetscSolverCahnHilliard2::initPreconditioner(PC pc)
{ {
FUNCNAME("PetscSolverCahnHilliard2::initPreconditioner()"); 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); PCSetType(pc, PCFIELDSPLIT);
PCShellSetApply(pc, pcChShell2); PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
PCShellSetContext(pc, &matShellContext); PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
const FiniteElemSpace *feSpace = componentSpaces[0]; createFieldSplit(pc, "ch_potential", chPotentialComponent);
createFieldSplit(pc, "ch_schur", chSchurComponent);
// === matrix M === 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); DOFMatrix massMatrix(feSpace, feSpace);
Operator massOp(feSpace, feSpace); Operator massOp(feSpace, feSpace);
Simple_ZOT zot; Simple_ZOT zot;
...@@ -144,8 +177,7 @@ namespace AMDiS { ...@@ -144,8 +177,7 @@ namespace AMDiS {
massMatrixSolver = createSubSolver(0, "mass_"); massMatrixSolver = createSubSolver(0, "mass_");
massMatrixSolver->fillPetscMatrix(&massMatrix); massMatrixSolver->fillPetscMatrix(&massMatrix);
// === Laplace matrix solver ===
// === matrix (M + eps*sqrt(delta)*K) ===
DOFMatrix laplaceMatrix(feSpace, feSpace); DOFMatrix laplaceMatrix(feSpace, feSpace);
Operator laplaceOp(feSpace, feSpace); Operator laplaceOp(feSpace, feSpace);
laplaceOp.addTerm(&zot); // M laplaceOp.addTerm(&zot); // M
...@@ -166,80 +198,47 @@ namespace AMDiS { ...@@ -166,80 +198,47 @@ namespace AMDiS {
// === Setup solver === // === Setup solver ===
matShellContext.kspMass = massMatrixSolver->getSolver();
matShellContext.kspMinusDeltaK = deltaKMatrixSolver->getSolver();
matShellContext.kspMplusK = laplaceMatrixSolver->getSolver(); matShellContext.kspMplusK = laplaceMatrixSolver->getSolver();
matShellContext.matM = massMatrixSolver->getMatInterior();
matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior(); matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior();
matShellContext.eps = eps; matShellContext.eps = eps;
matShellContext.delta = delta; matShellContext.delta = delta;
matShellContext.kspMass = massMatrixSolver->getSolver();
matShellContext.matMass = massMatrixSolver->getMatInterior();
matShellContext.mpiCommGlobal = &(meshDistributor->getMpiComm(0)); matShellContext.mpiCommGlobal = &(meshDistributor->getMpiComm(0));
petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 2); 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.kspMplusK, "MpK_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
// petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1); petsc_helper::setSolver(kspChPotential, "chPotential", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
// {
// PC pc; MSG("Setup of Cahn-Hilliard preconditioner needed %.5f seconds\n",
// KSPGetPC(matShellContext.kspLaplace, &pc); MPI::Wtime() - wtime);
// 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) void PetscSolverCahnHilliard2::exitPreconditioner(PC pc)
{ {
FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()"); FUNCNAME("PetscSolverCahnHilliard2::exitPreconditioner()");
massMatrixSolver->destroyMatrixData(); massMatrixSolver->destroyMatrixData();
massMatrixSolver->destroyVectorData();
laplaceMatrixSolver->destroyMatrixData(); laplaceMatrixSolver->destroyMatrixData();
laplaceMatrixSolver->destroyVectorData();
deltaKMatrixSolver->destroyMatrixData(); deltaKMatrixSolver->destroyMatrixData();
deltaKMatrixSolver->destroyVectorData();
massMatrixSolver->destroyVectorData(); massMatrixSolver->destroyVectorData();
laplaceMatrixSolver->destroyVectorData(); laplaceMatrixSolver->destroyVectorData();
deltaKMatrixSolver->destroyVectorData(); deltaKMatrixSolver->destroyVectorData();
delete massMatrixSolver; delete massMatrixSolver;
massMatrixSolver = NULL; massMatrixSolver = NULL;
delete laplaceMatrixSolver; delete laplaceMatrixSolver;
laplaceMatrixSolver = NULL; laplaceMatrixSolver = NULL;
delete deltaKMatrixSolver; delete deltaKMatrixSolver;
deltaKMatrixSolver = NULL; deltaKMatrixSolver = NULL;
}
}
} }
...@@ -23,7 +23,7 @@ ...@@ -23,7 +23,7 @@
#ifndef AMDIS_PETSC_SOLVER_CAHN_HILLIARD2_H #ifndef AMDIS_PETSC_SOLVER_CAHN_HILLIARD2_H
#define AMDIS_PETSC_SOLVER_CAHN_HILLIARD2_H #define AMDIS_PETSC_SOLVER_CAHN_HILLIARD2_H
#include "parallel/PetscSolverGlobalBlockMatrix.h" #include "parallel/PetscSolverGlobalMatrix.h"
namespace AMDiS { namespace AMDiS {
...@@ -31,12 +31,12 @@ namespace AMDiS { ...@@ -31,12 +31,12 @@ namespace AMDiS {
struct CahnHilliardData2 { struct CahnHilliardData2 {
KSP kspMass, kspMinusDeltaK, kspMplusK; KSP kspMass, kspMinusDeltaK, kspMplusK;
Mat matM, matMinusDeltaK; Mat matMass, matMinusDeltaK;
double *eps, *delta; double *eps, *delta;
MPI::Intracomm *mpiCommGlobal; MPI::Intracomm *mpiCommGlobal;
}; };
class PetscSolverCahnHilliard2 : public PetscSolverGlobalBlockMatrix class PetscSolverCahnHilliard2 : public PetscSolverGlobalMatrix
{ {
public: public:
/// Creator class /// Creator class
...@@ -51,33 +51,63 @@ namespace AMDiS { ...@@ -51,33 +51,63 @@ namespace AMDiS {
return new PetscSolverCahnHilliard2(this->name); 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) void setChData(double *epsPtr, double *deltaPtr)
{ {
eps = epsPtr; eps = epsPtr;
delta = deltaPtr; delta = deltaPtr;
} }
void setPhase(DOFVector<double> *d, double eP3=0)
{
phase = d;
epsPhase3 = eP3;
}
protected: protected:
void initSolver(KSP &ksp); void initSolver(KSP &ksp);
void initPreconditioner(PC pc); void initPreconditioner(PC pc);
void exitPreconditioner(PC pc); void exitPreconditioner(PC pc);
PetscSolver* createSubSolver(int component, string kspPrefix);
private: private:
int pressureComponent; 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; PetscSolver *massMatrixSolver, *laplaceMatrixSolver, *deltaKMatrixSolver;
CahnHilliardData2 matShellContext; CahnHilliardData2 matShellContext;
double *eps, *delta; double *eps, *delta;
double epsPhase3;
SystemVector* solution;
DOFVector<double>* phase;
}; };
} }
#endif // AMDIS_PETSC_SOLVER_CAHN_HILLIARD2_H #endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment