Commit 48245c1c authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Fixed memory issues in navier stokes solver.

parent 8f6cf14b
...@@ -657,12 +657,12 @@ namespace AMDiS { ...@@ -657,12 +657,12 @@ namespace AMDiS {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
vector<FixRefinementPatch::EdgeInEl> refineEdges; vector<FixRefinementPatch::EdgeInEl> refineEdges;
FixRefinementPatch::getOtherEl(stack, refineEdges); FixRefinementPatch::getOtherEl(stack, refineEdges);
if (refineEdges.size()) { // if (refineEdges.size()) {
MSG("FIX REFINEMENT PATH ON ELEMENT %d %d: %d additional edges\n", // MSG("FIX REFINEMENT PATH ON ELEMENT %d %d: %d additional edges\n",
elInfo->getElement()->getIndex(), // elInfo->getElement()->getIndex(),
elInfo->getMacroElement()->getIndex(), // elInfo->getMacroElement()->getIndex(),
refineEdges.size()); // refineEdges.size());
} // }
#endif #endif
// === Traverse and refine the refinement patch. ==== // === Traverse and refine the refinement patch. ====
...@@ -679,13 +679,6 @@ namespace AMDiS { ...@@ -679,13 +679,6 @@ namespace AMDiS {
for (int edgeIndex = 0; for (int edgeIndex = 0;
edgeIndex < static_cast<unsigned int>(refineEdges.size()); edgeIndex++) { edgeIndex < static_cast<unsigned int>(refineEdges.size()); edgeIndex++) {
MSG(" IN REF FIX: %d %d\n", refineEdges[edgeIndex].first->getIndex(), refineEdges[edgeIndex].second);
}
for (int edgeIndex = 0;
edgeIndex < static_cast<unsigned int>(refineEdges.size()); edgeIndex++) {
MSG(" MAKE -> %d\n", edgeIndex);
Element *otherEl = refineEdges[edgeIndex].first; Element *otherEl = refineEdges[edgeIndex].first;
TraverseStack stack2; TraverseStack stack2;
ElInfo *elInfo2 = ElInfo *elInfo2 =
......
...@@ -1049,7 +1049,6 @@ namespace AMDiS { ...@@ -1049,7 +1049,6 @@ namespace AMDiS {
// === Check the boundaries and adapt mesh if necessary. === // === Check the boundaries and adapt mesh if necessary. ===
MSG_DBG("Run checkAndAdaptBoundary ...\n"); MSG_DBG("Run checkAndAdaptBoundary ...\n");
double second = MPI::Wtime();
// Check for periodic boundaries within rank's subdomain. // Check for periodic boundaries within rank's subdomain.
for (InteriorBoundary::iterator it(intBoundary.getPeriodic()); for (InteriorBoundary::iterator it(intBoundary.getPeriodic());
!it.end(); ++it) { !it.end(); ++it) {
...@@ -1064,12 +1063,9 @@ namespace AMDiS { ...@@ -1064,12 +1063,9 @@ namespace AMDiS {
} }
} }
} }
MSG(" -> create mesh structure codes needed %.5f seconds\n", MPI::Wtime() - second);
second = MPI::Wtime();
meshChanged |= checkAndAdaptBoundary(allBound); meshChanged |= checkAndAdaptBoundary(allBound);
MSG(" -> checkAndAdaptBoundary needed %.5f seconds\n", MPI::Wtime() - second);
// === Check on all ranks if at least one rank's mesh has changed. === // === Check on all ranks if at least one rank's mesh has changed. ===
......
...@@ -249,10 +249,11 @@ namespace AMDiS { ...@@ -249,10 +249,11 @@ namespace AMDiS {
} }
void setSolver(KSP ksp, void setSolverWithLu(KSP ksp,
const char* kspPrefix, const char* kspPrefix,
KSPType kspType, KSPType kspType,
PCType pcType, PCType pcType,
MatSolverPackage matSolverPackage,
PetscReal rtol, PetscReal rtol,
PetscReal atol, PetscReal atol,
PetscInt maxIt) PetscInt maxIt)
...@@ -267,8 +268,23 @@ namespace AMDiS { ...@@ -267,8 +268,23 @@ namespace AMDiS {
PC pc; PC pc;
KSPGetPC(ksp, &pc); KSPGetPC(ksp, &pc);
PCSetType(pc, pcType); PCSetType(pc, pcType);
if (matSolverPackage != PETSC_NULL)
PCFactorSetMatSolverPackage(pc, matSolverPackage);
PCSetFromOptions(pc); PCSetFromOptions(pc);
} }
void setSolver(KSP ksp,
const char* kspPrefix,
KSPType kspType,
PCType pcType,
PetscReal rtol,
PetscReal atol,
PetscInt maxIt)
{
setSolverWithLu(ksp, kspPrefix, kspType, pcType, PETSC_NULL,
rtol, atol, maxIt);
}
} }
} }
...@@ -91,6 +91,15 @@ namespace AMDiS { ...@@ -91,6 +91,15 @@ namespace AMDiS {
*/ */
void matNestConvert(Mat matNest, Mat &mat); void matNestConvert(Mat matNest, Mat &mat);
void setSolverWithLu(KSP ksp,
const char* kspPrefix,
KSPType kspType,
PCType pcType,
MatSolverPackage matSolverPackage,
PetscReal rtol = PETSC_DEFAULT,
PetscReal atol = PETSC_DEFAULT,
PetscInt maxIt = PETSC_DEFAULT);
void setSolver(KSP ksp, void setSolver(KSP ksp,
const char* kspPrefix, const char* kspPrefix,
KSPType kspType, KSPType kspType,
......
...@@ -47,6 +47,9 @@ namespace AMDiS { ...@@ -47,6 +47,9 @@ namespace AMDiS {
pressureComponent(-1), pressureComponent(-1),
pressureNullSpace(true), pressureNullSpace(true),
useOldInitialGuess(false), useOldInitialGuess(false),
velocitySolutionMode(0),
massSolutionMode(0),
laplaceSolutionMode(0),
massMatrixSolver(NULL), massMatrixSolver(NULL),
laplaceMatrixSolver(NULL), laplaceMatrixSolver(NULL),
nu(NULL), nu(NULL),
...@@ -65,6 +68,15 @@ namespace AMDiS { ...@@ -65,6 +68,15 @@ namespace AMDiS {
Parameters::get(initFileStr + "->navierstokes->use old initial guess", Parameters::get(initFileStr + "->navierstokes->use old initial guess",
useOldInitialGuess); useOldInitialGuess);
Parameters::get(initFileStr + "->navierstokes->velocity solver",
velocitySolutionMode);
Parameters::get(initFileStr + "->navierstokes->mass solver",
massSolutionMode);
Parameters::get(initFileStr + "->navierstokes->laplace solver",
laplaceSolutionMode);
} }
...@@ -138,7 +150,19 @@ namespace AMDiS { ...@@ -138,7 +150,19 @@ namespace AMDiS {
KSP kspSchur = subKsp[1]; KSP kspSchur = subKsp[1];
PetscFree(subKsp); PetscFree(subKsp);
petsc_helper::setSolver(kspVelocity, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1); switch (velocitySolutionMode) {
case 0:
petsc_helper::setSolver(kspVelocity, "",
KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
break;
case 1:
petsc_helper::setSolverWithLu(kspVelocity, "", KSPPREONLY,
PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1);
break;
default:
ERROR_EXIT("No velocity solution mode %d available!\n", velocitySolutionMode);
}
KSPSetType(kspSchur, KSPPREONLY); KSPSetType(kspSchur, KSPPREONLY);
PC pcSub; PC pcSub;
...@@ -155,12 +179,17 @@ namespace AMDiS { ...@@ -155,12 +179,17 @@ namespace AMDiS {
const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent]; const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace); DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
{
Operator massOp(pressureFeSpace, pressureFeSpace); Operator massOp(pressureFeSpace, pressureFeSpace);
ZeroOrderTerm *massTerm = NULL;
if (!phase) if (!phase)
massOp.addTerm(new Simple_ZOT); massTerm = new Simple_ZOT;
else else
massOp.addTerm(new VecAtQP_ZOT(phase, &idFct)); massTerm = new VecAtQP_ZOT(phase, &idFct);
massOp.addTerm(massTerm);
massMatrix.assembleOperator(massOp); massMatrix.assembleOperator(massOp);
delete massTerm;
}
massMatrixSolver = createSubSolver(pressureComponent, "mass_"); massMatrixSolver = createSubSolver(pressureComponent, "mass_");
massMatrixSolver->fillPetscMatrix(&massMatrix); massMatrixSolver->fillPetscMatrix(&massMatrix);
...@@ -168,12 +197,17 @@ namespace AMDiS { ...@@ -168,12 +197,17 @@ namespace AMDiS {
// === Laplace matrix solver === // === Laplace matrix solver ===
DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace); DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
{
Operator laplaceOp(pressureFeSpace, pressureFeSpace); Operator laplaceOp(pressureFeSpace, pressureFeSpace);
SecondOrderTerm *laplaceTerm = NULL;
if (!phase) if (!phase)
laplaceOp.addTerm(new Simple_SOT); laplaceTerm = new Simple_SOT;
else else
laplaceOp.addTerm(new VecAtQP_SOT(phase, &idFct)); laplaceTerm = new VecAtQP_SOT(phase, &idFct);
laplaceOp.addTerm(laplaceTerm);
laplaceMatrix.assembleOperator(laplaceOp); laplaceMatrix.assembleOperator(laplaceOp);
delete laplaceTerm;
}
laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_")); laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_"));
laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix); laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
...@@ -190,47 +224,60 @@ namespace AMDiS { ...@@ -190,47 +224,60 @@ namespace AMDiS {
vz.interpol(solution->getDOFVector(2)); vz.interpol(solution->getDOFVector(2));
DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace); DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
{
Operator conDifOp(pressureFeSpace, pressureFeSpace); Operator conDifOp(pressureFeSpace, pressureFeSpace);
ZeroOrderTerm *conDif0 = NULL;
SecondOrderTerm *conDif1 = NULL;
FirstOrderTerm *conDif2 = NULL, *conDif3 = NULL, *conDif4 = NULL;
if (!phase) { if (!phase) {
MSG("INIT WITHOUT PHASE!\n"); MSG("INIT WITHOUT PHASE!\n");
Simple_ZOT *conDif0 = new Simple_ZOT(*invTau);
conDif0 = new Simple_ZOT(*invTau);
conDifOp.addTerm(conDif0); conDifOp.addTerm(conDif0);
Simple_SOT *conDif1 = new Simple_SOT(*nu); conDif1 = new Simple_SOT(*nu);
conDifOp.addTerm(conDif1); conDifOp.addTerm(conDif1);
VecAtQP_FOT *conDif2 = new VecAtQP_FOT(&vx, &idFct, 0); conDif2 = new VecAtQP_FOT(&vx, &idFct, 0);
conDifOp.addTerm(conDif2, GRD_PHI); conDifOp.addTerm(conDif2, GRD_PHI);
VecAtQP_FOT *conDif3 = new VecAtQP_FOT(&vy, &idFct, 1); conDif3 = new VecAtQP_FOT(&vy, &idFct, 1);
conDifOp.addTerm(conDif3, GRD_PHI); conDifOp.addTerm(conDif3, GRD_PHI);
if (dim == 3) { if (dim == 3) {
VecAtQP_FOT *conDif4 = new VecAtQP_FOT(&vz, &idFct, 2); conDif4 = new VecAtQP_FOT(&vz, &idFct, 2);
conDifOp.addTerm(conDif4, GRD_PHI); conDifOp.addTerm(conDif4, GRD_PHI);
} }
} else { } else {
MSG("INIT WITH PHASE!\n"); MSG("INIT WITH PHASE: %e %e\n", *nu, *invTau);
vp.interpol(phase); vp.interpol(phase);
VecAtQP_ZOT *conDif0 = new VecAtQP_ZOT(&vp, new MultConstFct(*invTau)); conDif0 = new VecAtQP_ZOT(&vp, new MultConstFct(*invTau));
conDifOp.addTerm(conDif0); conDifOp.addTerm(conDif0);
conDif1 = new VecAtQP_SOT(&vp, new MultConstFct(*nu));
VecAtQP_SOT *conDif1 = new VecAtQP_SOT(&vp, new MultConstFct(*nu));
conDifOp.addTerm(conDif1); conDifOp.addTerm(conDif1);
conDif2 = new Vec2AtQP_FOT(&vx, &vp, new Multiplier3(), 0);
Vec2AtQP_FOT *conDif2 = new Vec2AtQP_FOT(&vx, &vp, new Multiplier3(), 0);
conDifOp.addTerm(conDif2, GRD_PHI); conDifOp.addTerm(conDif2, GRD_PHI);
Vec2AtQP_FOT *conDif3 = new Vec2AtQP_FOT(&vy, &vp, new Multiplier3(), 1); conDif3 = new Vec2AtQP_FOT(&vy, &vp, new Multiplier3(), 1);
conDifOp.addTerm(conDif3, GRD_PHI); conDifOp.addTerm(conDif3, GRD_PHI);
if (dim == 3) { if (dim == 3) {
Vec2AtQP_FOT *conDif4 = new Vec2AtQP_FOT(&vz, &vp, new Multiplier3(), 2); conDif4 = new Vec2AtQP_FOT(&vz, &vp, new Multiplier3(), 2);
conDifOp.addTerm(conDif4, GRD_PHI); conDifOp.addTerm(conDif4, GRD_PHI);
} }
} }
conDifMatrix.assembleOperator(conDifOp); conDifMatrix.assembleOperator(conDifOp);
delete conDif0;
delete conDif1;
delete conDif2;
delete conDif3;
if (dim == 3)
delete conDif4;
}
conDifMatrixSolver = createSubSolver(pressureComponent, "condif_"); conDifMatrixSolver = createSubSolver(pressureComponent, "condif_");
conDifMatrixSolver->fillPetscMatrix(&conDifMatrix); conDifMatrixSolver->fillPetscMatrix(&conDifMatrix);
...@@ -241,11 +288,52 @@ namespace AMDiS { ...@@ -241,11 +288,52 @@ namespace AMDiS {
matShellContext.kspLaplace = laplaceMatrixSolver->getSolver(); matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
matShellContext.matConDif = conDifMatrixSolver->getMatInterior(); matShellContext.matConDif = conDifMatrixSolver->getMatInterior();
petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 2); switch (massSolutionMode) {
petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1); case 0:
petsc_helper::setSolver(matShellContext.kspMass, "mass_",
KSPCG, PCJACOBI, 0.0, 1e-14, 2);
break;
case 1:
petsc_helper::setSolverWithLu(matShellContext.kspMass, "mass_", KSPRICHARDSON,
PCLU, MATSOLVERMUMPS, 0.0, 1e-14, 1);
break;
default:
ERROR_EXIT("No mass solution mode %d available!\n", massSolutionMode);
}
switch (laplaceSolutionMode) {
case 0:
petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_",
KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
break;
case 1:
petsc_helper::setSolverWithLu(matShellContext.kspLaplace, "mass_",
KSPRICHARDSON, PCLU, MATSOLVERMUMPS,
0.0, 1e-14, 1);
break;
default:
ERROR_EXIT("No laplace solution mode %d available!\n", laplaceSolutionMode);
}
setConstantNullSpace(matShellContext.kspLaplace); setConstantNullSpace(matShellContext.kspLaplace);
} }
void PetscSolverNavierStokes::exitPreconditioner(PC pc)
{
FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");
massMatrixSolver->destroyMatrixData();
laplaceMatrixSolver->destroyMatrixData();
conDifMatrixSolver->destroyMatrixData();
delete massMatrixSolver;
massMatrixSolver = NULL;
delete laplaceMatrixSolver;
laplaceMatrixSolver = NULL;
delete conDifMatrixSolver;
conDifMatrixSolver = NULL;
}
} }
...@@ -142,6 +142,8 @@ namespace AMDiS { ...@@ -142,6 +142,8 @@ namespace AMDiS {
void initPreconditioner(PC pc); void initPreconditioner(PC pc);
void exitPreconditioner(PC pc);
private: private:
int pressureComponent; int pressureComponent;
...@@ -150,6 +152,15 @@ namespace AMDiS { ...@@ -150,6 +152,15 @@ namespace AMDiS {
/// If true, old solution is used for initial guess in solver phase. /// If true, old solution is used for initial guess in solver phase.
bool useOldInitialGuess; bool useOldInitialGuess;
/// 0: approximate solve 1: direct solver
int velocitySolutionMode;
/// 0: approximate solve 1: direct solver
int massSolutionMode;
/// 0: approximate solve 1: direct solver
int laplaceSolutionMode;
PetscSolver *massMatrixSolver, *laplaceMatrixSolver, *conDifMatrixSolver; PetscSolver *massMatrixSolver, *laplaceMatrixSolver, *conDifMatrixSolver;
NavierStokesSchurData matShellContext; NavierStokesSchurData matShellContext;
......
Supports Markdown
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