diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index 59bd92234e30807c46fa55bc0e963d1375998b8d..7c831b20e1fd050d702ed3bd62700f61ad8cf6d0 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -145,6 +145,8 @@
#include "parallel/Mtl4Solver.h"
#else
#include "parallel/PetscProblemStat.h"
+#include "parallel/PetscSolver.h"
+#include "parallel/PetscSolverNavierStokes.h"
#endif
#endif
diff --git a/AMDiS/src/AMDiS_fwd.h b/AMDiS/src/AMDiS_fwd.h
index a3ecfb0e8dcd4b7118aeecd7f15284b531a503e0..09f1ac32e38f4eeb5b0d9f192d4d36ab97ea7190 100644
--- a/AMDiS/src/AMDiS_fwd.h
+++ b/AMDiS/src/AMDiS_fwd.h
@@ -103,6 +103,7 @@ namespace AMDiS {
class FeSpaceDofMap;
class MatrixNnzStructure;
class MeshLevelData;
+ class PetscSolver;
class PetscSolverFeti;
class PetscSolverFetiDebug;
#endif
diff --git a/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.cc b/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.cc
index 3fbc4cc0d7635a264cefd8638663c11ca7435cc0..6728b3bb73b1bc70403c9b7b3a6894a3954a43f1 100644
--- a/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.cc
+++ b/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.cc
@@ -10,8 +10,8 @@
// See also license.opensource.txt in the distribution.
-#include "AMDiS.h"
#include "parallel/ParallelCoarseSpaceMatVec.h"
+#include "parallel/ParallelDofMapping.h"
#include "parallel/MatrixNnzStructure.h"
namespace AMDiS {
diff --git a/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h b/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h
index b4911552195c7006dd666d0ac533f59a66baa196..e1bcd2ce2ef4ad7f427063ae833d89478f9af9bb 100644
--- a/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h
+++ b/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h
@@ -23,10 +23,13 @@
#ifndef AMDIS_PARALLEL_COARSE_SPACE_MAT_VEC_H
#define AMDIS_PARALLEL_COARSE_SPACE_MAT_VEC_H
+#include <mpi.h>
#include <vector>
#include <map>
#include <petsc.h>
#include "AMDiS_fwd.h"
+#include "parallel/ParallelDofMapping.h"
+#include "parallel/MeshDistributor.h"
#include "parallel/MatrixNnzStructure.h"
namespace AMDiS {
diff --git a/AMDiS/src/parallel/PetscHelper.cc b/AMDiS/src/parallel/PetscHelper.cc
index 800709040db4ce7896f4671efabbd50191fe830d..2e560124afd7d585cc2f00ca4ff51f39d25dda3b 100644
--- a/AMDiS/src/parallel/PetscHelper.cc
+++ b/AMDiS/src/parallel/PetscHelper.cc
@@ -11,6 +11,7 @@
#include "parallel/PetscHelper.h"
+#include "parallel/PetscSolver.h"
#include "Global.h"
namespace AMDiS {
@@ -247,6 +248,25 @@ namespace AMDiS {
MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
}
+
+ void setSolver(KSP ksp,
+ string kspPrefix,
+ KSPType kspType,
+ PCType pcType,
+ PetscReal rtol,
+ PetscReal atol,
+ PetscInt maxIt)
+ {
+ KSPSetType(ksp, kspType);
+ KSPSetTolerances(ksp, rtol, atol, PETSC_DEFAULT, maxIt);
+ if (kspPrefix != "")
+ KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
+ KSPSetFromOptions(ksp);
+
+ PC pc;
+ KSPGetPC(ksp, &pc);
+ PCSetType(pc, pcType);
+ }
}
}
diff --git a/AMDiS/src/parallel/PetscHelper.h b/AMDiS/src/parallel/PetscHelper.h
index 9362aad873e912491ac7dce4ee25e7e6aaa3c0c0..e6ac70ea8ce6957d70d0139d6e55207080b9eec7 100644
--- a/AMDiS/src/parallel/PetscHelper.h
+++ b/AMDiS/src/parallel/PetscHelper.h
@@ -27,6 +27,7 @@
#include <map>
#include <vector>
#include <petsc.h>
+#include "AMDiS_fwd.h"
namespace AMDiS {
@@ -89,6 +90,14 @@ namespace AMDiS {
* \param[out] mat matrix of type MATAIJ, created inside this function.
*/
void matNestConvert(Mat matNest, Mat &mat);
+
+ void setSolver(KSP ksp,
+ string kspPrefix,
+ KSPType kspType,
+ PCType pcType,
+ PetscReal rtol = PETSC_DEFAULT,
+ PetscReal atol = PETSC_DEFAULT,
+ PetscInt maxIt = PETSC_DEFAULT);
}
}
diff --git a/AMDiS/src/parallel/PetscSolver.h b/AMDiS/src/parallel/PetscSolver.h
index 463c5ac1365ba1a5e645647faad7fcb7ca3bc56a..0e75538bfd665041519259d5ab2c271aeb9923ed 100644
--- a/AMDiS/src/parallel/PetscSolver.h
+++ b/AMDiS/src/parallel/PetscSolver.h
@@ -150,6 +150,11 @@ namespace AMDiS {
return dofMap;
}
+ vector<const FiniteElemSpace*>& getComponentSpaces()
+ {
+ return componentSpaces;
+ }
+
protected:
/** \brief
* Copies between to PETSc vectors by using different index sets for the
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 8eeb8a71a2d6a503c34d30e7761cc40dd12fb095..6597dc2c57e437a568280c0d33d73de7e86c4a3f 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -10,7 +10,6 @@
// See also license.opensource.txt in the distribution.
-#include "AMDiS.h"
#include "MatrixVector.h"
#include "parallel/PetscHelper.h"
#include "parallel/PetscSolverFeti.h"
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 7700d26c447132a9b17f8f4b1f59cf784c6fad66..18e75fa99b630f5bf47064a046da8c90b1825bfe 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -9,8 +9,8 @@
//
// See also license.opensource.txt in the distribution.
-
-#include "AMDiS.h"
+#include <mpi.h>
+#include "DirichletBC.h"
#include "parallel/PetscSolverGlobalMatrix.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
@@ -904,4 +904,39 @@ namespace AMDiS {
return subSolver;
}
+
+ void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp,
+ int constFeSpace,
+ bool test)
+ {
+ Vec nullSpaceBasis;
+ VecDuplicate(getVecSolInterior(), &nullSpaceBasis);
+
+ SystemVector basisVec("tmp", componentSpaces, componentSpaces.size(), true);
+ basisVec.set(0.0);
+ basisVec.getDOFVector(constFeSpace)->set(1.0);
+
+ setDofVector(nullSpaceBasis, basisVec, true);
+ VecAssemblyBegin(nullSpaceBasis);
+ VecAssemblyEnd(nullSpaceBasis);
+ VecNormalize(nullSpaceBasis, PETSC_NULL);
+
+ if (test) {
+ Vec tmp;
+ MatGetVecs(getMatInterior(), &tmp, PETSC_NULL);
+ MatMult(getMatInterior(), nullSpaceBasis, tmp);
+ PetscReal n;
+ VecNorm(tmp, NORM_2, &n);
+ MSG("NORM IS: %e\n", n);
+ VecDestroy(&tmp);
+ }
+
+
+ MatNullSpace matNullSpace;
+ MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullSpace);
+ KSPSetNullSpace(ksp, matNullSpace);
+ MatNullSpaceDestroy(&matNullSpace);
+ }
+
+
}
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
index e35a3bd951029ca8b567a40273095fc21cb68750..b3ff065d8b4a7812c570a05739c0cfdf226ae9ec 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
@@ -130,6 +130,8 @@ namespace AMDiS {
PetscSolver* createSubSolver(int component, string kspPrefix);
+ void setConstantNullSpace(KSP ksp, int constFeSpace, bool test = false);
+
protected:
bool zeroStartVector;
diff --git a/AMDiS/src/parallel/PetscSolverNavierStokes.cc b/AMDiS/src/parallel/PetscSolverNavierStokes.cc
index 070c7f471220eeefd1ef6be590e3a395c66b66d8..7bfd5d6565adc5a0ca513416fae8b96e1117fac5 100644
--- a/AMDiS/src/parallel/PetscSolverNavierStokes.cc
+++ b/AMDiS/src/parallel/PetscSolverNavierStokes.cc
@@ -12,98 +12,22 @@
#include "parallel/PetscSolverNavierStokes.h"
+#include "parallel/PetscHelper.h"
namespace AMDiS {
using namespace std;
- PetscErrorCode pcNs(PC pc, Vec x, Vec y)
+ PetscErrorCode pcSchurShell(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);
+ NavierStokesSchurData* data = static_cast<NavierStokesSchurData*>(ctx);
- 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;
- MatShellGetContext(mat, &ctx);
- MatShellNavierStokesSchur* data = static_cast<MatShellNavierStokesSchur*>(ctx);
-
- switch (data->solverMode) {
- case 0:
- {
- Vec vec0, vec1;
- VecDuplicate(x, &vec0);
- MatGetVecs(data->A00, &vec1, PETSC_NULL);
-
- MatMult(data->A11, x, y);
- MatMult(data->A01, x, vec1);
- KSPSolve(data->kspVelocity, vec1, vec1);
- MatMult(data->A10, vec1, y);
- VecAYPX(y, -1.0, vec0);
-
- VecDestroy(&vec0);
- VecDestroy(&vec1);
- }
- break;
- case 1:
- {
- Vec vec0, vec1;
- VecDuplicate(x, &vec0);
- VecDuplicate(x, &vec1);
-
- // KSPSolve(data->kspLaplace, x, vec0);
- // MatMult(data->matConDif, vec0, vec1);
- //KSPSolve(data->kspMass, vec0, y);
- VecSet(y, 0.0);
-
- VecDestroy(&vec0);
- VecDestroy(&vec1);
- }
- break;
- default:
- ERROR_EXIT("Wrong solver mode %d\n", data->solverMode);
- }
+ KSPSolve(data->kspLaplace, x, y);
+ MatMult(data->matConDif, y, x);
+ KSPSolve(data->kspMass, x, y);
}
@@ -111,7 +35,10 @@ namespace AMDiS {
: PetscSolverGlobalMatrix(name),
pressureComponent(-1),
massMatrixSolver(NULL),
- laplaceMatrixSolver(NULL)
+ laplaceMatrixSolver(NULL),
+ nu(NULL),
+ invTau(NULL),
+ solution(NULL)
{
Parameters::get(initFileStr + "->stokes->pressure component",
pressureComponent);
@@ -127,32 +54,14 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverNavierStokes::initSolver()");
+ // Create FGMRES based outer solver
KSPCreate(mpiCommGlobal, &ksp);
- KSPSetOperators(ksp, getMatInterior(), getMatInterior(),
- SAME_NONZERO_PATTERN);
- KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
- KSPSetType(ksp, KSPGMRES);
- KSPSetOptionsPrefix(ksp, "ns_");
- KSPSetFromOptions(ksp);
+ KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
+ petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 100);
-
- // === Create null space information. ===
-
- Vec nullSpaceBasis;
- VecDuplicate(getVecSolInterior(), &nullSpaceBasis);
-
- SystemVector basisVec("tmp", componentSpaces, componentSpaces.size(), true);
- basisVec.set(0.0);
- basisVec.getDOFVector(pressureComponent)->set(1.0);
-
- setDofVector(nullSpaceBasis, basisVec, true);
- VecAssemblyBegin(nullSpaceBasis);
- VecAssemblyEnd(nullSpaceBasis);
- VecNormalize(nullSpaceBasis, PETSC_NULL);
-
- MatNullSpace matNullSpace;
- MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullSpace);
+ // Create null space information.
+ setConstantNullSpace(ksp, pressureComponent);
}
@@ -160,35 +69,45 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverNavierStokes::initPreconditioner()");
+ vector<int> velocityComponents;
+ velocityComponents.push_back(0);
+ velocityComponents.push_back(1);
- 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);
+ 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);
+
+ petsc_helper::setSolver(kspVelocity, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
+
+ KSPSetType(kspSchur, KSPPREONLY);
PC pcSub;
- KSPGetPC(matShellContext.kspVelocity, &pcSub);
- PCSetType(pcSub, PCLU);
- PCFactorSetMatSolverPackage(pcSub, MATSOLVERMUMPS);
-
+ KSPGetPC(kspSchur, &pcSub);
+ PCSetType(pcSub, PCSHELL);
+ PCShellSetApply(pcSub, pcSchurShell);
+ PCShellSetContext(pcSub, &matShellContext);
- // === Mass matrix solver ===
+ MatNullSpace matNullSpace;
+ MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
+ KSPSetNullSpace(kspSchur, matNullSpace);
+ MatNullSpaceDestroy(&matNullSpace);
- MSG("Initialize mass matrix solver ...\n");
+
+ // === Mass matrix solver ===
const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
@@ -196,49 +115,37 @@ namespace AMDiS {
Simple_ZOT zot;
massOp.addTerm(&zot);
massMatrix.assembleOperator(massOp);
-
massMatrixSolver = createSubSolver(pressureComponent, "mass_");
massMatrixSolver->fillPetscMatrix(&massMatrix);
- MSG("... OK!\n");
-
// === Laplace matrix solver ===
- MSG("Initialize laplace matrix solver ...\n");
-
DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
Operator laplaceOp(pressureFeSpace, pressureFeSpace);
Simple_SOT sot;
laplaceOp.addTerm(&sot);
laplaceMatrix.assembleOperator(laplaceOp);
-
laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_");
laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
- MSG("... OK!\n");
-
// === Create convection-diffusion operator ===
- MSG("Initialize pressure convection-diffusion operator ...\n");
-
- double timestep = 1.0;
- Parameters::get("navierstokes->adapt->timestep", timestep);
- timestep = 1.0 / timestep;
-
- double nu = 0.01;
- Parameters::get("nu", nu);
+ DOFVector<double> vx(pressureFeSpace, "vx");
+ DOFVector<double> vy(pressureFeSpace, "vy");
+ vx.interpol(solution->getDOFVector(0));
+ vy.interpol(solution->getDOFVector(1));
DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
Operator conDifOp(pressureFeSpace, pressureFeSpace);
- Simple_ZOT conDif0(timestep);
+ Simple_ZOT conDif0(*invTau);
conDifOp.addTerm(&conDif0);
- Simple_SOT conDif1(nu);
+ Simple_SOT conDif1(*nu);
conDifOp.addTerm(&conDif1);
- Vector_FOT conDif2(0);
+ VecAtQP_FOT conDif2(&vx, &idFct, 0);
conDifOp.addTerm(&conDif2, GRD_PHI);
- Vector_FOT conDif3(1);
+ VecAtQP_FOT conDif3(&vy, &idFct, 1);
conDifOp.addTerm(&conDif3, GRD_PHI);
conDifMatrix.assembleOperator(conDifOp);
@@ -246,87 +153,16 @@ namespace AMDiS {
conDifMatrixSolver = createSubSolver(pressureComponent, "condif_");
conDifMatrixSolver->fillPetscMatrix(&conDifMatrix);
- MSG("... OK!\n");
+
+ // === Setup solver ===
matShellContext.kspMass = massMatrixSolver->getSolver();
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);
- PCSetType(pcSub, PCLU);
- PCFactorSetMatSolverPackage(pcSub, MATSOLVERMUMPS);
- 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
+ petsc_helper::setSolver(matShellContext.kspMass, "", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
+ petsc_helper::setSolver(matShellContext.kspLaplace, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
}
+
}
diff --git a/AMDiS/src/parallel/PetscSolverNavierStokes.h b/AMDiS/src/parallel/PetscSolverNavierStokes.h
index 302f2b7037130baf3dba7cb9498f0cdcd4d2e96c..354e235a9b6a89a59513367f4df9aaa73029ce53 100644
--- a/AMDiS/src/parallel/PetscSolverNavierStokes.h
+++ b/AMDiS/src/parallel/PetscSolverNavierStokes.h
@@ -29,32 +29,38 @@ namespace AMDiS {
using namespace std;
- struct MatShellNavierStokesSchur {
- int solverMode;
-
- // === Data for mode 0 (solve the schur complement directly) ===
-
- Mat A00, A01, A10, A11;
-
- KSP kspVelocity;
-
-
- // === Data for mode 1 ===
-
+ struct NavierStokesSchurData {
KSP kspMass, kspLaplace;
Mat matConDif;
-
-
-
- IS isVelocity, isPressure;
};
class PetscSolverNavierStokes : public PetscSolverGlobalMatrix
{
+ private:
+ class IdFct : public AbstractFunction<double, double>
+ {
+ public:
+ IdFct()
+ : AbstractFunction<double, double>(0)
+ {}
+
+ double operator()(const double& x) const
+ {
+ return x;
+ }
+ };
+
public:
PetscSolverNavierStokes(string name);
+ void setStokesData(double *nuPtr, double *invTauPtr, SystemVector *vec)
+ {
+ nu = nuPtr;
+ invTau = invTauPtr;
+ solution = vec;
+ }
+
protected:
void initSolver(KSP &ksp);
@@ -65,7 +71,13 @@ namespace AMDiS {
PetscSolver *massMatrixSolver, *laplaceMatrixSolver, *conDifMatrixSolver;
- MatShellNavierStokesSchur matShellContext;
+ NavierStokesSchurData matShellContext;
+
+ double *nu, *invTau;
+
+ SystemVector* solution;
+
+ IdFct idFct;
};
}