Skip to content

GitLab

Commit e77df268 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

FETI-DP first version for different FE space, still not tested.

parent ab968408
Hide whitespace changes
Inline Side-by-side
AMDiS/src/parallel/FeSpaceMapping.cc
View file @ e77df268
......@@ -37,20 +37,20 @@ namespace AMDiS {
mpi::getDofNumbering(*mpiComm, nRankDofs, rStartDofs, nOverallDofs);
if (needGlobalMapping)
addOffset(rStartDofs);
computeGlobalMapping(rStartDofs);
if (overlap)
if (hasNonLocalDofs)
computeNonLocalIndices();
}
void GlobalDofMap::addOffset(int offset)
void GlobalDofMap::computeGlobalMapping(int offset)
{
for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
it->second.global = it->second.local + offset;
// it->second.local += offset;
}
void GlobalDofMap::computeNonLocalIndices()
{
FUNCNAME("GlobalDofMap::computeNonLocalIndices()");
......
AMDiS/src/parallel/FeSpaceMapping.h
View file @ e77df268
......@@ -106,7 +106,7 @@ namespace AMDiS {
nLocalDofs(0),
nOverallDofs(0),
rStartDofs(0),
overlap(false)
hasNonLocalDofs(false)
{}
void clear();
......@@ -163,7 +163,7 @@ namespace AMDiS {
void update();
void addOffset(int offset);
void computeGlobalMapping(int offset);
void computeNonLocalIndices();
......@@ -174,9 +174,9 @@ namespace AMDiS {
feSpace = fe;
}
void setOverlap(bool b)
void setNonLocalDofs(bool b)
{
overlap = b;
hasNonLocalDofs = b;
}
void setNeedGlobalMapping(bool b)
......@@ -209,7 +209,10 @@ namespace AMDiS {
///
int nRankDofs, nLocalDofs, nOverallDofs, rStartDofs;
bool overlap;
/// Is true if there are DOFs in at least one subdomain that are not owned
/// by the rank. If the value is false, each rank contains only DOFs that
/// are also owned by this rank.
bool hasNonLocalDofs;
};
......@@ -221,7 +224,7 @@ namespace AMDiS {
: mpiComm(NULL),
sendDofs(NULL),
recvDofs(NULL),
overlap(false),
hasNonLocalDofs(false),
feSpaces(0),
nRankDofs(-1),
nOverallDofs(-1),
......@@ -328,17 +331,17 @@ namespace AMDiS {
void init(MPI::Intracomm *m,
vector<const FiniteElemSpace*> &fe,
bool needGlobalMapping,
bool ol)
bool bNonLocalDofs)
{
mpiComm = m;
feSpaces = fe;
overlap = ol;
hasNonLocalDofs = bNonLocalDofs;
for (unsigned int i = 0; i < feSpaces.size(); i++) {
feSpacesUnique.insert(feSpaces[i]);
addFeSpace(feSpaces[i]);
data[feSpaces[i]].setNeedGlobalMapping(needGlobalMapping);
data[feSpaces[i]].setOverlap(overlap);
data[feSpaces[i]].setNonLocalDofs(hasNonLocalDofs);
}
}
......@@ -370,14 +373,13 @@ namespace AMDiS {
for (ItType it = dofMap.begin(); it != dofMap.end(); ++it) {
if (data[feSpaces[i]].isRankDof(it->first)) {
int globalMatIndex = it->second.local + offset;
// int globalMatIndex =
// it->second.local - data[feSpaces[i]].rStartDofs + offset;
dofToMatIndex.add(i, it->first, globalMatIndex);
}
}
offset += data[feSpaces[i]].nRankDofs;
if (!overlap)
if (!hasNonLocalDofs)
continue;
TEST_EXIT_DBG(sendDofs != NULL && recvDofs != NULL)
......@@ -413,50 +415,9 @@ namespace AMDiS {
}
}
}
offset += data[feSpaces[i]].nRankDofs;
}
}
inline int mapLocal(int index, int ithFeSpace)
{
int result = 0;
for (int i = 0; i < ithFeSpace; i++)
result += data[feSpaces[i]].nRankDofs;
result += data[feSpaces[ithFeSpace]][index].local;
return result;
}
inline int mapLocal(int index, const FiniteElemSpace *feSpace)
{
for (unsigned int i = 0; i < feSpaces.size(); i++)
if (feSpaces[i] == feSpace)
return mapLocal(index, feSpace, i);
return -1;
}
inline int mapGlobal(int index, int ithFeSpace)
{
int result = rStartDofs;
for (int i = 0; i < ithFeSpace; i++)
result += data[feSpaces[i]].nRankDofs;
result += data[feSpaces[ithFeSpace]][index].local;
return result;
}
inline int mapGlobal(int index, const FiniteElemSpace *feSpace)
{
for (unsigned int i = 0; i < feSpaces.size(); i++)
if (feSpaces[i] == feSpace)
return mapGlobal(index, feSpace, i);
return -1;
}
void setDofComm(DofComm &pSend, DofComm &pRecv)
{
sendDofs = &pSend;
......@@ -469,11 +430,19 @@ namespace AMDiS {
inline int getMatIndex(int ithComponent, DegreeOfFreedom d)
{
FUNCNAME("FeSpaceData::getMatIndex()");
return dofToMatIndex.get(ithComponent, d);
}
inline int getLocalMatIndex(int ithComponent, DegreeOfFreedom d)
{
FUNCNAME("FeSpaceData::getLocalMatIndex()");
TEST_EXIT_DBG(data[feSpaces[ithComponent]].isRankDof(d))
("Should not happen!\n");
return dofToMatIndex.get(ithComponent, d) - rStartDofs;
}
private:
MPI::Intracomm* mpiComm;
......@@ -481,7 +450,10 @@ namespace AMDiS {
DofComm *recvDofs;
bool overlap;
/// Is true if there are DOFs in at least one subdomain that are not owned
/// by the rank. If the value is false, each rank contains only DOFs that
/// are also owned by this rank.
bool hasNonLocalDofs;
map<const FiniteElemSpace*, T> data;
......
AMDiS/src/parallel/PetscSolverFeti.cc
View file @ e77df268
......@@ -232,6 +232,8 @@ namespace AMDiS {
dualDofMap.update();
lagrangeMap.update();
localDofMap.update();
if (fetiPreconditioner == FETI_DIRICHLET)
interiorDofMap.update();
for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
......@@ -240,12 +242,10 @@ namespace AMDiS {
primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs);
MSG("nRankDuals = %d nOverallDuals = %d\n",
dualDofMap[feSpace].nRankDofs,
dualDofMap[feSpace].nOverallDofs);
dualDofMap[feSpace].nRankDofs, dualDofMap[feSpace].nOverallDofs);
MSG("nRankLagrange = %d nOverallLagrange = %d\n",
lagrangeMap[feSpace].nRankDofs,
lagrangeMap[feSpace].nOverallDofs);
lagrangeMap[feSpace].nRankDofs, lagrangeMap[feSpace].nOverallDofs);
TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() ==
static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
......@@ -306,14 +306,14 @@ namespace AMDiS {
// === Create for each dual node that is owned by the rank, the set ===
// === of ranks that contain this node (denoted by W(x_j)). ===
boundaryDofRanks.clear();
boundaryDofRanks[feSpace].clear();
for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
!it.end(); it.nextRank())
for (; !it.endDofIter(); it.nextDof()) {
if (!isPrimal(feSpace, it.getDofIndex())) {
boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);
boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
}
}
......@@ -327,7 +327,7 @@ namespace AMDiS {
!it.end(); it.nextRank())
for (; !it.endDofIter(); it.nextDof())
if (!isPrimal(feSpace, it.getDofIndex()))
stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
stdMpi.updateSendDataSize();
......@@ -352,7 +352,7 @@ namespace AMDiS {
int i = 0;
for (; !it.endDofIter(); it.nextDof())
if (!isPrimal(feSpace, it.getDofIndex()))
boundaryDofRanks[it.getDofIndex()] =
boundaryDofRanks[feSpace][it.getDofIndex()] =
stdMpi.getRecvData(it.getRank())[i++];
}
......@@ -365,7 +365,7 @@ namespace AMDiS {
for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
if (meshDistributor->getIsRankDof(feSpace, it->first)) {
lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
int degree = boundaryDofRanks[it->first].size();
int degree = boundaryDofRanks[feSpace][it->first].size();
nRankLagrange += (degree * (degree - 1)) / 2;
} else {
lagrangeMap[feSpace].insert(it->first);
......@@ -385,13 +385,17 @@ namespace AMDiS {
// === the global index of all B nodes, insert all interior nodes first, ===
// === without defining a correct index. ===
nLocalInterior[feSpace] = 0;
int nLocalInterior = 0;
for (int i = 0; i < admin->getUsedSize(); i++) {
if (admin->isDofFree(i) == false &&
isPrimal(feSpace, i) == false &&
dualDofMap[feSpace].isSet(i) == false) {
localDofMap[feSpace].insertRankDof(i, nLocalInterior[feSpace]);
nLocalInterior[feSpace]++;
localDofMap[feSpace].insertRankDof(i, nLocalInterior);
if (fetiPreconditioner == FETI_DIRICHLET)
interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
nLocalInterior++;
}
}
......@@ -428,25 +432,24 @@ namespace AMDiS {
for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin();
it != dualMap.end(); ++it) {
TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].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());
vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(),
boundaryDofRanks[feSpaces[k]][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);
int colIndex = localDofMap.getMatIndex(k, it->first);
double value = (W[i] == mpiRank ? 1.0 : -1.0);
MatSetValue(mat_lagrange, index, colIndex, value,
INSERT_VALUES);
MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
}
index++;
}
......@@ -831,7 +834,7 @@ namespace AMDiS {
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);
int petscIndex = localDofMap.getLocalMatIndex(i, it->first);
dofVec[it->first] = localSolB[petscIndex];
}
......@@ -849,65 +852,70 @@ namespace AMDiS {
void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
{
FUNCNAME("PetscSolverFeti::fillPetscMatrix()");
// === 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);
localDofMap.init(mpiComm, feSpaces, false, false);
if (fetiPreconditioner == FETI_DIRICHLET)
interiorDofMap.init(mpiComm, feSpaces, false, false);
updateDofData();
// === Create matrices for the FETI-DP method. ===
int nRowsRankB = localDofMap.getRankDofs();
int nRowsOverallB = localDofMap.getOverallDofs();
int nRowsRankPrimal = primalDofMap.getRankDofs();
int nRowsOverallPrimal = primalDofMap.getOverallDofs();
int nRowsDual = dualDofMap.getRankDofs();
int nRowsInterior = nRowsRankB - nRowsDual;
MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
&mat_b_b);
MatCreateMPIAIJ(PETSC_COMM_WORLD,
nRowsRankPrimal, nRowsRankPrimal,
nRowsOverallPrimal, nRowsOverallPrimal,
30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
MatCreateMPIAIJ(PETSC_COMM_WORLD,
nRowsRankB, nRowsRankPrimal,
nRowsOverallB, nRowsOverallPrimal,
30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
MatCreateMPIAIJ(PETSC_COMM_WORLD,
nRowsRankPrimal, nRowsRankB,
nRowsOverallPrimal, nRowsOverallB,
15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
// === Create matrices for FETI-DP preconditioner. ===
if (fetiPreconditioner != FETI_NONE) {
int nRowsDual = dualDofMap.getRankDofs();
if (fetiPreconditioner != FETI_NONE)
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsDual, nRowsDual, 30, PETSC_NULL,
&mat_duals_duals);
if (fetiPreconditioner == FETI_DIRICHLET) {
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsInterior, nRowsInterior, 30, PETSC_NULL,
&mat_interior_interior);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsInterior, nRowsDual, 30, PETSC_NULL,
&mat_interior_duals);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsDual, nRowsInterior, 30, PETSC_NULL,
&mat_duals_interior);
if (fetiPreconditioner == FETI_DIRICHLET) {
int nRowsInterior = interiorDofMap.getRankDofs();
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsInterior, nRowsInterior, 30, PETSC_NULL,
&mat_interior_interior);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsInterior, nRowsDual, 30, PETSC_NULL,
&mat_interior_duals);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsDual, nRowsInterior, 30, PETSC_NULL,
&mat_duals_interior);
}
}
......@@ -992,32 +1000,24 @@ namespace AMDiS {
// === For preconditioner ===
if (!rowPrimal && !colPrimal) {
int rowIndex = localDofMap[feSpaces[i]][*cursor].local;
int colIndex = localDofMap[feSpaces[j]][col(*icursor)].local;
if (rowIndex < nLocalInterior[feSpaces[i]]) {
if (colIndex < nLocalInterior[feSpaces[j]]) {
int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
colsLocal.push_back(colIndex2);
if (!isDual(feSpaces[i], *cursor)) {
if (!isDual(feSpaces[j], col(*icursor))) {
int colIndex = interiorDofMap.getLocalMatIndex(j, col(*icursor));
colsLocal.push_back(colIndex);
valuesLocal.push_back(value(*icursor));
} else {
int colIndex2 =
(localDofMap[feSpaces[j]][col(*icursor)].local - nLocalInterior[feSpaces[j]]) *
nComponents + j;
colsLocalOther.push_back(colIndex2);
int colIndex = dualDofMap.getLocalMatIndex(j, col(*icursor));
colsLocalOther.push_back(colIndex);
valuesLocalOther.push_back(value(*icursor));
}
} else {
if (colIndex < nLocalInterior[feSpaces[j]]) {
int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
colsLocalOther.push_back(colIndex2);
if (!isDual(feSpaces[j], col(*icursor))) {
int colIndex = interiorDofMap.getLocalMatIndex(j, col(*icursor));
colsLocalOther.push_back(colIndex);
valuesLocalOther.push_back(value(*icursor));
} else {
int colIndex2 = (localDofMap[feSpaces[j]][col(*icursor)].local - nLocalInterior[feSpaces[j]]) *
nComponents + j;
colsLocal.push_back(colIndex2);
int colIndex = dualDofMap.getLocalMatIndex(j, col(*icursor));
colsLocal.push_back(colIndex);
valuesLocal.push_back(value(*icursor));
}
}
......@@ -1039,22 +1039,22 @@ namespace AMDiS {
if (colsOther.size()) {
for (unsigned int k = 0; k < colsOther.size(); k++)
colsOther[k] = localDofMap.mapGlobal(colsOther[k], j);
colsOther[k] = localDofMap.getMatIndex(j, colsOther[k]);
MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
}
} else {
int localRowIndex = localDofMap.mapLocal(*cursor, i);
int localRowIndex = localDofMap.getLocalMatIndex(i, *cursor);
for (unsigned int k = 0; k < cols.size(); k++)
cols[k] = localDofMap.mapLocal(cols[k], j);
cols[k] = localDofMap.getLocalMatIndex(j, cols[k]);
MatSetValues(mat_b_b, 1, &localRowIndex, cols.size(),
&(cols[0]), &(values[0]), ADD_VALUES);
if (colsOther.size()) {
int globalRowIndex = localDofMap.mapGlobal(*cursor, i);
int globalRowIndex = localDofMap.getMatIndex(i, *cursor);
for (unsigned int k = 0; k < colsOther.size(); k++)
colsOther[k] = primalDofMap.getMatIndex(j, colsOther[k]);
......@@ -1068,44 +1068,28 @@ namespace AMDiS {
if (!rowPrimal) {
switch (fetiPreconditioner) {
case FETI_DIRICHLET:
{
int rowIndex = localDofMap[feSpaces[i]][*cursor].local;
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);
if (colsLocalOther.size())
MatSetValues(mat_interior_duals, 1, &rowIndex2, colsLocalOther.size(),
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
} else {
int rowIndex2 =
(localDofMap[feSpaces[i]][*cursor].local - nLocalInterior[feSpaces[i]]) * nComponents + i;
MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
if (colsLocalOther.size())
MatSetValues(mat_duals_interior, 1, &rowIndex2, colsLocalOther.size(),
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
}
if (!isDual(feSpaces[i], *cursor)) {
int rowIndex = interiorDofMap.getLocalMatIndex(i, *cursor);
MatSetValues(mat_interior_interior, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
if (colsLocalOther.size())
MatSetValues(mat_interior_duals, 1, &rowIndex, colsLocalOther.size(),
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
} else {
int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor);
MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
if (colsLocalOther.size())
MatSetValues(mat_duals_interior, 1, &rowIndex, colsLocalOther.size(),
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
}
break;
case FETI_LUMPED:
{
int rowIndex = localDofMap[feSpaces[i]][*cursor].local;
if (rowIndex >= nLocalInterior[feSpaces[i]]) {
int rowIndex2 =
(localDofMap[feSpaces[i]][*cursor].local - nLocalInterior[feSpaces[i]]) * nComponents + i;
MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
}
if (isDual(feSpaces[i], *cursor)) {
int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor);
MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
}
break;
default:
......@@ -1194,7 +1178,7 @@ namespace AMDiS {
index = primalDofMap.getMatIndex(i, index);
VecSetValue(f_primal, index, *dofIt, ADD_VALUES);
} else {
index = localDofMap.mapGlobal(index, i);
index = localDofMap.getMatIndex(i, index);
VecSetValue(f_b, index, *dofIt, INSERT_VALUES);
}
}
......
AMDiS/src/parallel/PetscSolverFeti.h
View file @ e77df268
......@@ -144,12 +144,18 @@ namespace AMDiS {
void printStatistics();
/// Checks whether a given DOF in a gifen FE space is a primal DOF.
/// Checks whether a given DOF in a given FE space is a primal DOF.
inline bool isPrimal(const FiniteElemSpace *feSpace, DegreeOfFreedom dof)
{
return primalDofMap[feSpace].isSet(dof);
}
/// Checks whether a given DOF in a give FE space is a dual DOF.
inline bool isDual(const FiniteElemSpace *feSpace, DegreeOfFreedom dof)
{
return dualDofMap[feSpace].isSet(dof);
}
protected:
/// Mapping from primal DOF indices to a global index of primals.
FeSpaceData<GlobalDofMap> primalDofMap;
......@@ -161,13 +167,16 @@ namespace AMDiS {
/// constraint that is assigned to this DOF.
FeSpaceData<GlobalDofMap> lagrangeMap;
/// Index for each non primal DOF to the global index of
/// B variables.
/// Index for each non primal DOF to the global index of B variables.
FeSpaceData<GlobalDofMap> localDofMap;
/// Stores to each dual boundary DOF the set of ranks in which the DOF
/// is contained in.
DofIndexToPartitions boundaryDofRanks;
/// Mapping of pure local DOF indices, thus no primal and no dual DOFs are
/// in this map. Is used for the Dirichlet preconditioner only.
FeSpaceData<GlobalDofMap> interiorDofMap;