Commit 7f81c570 authored by Thomas Witkowski's avatar Thomas Witkowski

Fixed some issues for removing periodic boundary conditions.

parent 35eda406
......@@ -505,6 +505,7 @@ namespace AMDiS {
{
FUNCNAME("ProblemStat::createSolver()");
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
// === create solver ===
string solverType("0");
string initFileStr = name + "->solver";
......@@ -515,6 +516,7 @@ namespace AMDiS {
solverCreator->setName(initFileStr);
solver = solverCreator->create();
solver->initParameters();
#endif
}
......
......@@ -48,9 +48,7 @@ namespace AMDiS {
DegreeOfFreedom **dof = macroInfo->dof;
// === read periodic data =================================
if (periodicFile != "") {
WARNING("Periodic boundaries may lead to errors in small meshes if element neighbours are not set!\n");
if (periodicFile != "") {
FILE *file = fopen(periodicFile.c_str(), "r");
TEST_EXIT(file)("can't open file %s\n", periodicFile.c_str());
......
......@@ -105,16 +105,15 @@ namespace AMDiS {
periodicEdgeAssoc[edge0].insert(edge1);
// Add both vertices of the edge to be periodic.
periodicVertices[make_pair(edge0.first, edge1.first)] = boundaryType;
periodicVertices[make_pair(edge0.second, edge1.second)] = boundaryType;
DegreeOfFreedom mappedDof0 =
mesh->getPeriodicAssociations(boundaryType)[edge0.first];
DegreeOfFreedom mappedDof1 =
mesh->getPeriodicAssociations(boundaryType)[edge0.second];
periodicVertices[make_pair(edge0.first, mappedDof0)] = boundaryType;
periodicVertices[make_pair(edge0.second, mappedDof1)] = boundaryType;
periodicDofAssoc[edge0.first].insert(boundaryType);
periodicDofAssoc[edge0.second].insert(boundaryType);
TEST_EXIT_DBG(edge0.first ==
mesh->getPeriodicAssociations(boundaryType)[edge1.first] &&
edge0.second ==
mesh->getPeriodicAssociations(boundaryType)[edge1.second])
("Should not happen!\n");
periodicDofAssoc[edge0.second].insert(boundaryType);
}
}
break;
......@@ -134,12 +133,16 @@ namespace AMDiS {
periodicFaces[make_pair(face0, face1)] = elInfo->getBoundary(i);
/// Add all three vertices of the face to be periodic.
periodicVertices[make_pair(face0.get<0>(), face1.get<0>())] =
boundaryType;
periodicVertices[make_pair(face0.get<1>(), face1.get<1>())] =
boundaryType;
periodicVertices[make_pair(face0.get<2>(), face1.get<2>())] =
boundaryType;
DegreeOfFreedom mappedDof0 =
mesh->getPeriodicAssociations(boundaryType)[face0.get<0>()];
DegreeOfFreedom mappedDof1 =
mesh->getPeriodicAssociations(boundaryType)[face0.get<1>()];
DegreeOfFreedom mappedDof2 =
mesh->getPeriodicAssociations(boundaryType)[face0.get<2>()];
periodicVertices[make_pair(face0.get<0>(), mappedDof0)] = boundaryType;
periodicVertices[make_pair(face0.get<1>(), mappedDof1)] = boundaryType;
periodicVertices[make_pair(face0.get<2>(), mappedDof2)] = boundaryType;
periodicDofAssoc[face0.get<0>()].insert(boundaryType);
periodicDofAssoc[face0.get<1>()].insert(boundaryType);
......
......@@ -309,24 +309,32 @@ namespace AMDiS {
return faceElements[face];
}
/// Returns a map that maps to each rank all macro elements in this rank that
/// have a given vertex DOF in common.
/// Returns, for a given vertex, a map that maps from rank numbers to
/// element data objects, which identify on the rank the element which
/// contains this vertex. If more than one element in one subdomain contains
/// the vertex, the element with the highest element index is given. If the
/// vertex is not contained in a rank's subdomain, it will not be considered
/// in this mapping.
map<int, ElementObjectData>& getElementsInRank(DegreeOfFreedom vertex)
{
return vertexInRank[vertex];
}
/// Returns a map that maps to each rank all macro elements in this rank that
/// have a given edge in common.
/// Returns, for a given edge, a map that maps from rank numbers to
/// element data objects, which identify on the rank the element which
/// contains this edge. If more than one element in one subdomain contains
/// the edge, the element with the highest element index is given. If the
/// edge is not contained in a rank's subdomain, it will not be considered
/// in this mapping.
map<int, ElementObjectData>& getElementsInRank(DofEdge edge)
{
return edgeInRank[edge];
}
/// Returns a map that maps to each rank all macro elements in this rank that
/// have a given face in common.
/// Returns, for a given face, a map that maps from rank numbers to
/// element data objects, which identify on the rank the element which
/// contains this face. If the face is not contained in a rank's subdomain,
/// it will not be considered in this mapping.
map<int, ElementObjectData>& getElementsInRank(DofFace face)
{
return faceInRank[face];
......
This diff is collapsed.
......@@ -86,12 +86,32 @@ namespace AMDiS {
void deserialize(istream &in, Mesh *mesh);
private:
/// In this function, we put some verification procedures to check for
/// consistency if the created interior boundary data. This function is
/// enabled even for optimized mode to check all possible meshes. Thus,
/// everything herein should be fast. If the code is stable, we can think
/// to remove this function.
void verifyBoundary();
AtomicBoundary& getNewOwn(int rank);
AtomicBoundary& getNewOther(int rank);
AtomicBoundary& getNewPeriodic(int rank);
/// Checks whether the given boundary is already a other boundary with
/// given rank number. Returns true, if the bound is already in the other
/// boundary database.
bool checkOther(AtomicBoundary& bound, int rank);
/// Removes the given boundary object from all owned boundaries. Returns
/// true, if at least one object has been removed, otherwise false.
bool removeOwn(BoundaryObject& bound);
/// Removes the given boundary object from all other boundaries. Returns
/// true, if at least one object has been removed, otherwise false.
bool removeOther(BoundaryObject& bound);
void serialize(ostream &out, RankToBoundMap& boundary);
void deserialize(istream &in,
......
......@@ -68,8 +68,6 @@ namespace AMDiS {
return (*dof1 < *dof2);
}
bool MeshDistributor::sebastianMode = false;
MeshDistributor::MeshDistributor()
: problemStat(0),
initialized(false),
......@@ -1696,7 +1694,7 @@ namespace AMDiS {
dofMaps[i]->update();
#if (DEBUG != 0)
dofMaps[i]->printInfo();
// dofMaps[i]->printInfo();
#endif
}
......
......@@ -290,7 +290,7 @@ namespace AMDiS {
void setBoundaryDofRequirement(Flag flag)
{
createBoundaryDofFlag = flag;
createBoundaryDofFlag |= flag;
}
BoundaryDofInfo& getBoundaryDofInfo(const FiniteElemSpace *feSpace,
......@@ -584,8 +584,6 @@ namespace AMDiS {
vector<ParallelDofMapping*> dofMaps;
public:
static bool sebastianMode;
/// The boundary DOFs are sorted by subobject entities, i.e., first all
/// face DOFs, edge DOFs and to the last vertex DOFs will be set to
/// communication structure vectors, \ref sendDofs and \ref recvDofs.
......
......@@ -47,6 +47,7 @@ namespace AMDiS {
*/
class ParallelCoarseSpaceMatVec {
public:
/// Constructor
ParallelCoarseSpaceMatVec();
/// Set parallel DOF mapping for the interior DOFs.
......@@ -54,6 +55,12 @@ namespace AMDiS {
{
interiorMap = interiorDofs;
}
/// Returns the parallel DOF mapping for the interior DOFs.
ParallelDofMapping* getDofMapping()
{
return interiorMap;
}
/** \brief
* Sets the coarse space for all or a specific component.
......
......@@ -40,8 +40,8 @@ namespace AMDiS {
virtual ~ParallelProblemStatBase() {}
void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
bool assembleMatrix,
bool assembleVector);
bool assembleMatrix = true,
bool assembleVector = true);
void initialize(Flag initFlag,
ProblemStatSeq *adoptProblem = NULL,
......
......@@ -33,15 +33,16 @@ namespace AMDiS {
FUNCNAME("PetscProblemStat::PetscProblemStat()");
string tmp("");
Parameters::get("parallel->solver", tmp);
string initFileStr = nameStr + "->solver";
Parameters::get(initFileStr, tmp);
if (tmp == "petsc-schur") {
petscSolver = new PetscSolverSchur();
} else if (tmp == "petsc-feti") {
petscSolver = new PetscSolverFeti();
petscSolver = new PetscSolverFeti(initFileStr);
} else if (tmp == "petsc-block") {
petscSolver = new PetscSolverGlobalBlockMatrix();
} else if (tmp == "petsc" || tmp == "") {
} else if (tmp == "petsc") {
petscSolver = new PetscSolverGlobalMatrix();
} else if (tmp == "bddcml") {
#ifdef HAVE_BDDC_ML
......
......@@ -28,8 +28,9 @@ namespace AMDiS {
using namespace std;
PetscSolverFeti::PetscSolverFeti()
PetscSolverFeti::PetscSolverFeti(string name)
: PetscSolver(),
initFileStr(name),
primalDofMap(COMPONENT_WISE),
dualDofMap(COMPONENT_WISE),
interfaceDofMap(COMPONENT_WISE),
......@@ -51,8 +52,8 @@ namespace AMDiS {
FUNCNAME("PetscSolverFeti::PetscSolverFeti()");
string preconditionerName = "";
Parameters::get("parallel->solver->precon", preconditionerName);
if (preconditionerName == "" || preconditionerName == "none") {
Parameters::get(initFileStr + "->left precon", preconditionerName);
if (preconditionerName == "" || preconditionerName == "no") {
MSG("Create FETI-DP solver with no preconditioner!\n");
fetiPreconditioner = FETI_NONE;
} else if (preconditionerName == "dirichlet") {
......@@ -66,20 +67,34 @@ namespace AMDiS {
preconditionerName.c_str());
}
Parameters::get("parallel->feti->schur primal solver", schurPrimalSolver);
preconditionerName = "";
Parameters::get(initFileStr + "->right precon", preconditionerName);
if (preconditionerName != "" && preconditionerName != "no") {
ERROR_EXIT("FETI-DP does not support right preconditioning! (parameter \"%s->right precon\" has value \"%s\")\n",
initFileStr.c_str(), preconditionerName.c_str());
}
Parameters::get(initFileStr + "->feti->schur primal solver", schurPrimalSolver);
TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1)
("Wrong solver \"%d\"for the Schur primal complement!\n",
schurPrimalSolver);
Parameters::get(initFileStr + "->feti->stokes mode", stokesMode);
if (stokesMode) {
Parameters::get(initFileStr + "->feti->pressure component", pressureComponent);
TEST_EXIT(pressureComponent >= 0)
("FETI-DP in Stokes mode, no pressure component defined!\n");
}
Parameters::get(initFileStr + "->feti->augmented lagrange", augmentedLagrange);
Parameters::get(initFileStr + "->feti->symmetric", isSymmetric);
Parameters::get("parallel->multi level test", multiLevelTest);
if (multiLevelTest)
meshLevel = 1;
Parameters::get("parallel->print timings", printTimings);
Parameters::get("parallel->feti->augmented lagrange", augmentedLagrange);
Parameters::get("parallel->feti->symmetric", isSymmetric);
}
......@@ -92,14 +107,6 @@ namespace AMDiS {
MeshLevelData& levelData = meshDistributor->getMeshLevelData();
stokesMode = false;
Parameters::get("parallel->feti->stokes mode", stokesMode);
if (stokesMode) {
Parameters::get("parallel->feti->pressure component", pressureComponent);
TEST_EXIT(pressureComponent >= 0)
("FETI-DP in Stokes mode, no pressure component defined!\n");
}
if (subdomain == NULL) {
subdomain = new PetscSolverGlobalMatrix();
subdomain->setSymmetric(isSymmetric);
......@@ -1149,8 +1156,19 @@ namespace AMDiS {
KSPSetType(ksp_feti, KSPGMRES);
KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
KSPSetFromOptions(ksp_feti);
if (stokesMode)
KSPMonitorSet(ksp_feti, KSPMonitorFetiStokes, &fetiKspData, PETSC_NULL);
// === Set KSP monitor. ===
bool monitor = false;
Parameters::get(initFileStr + "->feti->monitor", monitor);
if (monitor) {
if (stokesMode)
KSPMonitorSet(ksp_feti, KSPMonitorFetiStokes, &fetiKspData, PETSC_NULL);
else
KSPMonitorSet(ksp_feti, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
}
// === Create null space objects. ===
......@@ -1257,64 +1275,53 @@ namespace AMDiS {
}
if (stokesMode) {
// === Create H2 vec ===
// === Create mass matrix solver ===
const FiniteElemSpace *pressureFeSpace =
componentSpaces[pressureComponent];
DOFVector<double> tmpvec(pressureFeSpace, "tmpvec");
createH2vec(tmpvec);
interfaceDofMap.createVec(fetiInterfaceLumpedPreconData.h2vec);
DofMap& m = interfaceDofMap[pressureComponent].getMap();
for (DofMap::iterator it = m.begin(); it != m.end(); ++it) {
if (dofMap[pressureFeSpace].isRankDof(it->first)) {
int index = interfaceDofMap.getMatIndex(pressureComponent, it->first);
VecSetValue(fetiInterfaceLumpedPreconData.h2vec,
index, tmpvec[it->first], INSERT_VALUES);
}
}
VecAssemblyBegin(fetiInterfaceLumpedPreconData.h2vec);
VecAssemblyEnd(fetiInterfaceLumpedPreconData.h2vec);
// === Create mass matrix solver ===
if (!massMatrixSolver) {
ParallelDofMapping *massMapping = new ParallelDofMapping(COMPONENT_WISE);
// Create parallel DOF mapping in pressure space.
ParallelDofMapping *massMapping = NULL;
if (massMatrixSolver) {
massMapping = massMatrixSolver->getDofMapping();
} else {
massMapping =
new ParallelDofMapping(COMPONENT_WISE);
massMapping->init(meshDistributor->getMeshLevelData(),
pressureFeSpace, pressureFeSpace);
massMapping->setDofComm(meshDistributor->getDofComm());
massMapping->setMpiComm(meshDistributor->getMeshLevelData().getMpiComm(0), 0);
massMapping->setComputeMatIndex(true);
(*massMapping)[0] = interfaceDofMap[pressureComponent];
massMapping->update();
DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
Operator op(pressureFeSpace, pressureFeSpace);
Simple_ZOT zot;
op.addTerm(&zot);
massMatrix.assembleOperator(op);
}
(*massMapping)[0] = interfaceDofMap[pressureComponent];
massMapping->update();
DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
Operator op(pressureFeSpace, pressureFeSpace);
Simple_ZOT zot;
op.addTerm(&zot);
massMatrix.assembleOperator(op);
if (!massMatrixSolver) {
massMatrixSolver = new PetscSolverGlobalMatrix;
massMatrixSolver->setKspPrefix("mass_");
massMatrixSolver->setMeshDistributor(meshDistributor,
mpiCommGlobal,
mpiCommLocal);
massMatrixSolver->setDofMapping(massMapping);
massMatrixSolver->fillPetscMatrix(&massMatrix);
int r, c;
MatGetSize(massMatrixSolver->getMatInterior(), &r, &c);
MatInfo info;
MatGetInfo(massMatrixSolver->getMatInterior(), MAT_GLOBAL_SUM, &info);
MSG("MASS MAT INFO: size = %d x %d nnz = %d\n",
r, c,
static_cast<int>(info.nz_used));
fetiInterfaceLumpedPreconData.ksp_mass =
massMatrixSolver->getSolver();
}
massMatrixSolver->fillPetscMatrix(&massMatrix);
int r, c;
MatGetSize(massMatrixSolver->getMatInterior(), &r, &c);
MatInfo info;
MatGetInfo(massMatrixSolver->getMatInterior(), MAT_GLOBAL_SUM, &info);
MSG("MASS MAT INFO: size = %d x %d nnz = %d\n",
r, c, static_cast<int>(info.nz_used));
fetiInterfaceLumpedPreconData.ksp_mass = massMatrixSolver->getSolver();
// === Create tmp vectors ===
......@@ -2374,59 +2381,4 @@ namespace AMDiS {
MeshDistributor::globalMeshDistributor->synchVector(vec);
}
void PetscSolverFeti::createH2vec(DOFVector<double> &vec)
{
FUNCNAME("PetscSolverFeti::createH2vec()");
vec.set(0.0);
Mesh *mesh = vec.getFeSpace()->getMesh();
TEST_EXIT(mesh->getDim() == 2)("This works only in 2D!\n");
int n0 = vec.getFeSpace()->getAdmin()->getNumberOfPreDofs(VERTEX);
TraverseStack stack;
ElInfo *elInfo =
stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
while (elInfo) {
Element *el = elInfo->getElement();
std::vector<double> length(3, 0.0);
double sumLength = 0.0;
for (int i = 0; i < 3; i++) {
int idx0 = el->getVertexOfEdge(i, 0);
int idx1 = el->getVertexOfEdge(i, 1);
DegreeOfFreedom dof0 = el->getDof(idx0, n0);
DegreeOfFreedom dof1 = el->getDof(idx1, n0);
WorldVector<double> c0 = elInfo->getCoord(idx0);
WorldVector<double> c1 = elInfo->getCoord(idx1);
c0 -= c1;
length[i] = norm(c0);
sumLength += length[i];
}
sumLength *= 0.5;
double area =
sqrt(sumLength * (sumLength - length[0]) * (sumLength - length[1]) * (sumLength - length[2]));
for (int i = 0; i < 3; i++) {
DegreeOfFreedom d = el->getDof(i, n0);
vec[d] += area;
}
elInfo = stack.traverseNext(elInfo);
}
meshDistributor->synchAddVector(vec);
double scalePower = 2.0;
Parameters::get("lp->scale power", scalePower);
DOFIterator<double> dofIt(&vec, USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
DegreeOfFreedom d = dofIt.getDOFIndex();
vec[d] = 1.0 / pow(vec[d], scalePower);
}
}
}
......@@ -41,7 +41,8 @@ namespace AMDiS {
class PetscSolverFeti : public PetscSolver
{
public:
PetscSolverFeti();
/// Constructor of FETI-DP solver class.
PetscSolverFeti(string name);
/// Assemble the sequentially created matrices to the global matrices
/// required by the FETI-DP method.
......@@ -228,11 +229,10 @@ namespace AMDiS {
return false;
}
/// Traverse the mesh and create a DOF vector where each DOF is in the
/// order of: 1/h^2
void createH2vec(DOFVector<double> &vec);
protected:
/// Prefix string for parameters in init file.
string initFileStr;
/// Mapping from primal DOF indices to a global index of primals.
ParallelDofMapping primalDofMap;
......
......@@ -35,7 +35,7 @@ namespace AMDiS {
VecDestroy(&v);
VecDestroy(&w);
PetscPrintf(PETSC_COMM_WORLD, "%3D KSP residual norm %1.12e [%1.12e %1.12e] and preconditioned norm [%1.12e]\n",
PetscPrintf(PETSC_COMM_WORLD, "%3D KSP residual norm %1.12e [ %1.12e %1.12e ] and preconditioned norm [%1.12e]\n",
n, norm, norm0, norm1, rnorm);
return 0;
......
......@@ -466,21 +466,7 @@ namespace AMDiS {
VecNestGetSubVec(yvec, 0, &y_interface);
VecNestGetSubVec(yvec, 1, &y_lagrange);
bool useMassMatrix = false;
Parameters::get("lp->mass matrix", useMassMatrix);
if (useMassMatrix) {
KSPSolve(data->ksp_mass, x_interface, y_interface);
} else {
VecCopy(x_interface, y_interface);
double scaleFactor = 1.0;
Parameters::get("lp->scal", scaleFactor);
bool mult = false;
Parameters::get("lp->mult", mult);
if (mult)
VecPointwiseMult(y_interface, x_interface, data->h2vec);
if (scaleFactor != 0.0)
VecScale(y_interface, scaleFactor);
}
KSPSolve(data->ksp_mass, x_interface, y_interface);
MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b0);
......
......@@ -144,8 +144,6 @@ namespace AMDiS {
PetscSolver* subSolver;
Vec h2vec;
Vec tmp_primal;
KSP ksp_mass;
......
......@@ -487,8 +487,8 @@ namespace AMDiS {
if (bIt->second && bIt->second->isDirichlet()) {
dirichletRow = true;
if ((dynamic_cast<DirichletBC*>(bIt->second))->applyBoundaryCondition()) {
TEST_EXIT_DBG(dirichletMainCol == -1)("Should not happen!\n");
dirichletMainCol = j;
break;
}
}
}
......
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