Commit 76f3f1db authored by Thomas Witkowski's avatar Thomas Witkowski

Fixed some problem with higher order DOFs in parallel computations.

parent 7b91a0b8
......@@ -414,49 +414,86 @@ namespace AMDiS {
break;
case FACE:
{
BoundaryObject nextBound = bound;
nextBound.elType = (bound.elType + 1) % 3;
if (child[0]) {
int childFace0 = sideOfChild[bound.elType][0][bound.ithObj];
int childFace1 = sideOfChild[bound.elType][1][bound.ithObj];
TEST_EXIT(childFace0 != -1 || childFace1 != -1)
BoundaryObject nextBound0 = bound, nextBound1 = bound;
prepareNextBound(nextBound0, 0);
prepareNextBound(nextBound1, 1);
TEST_EXIT_DBG(nextBound0.ithObj != -1 || nextBound1.ithObj != 1)
("No new face for child elements!\n");
TEST_EXIT_DBG(!bound.reverseMode)("Not yet implemented!\n");
if (childFace0 != -1) {
nextBound.ithObj = childFace0;
child[0]->getHigherOrderDofs(feSpace, nextBound, dofs,
if (nextBound0.ithObj != -1)
child[0]->getHigherOrderDofs(feSpace, nextBound0, dofs,
baseDofPtr, dofGeoIndex);
}
if (childFace1 != -1) {
nextBound.ithObj = childFace1;
child[1]->getHigherOrderDofs(feSpace, nextBound, dofs,
if (nextBound1.ithObj != -1)
child[1]->getHigherOrderDofs(feSpace, nextBound1, dofs,
baseDofPtr, dofGeoIndex);
}
} else {
// === On leaf elements iterate over all DOFs and collect the ===
// === right ones. ===
ElementDofIterator elDofIter(feSpace, true);
elDofIter.reset(this);
if (baseDofPtr) {
do {
if (elDofIter.getCurrentPos() >= 1 &&
elDofIter.getCurrentElementPos() == bound.ithObj)
bool next = true;
do {
bool addDof = false;
switch (elDofIter.getCurrentPos()) {
case 0:
// VERTEX is never a higher order DOF
addDof = false;
break;
case 1:
// On EDGE, first check whether this edge is part of the FACE and
// then check if this EDGE is not inside of excludedSubstructures.
{
int localEdge = elDofIter.getCurrentElementPos();
if (edgeOfFace[bound.ithObj][0] == localEdge ||
edgeOfFace[bound.ithObj][1] == localEdge ||
edgeOfFace[bound.ithObj][2] == localEdge) {
addDof = true;
for (ExcludeList::iterator it = bound.excludedSubstructures.begin();
it != bound.excludedSubstructures.end(); ++it)
if (it->first == EDGE && it->second == localEdge)
addDof = false;
} else {
addDof = false;
}
}
break;
case 2:
// On FACE, check for the right one.
if (elDofIter.getCurrentElementPos() == bound.ithObj)
addDof = true;
else
addDof = false;
break;
default:
ERROR_EXIT("Should not happen!\n");
}
if (addDof) {
if (baseDofPtr)
dofs.push_back(elDofIter.getBaseDof());
if (dofGeoIndex != NULL)
dofGeoIndex->push_back(elDofIter.getPosIndex());
} while (elDofIter.nextStrict());
} else {
do {
if (elDofIter.getCurrentPos() >= 1 &&
elDofIter.getCurrentElementPos() == bound.ithObj)
else
dofs.push_back(elDofIter.getDofPtr());
if (dofGeoIndex != NULL)
dofGeoIndex->push_back(elDofIter.getPosIndex());
} while (elDofIter.next());
}
if (dofGeoIndex != NULL)
dofGeoIndex->push_back(elDofIter.getPosIndex());
}
next = (baseDofPtr ? elDofIter.nextStrict() : elDofIter.next());
} while (next);
}
}
break;
......@@ -468,8 +505,7 @@ namespace AMDiS {
void Tetrahedron::getSubBoundary(BoundaryObject bound,
vector<BoundaryObject> &subBound) const
{
}
{}
void Tetrahedron::prepareNextBound(BoundaryObject &bound, int ithChild) const
......@@ -479,10 +515,9 @@ namespace AMDiS {
switch (bound.subObj) {
case FACE:
for (ExcludeList::iterator it = bound.excludedSubstructures.begin();
it != bound.excludedSubstructures.end(); ++it) {
it != bound.excludedSubstructures.end(); ++it)
if (it->first == EDGE && it->second != -1)
it->second = edgeOfChild[bound.elType][ithChild][it->second];
}
bound.ithObj = sideOfChild[bound.elType][ithChild][bound.ithObj];
bound.elType = (bound.elType + 1) % 3;
......
......@@ -161,8 +161,6 @@ namespace AMDiS {
{
FUNCNAME("VtkWriter::writeFile()");
return;
DataCollector<> dc(values->getFeSpace(), values);
vector<DataCollector<>*> dcList(0);
dcList.push_back(&dc);
......
......@@ -324,7 +324,7 @@ namespace AMDiS {
if (globalRefinement > 0) {
refineManager->globalRefine(mesh, globalRefinement);
updateLocalGlobalNumbering();
#if (DEBUG != 0)
......
......@@ -482,8 +482,10 @@ namespace AMDiS {
DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
for (it.reset(); !it.end(); ++it) {
coordsToIndex[(*it)] = (*(pdb.dofMaps[0]))[feSpace][it.getDOFIndex()].global;
// MSG(" CHECK FOR DOF %d AT COORDS %f %f %f\n",
// coordsToIndex[(*it)], (*it)[0], (*it)[1], (*it)[2]);
// MSG(" CHECK FOR DOF %d/%d AT COORDS %f %f %f\n",
// it.getDOFIndex(),
// coordsToIndex[(*it)],
// (*it)[0], (*it)[1], (pdb.mesh->getDim() == 3 ? (*it)[2] : 0.0));
}
StdMpi<CoordsIndexMap> stdMpi(pdb.mpiComm, true);
......
......@@ -26,6 +26,16 @@ namespace AMDiS {
PCShellGetContext(pc, &ctx);
NavierStokesSchurData* data = static_cast<NavierStokesSchurData*>(ctx);
// Project out constant null space
{
int vecSize;
VecGetSize(x, &vecSize);
PetscScalar vecSum;
VecSum(x, &vecSum);
vecSum = vecSum / static_cast<PetscScalar>(-1.0 * vecSize);
VecShift(x, vecSum);
}
KSPSolve(data->kspLaplace, x, y);
MatMult(data->matConDif, y, x);
KSPSolve(data->kspMass, x, y);
......@@ -35,6 +45,7 @@ namespace AMDiS {
PetscSolverNavierStokes::PetscSolverNavierStokes(string name)
: PetscSolverGlobalMatrix(name),
pressureComponent(-1),
pressureNullSpace(true),
massMatrixSolver(NULL),
laplaceMatrixSolver(NULL),
nu(NULL),
......@@ -44,11 +55,13 @@ namespace AMDiS {
{
Parameters::get(initFileStr + "->navierstokes->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");
Parameters::get(initFileStr + "->navierstokes->pressure null space",
pressureNullSpace);
TEST_EXIT(pressureNullSpace)("This is not yet tested, may be wrong!\n");
}
......@@ -60,10 +73,11 @@ namespace AMDiS {
KSPCreate(mpiCommGlobal, &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);
petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 10000);
// Create null space information.
setConstantNullSpace(ksp, pressureComponent, true);
if (pressureNullSpace)
setConstantNullSpace(ksp, pressureComponent, true);
}
......@@ -98,7 +112,7 @@ namespace AMDiS {
KSP kspSchur = subKsp[1];
PetscFree(subKsp);
petsc_helper::setSolver(kspVelocity, "velocity_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
petsc_helper::setSolver(kspVelocity, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
KSPSetType(kspSchur, KSPPREONLY);
PC pcSub;
......@@ -107,7 +121,8 @@ namespace AMDiS {
PCShellSetApply(pcSub, pcSchurShell);
PCShellSetContext(pcSub, &matShellContext);
setConstantNullSpace(kspSchur);
if (pressureNullSpace)
setConstantNullSpace(kspSchur);
// === Mass matrix solver ===
......@@ -115,10 +130,10 @@ namespace AMDiS {
const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
Operator massOp(pressureFeSpace, pressureFeSpace);
// if (!phase)
if (!phase)
massOp.addTerm(new Simple_ZOT);
// else
// massOp.addTerm(new VecAtQP_ZOT(phase, &idFct));
else
massOp.addTerm(new VecAtQP_ZOT(phase, &idFct));
massMatrix.assembleOperator(massOp);
massMatrixSolver = createSubSolver(pressureComponent, "mass_");
massMatrixSolver->fillPetscMatrix(&massMatrix);
......@@ -190,29 +205,9 @@ namespace AMDiS {
matShellContext.matConDif = conDifMatrixSolver->getMatInterior();
petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
// petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
KSPSetType(matShellContext.kspLaplace, KSPRICHARDSON);
KSPSetTolerances(matShellContext.kspLaplace, 0.0, 1e-14, PETSC_DEFAULT, 1);
KSPSetOptionsPrefix(matShellContext.kspLaplace, "laplace_");
KSPSetFromOptions(matShellContext.kspLaplace);
MatNullSpace matNullSpace;
MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
KSPSetNullSpace(matShellContext.kspLaplace, matNullSpace);
MatNullSpaceDestroy(&matNullSpace);
{
PC pc;
KSPGetPC(matShellContext.kspLaplace, &pc);
PCSetType(pc, PCHYPRE);
PCSetFromOptions(pc);
}
// setConstantNullSpace(matShellContext.kspLaplace);
petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
setConstantNullSpace(matShellContext.kspLaplace);
}
......
......@@ -30,7 +30,9 @@ namespace AMDiS {
using namespace std;
struct NavierStokesSchurData {
KSP kspMass, kspLaplace;
KSP kspMass;
KSP kspLaplace;
Mat matConDif;
};
......@@ -128,6 +130,8 @@ namespace AMDiS {
private:
int pressureComponent;
bool pressureNullSpace;
PetscSolver *massMatrixSolver, *laplaceMatrixSolver, *conDifMatrixSolver;
NavierStokesSchurData matShellContext;
......
......@@ -109,12 +109,12 @@ public:
// ---> for test purposes: print result of update counting
void printUpdateCntr()
{
cout << "\n";
cout << "\tUpdate statistic: \n";
cout << "\t-----------------\n";
cout << "\tcalculated updates: " << calcUpdate_Cntr << "\n";
cout << "\tset updates: " << setUpdate_Cntr << "\n";
cout << "\n";
/* cout << "\n"; */
/* cout << "\tUpdate statistic: \n"; */
/* cout << "\t-----------------\n"; */
/* cout << "\tcalculated updates: " << calcUpdate_Cntr << "\n"; */
/* cout << "\tset updates: " << setUpdate_Cntr << "\n"; */
/* cout << "\n"; */
}
// ---> end: for test purposes
......
Markdown is supported
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