From a41536c3986d075d89c21f0ebe1e509567c39ca7 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 26 Oct 2012 07:30:37 +0000
Subject: [PATCH] Merge to juropa
---
AMDiS/src/parallel/PetscSolverNavierStokes.cc | 181 +++++++++++++-----
AMDiS/src/parallel/PetscSolverNavierStokes.h | 4 +
2 files changed, 133 insertions(+), 52 deletions(-)
diff --git a/AMDiS/src/parallel/PetscSolverNavierStokes.cc b/AMDiS/src/parallel/PetscSolverNavierStokes.cc
index 6b2cf099..070c7f47 100644
--- a/AMDiS/src/parallel/PetscSolverNavierStokes.cc
+++ b/AMDiS/src/parallel/PetscSolverNavierStokes.cc
@@ -18,6 +18,51 @@ namespace AMDiS {
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)
{
void *ctx;
@@ -47,9 +92,10 @@ namespace AMDiS {
VecDuplicate(x, &vec0);
VecDuplicate(x, &vec1);
- KSPSolve(data->kspLaplace, x, vec0);
- MatMult(data->matConDif, vec0, vec1);
- KSPSolve(data->kspMass, vec1, y);
+ // KSPSolve(data->kspLaplace, x, vec0);
+ // MatMult(data->matConDif, vec0, vec1);
+ //KSPSolve(data->kspMass, vec0, y);
+ VecSet(y, 0.0);
VecDestroy(&vec0);
VecDestroy(&vec1);
@@ -114,57 +160,31 @@ namespace AMDiS {
{
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;
- 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);
- MatSetUp(matShell);
-
- KSPSetType(kspVelocity, KSPPREONLY);
+ PCSetType(pc, PCSHELL);
+ PCShellSetApply(pc, pcNs);
+ PCShellSetContext(pc, &matShellContext);
+
+ interiorMap->createIndexSet(matShellContext.isVelocity, 0, 2);
+ interiorMap->createIndexSet(matShellContext.isPressure, 2, 1);
+ MatGetSubMatrix(getMatInterior(), matShellContext.isVelocity,
+ matShellContext.isVelocity, MAT_INITIAL_MATRIX, &matShellContext.A00);
+ MatGetSubMatrix(getMatInterior(), matShellContext.isVelocity,
+ matShellContext.isPressure, MAT_INITIAL_MATRIX, &matShellContext.A01);
+ MatGetSubMatrix(getMatInterior(), matShellContext.isPressure,
+ matShellContext.isVelocity, MAT_INITIAL_MATRIX, &matShellContext.A10);
+ MatGetSubMatrix(getMatInterior(), matShellContext.isPressure,
+ matShellContext.isPressure, MAT_INITIAL_MATRIX, &matShellContext.A11);
+
+ 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;
- KSPGetPC(kspVelocity, &pcSub);
+ KSPGetPC(matShellContext.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);
-
+
// === Mass matrix solver ===
@@ -217,9 +237,9 @@ namespace AMDiS {
Simple_SOT conDif1(nu);
conDifOp.addTerm(&conDif1);
Vector_FOT conDif2(0);
- // conDifOp.addTerm(&conDif2, GRD_PHI);
+ conDifOp.addTerm(&conDif2, GRD_PHI);
Vector_FOT conDif3(1);
- // conDifOp.addTerm(&conDif3, GRD_PHI);
+ conDifOp.addTerm(&conDif3, GRD_PHI);
conDifMatrix.assembleOperator(conDifOp);
@@ -232,6 +252,7 @@ namespace AMDiS {
matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
matShellContext.matConDif = conDifMatrixSolver->getMatInterior();
+
KSPSetType(matShellContext.kspMass, KSPRICHARDSON);
KSPSetTolerances(matShellContext.kspMass, 0, 1e-13, 1e+3, 1);
KSPGetPC(matShellContext.kspMass, &pcSub);
@@ -240,16 +261,72 @@ namespace AMDiS {
KSPSetFromOptions(matShellContext.kspMass);
KSPSetType(matShellContext.kspLaplace, KSPCG);
+ KSPSetInitialGuessNonzero(matShellContext.kspLaplace, PETSC_TRUE);
KSPSetTolerances(matShellContext.kspLaplace, 0, 1e-13, 1e+3, 10000);
KSPGetPC(matShellContext.kspLaplace, &pcSub);
PCSetType(pcSub, PCJACOBI);
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;
KSPSetFromOptions(kspVelocity);
KSPSetFromOptions(kspSchur);
+#endif
}
}
diff --git a/AMDiS/src/parallel/PetscSolverNavierStokes.h b/AMDiS/src/parallel/PetscSolverNavierStokes.h
index 6f97c66c..302f2b70 100644
--- a/AMDiS/src/parallel/PetscSolverNavierStokes.h
+++ b/AMDiS/src/parallel/PetscSolverNavierStokes.h
@@ -44,6 +44,10 @@ namespace AMDiS {
KSP kspMass, kspLaplace;
Mat matConDif;
+
+
+
+ IS isVelocity, isPressure;
};
class PetscSolverNavierStokes : public PetscSolverGlobalMatrix
--
GitLab