Commit ace7781a authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Navier stokes solver ready to be used.

parent a41536c3
...@@ -145,6 +145,8 @@ ...@@ -145,6 +145,8 @@
#include "parallel/Mtl4Solver.h" #include "parallel/Mtl4Solver.h"
#else #else
#include "parallel/PetscProblemStat.h" #include "parallel/PetscProblemStat.h"
#include "parallel/PetscSolver.h"
#include "parallel/PetscSolverNavierStokes.h"
#endif #endif
#endif #endif
......
...@@ -103,6 +103,7 @@ namespace AMDiS { ...@@ -103,6 +103,7 @@ namespace AMDiS {
class FeSpaceDofMap; class FeSpaceDofMap;
class MatrixNnzStructure; class MatrixNnzStructure;
class MeshLevelData; class MeshLevelData;
class PetscSolver;
class PetscSolverFeti; class PetscSolverFeti;
class PetscSolverFetiDebug; class PetscSolverFetiDebug;
#endif #endif
......
...@@ -10,8 +10,8 @@ ...@@ -10,8 +10,8 @@
// See also license.opensource.txt in the distribution. // See also license.opensource.txt in the distribution.
#include "AMDiS.h"
#include "parallel/ParallelCoarseSpaceMatVec.h" #include "parallel/ParallelCoarseSpaceMatVec.h"
#include "parallel/ParallelDofMapping.h"
#include "parallel/MatrixNnzStructure.h" #include "parallel/MatrixNnzStructure.h"
namespace AMDiS { namespace AMDiS {
......
...@@ -23,10 +23,13 @@ ...@@ -23,10 +23,13 @@
#ifndef AMDIS_PARALLEL_COARSE_SPACE_MAT_VEC_H #ifndef AMDIS_PARALLEL_COARSE_SPACE_MAT_VEC_H
#define AMDIS_PARALLEL_COARSE_SPACE_MAT_VEC_H #define AMDIS_PARALLEL_COARSE_SPACE_MAT_VEC_H
#include <mpi.h>
#include <vector> #include <vector>
#include <map> #include <map>
#include <petsc.h> #include <petsc.h>
#include "AMDiS_fwd.h" #include "AMDiS_fwd.h"
#include "parallel/ParallelDofMapping.h"
#include "parallel/MeshDistributor.h"
#include "parallel/MatrixNnzStructure.h" #include "parallel/MatrixNnzStructure.h"
namespace AMDiS { namespace AMDiS {
......
...@@ -11,6 +11,7 @@ ...@@ -11,6 +11,7 @@
#include "parallel/PetscHelper.h" #include "parallel/PetscHelper.h"
#include "parallel/PetscSolver.h"
#include "Global.h" #include "Global.h"
namespace AMDiS { namespace AMDiS {
...@@ -247,6 +248,25 @@ namespace AMDiS { ...@@ -247,6 +248,25 @@ namespace AMDiS {
MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY); 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);
}
} }
} }
...@@ -27,6 +27,7 @@ ...@@ -27,6 +27,7 @@
#include <map> #include <map>
#include <vector> #include <vector>
#include <petsc.h> #include <petsc.h>
#include "AMDiS_fwd.h"
namespace AMDiS { namespace AMDiS {
...@@ -89,6 +90,14 @@ namespace AMDiS { ...@@ -89,6 +90,14 @@ namespace AMDiS {
* \param[out] mat matrix of type MATAIJ, created inside this function. * \param[out] mat matrix of type MATAIJ, created inside this function.
*/ */
void matNestConvert(Mat matNest, Mat &mat); 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);
} }
} }
......
...@@ -150,6 +150,11 @@ namespace AMDiS { ...@@ -150,6 +150,11 @@ namespace AMDiS {
return dofMap; return dofMap;
} }
vector<const FiniteElemSpace*>& getComponentSpaces()
{
return componentSpaces;
}
protected: protected:
/** \brief /** \brief
* Copies between to PETSc vectors by using different index sets for the * Copies between to PETSc vectors by using different index sets for the
......
...@@ -10,7 +10,6 @@ ...@@ -10,7 +10,6 @@
// See also license.opensource.txt in the distribution. // See also license.opensource.txt in the distribution.
#include "AMDiS.h"
#include "MatrixVector.h" #include "MatrixVector.h"
#include "parallel/PetscHelper.h" #include "parallel/PetscHelper.h"
#include "parallel/PetscSolverFeti.h" #include "parallel/PetscSolverFeti.h"
......
...@@ -9,8 +9,8 @@ ...@@ -9,8 +9,8 @@
// //
// See also license.opensource.txt in the distribution. // See also license.opensource.txt in the distribution.
#include <mpi.h>
#include "AMDiS.h" #include "DirichletBC.h"
#include "parallel/PetscSolverGlobalMatrix.h" #include "parallel/PetscSolverGlobalMatrix.h"
#include "parallel/StdMpi.h" #include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h" #include "parallel/MpiHelper.h"
...@@ -904,4 +904,39 @@ namespace AMDiS { ...@@ -904,4 +904,39 @@ namespace AMDiS {
return subSolver; 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);
}
} }
...@@ -130,6 +130,8 @@ namespace AMDiS { ...@@ -130,6 +130,8 @@ namespace AMDiS {
PetscSolver* createSubSolver(int component, string kspPrefix); PetscSolver* createSubSolver(int component, string kspPrefix);
void setConstantNullSpace(KSP ksp, int constFeSpace, bool test = false);
protected: protected:
bool zeroStartVector; bool zeroStartVector;
......
...@@ -12,98 +12,22 @@ ...@@ -12,98 +12,22 @@
#include "parallel/PetscSolverNavierStokes.h" #include "parallel/PetscSolverNavierStokes.h"
#include "parallel/PetscHelper.h"
namespace AMDiS { namespace AMDiS {
using namespace std; using namespace std;
PetscErrorCode pcNs(PC pc, Vec x, Vec y) PetscErrorCode pcSchurShell(PC pc, Vec x, Vec y)
{ {
void *ctx; void *ctx;
PCShellGetContext(pc, &ctx); PCShellGetContext(pc, &ctx);
MatShellNavierStokesSchur* data = static_cast<MatShellNavierStokesSchur*>(ctx); NavierStokesSchurData* data = static_cast<NavierStokesSchurData*>(ctx);
Vec velocity, pressure; KSPSolve(data->kspLaplace, x, y);
VecGetSubVector(x, data->isVelocity, &velocity); MatMult(data->matConDif, y, x);
VecGetSubVector(x, data->isPressure, &pressure); KSPSolve(data->kspMass, x, y);
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);
}
} }
...@@ -111,7 +35,10 @@ namespace AMDiS { ...@@ -111,7 +35,10 @@ namespace AMDiS {
: PetscSolverGlobalMatrix(name), : PetscSolverGlobalMatrix(name),
pressureComponent(-1), pressureComponent(-1),
massMatrixSolver(NULL), massMatrixSolver(NULL),
laplaceMatrixSolver(NULL) laplaceMatrixSolver(NULL),
nu(NULL),
invTau(NULL),
solution(NULL)
{ {
Parameters::get(initFileStr + "->stokes->pressure component", Parameters::get(initFileStr + "->stokes->pressure component",
pressureComponent); pressureComponent);
...@@ -127,68 +54,60 @@ namespace AMDiS { ...@@ -127,68 +54,60 @@ namespace AMDiS {
{ {
FUNCNAME("PetscSolverNavierStokes::initSolver()"); FUNCNAME("PetscSolverNavierStokes::initSolver()");
// Create FGMRES based outer solver
KSPCreate(mpiCommGlobal, &ksp); KSPCreate(mpiCommGlobal, &ksp);
KSPSetOperators(ksp, getMatInterior(), getMatInterior(), KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
SAME_NONZERO_PATTERN);
KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
KSPSetType(ksp, KSPGMRES);
KSPSetOptionsPrefix(ksp, "ns_");
KSPSetFromOptions(ksp);
KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL); KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 100);
// Create null space information.
setConstantNullSpace(ksp, pressureComponent);
}
// === Create null space information. ===
Vec nullSpaceBasis; void PetscSolverNavierStokes::initPreconditioner(PC pc)
VecDuplicate(getVecSolInterior(), &nullSpaceBasis); {
FUNCNAME("PetscSolverNavierStokes::initPreconditioner()");
SystemVector basisVec("tmp", componentSpaces, componentSpaces.size(), true); vector<int> velocityComponents;
basisVec.set(0.0); velocityComponents.push_back(0);
basisVec.getDOFVector(pressureComponent)->set(1.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);
setDofVector(nullSpaceBasis, basisVec, true); KSPSetUp(kspInterior);
VecAssemblyBegin(nullSpaceBasis);
VecAssemblyEnd(nullSpaceBasis);
VecNormalize(nullSpaceBasis, PETSC_NULL);
MatNullSpace matNullSpace; KSP *subKsp;
MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullSpace); int nSubKsp;
} PCFieldSplitGetSubKSP(pc, &nSubKsp, &subKsp);
TEST_EXIT(nSubKsp == 2)
("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");
void PetscSolverNavierStokes::initPreconditioner(PC pc) KSP kspVelocity = subKsp[0];
{ KSP kspSchur = subKsp[1];
FUNCNAME("PetscSolverNavierStokes::initPreconditioner()"); PetscFree(subKsp);
petsc_helper::setSolver(kspVelocity, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
PCSetType(pc, PCSHELL); KSPSetType(kspSchur, KSPPREONLY);
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; PC pcSub;
KSPGetPC(matShellContext.kspVelocity, &pcSub); KSPGetPC(kspSchur, &pcSub);
PCSetType(pcSub, PCLU); PCSetType(pcSub, PCSHELL);
PCFactorSetMatSolverPackage(pcSub, MATSOLVERMUMPS); PCShellSetApply(pcSub, pcSchurShell);
PCShellSetContext(pcSub, &matShellContext);
MatNullSpace matNullSpace;
MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
KSPSetNullSpace(kspSchur, matNullSpace);
MatNullSpaceDestroy(&matNullSpace);
// === Mass matrix solver ===
MSG("Initialize mass matrix solver ...\n"); // === Mass matrix solver ===
const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent]; const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace); DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
...@@ -196,49 +115,37 @@ namespace AMDiS { ...@@ -196,49 +115,37 @@ namespace AMDiS {
Simple_ZOT zot; Simple_ZOT zot;
massOp.addTerm(&zot); massOp.addTerm(&zot);
massMatrix.assembleOperator(massOp); massMatrix.assembleOperator(massOp);
massMatrixSolver = createSubSolver(pressureComponent, "mass_"); massMatrixSolver = createSubSolver(pressureComponent, "mass_");
massMatrixSolver->fillPetscMatrix(&massMatrix); massMatrixSolver->fillPetscMatrix(&massMatrix);
MSG("... OK!\n");
// === Laplace matrix solver === // === Laplace matrix solver ===
MSG("Initialize laplace matrix solver ...\n");
DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace); DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
Operator laplaceOp(pressureFeSpace, pressureFeSpace); Operator laplaceOp(pressureFeSpace, pressureFeSpace);
Simple_SOT sot; Simple_SOT sot;
laplaceOp.addTerm(&sot); laplaceOp.addTerm(&sot);
laplaceMatrix.assembleOperator(laplaceOp); laplaceMatrix.assembleOperator(laplaceOp);
laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_"); laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_");
laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix); laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
MSG("... OK!\n");
// === Create convection-diffusion operator === // === Create convection-diffusion operator ===
MSG("Initialize pressure convection-diffusion operator ...\n"); DOFVector<double> vx(pressureFeSpace, "vx");
DOFVector<double> vy(pressureFeSpace, "vy");
double timestep = 1.0; vx.interpol(solution->getDOFVector(0));
Parameters::get("navierstokes->adapt->timestep", timestep); vy.interpol(solution->getDOFVector(1));
timestep = 1.0 / timestep;
double nu = 0.01;
Parameters::get("nu", nu);
DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace); DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
Operator conDifOp(pressureFeSpace, pressureFeSpace); Operator conDifOp(pressureFeSpace, pressureFeSpace);
Simple_ZOT conDif0(timestep); Simple_ZOT conDif0(*invTau);
conDifOp.addTerm(&conDif0); conDifOp.addTerm(&conDif0);
Simple_SOT conDif1(nu); Simple_SOT conDif1(*nu);
conDifOp.addTerm(&conDif1); conDifOp.addTerm(&conDif1);
Vector_FOT conDif2(0); VecAtQP_FOT conDif2(&vx, &idFct, 0);
conDifOp.addTerm(&conDif2, GRD_PHI); conDifOp.addTerm(&conDif2, GRD_PHI);
Vector_FOT conDif3(1); VecAtQP_FOT conDif3(&vy, &idFct, 1);
conDifOp.addTerm(&conDif3, GRD_PHI); conDifOp.addTerm(&conDif3, GRD_PHI);
conDifMatrix.assembleOperator(conDifOp); conDifMatrix.assembleOperator(conDifOp);
...@@ -246,87 +153,16 @@ namespace AMDiS { ...@@ -246,87 +153,16 @@ namespace AMDiS {
conDifMatrixSolver = createSubSolver(pressureComponent, "condif_"); conDifMatrixSolver = createSubSolver(pressureComponent, "condif_");
conDifMatrixSolver->fillPetscMatrix(&conDifMatrix); conDifMatrixSolver->fillPetscMatrix(&conDifMatrix);
MSG("... OK!\n");
// === Setup solver ===
matShellContext.kspMass = massMatrixSolver->getSolver(); matShellContext.kspMass = massMatrixSolver->getSolver();
matShellContext.kspLaplace = laplaceMatrixSolver->getSolver(); matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
matShellContext.matConDif = conDifMatrixSolver->getMatInterior(); matShellContext.matConDif = conDifMatrixSolver->getMatInterior();
petsc_helper::setSolver(matShellContext.kspMass, "", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
KSPSetType(matShellContext.kspMass, KSPRICHARDSON); petsc_helper::setSolver(matShellContext.kspLaplace, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
KSPSetTolerances(matShellContext.kspMass, 0, 1e-13, 1e+3, 1);