Commit 0ca6954f authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Work on this totaly bull shit navier stokes what ever I have never heard about...

Work on this totaly bull shit navier stokes what ever I have never heard about solver, i will go home and eat one or two cakes.
parent 9cd0a12c
...@@ -30,7 +30,7 @@ namespace AMDiS { ...@@ -30,7 +30,7 @@ namespace AMDiS {
SerUtil::deserialize(in, size); SerUtil::deserialize(in, size);
vectors.resize(size); vectors.resize(size);
for (int i = oldSize; i < size; i++) for (int i = oldSize; i < size; i++)
vectors[i] = new DOFVector<double>(feSpace[i], ""); vectors[i] = new DOFVector<double>(componentSpaces[i], "");
for (int i = 0; i < size; i++) for (int i = 0; i < size; i++)
vectors[i]->deserialize(in); vectors[i]->deserialize(in);
} }
......
...@@ -37,16 +37,21 @@ namespace AMDiS { ...@@ -37,16 +37,21 @@ namespace AMDiS {
/// Constructor. /// Constructor.
SystemVector(std::string name_, SystemVector(std::string name_,
std::vector<const FiniteElemSpace*> feSpace_, std::vector<const FiniteElemSpace*> feSpace_,
int size) int size,
bool createVec = false)
: name(name_), : name(name_),
feSpace(feSpace_), componentSpaces(feSpace_),
vectors(size) vectors(size)
{} {
if (createVec)
for (int i = 0; i < size; i++)
vectors[i] = new DOFVector<double>(componentSpaces[i], "tmp");
}
/// Copy Constructor. /// Copy Constructor.
SystemVector(const SystemVector& rhs) SystemVector(const SystemVector& rhs)
: name(rhs.getName()), : name(rhs.getName()),
feSpace(rhs.getFeSpaces()), componentSpaces(rhs.getFeSpaces()),
vectors(rhs.getNumVectors()) vectors(rhs.getNumVectors())
{ {
for (unsigned int i = 0; i < vectors.size(); i++) for (unsigned int i = 0; i < vectors.size(); i++)
...@@ -111,13 +116,13 @@ namespace AMDiS { ...@@ -111,13 +116,13 @@ namespace AMDiS {
/// Returns the fe space for a given component. /// Returns the fe space for a given component.
inline const FiniteElemSpace *getFeSpace(int i) const inline const FiniteElemSpace *getFeSpace(int i) const
{ {
return feSpace[i]; return componentSpaces[i];
} }
/// Returns the fe spaces for all components. /// Returns the fe spaces for all components.
inline std::vector<const FiniteElemSpace*> getFeSpaces() const inline std::vector<const FiniteElemSpace*> getFeSpaces() const
{ {
return feSpace; return componentSpaces;
} }
/// Here the system vector is interpreted as one large vector. The given /// Here the system vector is interpreted as one large vector. The given
...@@ -228,7 +233,7 @@ namespace AMDiS { ...@@ -228,7 +233,7 @@ namespace AMDiS {
std::string name; std::string name;
/// Finite element space. /// Finite element space.
std::vector<const FiniteElemSpace*> feSpace; std::vector<const FiniteElemSpace*> componentSpaces;
/// Local dof vectors. /// Local dof vectors.
std::vector<DOFVector<double>*> vectors; std::vector<DOFVector<double>*> vectors;
......
...@@ -35,8 +35,8 @@ namespace AMDiS { ...@@ -35,8 +35,8 @@ namespace AMDiS {
class BddcMlSolver : public PetscSolver class BddcMlSolver : public PetscSolver
{ {
public: public:
BddcMlSolver() BddcMlSolver(string name)
: PetscSolver(), : PetscSolver(name),
rhsVec(NULL), rhsVec(NULL),
mat(NULL) mat(NULL)
{} {}
......
...@@ -37,21 +37,21 @@ namespace AMDiS { ...@@ -37,21 +37,21 @@ namespace AMDiS {
Parameters::get(initFileStr, tmp); Parameters::get(initFileStr, tmp);
if (tmp == "petsc-schur") { if (tmp == "petsc-schur") {
petscSolver = new PetscSolverSchur(); petscSolver = new PetscSolverSchur(initFileStr);
} else if (tmp == "petsc-feti") { } else if (tmp == "petsc-feti") {
petscSolver = new PetscSolverFeti(initFileStr); petscSolver = new PetscSolverFeti(initFileStr);
} else if (tmp == "petsc-block") { } else if (tmp == "petsc-block") {
petscSolver = new PetscSolverGlobalBlockMatrix(); petscSolver = new PetscSolverGlobalBlockMatrix(initFileStr);
} else if (tmp == "petsc") { } else if (tmp == "petsc") {
petscSolver = new PetscSolverGlobalMatrix(); petscSolver = new PetscSolverGlobalMatrix(initFileStr);
} else if (tmp == "bddcml") { } else if (tmp == "bddcml") {
#ifdef HAVE_BDDC_ML #ifdef HAVE_BDDC_ML
petscSolver = new BddcMlSolver(); petscSolver = new BddcMlSolver(initFileStr);
#else #else
ERROR_EXIT("AMDiS was compiled without BDDC-ML support!\n"); ERROR_EXIT("AMDiS was compiled without BDDC-ML support!\n");
#endif #endif
} else if (tmp == "petsc-stokes") { } else if (tmp == "petsc-stokes") {
petscSolver = new PetscSolverNavierStokes(); petscSolver = new PetscSolverNavierStokes(initFileStr);
} else { } else {
ERROR_EXIT("No parallel solver %s available!\n", tmp.c_str()); ERROR_EXIT("No parallel solver %s available!\n", tmp.c_str());
} }
......
...@@ -20,8 +20,9 @@ namespace AMDiS { ...@@ -20,8 +20,9 @@ namespace AMDiS {
using namespace std; using namespace std;
PetscSolver::PetscSolver() PetscSolver::PetscSolver(string name)
: ParallelCoarseSpaceMatVec(), : ParallelCoarseSpaceMatVec(),
initFileStr(name),
dofMap(FESPACE_WISE, true), dofMap(FESPACE_WISE, true),
dofMapSd(FESPACE_WISE, true), dofMapSd(FESPACE_WISE, true),
kspPrefix(""), kspPrefix(""),
......
...@@ -50,7 +50,7 @@ namespace AMDiS { ...@@ -50,7 +50,7 @@ namespace AMDiS {
class PetscSolver : public ParallelCoarseSpaceMatVec class PetscSolver : public ParallelCoarseSpaceMatVec
{ {
public: public:
PetscSolver(); PetscSolver(string name);
virtual ~PetscSolver() {} virtual ~PetscSolver() {}
...@@ -145,6 +145,11 @@ namespace AMDiS { ...@@ -145,6 +145,11 @@ namespace AMDiS {
handleDirichletRows = b; handleDirichletRows = b;
} }
ParallelDofMapping& getDofMap()
{
return dofMap;
}
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
...@@ -163,6 +168,9 @@ namespace AMDiS { ...@@ -163,6 +168,9 @@ namespace AMDiS {
/// Run test, if matrix is symmetric. /// Run test, if matrix is symmetric.
bool testMatrixSymmetric(Mat mat, bool advancedTest = false); bool testMatrixSymmetric(Mat mat, bool advancedTest = false);
protected: protected:
/// Prefix string for parameters in init file.
string initFileStr;
/// FE spaces of all components for the stationary problem the specific /// FE spaces of all components for the stationary problem the specific
/// solver object is registered to. /// solver object is registered to.
vector<const FiniteElemSpace*> componentSpaces; vector<const FiniteElemSpace*> componentSpaces;
......
...@@ -29,8 +29,7 @@ namespace AMDiS { ...@@ -29,8 +29,7 @@ namespace AMDiS {
using namespace std; using namespace std;
PetscSolverFeti::PetscSolverFeti(string name) PetscSolverFeti::PetscSolverFeti(string name)
: PetscSolver(), : PetscSolver(name),
initFileStr(name),
primalDofMap(COMPONENT_WISE), primalDofMap(COMPONENT_WISE),
dualDofMap(COMPONENT_WISE), dualDofMap(COMPONENT_WISE),
interfaceDofMap(COMPONENT_WISE), interfaceDofMap(COMPONENT_WISE),
...@@ -108,7 +107,7 @@ namespace AMDiS { ...@@ -108,7 +107,7 @@ namespace AMDiS {
MeshLevelData& levelData = meshDistributor->getMeshLevelData(); MeshLevelData& levelData = meshDistributor->getMeshLevelData();
if (subdomain == NULL) { if (subdomain == NULL) {
subdomain = new PetscSolverGlobalMatrix(); subdomain = new PetscSolverGlobalMatrix("");
subdomain->setSymmetric(isSymmetric); subdomain->setSymmetric(isSymmetric);
subdomain->setHandleDirichletRows(false); subdomain->setHandleDirichletRows(false);
...@@ -1302,7 +1301,7 @@ namespace AMDiS { ...@@ -1302,7 +1301,7 @@ namespace AMDiS {
massMatrix.assembleOperator(op); massMatrix.assembleOperator(op);
if (!massMatrixSolver) { if (!massMatrixSolver) {
massMatrixSolver = new PetscSolverGlobalMatrix; massMatrixSolver = new PetscSolverGlobalMatrix("");
massMatrixSolver->setKspPrefix("mass_"); massMatrixSolver->setKspPrefix("mass_");
massMatrixSolver->setMeshDistributor(meshDistributor, massMatrixSolver->setMeshDistributor(meshDistributor,
mpiCommGlobal, mpiCommGlobal,
......
...@@ -230,9 +230,6 @@ namespace AMDiS { ...@@ -230,9 +230,6 @@ namespace AMDiS {
} }
protected: protected:
/// Prefix string for parameters in init file.
string initFileStr;
/// Mapping from primal DOF indices to a global index of primals. /// Mapping from primal DOF indices to a global index of primals.
ParallelDofMapping primalDofMap; ParallelDofMapping primalDofMap;
......
...@@ -34,8 +34,8 @@ namespace AMDiS { ...@@ -34,8 +34,8 @@ namespace AMDiS {
class PetscSolverGlobalBlockMatrix : public PetscSolver class PetscSolverGlobalBlockMatrix : public PetscSolver
{ {
public: public:
PetscSolverGlobalBlockMatrix() PetscSolverGlobalBlockMatrix(string name)
: PetscSolver(), : PetscSolver(name),
nComponents(0), nComponents(0),
nBlocks(-1) nBlocks(-1)
{} {}
......
...@@ -321,8 +321,7 @@ namespace AMDiS { ...@@ -321,8 +321,7 @@ namespace AMDiS {
if (nullspace.size() > 0) { if (nullspace.size() > 0) {
VecDuplicate(getVecSolInterior(), &nullspaceBasis); VecDuplicate(getVecSolInterior(), &nullspaceBasis);
for (int i = 0; i < nComponents; i++) setDofVector(nullspaceBasis, *(nullspace[0]), true);
setDofVector(nullspaceBasis, nullspace[0]->getDOFVector(i), i, true);
VecAssemblyBegin(nullspaceBasis); VecAssemblyBegin(nullspaceBasis);
VecAssemblyEnd(nullspaceBasis); VecAssemblyEnd(nullspaceBasis);
...@@ -882,4 +881,27 @@ namespace AMDiS { ...@@ -882,4 +881,27 @@ namespace AMDiS {
} }
} }
PetscSolver* PetscSolverGlobalMatrix::createSubSolver(int component,
string kspPrefix)
{
FUNCNAME("PetscSolverGlobalMatrix::createSubSolver()");
vector<const FiniteElemSpace*> fe;
fe.push_back(componentSpaces[component]);
PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
subSolver->setKspPrefix(kspPrefix.c_str());
subSolver->setMeshDistributor(meshDistributor,
mpiCommGlobal,
mpiCommLocal);
subSolver->init(fe, fe);
ParallelDofMapping &subDofMap = subSolver->getDofMap();
subDofMap[0] = dofMap[component];
subDofMap.update();
return subSolver;
}
} }
...@@ -39,8 +39,8 @@ namespace AMDiS { ...@@ -39,8 +39,8 @@ namespace AMDiS {
class PetscSolverGlobalMatrix : public PetscSolver class PetscSolverGlobalMatrix : public PetscSolver
{ {
public: public:
PetscSolverGlobalMatrix() PetscSolverGlobalMatrix(string name)
: PetscSolver(), : PetscSolver(name),
zeroStartVector(false), zeroStartVector(false),
printMatInfo(false) printMatInfo(false)
{ {
...@@ -88,6 +88,14 @@ namespace AMDiS { ...@@ -88,6 +88,14 @@ namespace AMDiS {
*/ */
void createFieldSplit(PC pc, string splitName, vector<int> &components); void createFieldSplit(PC pc, string splitName, vector<int> &components);
/// Wrapper to create field split from only one component.
void createFieldSplit(PC pc, string splitName, int component)
{
vector<int> components;
components.push_back(component);
createFieldSplit(pc, splitName, components);
}
virtual void initSolver(KSP &ksp); virtual void initSolver(KSP &ksp);
virtual void exitSolver(KSP ksp); virtual void exitSolver(KSP ksp);
...@@ -112,6 +120,15 @@ namespace AMDiS { ...@@ -112,6 +120,15 @@ namespace AMDiS {
setDofVector(vecInterior, PETSC_NULL, vec, nRowVec, rankOnly); setDofVector(vecInterior, PETSC_NULL, vec, nRowVec, rankOnly);
} }
inline void setDofVector(Vec vecInterior,
SystemVector &vec,
bool rankOnly = false)
{
for (int i = 0; i < vec.getSize(); i++)
setDofVector(vecInterior, PETSC_NULL, vec.getDOFVector(i), i, rankOnly);
}
PetscSolver* createSubSolver(int component, string kspPrefix);
protected: protected:
bool zeroStartVector; bool zeroStartVector;
......
...@@ -17,13 +17,70 @@ namespace AMDiS { ...@@ -17,13 +17,70 @@ namespace AMDiS {
using namespace std; using namespace std;
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, vec1, y);
VecDestroy(&vec0);
VecDestroy(&vec1);
}
break;
default:
ERROR_EXIT("Wrong solver mode %d\n", data->solverMode);
}
}
PetscSolverNavierStokes::PetscSolverNavierStokes(string name)
: PetscSolverGlobalMatrix(name),
pressureComponent(-1),
massMatrixSolver(NULL),
laplaceMatrixSolver(NULL)
{
Parameters::get(initFileStr + "->stokes->pressure component",
pressureComponent);
TEST_EXIT(pressureComponent >= 0)
("For using PETSc stokes solver you must define a pressure component!\n");
TEST_EXIT(pressureComponent == 2)("Fixed for pressure component = 2\n");
}
void PetscSolverNavierStokes::initSolver(KSP &ksp) void PetscSolverNavierStokes::initSolver(KSP &ksp)
{ {
FUNCNAME("PetscSolverNavierStokes::initSolver()"); FUNCNAME("PetscSolverNavierStokes::initSolver()");
MSG("RUN NAVIER STOKES SOLVER INIT!\n");
KSPCreate(mpiCommGlobal, &ksp); KSPCreate(mpiCommGlobal, &ksp);
KSPSetOperators(ksp, getMatInterior(), getMatInterior(), KSPSetOperators(ksp, getMatInterior(), getMatInterior(),
SAME_NONZERO_PATTERN); SAME_NONZERO_PATTERN);
...@@ -31,7 +88,25 @@ namespace AMDiS { ...@@ -31,7 +88,25 @@ namespace AMDiS {
KSPSetType(ksp, KSPGMRES); KSPSetType(ksp, KSPGMRES);
KSPSetOptionsPrefix(ksp, "ns_"); KSPSetOptionsPrefix(ksp, "ns_");
KSPSetFromOptions(ksp); KSPSetFromOptions(ksp);
KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL); KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
// === 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);
} }
...@@ -39,13 +114,9 @@ namespace AMDiS { ...@@ -39,13 +114,9 @@ namespace AMDiS {
{ {
FUNCNAME("PetscSolverNavierStokes::initPreconditioner()"); FUNCNAME("PetscSolverNavierStokes::initPreconditioner()");
MSG("RUN NAVIER STOKES PRECONDITIONER INIT!\n");
vector<int> velocityComponents; vector<int> velocityComponents;
velocityComponents.push_back(0); velocityComponents.push_back(0);
velocityComponents.push_back(1); velocityComponents.push_back(1);
vector<int> pressureComponent;
pressureComponent.push_back(2);
PCSetType(pc, PCFIELDSPLIT); PCSetType(pc, PCFIELDSPLIT);
PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR); PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
...@@ -59,42 +130,126 @@ namespace AMDiS { ...@@ -59,42 +130,126 @@ namespace AMDiS {
int nSubKsp; int nSubKsp;
PCFieldSplitGetSubKSP(pc, &nSubKsp, &subKsp); PCFieldSplitGetSubKSP(pc, &nSubKsp, &subKsp);
TEST_EXIT(nSubKsp == 2)("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n"); TEST_EXIT(nSubKsp == 2)
("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");
KSP velocityKsp = subKsp[0]; KSP kspVelocity = subKsp[0];
KSP schurKsp = subKsp[1]; KSP kspSchur = subKsp[1];
PetscFree(subKsp); PetscFree(subKsp);
Mat A00, A01, A10, A11; Mat A00, A01, A10, A11;
PCFieldSplitGetSchurBlocks(pc, &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);
PC pcSub;
KSPGetPC(kspVelocity, &pcSub);
PCSetType(pcSub, PCLU);
PCFactorSetMatSolverPackage(pcSub, MATSOLVERMUMPS);
MatInfo minfo;
PetscInt nRow, nCol;
MatGetSize(A00, &nRow, &nCol);
MatGetInfo(A00, MAT_GLOBAL_SUM, &minfo);
MSG("A00: %d x %d with %d nnz\n", nRow, nCol, static_cast<int>(minfo.nz_used));
MatGetSize(A01, &nRow, &nCol); KSPSetType(kspSchur, KSPGMRES);
MatGetInfo(A01, MAT_GLOBAL_SUM, &minfo); KSPSetTolerances(kspSchur, 0, 1e-14, 1e+3, 1000);
MSG("A01: %d x %d with %d nnz\n", nRow, nCol, static_cast<int>(minfo.nz_used)); KSPSetOperators(kspSchur, matShell, matShell, SAME_NONZERO_PATTERN);
KSPGetPC(kspSchur, &pcSub);
PCSetType(pcSub, PCNONE);
MatGetSize(A10, &nRow, &nCol);
MatGetInfo(A10, MAT_GLOBAL_SUM, &minfo);
MSG("A10: %d x %d with %d nnz\n", nRow, nCol, static_cast<int>(minfo.nz_used));
MatGetSize(A11, &nRow, &nCol); // === Mass matrix solver ===
MatGetInfo(A11, MAT_GLOBAL_SUM, &minfo);
MSG("A11: %d x %d with %d nnz\n", nRow, nCol, static_cast<int>(minfo.nz_used));
KSPSetType(velocityKsp, KSPPREONLY); MSG("Initialize mass matrix solver ...\n");
PC pcSub;
KSPGetPC(velocityKsp, &pcSub); const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
Operator massOp(pressureFeSpace, pressureFeSpace);
Simple_ZOT zot;
massOp.addTerm(&zot);
massMatrix.assembleOperator(massOp);
massMatrixSolver = createSubSolver(pressureComponent, "mass_");
massMatrixSolver->fillPetscMatrix(&massMatrix);
MSG("... OK!\n");