From 9ce826f5164ad5a1cf4e424f34b527416a69460d Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 13 Mar 2012 15:47:11 +0000
Subject: [PATCH] Mh, seems to work.
---
AMDiS/src/parallel/FeSpaceMapping.cc | 7 +-
AMDiS/src/parallel/FeSpaceMapping.h | 15 ++-
AMDiS/src/parallel/PetscSolverFeti.cc | 147 ++++++++++++--------------
AMDiS/src/parallel/PetscSolverFeti.h | 3 -
4 files changed, 84 insertions(+), 88 deletions(-)
diff --git a/AMDiS/src/parallel/FeSpaceMapping.cc b/AMDiS/src/parallel/FeSpaceMapping.cc
index 0d582345..36d1d7f4 100644
--- a/AMDiS/src/parallel/FeSpaceMapping.cc
+++ b/AMDiS/src/parallel/FeSpaceMapping.cc
@@ -47,7 +47,8 @@ namespace AMDiS {
void GlobalDofMap::addOffset(int offset)
{
for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
- it->second.local += offset;
+ it->second.global = it->second.local + offset;
+ // it->second.local += offset;
}
void GlobalDofMap::computeNonLocalIndices()
@@ -60,7 +61,7 @@ namespace AMDiS {
for (DofComm::Iterator it(*sendDofs, feSpace); !it.end(); it.nextRank())
for (; !it.endDofIter(); it.nextDof())
if (dofMap.count(it.getDofIndex()) && !nonRankDofs.count(it.getDofIndex()))
- stdMpi.getSendData(it.getRank()).push_back(dofMap[it.getDofIndex()].local);
+ stdMpi.getSendData(it.getRank()).push_back(dofMap[it.getDofIndex()].global);
stdMpi.updateSendDataSize();
@@ -84,7 +85,7 @@ namespace AMDiS {
int i = 0;
for (; !it.endDofIter(); it.nextDof())
if (nonRankDofs.count(it.getDofIndex()))
- dofMap[it.getDofIndex()].local = stdMpi.getRecvData(it.getRank())[i++];
+ dofMap[it.getDofIndex()].global = stdMpi.getRecvData(it.getRank())[i++];
}
}
diff --git a/AMDiS/src/parallel/FeSpaceMapping.h b/AMDiS/src/parallel/FeSpaceMapping.h
index 81ee57c0..b2832527 100644
--- a/AMDiS/src/parallel/FeSpaceMapping.h
+++ b/AMDiS/src/parallel/FeSpaceMapping.h
@@ -369,8 +369,10 @@ namespace AMDiS {
typedef map<DegreeOfFreedom, MultiIndex>::iterator ItType;
for (ItType it = dofMap.begin(); it != dofMap.end(); ++it) {
if (data[feSpaces[i]].isRankDof(it->first)) {
- int globalMatIndex =
- it->second.local - data[feSpaces[i]].rStartDofs + offset;
+ int globalMatIndex = it->second.local + offset;
+
+// int globalMatIndex =
+// it->second.local - data[feSpaces[i]].rStartDofs + offset;
dofToMatIndex.add(i, it->first, globalMatIndex);
}
}
@@ -400,9 +402,9 @@ namespace AMDiS {
stdMpi.startCommunication();
{
- int counter = 0;
for (DofComm::Iterator it(*recvDofs, feSpaces[i]);
!it.end(); it.nextRank()) {
+ int counter = 0;
for (; !it.endDofIter(); it.nextDof()) {
if (dofMap.count(it.getDofIndex())) {
DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
@@ -465,6 +467,13 @@ namespace AMDiS {
data[*it].setDofComm(pSend, pRecv);
}
+ inline int getMatIndex(int ithComponent, DegreeOfFreedom d)
+ {
+ FUNCNAME("FeSpaceData::getMatIndex()");
+
+ return dofToMatIndex.get(ithComponent, d);
+ }
+
private:
MPI::Intracomm* mpiComm;
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 4fbe4262..6201c4f9 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -181,7 +181,6 @@ namespace AMDiS {
PetscSolverFeti::PetscSolverFeti()
: PetscSolver(),
- nComponents(-1),
schurPrimalSolver(0)
{
FUNCNAME("PetscSolverFeti::PetscSolverFeti()");
@@ -408,8 +407,6 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::createMatLagrange()");
- const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
-
// === Create distributed matrix for Lagrange constraints. ===
MatCreateMPIAIJ(PETSC_COMM_WORLD,
@@ -425,35 +422,34 @@ namespace AMDiS {
// === m == r, than the rank sets -1.0 for the corresponding ===
// === constraint. ===
- map<DegreeOfFreedom, MultiIndex> &dualMap = dualDofMap[feSpace].getMap();
- for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
- TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
- ("Should not happen!\n");
+ for (unsigned int k = 0; k < feSpaces.size(); k++) {
+ map<DegreeOfFreedom, MultiIndex> &dualMap =
+ dualDofMap[feSpaces[k]].getMap();
- // Global index of the first Lagrange constriant for this node.
- int index = lagrangeMap[feSpace][it->first].local;
- // Copy set of all ranks that contain this dual node.
- vector<int> W(boundaryDofRanks[it->first].begin(),
- boundaryDofRanks[it->first].end());
- // Number of ranks that contain this dual node.
- int degree = W.size();
-
- for (int i = 0; i < degree; i++) {
- for (int j = i + 1; j < degree; j++) {
- if (W[i] == mpiRank || W[j] == mpiRank) {
- // Set the constraint for all components of the system.
- for (int k = 0; k < nComponents; k++) {
- int rowIndex = index * nComponents + k;
- int colIndex =
- (localDofMap[feSpace].rStartDofs +
- localDofMap[feSpace][it->first].local) * nComponents + k;
+ for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin();
+ it != dualMap.end(); ++it) {
+ TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
+ ("Should not happen!\n");
+
+ // Global index of the first Lagrange constriant for this node.
+ int index = lagrangeMap.getMatIndex(k, it->first);
+
+ // Copy set of all ranks that contain this dual node.
+ vector<int> W(boundaryDofRanks[it->first].begin(),
+ boundaryDofRanks[it->first].end());
+ // Number of ranks that contain this dual node.
+ int degree = W.size();
+
+ for (int i = 0; i < degree; i++) {
+ for (int j = i + 1; j < degree; j++) {
+ if (W[i] == mpiRank || W[j] == mpiRank) {
+ int colIndex = localDofMap.mapGlobal(it->first, k);
double value = (W[i] == mpiRank ? 1.0 : -1.0);
- MatSetValue(mat_lagrange, rowIndex, colIndex, value,
+ MatSetValue(mat_lagrange, index, colIndex, value,
INSERT_VALUES);
}
+ index++;
}
-
- index++;
}
}
}
@@ -767,8 +763,6 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::recoverSolution()");
- const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
-
// === Get local part of the solution for B variables. ===
PetscScalar *localSolB;
@@ -778,19 +772,25 @@ namespace AMDiS {
// === Create scatter to get solutions of all primal nodes that are ===
// === contained in rank's domain. ===
+ vector<const FiniteElemSpace*> feSpaces = getFeSpaces(&vec);
+
vector<PetscInt> globalIsIndex, localIsIndex;
globalIsIndex.reserve(primalDofMap.getLocalDofs());
localIsIndex.reserve(primalDofMap.getLocalDofs());
{
- int counter = 0;
- for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
- it != primalDofMap[feSpace].getMap().end(); ++it) {
- for (int i = 0; i < nComponents; i++) {
- globalIsIndex.push_back(it->second.local * nComponents + i);
- localIsIndex.push_back(counter++);
+ int cnt = 0;
+ for (unsigned int i = 0; i < feSpaces.size(); i++) {
+ map<DegreeOfFreedom, MultiIndex>& dofMap =
+ primalDofMap[feSpaces[i]].getMap();
+ for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin();
+ it != dofMap.end(); ++it) {
+ globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first));
+ localIsIndex.push_back(cnt++);
}
}
+
+ TEST_EXIT_DBG(cnt == primalDofMap.getLocalDofs())("Should not happen!\n");
}
IS globalIs, localIs;
@@ -825,21 +825,19 @@ namespace AMDiS {
// === And copy from PETSc local vectors to the DOF vectors. ===
- for (int i = 0; i < nComponents; i++) {
+ int cnt = 0;
+ for (int i = 0; i < vec.getSize(); i++) {
DOFVector<double>& dofVec = *(vec.getDOFVector(i));
- for (map<DegreeOfFreedom, MultiIndex>::iterator it = localDofMap[feSpace].getMap().begin();
- it != localDofMap[feSpace].getMap().end(); ++it) {
+ for (map<DegreeOfFreedom, MultiIndex>::iterator it = localDofMap[feSpaces[i]].getMap().begin();
+ it != localDofMap[feSpaces[i]].getMap().end(); ++it) {
int petscIndex = localDofMap.mapLocal(it->first, i);
dofVec[it->first] = localSolB[petscIndex];
}
- int counter = 0;
- for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
- it != primalDofMap[feSpace].getMap().end(); ++it) {
- dofVec[it->first] = localSolPrimal[counter * nComponents + i];
- counter++;
- }
+ for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpaces[i]].getMap().begin();
+ it != primalDofMap[feSpaces[i]].getMap().end(); ++it)
+ dofVec[it->first] = localSolPrimal[cnt++];
}
VecRestoreArray(vec_sol_b, &localSolB);
@@ -852,13 +850,10 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::fillPetscMatrix()");
- vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
-
- const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
- nComponents = mat->getSize();
-
// === Create all sets and indices. ===
+ vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
+
primalDofMap.init(mpiComm, feSpaces, true, true);
dualDofMap.init(mpiComm, feSpaces, false, false);
lagrangeMap.init(mpiComm, feSpaces, true, true);
@@ -943,6 +938,7 @@ namespace AMDiS {
// === Traverse all sequentially created matrices and add the values to ===
// === the global PETSc matrices. ===
+ int nComponents = mat->getSize();
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
if (!(*mat)[i][j])
@@ -955,7 +951,7 @@ namespace AMDiS {
for (cursor_type cursor = begin<row>((*mat)[i][j]->getBaseMatrix()),
cend = end<row>((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
- bool rowPrimal = isPrimal(feSpace, *cursor);
+ bool rowPrimal = isPrimal(feSpaces[i], *cursor);
cols.clear();
colsOther.clear();
@@ -972,7 +968,7 @@ namespace AMDiS {
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
icursor != icend; ++icursor) {
- bool colPrimal = primalDofMap[feSpace].isSet(col(*icursor));
+ bool colPrimal = isPrimal(feSpaces[j], col(*icursor));
if (colPrimal) {
if (rowPrimal) {
@@ -996,29 +992,29 @@ namespace AMDiS {
// === For preconditioner ===
if (!rowPrimal && !colPrimal) {
- int rowIndex = localDofMap[feSpace][*cursor].local;
- int colIndex = localDofMap[feSpace][col(*icursor)].local;
+ int rowIndex = localDofMap[feSpaces[i]][*cursor].local;
+ int colIndex = localDofMap[feSpaces[j]][col(*icursor)].local;
- if (rowIndex < nLocalInterior[feSpace]) {
- if (colIndex < nLocalInterior[feSpace]) {
+ if (rowIndex < nLocalInterior[feSpaces[i]]) {
+ if (colIndex < nLocalInterior[feSpaces[j]]) {
int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
colsLocal.push_back(colIndex2);
valuesLocal.push_back(value(*icursor));
} else {
int colIndex2 =
- (localDofMap[feSpace][col(*icursor)].local - nLocalInterior[feSpace]) *
+ (localDofMap[feSpaces[j]][col(*icursor)].local - nLocalInterior[feSpaces[j]]) *
nComponents + j;
colsLocalOther.push_back(colIndex2);
valuesLocalOther.push_back(value(*icursor));
}
} else {
- if (colIndex < nLocalInterior[feSpace]) {
+ if (colIndex < nLocalInterior[feSpaces[j]]) {
int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
colsLocalOther.push_back(colIndex2);
valuesLocalOther.push_back(value(*icursor));
} else {
- int colIndex2 = (localDofMap[feSpace][col(*icursor)].local - nLocalInterior[feSpace]) *
+ int colIndex2 = (localDofMap[feSpaces[j]][col(*icursor)].local - nLocalInterior[feSpaces[j]]) *
nComponents + j;
colsLocal.push_back(colIndex2);
@@ -1034,11 +1030,9 @@ namespace AMDiS {
// === Set matrix values. ===
if (rowPrimal) {
- int rowIndex =
- primalDofMap[feSpace][*cursor].local * nComponents + i;
-
+ int rowIndex = primalDofMap.getMatIndex(i, *cursor);
for (unsigned int k = 0; k < cols.size(); k++)
- cols[k] = primalDofMap[feSpace][cols[k]].local * nComponents + j;
+ cols[k] = primalDofMap.getMatIndex(j, cols[k]);
MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(),
&(cols[0]), &(values[0]), ADD_VALUES);
@@ -1061,10 +1055,8 @@ namespace AMDiS {
if (colsOther.size()) {
int globalRowIndex = localDofMap.mapGlobal(*cursor, i);
-
for (unsigned int k = 0; k < colsOther.size(); k++)
- colsOther[k] =
- primalDofMap[feSpace][colsOther[k]].local * nComponents + j;
+ colsOther[k] = primalDofMap.getMatIndex(j, colsOther[k]);
MatSetValues(mat_b_primal, 1, &globalRowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
@@ -1077,10 +1069,10 @@ namespace AMDiS {
switch (fetiPreconditioner) {
case FETI_DIRICHLET:
{
- int rowIndex = localDofMap[feSpace][*cursor].local;
+ int rowIndex = localDofMap[feSpaces[i]][*cursor].local;
- if (rowIndex < nLocalInterior[feSpace]) {
- int rowIndex2 = localDofMap[feSpace][*cursor].local * nComponents + i;
+ if (rowIndex < nLocalInterior[feSpaces[i]]) {
+ int rowIndex2 = localDofMap[feSpaces[i]][*cursor].local * nComponents + i;
MatSetValues(mat_interior_interior, 1, &rowIndex2, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
@@ -1090,7 +1082,7 @@ namespace AMDiS {
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
} else {
int rowIndex2 =
- (localDofMap[feSpace][*cursor].local - nLocalInterior[feSpace]) * nComponents + i;
+ (localDofMap[feSpaces[i]][*cursor].local - nLocalInterior[feSpaces[i]]) * nComponents + i;
MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
@@ -1105,11 +1097,11 @@ namespace AMDiS {
case FETI_LUMPED:
{
- int rowIndex = localDofMap[feSpace][*cursor].local;
+ int rowIndex = localDofMap[feSpaces[i]][*cursor].local;
- if (rowIndex >= nLocalInterior[feSpace]) {
+ if (rowIndex >= nLocalInterior[feSpaces[i]]) {
int rowIndex2 =
- (localDofMap[feSpace][*cursor].local - nLocalInterior[feSpace]) * nComponents + i;
+ (localDofMap[feSpaces[i]][*cursor].local - nLocalInterior[feSpaces[i]]) * nComponents + i;
MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
@@ -1187,8 +1179,6 @@ namespace AMDiS {
FUNCNAME("PetscSolverFeti::fillPetscRhs()");
vector<const FiniteElemSpace*> feSpaces = getFeSpaces(vec);
- const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
- int nComponents = vec->getSize();
VecCreateMPI(PETSC_COMM_WORLD,
localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &f_b);
@@ -1196,14 +1186,13 @@ namespace AMDiS {
primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
&f_primal);
- for (int i = 0; i < nComponents; i++) {
+ for (unsigned int i = 0; i < feSpaces.size(); i++) {
DOFVector<double>::Iterator dofIt(vec->getDOFVector(i), USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
int index = dofIt.getDOFIndex();
- if (isPrimal(feSpace, index)) {
- index = primalDofMap[feSpace][index].local * nComponents + i;
- double value = *dofIt;
- VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
+ if (isPrimal(feSpaces[i], index)) {
+ index = primalDofMap.getMatIndex(i, index);
+ VecSetValue(f_primal, index, *dofIt, ADD_VALUES);
} else {
index = localDofMap.mapGlobal(index, i);
VecSetValue(f_b, index, *dofIt, INSERT_VALUES);
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index b9d91767..a2dcadd2 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -151,9 +151,6 @@ namespace AMDiS {
}
protected:
- /// Number of components in the PDE to be solved.
- int nComponents;
-
/// Mapping from primal DOF indices to a global index of primals.
FeSpaceData<GlobalDofMap> primalDofMap;
--
GitLab