From 2f5be996fced70b0fa8d2e41f67e45c0204ce617 Mon Sep 17 00:00:00 2001
From: Sebastian Aland <sebastian.aland@yahoo.com>
Date: Tue, 5 Feb 2013 12:32:40 +0000
Subject: [PATCH] PetscSolverCahnHilliard geht jetzt auch mit diffuse domain
---
AMDiS/src/parallel/PetscSolverCahnHilliard.cc | 272 +++++++++++-------
AMDiS/src/parallel/PetscSolverCahnHilliard.h | 5 +-
2 files changed, 165 insertions(+), 112 deletions(-)
diff --git a/AMDiS/src/parallel/PetscSolverCahnHilliard.cc b/AMDiS/src/parallel/PetscSolverCahnHilliard.cc
index 1ee3f22c..3b3e0c80 100644
--- a/AMDiS/src/parallel/PetscSolverCahnHilliard.cc
+++ b/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;
}
-
-
+
+
}
diff --git a/AMDiS/src/parallel/PetscSolverCahnHilliard.h b/AMDiS/src/parallel/PetscSolverCahnHilliard.h
index 9ff119c7..2ebcf8be 100644
--- a/AMDiS/src/parallel/PetscSolverCahnHilliard.h
+++ b/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;
};
--
GitLab