From 2bec218fb2cd21b726a4b9bcc1a616c46f5a682c Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 22 Jun 2012 07:10:29 +0000
Subject: [PATCH] Added more nullspace support.
---
AMDiS/src/Mesh.h | 14 +-
AMDiS/src/ProblemStat.h | 123 +++++++-----------
AMDiS/src/parallel/ElementObjectDatabase.cc | 1 +
AMDiS/src/parallel/MeshDistributor.h | 16 +--
AMDiS/src/parallel/PetscSolver.cc | 2 +
AMDiS/src/parallel/PetscSolver.h | 2 +
AMDiS/src/parallel/PetscSolverGlobalMatrix.cc | 56 ++++++--
AMDiS/src/parallel/PetscSolverGlobalMatrix.h | 10 +-
8 files changed, 120 insertions(+), 104 deletions(-)
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index e5622c4b..e29b0b07 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -402,11 +402,9 @@ namespace AMDiS {
/// Adds a macro element to the mesh
void addMacroElement(MacroElement* me);
- /* \brief
- * Removes a set of macro elements from the mesh. This works only for the
- * case, that there are no global or local refinements, i.e., all macro
- * elements have no children.
- */
+ /// Removes a set of macro elements from the mesh. This works only for the
+ /// case, that there are no global or local refinements, i.e., all macro
+ /// elements have no children.
void removeMacroElements(std::set<MacroElement*>& macros,
vector<const FiniteElemSpace*>& feSpaces);
@@ -484,10 +482,8 @@ namespace AMDiS {
}
- /** \brief
- * This function is equal to \ref getDofIndexCoords as defined above, but takes
- * a DOF index instead of a DOF pointer.
- */
+ /// This function is equal to \ref getDofIndexCoords as defined above, but
+ /// takes a DOF index instead of a DOF pointer.
bool getDofIndexCoords(DegreeOfFreedom dof,
const FiniteElemSpace* feSpace,
WorldVector<double>& coords);
diff --git a/AMDiS/src/ProblemStat.h b/AMDiS/src/ProblemStat.h
index 5e9439a1..faa4f8cb 100644
--- a/AMDiS/src/ProblemStat.h
+++ b/AMDiS/src/ProblemStat.h
@@ -49,11 +49,9 @@ namespace AMDiS {
};
- /** \brief
- * This class defines the stationary problem definition in sequential
- * computations. For parallel computations, see
- * \ref ParallelProblemStatBase.
- */
+ /// This class defines the stationary problem definition in sequential
+ /// computations. For parallel computations, see
+ /// \ref ParallelProblemStatBase.
class ProblemStatSeq : public ProblemStatBase,
public StandardProblemIteration
{
@@ -149,11 +147,9 @@ namespace AMDiS {
void dualAssemble(AdaptInfo *adaptInfo, Flag flag,
bool asmMatrix = true, bool asmVector = true);
- /** \brief
- * Determines the execution order of the single adaption steps. If adapt is
- * true, mesh adaption will be performed. This allows to avoid mesh adaption,
- * e.g. in timestep adaption loops of timestep adaptive strategies.
- */
+ /// Determines the execution order of the single adaption steps. If adapt is
+ /// true, mesh adaption will be performed. This allows to avoid mesh adaption,
+ /// e.g. in timestep adaption loops of timestep adaptive strategies.
virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION);
/// Returns number of managed problems
@@ -168,10 +164,8 @@ namespace AMDiS {
return nComponents;
}
- /** \brief
- * Returns the problem with the given number. If only one problem
- * is managed by this master problem, the number hasn't to be given.
- */
+ /// Returns the problem with the given number. If only one problem
+ /// is managed by this master problem, the number hasn't to be given.
virtual ProblemStatBase *getProblem(int number = 0)
{
return this;
@@ -260,10 +254,8 @@ namespace AMDiS {
addBoundaryVectorOperator(type, &op, row);
}
- /** \brief
- * This function assembles a DOFMatrix and a DOFVector for the case,
- * the meshes from row and col FE-space are equal.
- */
+ /// This function assembles a DOFMatrix and a DOFVector for the case,
+ /// the meshes from row and col FE-space are equal.
void assembleOnOneMesh(FiniteElemSpace *feSpace,
Flag assembleFlag,
DOFMatrix *matrix, DOFVector<double> *vector);
@@ -533,10 +525,8 @@ namespace AMDiS {
}
/** \} */
- /** \brief
- * Outputs the mesh of the given component, but the values are taken from the
- * residual error estimator.
- */
+ /// Outputs the mesh of the given component, but the values are taken from
+ /// the residual error estimator.
void writeResidualMesh(int comp, AdaptInfo *adaptInfo, string name);
/// Function that implements the serialization procedure.
@@ -553,10 +543,8 @@ namespace AMDiS {
}
protected:
- /** \brief
- * If the exact solution is known, the problem can compute the exact
- * error instead of the error estimation. This is done in this function.
- */
+ /// If the exact solution is known, the problem can compute the exact
+ /// error instead of the error estimation. This is done in this function.
void computeError(AdaptInfo *adaptInfo);
protected:
@@ -571,11 +559,9 @@ namespace AMDiS {
/// vectors, \ref solution.
vector<string> componentNames;
- /** \brief
- * Number of problem meshes. If all components are defined on the same mesh,
- * this number is 1. Otherwise, this variable is the number of different meshes
- * within the problem.
- */
+ /// Number of problem meshes. If all components are defined on the same mesh,
+ /// this number is 1. Otherwise, this variable is the number of different meshes
+ /// within the problem.
int nMeshes;
/// FE spaces of this problem.
@@ -590,11 +576,9 @@ namespace AMDiS {
/// Pointer to the meshes for the different problem components
vector<Mesh*> componentMeshes;
- /** \brief
- * Stores information about which meshes must be traversed to assemble the
- * specific components. I.e., it was implemented to make use of different
- * meshes for different components.
- */
+ /// Stores information about which meshes must be traversed to assemble the
+ /// specific components. I.e., it was implemented to make use of different
+ /// meshes for different components.
ComponentTraverseInfo traverseInfo;
/// Responsible for element marking.
@@ -618,21 +602,17 @@ namespace AMDiS {
/// Composed system matrix
SolverMatrix<Matrix<DOFMatrix*> > solverMatrix;
- /** \brief
- * Some DOFMatrices of the systemMatrix may be assembled only once (for
- * example if they are independent of the time or older solutions). If
- * [i][j] of this field is set to true, the corresponding DOFMatrix will
- * be assembled only once. All other matrices will be assembled at every
- * time step.
- */
+ /// Some DOFMatrices of the systemMatrix may be assembled only once (for
+ /// example if they are independent of the time or older solutions). If
+ /// [i][j] of this field is set to true, the corresponding DOFMatrix will
+ /// be assembled only once. All other matrices will be assembled at every
+ /// time step.
vector<vector<bool> > assembleMatrixOnlyOnce;
- /** \brief
- * If [i][j] of this field is set to true, the corresponding DOFMatrix of
- * the systemMatrix has been assembled at least once. This field is used
- * to determine, if assembling of a matrix can be ommitted, if it is set
- * to be assembled only once.
- */
+ /// If [i][j] of this field is set to true, the corresponding DOFMatrix of
+ /// the systemMatrix has been assembled at least once. This field is used
+ /// to determine, if assembling of a matrix can be ommitted, if it is set
+ /// to be assembled only once.
vector<vector<bool> > assembledMatrix;
/// Determines whether domain boundaries should be considered at assembling.
@@ -641,47 +621,40 @@ namespace AMDiS {
/// Writes the meshes and solution after the adaption loop.
vector<FileWriterInterface*> fileWriters;
- /** \brief
- * All actions of mesh refinement are performed by refinementManager.
- * If new refinement algorithms should be realized, one has to override
- * RefinementManager and give one instance of it to AdaptStationary.
- */
+ /// All actions of mesh refinement are performed by refinementManager.
+ /// If new refinement algorithms should be realized, one has to override
+ /// RefinementManager and give one instance of it to AdaptStationary.
RefinementManager *refinementManager;
- /** \brief
- * All actions of mesh coarsening are performed by coarseningManager.
- * If new coarsening algorithms should be realized, one has to override
- * CoarseningManager and give one instance of it to AdaptStationary.
- */
+ /// All actions of mesh coarsening are performed by coarseningManager.
+ /// If new coarsening algorithms should be realized, one has to override
+ /// CoarseningManager and give one instance of it to AdaptStationary.
CoarseningManager *coarseningManager;
/// Info level.
int info;
- /// If true, the stationary problem was deserialized from some serialization file.
+ /// If true, the stationary problem was deserialized from some serialization
+ /// file.
bool deserialized;
- /** \brief
- * This vectors stores pointers to functions defining the exact solution of
- * the problem. This may be used to compute the real error of the computed
- * solution.
- */
+ /// This vectors stores pointers to functions defining the exact solution of
+ /// the problem. This may be used to compute the real error of the computed
+ /// solution.
vector<AbstractFunction<double, WorldVector<double> >*> exactSolutionFcts;
- /** \brief
- * If true, the error is not estimated but computed from the exact solution
- * defined by \ref exactSolutionFcts.
- */
+ /// If true, the error is not estimated but computed from the exact solution
+ /// defined by \ref exactSolutionFcts.
bool computeExactError;
- /** \brief
- * If at least on boundary condition is set, this variable is true. It is used
- * to ensure that no operators are added after boundary condition were set. If
- * this would happen, boundary conditions could set wrong on off diagonal matrices.
- */
+ /// If at least on boundary condition is set, this variable is true. It is
+ /// used to ensure that no operators are added after boundary condition were
+ /// set. If this would happen, boundary conditions could set wrong on off
+ /// diagonal matrices.
bool boundaryConditionSet;
- /// If true, AMDiS prints information about the assembling procedure to the screen.
+ /// If true, AMDiS prints information about the assembling procedure to
+ /// the screen.
bool writeAsmInfo;
map<Operator*, vector<OperatorPos> > operators;
diff --git a/AMDiS/src/parallel/ElementObjectDatabase.cc b/AMDiS/src/parallel/ElementObjectDatabase.cc
index fea60559..12ffab84 100644
--- a/AMDiS/src/parallel/ElementObjectDatabase.cc
+++ b/AMDiS/src/parallel/ElementObjectDatabase.cc
@@ -243,6 +243,7 @@ namespace AMDiS {
("Should not happen (%d)!\n", multPeriodicDof2.size());
TEST_EXIT_DBG(multPeriodicDof3.size() == 0)("Should not happen!\n");
}
+
if (mesh->getDim() == 3) {
TEST_EXIT_DBG(multPeriodicDof3.size() == 0 ||
multPeriodicDof3.size() == 8)
diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h
index 75c412dd..953cbec6 100644
--- a/AMDiS/src/parallel/MeshDistributor.h
+++ b/AMDiS/src/parallel/MeshDistributor.h
@@ -295,6 +295,14 @@ namespace AMDiS {
return levelData;
}
+ void updateLocalGlobalNumbering();
+
+ /// Updates the local and global DOF numbering after the mesh has been
+ /// changed.
+ void updateLocalGlobalNumbering(ParallelDofMapping &dmap,
+ DofComm &dcom,
+ const FiniteElemSpace *feSpace);
+
protected:
void addProblemStat(ProblemStatSeq *probStat);
@@ -308,14 +316,6 @@ namespace AMDiS {
/// partition.
void removeMacroElements();
- void updateLocalGlobalNumbering();
-
- /// Updates the local and global DOF numbering after the mesh has been
- /// changed.
- void updateLocalGlobalNumbering(ParallelDofMapping &dmap,
- DofComm &dcom,
- const FiniteElemSpace *feSpace);
-
/// Calls \ref createPeriodicMap(feSpace) for all FE spaces that are
/// handled by the mesh distributor.
void createPeriodicMap();
diff --git a/AMDiS/src/parallel/PetscSolver.cc b/AMDiS/src/parallel/PetscSolver.cc
index db9e5f58..9ce57756 100644
--- a/AMDiS/src/parallel/PetscSolver.cc
+++ b/AMDiS/src/parallel/PetscSolver.cc
@@ -37,6 +37,8 @@ namespace AMDiS {
Parameters::get("parallel->remove rhs null space", removeRhsNullspace);
Parameters::get("parallel->has constant null space", hasConstantNullspace);
+ Parameters::get("parallel->nullspace->const in comp",
+ constNullspaceComponent);
}
diff --git a/AMDiS/src/parallel/PetscSolver.h b/AMDiS/src/parallel/PetscSolver.h
index 60279806..dabfeb7c 100644
--- a/AMDiS/src/parallel/PetscSolver.h
+++ b/AMDiS/src/parallel/PetscSolver.h
@@ -256,6 +256,8 @@ namespace AMDiS {
bool removeRhsNullspace;
bool hasConstantNullspace;
+
+ vector<int> constNullspaceComponent;
};
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 27e48d03..578fd3b2 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -119,9 +119,8 @@ namespace AMDiS {
KSPSetType(kspInterior, KSPBCGS);
KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str());
KSPSetFromOptions(kspInterior);
- PCSetFromOptions(pcInterior);
- createFieldSplit(pcInterior);
+ initPreconditioner(pcInterior);
// Do not delete the solution vector, use it for the initial guess.
if (!zeroStartVector)
@@ -131,6 +130,14 @@ namespace AMDiS {
}
+ void PetscSolverGlobalMatrix::fillPetscMatrix(DOFMatrix *mat)
+ {
+ Matrix<DOFMatrix*> m(1, 1);
+ m[0][0] = mat;
+ fillPetscMatrix(&m);
+ }
+
+
void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *mat)
{
FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
@@ -317,7 +324,6 @@ namespace AMDiS {
KSPSetOperators(kspInterior, matIntInt, matIntInt, SAME_NONZERO_PATTERN);
KSPSetOptionsPrefix(kspInterior, "interior_");
KSPSetType(kspInterior, KSPPREONLY);
- PC pcInterior;
KSPGetPC(kspInterior, &pcInterior);
PCSetType(pcInterior, PCLU);
if (subdomainLevel == 0)
@@ -393,9 +399,21 @@ namespace AMDiS {
MatNullSpace matNullspace;
Vec nullspaceBasis;
- if (nullspace.size() > 0 || hasConstantNullspace) {
+ if (nullspace.size() > 0 ||
+ hasConstantNullspace ||
+ constNullspaceComponent.size() > 0) {
TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
+ if (constNullspaceComponent.size() > 0) {
+ nullspace.clear();
+ SystemVector *basisVec = new SystemVector(vec);
+ basisVec->set(0.0);
+ for (unsigned int i = 0; i < constNullspaceComponent.size(); i++)
+ basisVec->getDOFVector(constNullspaceComponent[i])->set(1.0);
+
+ nullspace.push_back(basisVec);
+ }
+
if (nullspace.size() > 0) {
VecDuplicate(petscSolVec, &nullspaceBasis);
for (int i = 0; i < nComponents; i++)
@@ -403,6 +421,8 @@ namespace AMDiS {
VecAssemblyBegin(nullspaceBasis);
VecAssemblyEnd(nullspaceBasis);
+
+ VecNormalize(nullspaceBasis, PETSC_NULL);
MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE),
1, &nullspaceBasis, &matNullspace);
@@ -415,6 +435,7 @@ namespace AMDiS {
MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
}
+ MatSetNullSpace(matIntInt, matNullspace);
KSPSetNullSpace(kspInterior, matNullspace);
// === Remove null space, if requested. ===
@@ -455,11 +476,6 @@ namespace AMDiS {
}
VecRestoreArray(petscSolVec, &vecPointer);
-
-
- // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to ===
- // === more than one partition. ===
- meshDistributor->synchVector(vec);
}
@@ -517,6 +533,8 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");
+ exitPreconditioner(pcInterior);
+
MatDestroy(&matIntInt);
KSPDestroy(&kspInterior);
@@ -580,6 +598,21 @@ namespace AMDiS {
}
+ void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
+ {
+ FUNCNAME("PetscSolverGlobalMatrix::initPreconditioner()");
+
+ PCSetFromOptions(pc);
+ createFieldSplit(pc);
+ }
+
+
+ void PetscSolverGlobalMatrix::exitPreconditioner(PC pc)
+ {
+ FUNCNAME("PetscSolverGlobalMatrix::exitPreconditioner()");
+ }
+
+
void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* mat,
int nRowMat, int nColMat)
{
@@ -766,7 +799,6 @@ namespace AMDiS {
else
perColDof = entry[i].second;
-
entry.push_back(make_pair(perRowDof, perColDof));
}
}
@@ -815,6 +847,7 @@ namespace AMDiS {
// Traverse all used DOFs in the dof vector.
DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
+
if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
continue;
@@ -836,7 +869,7 @@ namespace AMDiS {
else
index =
interiorMap->getMatIndex(nRowVec, dofIt.getDOFIndex()) + rStartInterior;
-
+
if (perMap.isPeriodic(feSpace, globalRowDof)) {
std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
double value = *dofIt / (perAsc.size() + 1.0);
@@ -846,6 +879,7 @@ namespace AMDiS {
perIt != perAsc.end(); ++perIt) {
int mappedDof = perMap.map(feSpace, *perIt, globalRowDof);
int mappedIndex = interiorMap->getMatIndex(nRowVec, mappedDof);
+
VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
}
} else {
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
index 0c9457d4..e3efdf02 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
@@ -50,6 +50,10 @@ namespace AMDiS {
void fillPetscMatrix(Matrix<DOFMatrix*> *mat);
+ /// This function is just a small wrapper that creates a 1x1 matrix that
+ /// contains exactly one DOFMatrix and than calls \ref fillPetscMatrix
+ void fillPetscMatrix(DOFMatrix *mat);
+
void fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *mat);
void fillPetscRhs(SystemVector *vec);
@@ -62,9 +66,13 @@ namespace AMDiS {
void destroyVectorData();
+ protected:
void createFieldSplit(PC pc);
- protected:
+ virtual void initPreconditioner(PC pc);
+
+ virtual void exitPreconditioner(PC pc);
+
/// Creates a new non zero pattern structure for the PETSc matrix.
void createPetscNnzStructure(Matrix<DOFMatrix*> *mat);
--
GitLab