Commit a41536c3 authored by Thomas Witkowski's avatar Thomas Witkowski

Merge to juropa

parent 0ca6954f
...@@ -18,6 +18,51 @@ namespace AMDiS { ...@@ -18,6 +18,51 @@ namespace AMDiS {
using namespace std; using namespace std;
PetscErrorCode pcNs(PC pc, Vec x, Vec y)
{
void *ctx;
PCShellGetContext(pc, &ctx);
MatShellNavierStokesSchur* data = static_cast<MatShellNavierStokesSchur*>(ctx);
Vec velocity, pressure;
VecGetSubVector(x, data->isVelocity, &velocity);
VecGetSubVector(x, data->isPressure, &pressure);
Vec velocityTmp;
VecDuplicate(velocity, &velocityTmp);
Vec pressureTmp;
VecDuplicate(pressure, &pressureTmp);
#if 0
MatMult(data->A01, pressure, velocityTmp);
VecAXPY(velocity, 1.0, velocityTmp);
KSPSolve(data->kspVelocity, velocity, velocity);
VecScale(pressure, -1.0);
#else
KSPSolve(data->kspVelocity, velocity, velocityTmp);
MatMult(data->A10, velocityTmp, pressureTmp);
VecAXPY(pressure, -1.0, pressureTmp);
KSPSolve(data->kspLaplace, pressure, pressure);
MatMult(data->matConDif, pressure, pressureTmp);
KSPSolve(data->kspMass, pressureTmp, pressure);
MatMult(data->A01, pressure, velocityTmp);
VecAXPY(velocity, -1.0, velocityTmp);
KSPSolve(data->kspVelocity, velocity, velocity);
#endif
VecRestoreSubVector(x, data->isVelocity, &velocity);
VecRestoreSubVector(x, data->isPressure, &pressure);
VecCopy(x, y);
VecDestroy(&velocityTmp);
VecDestroy(&pressureTmp);
}
int petscMultNavierStokesSchur(Mat mat, Vec x, Vec y) int petscMultNavierStokesSchur(Mat mat, Vec x, Vec y)
{ {
void *ctx; void *ctx;
...@@ -47,9 +92,10 @@ namespace AMDiS { ...@@ -47,9 +92,10 @@ namespace AMDiS {
VecDuplicate(x, &vec0); VecDuplicate(x, &vec0);
VecDuplicate(x, &vec1); VecDuplicate(x, &vec1);
KSPSolve(data->kspLaplace, x, vec0); // KSPSolve(data->kspLaplace, x, vec0);
MatMult(data->matConDif, vec0, vec1); // MatMult(data->matConDif, vec0, vec1);
KSPSolve(data->kspMass, vec1, y); //KSPSolve(data->kspMass, vec0, y);
VecSet(y, 0.0);
VecDestroy(&vec0); VecDestroy(&vec0);
VecDestroy(&vec1); VecDestroy(&vec1);
...@@ -114,57 +160,31 @@ namespace AMDiS { ...@@ -114,57 +160,31 @@ namespace AMDiS {
{ {
FUNCNAME("PetscSolverNavierStokes::initPreconditioner()"); FUNCNAME("PetscSolverNavierStokes::initPreconditioner()");
vector<int> velocityComponents;
velocityComponents.push_back(0);
velocityComponents.push_back(1);
PCSetType(pc, PCFIELDSPLIT);
PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
createFieldSplit(pc, "velocity", velocityComponents);
createFieldSplit(pc, "pressure", pressureComponent);
KSPSetUp(kspInterior);
KSP *subKsp;
int nSubKsp;
PCFieldSplitGetSubKSP(pc, &nSubKsp, &subKsp);
TEST_EXIT(nSubKsp == 2)
("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");
KSP kspVelocity = subKsp[0];
KSP kspSchur = subKsp[1];
PetscFree(subKsp);
Mat A00, A01, A10, A11; PCSetType(pc, PCSHELL);
PCFieldSplitGetSchurBlocks(pc, &A00, &A01, &A10, &A11); PCShellSetApply(pc, pcNs);
matShellContext.A00 = A00; PCShellSetContext(pc, &matShellContext);
matShellContext.A01 = A01;
matShellContext.A10 = A10; interiorMap->createIndexSet(matShellContext.isVelocity, 0, 2);
matShellContext.A11 = A11; interiorMap->createIndexSet(matShellContext.isPressure, 2, 1);
matShellContext.kspVelocity = kspVelocity; MatGetSubMatrix(getMatInterior(), matShellContext.isVelocity,
matShellContext.isVelocity, MAT_INITIAL_MATRIX, &matShellContext.A00);
Mat matShell; MatGetSubMatrix(getMatInterior(), matShellContext.isVelocity,
MatDuplicate(A11, MAT_DO_NOT_COPY_VALUES, &matShell); matShellContext.isPressure, MAT_INITIAL_MATRIX, &matShellContext.A01);
MatSetType(matShell, MATSHELL); MatGetSubMatrix(getMatInterior(), matShellContext.isPressure,
MatShellSetContext(matShell, &matShellContext); matShellContext.isVelocity, MAT_INITIAL_MATRIX, &matShellContext.A10);
MatShellSetOperation(matShell, MATOP_MULT, (void(*)(void))petscMultNavierStokesSchur); MatGetSubMatrix(getMatInterior(), matShellContext.isPressure,
MatSetUp(matShell); matShellContext.isPressure, MAT_INITIAL_MATRIX, &matShellContext.A11);
KSPSetType(kspVelocity, KSPPREONLY); KSPCreate(mpiCommGlobal, &(matShellContext.kspVelocity));
KSPSetOperators(matShellContext.kspVelocity, matShellContext.A00, matShellContext.A00, SAME_NONZERO_PATTERN);
// KSPSetType(matShellContext.kspVelocity, KSPPREONLY);
KSPSetInitialGuessNonzero(matShellContext.kspVelocity, PETSC_TRUE);
PC pcSub; PC pcSub;
KSPGetPC(kspVelocity, &pcSub); KSPGetPC(matShellContext.kspVelocity, &pcSub);
PCSetType(pcSub, PCLU); PCSetType(pcSub, PCLU);
PCFactorSetMatSolverPackage(pcSub, MATSOLVERMUMPS); PCFactorSetMatSolverPackage(pcSub, MATSOLVERMUMPS);
KSPSetType(kspSchur, KSPGMRES);
KSPSetTolerances(kspSchur, 0, 1e-14, 1e+3, 1000);
KSPSetOperators(kspSchur, matShell, matShell, SAME_NONZERO_PATTERN);
KSPGetPC(kspSchur, &pcSub);
PCSetType(pcSub, PCNONE);
// === Mass matrix solver === // === Mass matrix solver ===
...@@ -217,9 +237,9 @@ namespace AMDiS { ...@@ -217,9 +237,9 @@ namespace AMDiS {
Simple_SOT conDif1(nu); Simple_SOT conDif1(nu);
conDifOp.addTerm(&conDif1); conDifOp.addTerm(&conDif1);
Vector_FOT conDif2(0); Vector_FOT conDif2(0);
// conDifOp.addTerm(&conDif2, GRD_PHI); conDifOp.addTerm(&conDif2, GRD_PHI);
Vector_FOT conDif3(1); Vector_FOT conDif3(1);
// conDifOp.addTerm(&conDif3, GRD_PHI); conDifOp.addTerm(&conDif3, GRD_PHI);
conDifMatrix.assembleOperator(conDifOp); conDifMatrix.assembleOperator(conDifOp);
...@@ -232,6 +252,7 @@ namespace AMDiS { ...@@ -232,6 +252,7 @@ namespace AMDiS {
matShellContext.kspLaplace = laplaceMatrixSolver->getSolver(); matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
matShellContext.matConDif = conDifMatrixSolver->getMatInterior(); matShellContext.matConDif = conDifMatrixSolver->getMatInterior();
KSPSetType(matShellContext.kspMass, KSPRICHARDSON); KSPSetType(matShellContext.kspMass, KSPRICHARDSON);
KSPSetTolerances(matShellContext.kspMass, 0, 1e-13, 1e+3, 1); KSPSetTolerances(matShellContext.kspMass, 0, 1e-13, 1e+3, 1);
KSPGetPC(matShellContext.kspMass, &pcSub); KSPGetPC(matShellContext.kspMass, &pcSub);
...@@ -240,16 +261,72 @@ namespace AMDiS { ...@@ -240,16 +261,72 @@ namespace AMDiS {
KSPSetFromOptions(matShellContext.kspMass); KSPSetFromOptions(matShellContext.kspMass);
KSPSetType(matShellContext.kspLaplace, KSPCG); KSPSetType(matShellContext.kspLaplace, KSPCG);
KSPSetInitialGuessNonzero(matShellContext.kspLaplace, PETSC_TRUE);
KSPSetTolerances(matShellContext.kspLaplace, 0, 1e-13, 1e+3, 10000); KSPSetTolerances(matShellContext.kspLaplace, 0, 1e-13, 1e+3, 10000);
KSPGetPC(matShellContext.kspLaplace, &pcSub); KSPGetPC(matShellContext.kspLaplace, &pcSub);
PCSetType(pcSub, PCJACOBI); PCSetType(pcSub, PCJACOBI);
KSPSetFromOptions(matShellContext.kspLaplace); KSPSetFromOptions(matShellContext.kspLaplace);
#if 0
vector<int> velocityComponents;
velocityComponents.push_back(0);
velocityComponents.push_back(1);
PCSetType(pc, PCFIELDSPLIT);
PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
createFieldSplit(pc, "velocity", velocityComponents);
createFieldSplit(pc, "pressure", pressureComponent);
KSPSetUp(kspInterior);
KSP *subKsp;
int nSubKsp;
PCFieldSplitGetSubKSP(pc, &nSubKsp, &subKsp);
TEST_EXIT(nSubKsp == 2)
("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");
KSP kspVelocity = subKsp[0];
KSP kspSchur = subKsp[1];
PetscFree(subKsp);
Mat A00, A01, A10, A11;
PCFieldSplitGetSchurBlocks(pc, &A00, &A01, &A10, &A11);
matShellContext.A00 = A00;
matShellContext.A01 = A01;
matShellContext.A10 = A10;
matShellContext.A11 = A11;
matShellContext.kspVelocity = kspVelocity;
Mat matShell;
MatDuplicate(A11, MAT_DO_NOT_COPY_VALUES, &matShell);
MatSetType(matShell, MATSHELL);
MatShellSetContext(matShell, &matShellContext);
MatShellSetOperation(matShell, MATOP_MULT, (void(*)(void))petscMultNavierStokesSchur);
KSPSetType(kspVelocity, KSPPREONLY);
PC pcSub;
KSPGetPC(kspVelocity, &pcSub);
PCSetType(pcSub, PCLU);
PCFactorSetMatSolverPackage(pcSub, MATSOLVERMUMPS);
KSPSetType(kspSchur, KSPGMRES);
KSPSetTolerances(kspSchur, 0, 1e-14, 1e+3, 1000);
KSPSetOperators(kspSchur, matShell, matShell, SAME_NONZERO_PATTERN);
KSPGetPC(kspSchur, &pcSub);
PCSetType(pcSub, PCNONE);
matShellContext.solverMode = 1; matShellContext.solverMode = 1;
KSPSetFromOptions(kspVelocity); KSPSetFromOptions(kspVelocity);
KSPSetFromOptions(kspSchur); KSPSetFromOptions(kspSchur);
#endif
} }
} }
...@@ -44,6 +44,10 @@ namespace AMDiS { ...@@ -44,6 +44,10 @@ namespace AMDiS {
KSP kspMass, kspLaplace; KSP kspMass, kspLaplace;
Mat matConDif; Mat matConDif;
IS isVelocity, isPressure;
}; };
class PetscSolverNavierStokes : public PetscSolverGlobalMatrix class PetscSolverNavierStokes : public PetscSolverGlobalMatrix
......
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