From f042083cad0a2ea7f11bc9c6ca046e415687fe7a Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Wed, 15 Aug 2012 09:05:31 +0000
Subject: [PATCH] Better matrix nnz computation when using periodic boundary
conditions.
---
AMDiS/src/parallel/MatrixNnzStructure.cc | 54 +++++++++++------
AMDiS/src/parallel/MatrixNnzStructure.h | 8 ++-
AMDiS/src/parallel/PeriodicMap.cc | 39 ++++++++++++
AMDiS/src/parallel/PeriodicMap.h | 34 +++++++++++
AMDiS/src/parallel/PetscSolver.h | 1 -
AMDiS/src/parallel/PetscSolverGlobalMatrix.cc | 60 +++++--------------
6 files changed, 129 insertions(+), 67 deletions(-)
diff --git a/AMDiS/src/parallel/MatrixNnzStructure.cc b/AMDiS/src/parallel/MatrixNnzStructure.cc
index c03c8c56..e8484c33 100644
--- a/AMDiS/src/parallel/MatrixNnzStructure.cc
+++ b/AMDiS/src/parallel/MatrixNnzStructure.cc
@@ -37,6 +37,7 @@ namespace AMDiS {
ParallelDofMapping &rowDofMap,
ParallelDofMapping &colDofMap,
PeriodicMap *perMap,
+ const ElementObjectDatabase& elObjDb,
bool localMatrix)
{
FUNCNAME("MatrixNnzStructure::create()");
@@ -147,8 +148,12 @@ namespace AMDiS {
petscRowIdx = rowDofMap.getMatIndex(rowComp, *cursor);
}
- bool periodicRow =
- (perMap && perMap->isPeriodic(rowFeSpace, rowIt->second.global));
+
+ // Set of periodic row associations (is empty, if row DOF is not
+ // periodic.
+ std::set<int> perRowAsc;
+ if (perMap)
+ perMap->fillAssociations(rowFeSpace, rowIt->second.global, elObjDb, perRowAsc);
if (localMatrix || rowDofMap[rowFeSpace].isRankDof(*cursor)) {
@@ -186,24 +191,35 @@ namespace AMDiS {
icend = end<nz>(cursor); icursor != icend; ++icursor) {
if (colDofMap[colFeSpace].find(col(*icursor), colDofIndex) == false)
continue;
-
- int petscColIdx = (colDofMap.isMatIndexFromGlobal() ?
- colDofMap.getMatIndex(colComp, colDofIndex.global) :
- colDofMap.getMatIndex(colComp, col(*icursor)));
-
- // The row DOF is a rank DOF, if also the column is a rank DOF,
- // increment the diagNnz values for this row,
- // otherwise the offdiagNnz value.
- if (petscColIdx >= rankStartColIndex &&
- petscColIdx < rankStartColIndex + nRankCols)
- dnnz[localPetscRowIdx]++;
- else
- onnz[localPetscRowIdx]++;
- }
- if (periodicRow) {
- dnnz[localPetscRowIdx] *= 4;
- onnz[localPetscRowIdx] *= 4;
+ // Set of periodic row associations (is empty, if row DOF is not
+ // periodic.
+ std::set<int> perColAsc;
+ if (perMap)
+ perMap->fillAssociations(colFeSpace, colDofIndex.global, elObjDb, perColAsc);
+
+ if (perColAsc.empty()) {
+ if (colDofMap[colFeSpace].isRankDof(col(*icursor)))
+ dnnz[localPetscRowIdx]++;
+ else
+ onnz[localPetscRowIdx]++;
+ } else {
+ vector<int> newCols;
+ perMap->mapDof(colFeSpace, colDofIndex.global, perColAsc, newCols);
+
+ for (int aa = 0; aa < newCols.size(); aa++) {
+ int petscColIdx = colDofMap.getMatIndex(colComp, newCols[aa]);
+
+ // The row DOF is a rank DOF, if also the column is a rank DOF,
+ // increment the diagNnz values for this row,
+ // otherwise the offdiagNnz value.
+ if (petscColIdx >= rankStartColIndex &&
+ petscColIdx < rankStartColIndex + nRankCols)
+ dnnz[localPetscRowIdx]++;
+ else
+ onnz[localPetscRowIdx]++;
+ }
+ }
}
}
} else {
diff --git a/AMDiS/src/parallel/MatrixNnzStructure.h b/AMDiS/src/parallel/MatrixNnzStructure.h
index de266459..defb21de 100644
--- a/AMDiS/src/parallel/MatrixNnzStructure.h
+++ b/AMDiS/src/parallel/MatrixNnzStructure.h
@@ -44,16 +44,18 @@ namespace AMDiS {
MPI::Intracomm mpiComm,
ParallelDofMapping &rowDofMap,
ParallelDofMapping &colDofMap,
- PeriodicMap *perMap = NULL,
+ PeriodicMap *perMap,
+ const ElementObjectDatabase& elObjDb,
bool localMatrix = false);
void create(Matrix<DOFMatrix*> *mat,
MPI::Intracomm mpiComm,
ParallelDofMapping &dofMap,
- PeriodicMap *perMap = NULL,
+ PeriodicMap *perMap,
+ const ElementObjectDatabase& elObjDb,
bool localMatrix = false)
{
- create(mat, mpiComm, dofMap, dofMap, perMap, localMatrix);
+ create(mat, mpiComm, dofMap, dofMap, perMap, elObjDb, localMatrix);
}
protected:
diff --git a/AMDiS/src/parallel/PeriodicMap.cc b/AMDiS/src/parallel/PeriodicMap.cc
index 558a46ce..1fe38a99 100644
--- a/AMDiS/src/parallel/PeriodicMap.cc
+++ b/AMDiS/src/parallel/PeriodicMap.cc
@@ -13,6 +13,7 @@
#include <map>
#include <fstream>
#include "Global.h"
+#include "parallel/ElementObjectDatabase.h"
#include "parallel/PeriodicMap.h"
#include "parallel/ParallelTypes.h"
@@ -33,6 +34,44 @@ namespace AMDiS {
}
+ void PeriodicMap::fillAssociations(const FiniteElemSpace* feSpace,
+ int globalDofIndex,
+ const ElementObjectDatabase &elObjDb,
+ std::set<int>& perAsc)
+ {
+ if (isPeriodic(feSpace, globalDofIndex) == false)
+ return;
+
+ std::set<int>& tmpAsc = getAssociations(feSpace, globalDofIndex);
+ for (std::set<int>::iterator it = tmpAsc.begin(); it != tmpAsc.end(); ++it)
+ if (elObjDb.isValidPeriodicType(*it))
+ perAsc.insert(*it);
+ }
+
+
+ void PeriodicMap::mapDof(const FiniteElemSpace* feSpace,
+ int globalDofIndex,
+ const std::set<int>& perAsc,
+ vector<int>& mappedDofs)
+ {
+ mappedDofs.clear();
+ mappedDofs.push_back(globalDofIndex);
+
+ for (std::set<int>::const_iterator it = perAsc.begin();
+ it != perAsc.end(); ++it) {
+ int nDofs = static_cast<int>(mappedDofs.size());
+
+ for (int i = 0; i < nDofs; i++) {
+ TEST_EXIT_DBG(isPeriodic(feSpace, *it, mappedDofs[i]))
+ ("Wrong periodic DOF associations at boundary %d with DOF %d!\n",
+ *it, mappedDofs[i]);
+
+ mappedDofs.push_back(map(feSpace, *it, mappedDofs[i]));
+ }
+ }
+ }
+
+
void PeriodicMap::serialize(ostream &out,
vector<const FiniteElemSpace*> feSpaces)
{
diff --git a/AMDiS/src/parallel/PeriodicMap.h b/AMDiS/src/parallel/PeriodicMap.h
index 8ed5e2e7..0670eeb9 100644
--- a/AMDiS/src/parallel/PeriodicMap.h
+++ b/AMDiS/src/parallel/PeriodicMap.h
@@ -26,6 +26,7 @@
#include <map>
#include <set>
#include <vector>
+#include "AMDiS_fwd.h"
#include "parallel/ParallelTypes.h"
#include "Boundary.h"
#include "Serializer.h"
@@ -126,6 +127,39 @@ namespace AMDiS {
}
+ /** \brief
+ * Given a global DOF index in some given finite element space, this
+ * function creates all valid periodic associations of this DOF. If the
+ * DOF is not periodic, this function does not fail but has no effect on
+ * the output data.
+ *
+ * \param[in] feSpace feSpace on which the function should work on.
+ * \param[in] globalDofIndex global index of a DOF.
+ * \param[in] elObjDb Element object database that is used to check
+ * if a given periodic index is valid.
+ * \param[out] perAsc set of periodic associations, is not deleted
+ * by this function
+ */
+ void fillAssociations(const FiniteElemSpace* feSpace,
+ int globalDofIndex,
+ const ElementObjectDatabase &elObjDb,
+ std::set<int>& perAsc);
+
+
+ /** \brief
+ * Maps a given DOF index for all given periodic DOF associations.
+ *
+ * \param[in] feSpace feSpace on which the function should work on.
+ * \param[in] globalDofIndex global index of a DOF.
+ * \param[in] perAsc set of periodic associations.
+ * \param[out] mappedDofs set of global DOF indices.
+ */
+ void mapDof(const FiniteElemSpace* feSpace,
+ int globalDofIndex,
+ const std::set<int>& perAsc,
+ vector<int>& mappedDofs);
+
+
/// Returns true, if the DOF (global index) is a periodic DOF.
inline bool isPeriodic(const FiniteElemSpace *feSpace, int globalDofIndex)
{
diff --git a/AMDiS/src/parallel/PetscSolver.h b/AMDiS/src/parallel/PetscSolver.h
index 78fd2e55..57f78d28 100644
--- a/AMDiS/src/parallel/PetscSolver.h
+++ b/AMDiS/src/parallel/PetscSolver.h
@@ -113,7 +113,6 @@ namespace AMDiS {
void setKspPrefix(std::string s)
{
- MSG("======>>>>>> SET TO %s\n", s.c_str());
kspPrefix = s;
}
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index c1deae93..70b1c65e 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -37,7 +37,8 @@ namespace AMDiS {
if (checkMeshChange(mat))
nnzInterior.create(mat, mpiCommGlobal, *interiorMap,
- &(meshDistributor->getPeriodicMap()));
+ &(meshDistributor->getPeriodicMap()),
+ meshDistributor->getElementObjectDb());
// === Create PETSc vector (solution and a temporary vector). ===
@@ -109,7 +110,6 @@ namespace AMDiS {
KSPSetOperators(kspInterior, matIntInt, matIntInt, SAME_NONZERO_PATTERN);
KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
KSPSetType(kspInterior, KSPBCGS);
- MSG("SET PREFIX TOOO: %s\n", kspPrefix.c_str());
KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str());
KSPSetFromOptions(kspInterior);
@@ -148,12 +148,14 @@ namespace AMDiS {
bool localMatrix = (subdomainLevel == 0);
if (checkMeshChange(mat, localMatrix)) {
- nnzInterior.create(mat, mpiCommGlobal, *interiorMap, NULL, localMatrix);
+ nnzInterior.create(mat, mpiCommGlobal, *interiorMap, NULL,
+ meshDistributor->getElementObjectDb(),
+ localMatrix);
if (coarseSpaceMap) {
- nnzCoarse.create(mat, mpiCommGlobal, *coarseSpaceMap);
- nnzCoarseInt.create(mat, mpiCommGlobal, *coarseSpaceMap, *interiorMap);
- nnzIntCoarse.create(mat, mpiCommGlobal, *interiorMap, *coarseSpaceMap);
+ nnzCoarse.create(mat, mpiCommGlobal, *coarseSpaceMap, NULL, meshDistributor->getElementObjectDb());
+ nnzCoarseInt.create(mat, mpiCommGlobal, *coarseSpaceMap, *interiorMap, NULL, meshDistributor->getElementObjectDb());
+ nnzIntCoarse.create(mat, mpiCommGlobal, *interiorMap, *coarseSpaceMap, NULL, meshDistributor->getElementObjectDb());
}
}
@@ -743,12 +745,9 @@ namespace AMDiS {
// Create set of all periodic associations of the column index.
std::set<int> perAsc;
- std::set<int>& perColAsc = perMap.getAssociations(colFe, globalColDof);
- for (std::set<int>::iterator it = perColAsc.begin();
- it != perColAsc.end(); ++it)
- if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it))
- perAsc.insert(*it);
-
+ perMap.fillAssociations(colFe, globalColDof,
+ meshDistributor->getElementObjectDb(), perAsc);
+
// Scale value to the number of periodic associations of the column index.
double scaledValue =
value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));
@@ -758,24 +757,7 @@ namespace AMDiS {
// === associations of the column DOF index. ===
vector<int> newCols;
-
- // First, add the original matrix index.
- newCols.push_back(globalColDof);
-
- // And add all periodic matrix indices.
- for (std::set<int>::iterator it = perAsc.begin();
- it != perAsc.end(); ++it) {
- int nCols = static_cast<int>(newCols.size());
-
- for (int i = 0; i < nCols; i++) {
- TEST_EXIT_DBG(perMap.isPeriodic(colFe, *it, newCols[i]))
- ("Wrong periodic DOF associations at boundary %d with DOF %d!\n",
- *it, newCols[i]);
-
- newCols.push_back(perMap.map(colFe, *it, newCols[i]));
- }
- }
-
+ perMap.mapDof(colFe, globalColDof, perAsc, newCols);
for (unsigned int i = 0; i < newCols.size(); i++) {
cols.push_back(interiorMap->getMatIndex(nColMat, newCols[i]));
values.push_back(scaledValue);
@@ -809,20 +791,10 @@ namespace AMDiS {
// === indices to the set perAsc. ===
std::set<int> perAsc;
-
- if (perMap.isPeriodic(colFe, globalColDof)) {
- std::set<int>& perColAsc = perMap.getAssociations(colFe, globalColDof);
- for (std::set<int>::iterator it = perColAsc.begin();
- it != perColAsc.end(); ++it)
- if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it))
- perAsc.insert(*it);
- }
-
- std::set<int>& perRowAsc = perMap.getAssociations(rowFe, globalRowDof);
- for (std::set<int>::iterator it = perRowAsc.begin();
- it != perRowAsc.end(); ++it)
- if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it))
- perAsc.insert(*it);
+ perMap.fillAssociations(colFe, globalColDof,
+ meshDistributor->getElementObjectDb(), perAsc);
+ perMap.fillAssociations(rowFe, globalRowDof,
+ meshDistributor->getElementObjectDb(), perAsc);
// Scale the value with respect to the number of periodic associations.
double scaledValue =
--
GitLab