Skip to content
GitLab
Commit 014e5191 authored by Sebastian Aland's avatar Sebastian Aland
Browse files

PetscSolverNSCH added

parent 02638a26
Hide whitespace changes
Inline Side-by-side
AMDiS/src/parallel/PetscSolverNSCH.cc 0 → 100644
View file @ 014e5191
//
// 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 "parallel/PetscSolverNSCH.h"
#include "parallel/PetscHelper.h"
#include "TransformDOF.h"
namespace AMDiS {
using namespace std;
PetscErrorCode pcShell(PC pc, Vec b, Vec x) // solve Px=b
{FUNCNAME("PCApply()");
void *ctx;
PCShellGetContext(pc, &ctx);
CahnHilliardData2* data = static_cast<CahnHilliardData2*>(ctx);
/// extract vectors
Vec b1, b2, b34, b5, x1, x2, x34, x5, b3,b4;
data->globalMatrixSolver->extractVectorComponent(b, 0+3, &b1);
data->globalMatrixSolver->extractVectorComponent(b, 1+3, &b2);
data->globalMatrixSolver->extractVectorComponent(b, 0, &b34, 2);
data->globalMatrixSolver->extractVectorComponent(b, 2, &b5);
data->globalMatrixSolver->extractVectorComponent(x, 0+3, &x1);
data->globalMatrixSolver->extractVectorComponent(x, 1+3, &x2);
data->globalMatrixSolver->extractVectorComponent(x, 0, &x34, 2);
data->globalMatrixSolver->extractVectorComponent(x, 2, &x5);
data->globalMatrixSolver->extractVectorComponent(b, 0, &b3, 1);
data->globalMatrixSolver->extractVectorComponent(b, 1, &b4, 1);
/// solve Cahn-Hilliard Preconditioner
Vec y1, y2;
VecDuplicate(b1, &y1);
VecDuplicate(b2, &y2);
KSPSolve(data->kspMassCH, b1, y1); // M*y1 = b1
MatMultAdd(data->matMinusDeltaK, y1, b2, x1); // -> x1 := b2-delta*K*y1
KSPSolve(data->kspLaplaceCH, x1, y2); // (M+eps*sqrt(delta))*y2 = x1
MatMult(data->matMassCH, y2, x1); // x1 := M*y2
KSPSolve(data->kspLaplaceCH, 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);
/**/
/// solve Navier-Stokes Preconditioner
Vec tmp34, tmp5, tmp5_2, tmp34_2;
VecDuplicate(b34, &tmp34);
VecDuplicate(b34, &tmp34_2);
VecDuplicate(b5, &tmp5);
VecDuplicate(b5, &tmp5_2);
KSPSolve(data->kspVelocity, b34, tmp34);
VecScale(tmp34, -1.0);
MatMultAdd(data->matDiv, tmp34, b5, tmp5);
/// approximierte Schur Komplement
KSPSolve(data->kspLaplace, tmp5, x5);
{ //project out constant Null-space
int vecSize;
VecGetSize(x5, &vecSize);
PetscScalar vecSum;
VecSum(x5, &vecSum);
vecSum = vecSum / static_cast<PetscScalar>(-1.0 * vecSize);
VecShift(x5, vecSum);
//VecView(y, PETSC_VIEWER_STDOUT_WORLD);
}
MatMult(data->matConDif, x5, tmp5);
KSPSolve(data->kspMass, tmp5, x5);
VecScale(x5,-1.0);
MatMultAdd(data->matGrad, x5, b34, tmp34);
KSPSolve(data->kspVelocity, tmp34, x34);
VecScale(x5,-1.0);
/**/
VecDestroy(&tmp34);
VecDestroy(&tmp5);
VecDestroy(&tmp5_2);
VecDestroy(&tmp34_2);
VecDestroy(&b1);
VecDestroy(&b2);
VecDestroy(&b34);
VecDestroy(&b5);
VecDestroy(&x1);
VecDestroy(&x2);
VecDestroy(&x34);
VecDestroy(&x5);
}
PetscSolverNSCH::PetscSolverNSCH(string name)
: PetscSolverGlobalMatrix(name),
useOldInitialGuess(false),
laplaceSolutionMode(0),
massMatrixSolverCH(NULL),
laplaceMatrixSolverCH(NULL),
deltaKMatrixSolver(NULL),
pressureNullSpace(true),
velocitySolutionMode(0),
regularizeLaplace(0),
massSolutionMode(0),
massMatrixSolver(NULL),
laplaceMatrixSolver(NULL),
nu(NULL),
invTau(NULL),
solution(NULL),
phase(NULL)
{
Parameters::get(initFileStr + "->use old initial guess",
useOldInitialGuess);
Parameters::get(initFileStr + "->navierstokes->velocity solver",
velocitySolutionMode);
Parameters::get(initFileStr + "->navierstokes->mass solver",
massSolutionMode);
Parameters::get(initFileStr + "->navierstokes->laplace solver",
laplaceSolutionMode);
Parameters::get(initFileStr + "->navierstokes->regularize laplace",
regularizeLaplace);
}
void PetscSolverNSCH::solvePetscMatrix(SystemVector &vec,
AdaptInfo *adaptInfo)
{
FUNCNAME("PetscSolverNSCH::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 PetscSolverNSCH::initSolver(KSP &ksp)
{
FUNCNAME("PetscSolverNSCH::initSolver()");
// Create FGMRES based outer solver
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, PCSHELL, 1e-6, 1e-8, 10000);
setConstantNullSpace(ksp, componentSpaces[0]->getMesh()->getDim() , true);
}
void PetscSolverNSCH::initPreconditioner(PC pc)
{
FUNCNAME("PetscSolverNSCH::initPreconditioner()");
MPI::COMM_WORLD.Barrier();
double wtime = MPI::Wtime();
int dim = componentSpaces[0]->getMesh()->getDim();
pressureComponent=dim;
const FiniteElemSpace *cahnHilliardFeSpace = componentSpaces[0+3];
const FiniteElemSpace *velocityFeSpace= componentSpaces[0];
const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
PCSetType(pc, PCSHELL);
PCShellSetApply(pc, pcShell);
PCShellSetContext(pc, &matShellContext);
matShellContext.globalMatrixSolver = (this);
matShellContext.mpiCommGlobal= &(meshDistributor->getMpiComm(0));
/// Init Cahn-Hilliard part
DOFMatrix laplaceMatrixCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
Operator laplaceOpCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
DOFMatrix massMatrixCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
Operator massOpCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
DOFMatrix deltaKMatrix(cahnHilliardFeSpace, cahnHilliardFeSpace);
Operator laplaceOp2(cahnHilliardFeSpace, cahnHilliardFeSpace);
{
Simple_ZOT zot;
massOpCH.addTerm(&zot);
laplaceOpCH.addTerm(&zot); // M
Simple_SOT sot2((*eps)*sqrt(*delta));
laplaceOpCH.addTerm(&sot2); // eps*sqrt(delta)*K
Simple_SOT sot(-(*delta));
laplaceOp2.addTerm(&sot); // -delta*K
massMatrixCH.assembleOperator(massOpCH);
massMatrixSolverCH = createSubSolver(0+3, "mass_");
massMatrixSolverCH->fillPetscMatrix(&massMatrixCH);
// === matrix (M + eps*sqrt(delta)*K) ===
laplaceMatrixCH.assembleOperator(laplaceOpCH);
laplaceMatrixSolverCH = createSubSolver(0+3, "laplace_");
laplaceMatrixSolverCH->fillPetscMatrix(&laplaceMatrixCH);
// === matrix (-delta*K) ===
deltaKMatrix.assembleOperator(laplaceOp2);
deltaKMatrixSolver = createSubSolver(0+3, "laplace2_");
deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
}
// === Setup solver ===
matShellContext.kspMassCH = massMatrixSolverCH->getSolver();
matShellContext.kspLaplaceCH = laplaceMatrixSolverCH->getSolver();
matShellContext.matMassCH = massMatrixSolverCH->getMatInterior();
matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior();
matShellContext.eps = eps;
matShellContext.delta = delta;
petsc_helper::setSolver(matShellContext.kspMassCH, "", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
petsc_helper::setSolver(matShellContext.kspLaplaceCH, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
/// Init Navier-Stokes part
DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
Operator massOp(pressureFeSpace, pressureFeSpace);
ZeroOrderTerm *massTerm = new Simple_ZOT;
massOp.addTerm(massTerm);
massMatrix.assembleOperator(massOp);
delete massTerm;
massMatrixSolver = createSubSolver(pressureComponent, "mass_");
massMatrixSolver->fillPetscMatrix(&massMatrix);
matShellContext.kspMass = massMatrixSolver->getSolver();
DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
Operator laplaceOp(pressureFeSpace, pressureFeSpace);
SecondOrderTerm *laplaceTerm = new Simple_SOT;
laplaceOp.addTerm(laplaceTerm);
laplaceMatrix.assembleOperator(laplaceOp);
delete laplaceTerm;
laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_"));
laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
// === Create convection-diffusion operator ===
DOFVector<double> vx(pressureFeSpace, "vx");
DOFVector<double> vy(pressureFeSpace, "vy");
DOFVector<double> vz(pressureFeSpace, "vz");
DOFVector<double> vp(pressureFeSpace, "vp");
vx.interpol(solution->getDOFVector(0));
vy.interpol(solution->getDOFVector(1));
if (dim == 3) vz.interpol(solution->getDOFVector(2));
DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
{
Operator conDifOp(pressureFeSpace, pressureFeSpace);
ZeroOrderTerm *conDif0 = NULL;
SecondOrderTerm *conDif1 = NULL;
FirstOrderTerm *conDif2 = NULL, *conDif3 = NULL, *conDif4 = NULL;
vp.interpol(solution->getDOFVector(dim+2));
densityFunctionTau = new LinearInterpolation(*rho1,*rho2,*invTau);
viscosityFunction = new LinearInterpolation(*nu1,*nu2);
densityFunction = new LinearInterpolation2(*rho1,*rho2);
conDif0 = new VecAtQP_ZOT(&vp, densityFunctionTau);
conDifOp.addTerm(conDif0);
conDif1 = new VecAtQP_SOT(&vp, viscosityFunction );
conDifOp.addTerm(conDif1);
conDif2 = new Vec2AtQP_FOT(&vx, &vp, densityFunction, 0);
conDifOp.addTerm(conDif2, GRD_PHI);
conDif3 = new Vec2AtQP_FOT(&vy, &vp, densityFunction, 1);
conDifOp.addTerm(conDif3, GRD_PHI);
if (dim == 3) {
conDif4 = new Vec2AtQP_FOT(&vz, &vp, densityFunction, 2);
conDifOp.addTerm(conDif4, GRD_PHI);
}
conDifMatrix.assembleOperator(conDifOp);
delete conDif0;
delete conDif1;
delete conDif2;
delete conDif3;
if (dim == 3) delete conDif4;
}
conDifMatrixSolver = createSubSolver(pressureComponent, "condif_");
conDifMatrixSolver->fillPetscMatrix(&conDifMatrix);
matShellContext.matConDif = conDifMatrixSolver->getMatInterior();
extractMatrixComponent(mat[0][0], 2, 1, 0, 2, &(matShellContext.matDiv));
extractMatrixComponent(mat[0][0], 0, 2, 2, 1, &(matShellContext.matGrad));
extractMatrixComponent(mat[0][0], 0, 2, 0, 2, &(matShellContext.velocityMat));
///erstelle kspVelocity
KSPCreate((meshDistributor->getMpiComm(0)), &(matShellContext.kspVelocity));
KSPSetOperators(matShellContext.kspVelocity, matShellContext.velocityMat, matShellContext.velocityMat, SAME_NONZERO_PATTERN);
///regularisiere LaplaceMatrix
if (regularizeLaplace)
{
PetscInt rows[1];
rows[0]=0;
MatZeroRows(laplaceMatrixSolver->getMatInterior(), 1, rows, 0, PETSC_NULL, PETSC_NULL);
KSPCreate((meshDistributor->getMpiComm(0)), &(matShellContext.kspLaplace));
KSPSetOperators(matShellContext.kspLaplace, laplaceMatrixSolver->getMatInterior(), laplaceMatrixSolver->getMatInterior(), SAME_NONZERO_PATTERN);
}
else
{ matShellContext.kspLaplace=laplaceMatrixSolver->getSolver();
setConstantNullSpace(matShellContext.kspLaplace);
}
// === Setup solver ===
switch (massSolutionMode) {
case 0:
petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
break;
case 1:
petsc_helper::setSolverWithLu(matShellContext.kspMass, "mass_", KSPRICHARDSON, PCLU, MATSOLVERMUMPS, 0.0, 1e-14, 1);
break;
default:
ERROR_EXIT("No mass solution mode %d available!\n", massSolutionMode);
}
switch (laplaceSolutionMode) {
case 0:
petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
break;
case 1:
petsc_helper::setSolverWithLu(matShellContext.kspLaplace, "mass_", KSPRICHARDSON, PCLU, MATSOLVERMUMPS, 0.0, 1e-14, 1);
break;
default:
ERROR_EXIT("No laplace solution mode %d available!\n", laplaceSolutionMode);
}
switch (velocitySolutionMode) {
case 0:
petsc_helper::setSolver(matShellContext.kspVelocity, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
break;
case 1:
petsc_helper::setSolverWithLu(matShellContext.kspVelocity, "", KSPPREONLY, PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1);
break;
default:
ERROR_EXIT("No velocity solution mode %d available!\n", velocitySolutionMode);
}
PCSetFromOptions(pc);
MSG("Setup of Cahn-Hilliard preconditioner needed %.5f seconds\n",
MPI::Wtime() - wtime);
}
void PetscSolverNSCH::exitPreconditioner(PC pc)
{
FUNCNAME("PetscSolverNSCH::exitPreconditioner()");
massMatrixSolverCH->destroyMatrixData();
massMatrixSolverCH->destroyVectorData();
laplaceMatrixSolverCH->destroyMatrixData();
laplaceMatrixSolverCH->destroyVectorData();
deltaKMatrixSolver->destroyMatrixData();
deltaKMatrixSolver->destroyVectorData();
massMatrixSolverCH->destroyVectorData();
laplaceMatrixSolverCH->destroyVectorData();
deltaKMatrixSolver->destroyVectorData();
delete massMatrixSolverCH;
massMatrixSolverCH = NULL;
delete laplaceMatrixSolverCH;
laplaceMatrixSolverCH = NULL;
delete deltaKMatrixSolver;
deltaKMatrixSolver = NULL;
massMatrixSolver->destroyMatrixData();
massMatrixSolver->destroyVectorData();
laplaceMatrixSolver->destroyMatrixData();
laplaceMatrixSolver->destroyVectorData();
conDifMatrixSolver->destroyMatrixData();
conDifMatrixSolver->destroyVectorData();
massMatrixSolver->destroyVectorData();
laplaceMatrixSolver->destroyVectorData();
conDifMatrixSolver->destroyVectorData();
delete massMatrixSolver;
massMatrixSolver = NULL;
delete laplaceMatrixSolver;
laplaceMatrixSolver = NULL;
delete conDifMatrixSolver;
conDifMatrixSolver = NULL;
KSPDestroy(&(matShellContext.kspVelocity));
if (regularizeLaplace)
KSPDestroy(&(matShellContext.kspLaplace));
MatDestroy(&(matShellContext.matGrad));
MatDestroy(&(matShellContext.matDiv));
MatDestroy(&(matShellContext.velocityMat));
delete densityFunction;
delete densityFunctionTau;
delete viscosityFunction;
}
}
AMDiS/src/parallel/PetscSolverNSCH.h 0 → 100644
View file @ 014e5191
// ============================================================================
// == ==
// == 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 PetscSolverNSCH.h */
#include "parallel/PetscSolverGlobalMatrix.h"
namespace AMDiS {
using namespace std;
struct CahnHilliardData2 {
KSP kspMassCH, kspLaplaceCH, kspVelocity, kspLaplace, kspMass;
Mat matMassCH, matMinusDeltaK, matGrad, matDiv, matConDif, matSchur, velocityMat;
PetscSolverGlobalMatrix *globalMatrixSolver;
double *eps, *delta;
MPI::Intracomm *mpiCommGlobal;
};
class PetscSolverNSCH : public PetscSolverGlobalMatrix
{
public:
/// Creator class
class IdFct : public AbstractFunction<double, double>
{
public:
IdFct()
: AbstractFunction<double, double>(0)
{}
double operator()(const double& x) const
{
return x;
}
};
class MultConstFct : public AbstractFunction<double, double>
{
public:
MultConstFct(double c)
: AbstractFunction<double, double>(1),
mConst(c)
{}
double operator()(const double& x) const
{
return mConst * x;
}
private:
double mConst;
};
class LinearInterpolation : public AbstractFunction<double, double>
{
public:
LinearInterpolation(double c1, double c2, double factor=1.0)
: AbstractFunction<double, double>(0)
{ a = (c1-c2)/2.0*factor; b = (c1+c2)/2.0*factor; cmin=std::min(c1,c2)*factor; cmax=std::max(c1,c2)*factor;}
double operator()(const double& x) const
{ double result = b+a*x;
if (result<cmin) result = cmin;
if (result>cmax) result = cmax;
return result;
}
private:
double a,b,cmin,cmax;
};
class LinearInterpolation2 : public BinaryAbstractFunction<double, double, double>
{
public:
LinearInterpolation2(double c1, double c2, double factor=1.0)
: BinaryAbstractFunction<double, double, double>(0)
{ a = (c1-c2)/2.0*factor; b = (c1+c2)/2.0*factor; cmin=std::min(c1,c2)*factor; cmax=std::max(c1,c2)*factor;}
double operator()(const double& u, const double& x) const
{
double result = b+a*x;
if (result<cmin) result = cmin;
if (result>cmax) result = cmax;
return result * u;
}
private:
double a,b,cmin,cmax;
};
class Multiplier3 : public BinaryAbstractFunction<double, double, double>
{
public:
Multiplier3()
: BinaryAbstractFunction<double, double, double >(0)
{}
double operator()(const double& phi, const double& phase) const
{
return phase * phi;
}
};
class EinsMinus : public AbstractFunction<double, double>
{
public:
EinsMinus(double d)
: AbstractFunction<double, double>(0),
c(d)
{}
double operator()(const double& x) const
{
return c * std::max(1.0-x,0.000001);
}
private:
double c;
};
struct Multiplication : public BinaryAbstractFunction<double, double, double>{
double operator()(const double &v1, const double &v2) const{
return v2 / v1;
}
};
class Creator : public OEMSolverCreator
{
public:
virtual ~Creator() {}
/// Returns a new PetscSolver object.
OEMSolver* create()