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

Fixed some bugs in FETI-DP Stokes mode implementation.

parent 3acdc7cb
...@@ -186,7 +186,6 @@ namespace AMDiS { ...@@ -186,7 +186,6 @@ namespace AMDiS {
int nRankCoarseRows = cMap->getRankDofs(); int nRankCoarseRows = cMap->getRankDofs();
int nOverallCoarseRows = cMap->getOverallDofs(); int nOverallCoarseRows = cMap->getOverallDofs();
MSG("COARSE MAT %d SIZE %d %d\n", i + 1, nOverallCoarseRows, nOverallCoarseRows);
MatCreateAIJ(mpiCommGlobal, MatCreateAIJ(mpiCommGlobal,
nRankCoarseRows, nRankCoarseRows, nRankCoarseRows, nRankCoarseRows,
nOverallCoarseRows, nOverallCoarseRows, nOverallCoarseRows, nOverallCoarseRows,
...@@ -212,7 +211,6 @@ namespace AMDiS { ...@@ -212,7 +211,6 @@ namespace AMDiS {
int nColsRankMat = (j == 0 ? nRankInteriorRows : uniqueCoarseMap[j - 1]->getRankDofs()); int nColsRankMat = (j == 0 ? nRankInteriorRows : uniqueCoarseMap[j - 1]->getRankDofs());
int nColsOverallMat = (j == 0 ? nOverallInteriorRows : uniqueCoarseMap[j - 1]->getOverallDofs()); int nColsOverallMat = (j == 0 ? nOverallInteriorRows : uniqueCoarseMap[j - 1]->getOverallDofs());
MSG("COARSE MAT %d %d SIZE %d %d\n", i, j, nRowsOverallMat, nColsOverallMat);
MatCreateAIJ(mpiCommGlobal, MatCreateAIJ(mpiCommGlobal,
nRowsRankMat, nColsRankMat, nRowsRankMat, nColsRankMat,
nRowsOverallMat, nColsOverallMat, nRowsOverallMat, nColsOverallMat,
...@@ -221,7 +219,6 @@ namespace AMDiS { ...@@ -221,7 +219,6 @@ namespace AMDiS {
MSG("REMOVE THIS!\n"); MSG("REMOVE THIS!\n");
MatSetOption(mat[i][j], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); MatSetOption(mat[i][j], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
MSG("COARSE MAT %d %d SIZE %d %d\n", j, i, nRowsOverallMat, nColsOverallMat);
MatCreateAIJ(mpiCommGlobal, MatCreateAIJ(mpiCommGlobal,
nColsRankMat, nRowsRankMat, nColsRankMat, nRowsRankMat,
nColsOverallMat, nRowsOverallMat, nColsOverallMat, nRowsOverallMat,
......
...@@ -167,6 +167,13 @@ namespace AMDiS { ...@@ -167,6 +167,13 @@ namespace AMDiS {
return vecRhs[coarseSpace + 1]; return vecRhs[coarseSpace + 1];
} }
/// Get the RHS vector of the coarse space of a system component.
inline Vec& getVecRhsCoarseByComponent(int comp)
{
int vecIndex = componentIthCoarseMap[comp] + 1;
return vecRhs[vecIndex];
}
/// Get the solution vector of the interior. /// Get the solution vector of the interior.
inline Vec& getVecSolInterior() inline Vec& getVecSolInterior()
{ {
...@@ -179,6 +186,13 @@ namespace AMDiS { ...@@ -179,6 +186,13 @@ namespace AMDiS {
return vecSol[coarseSpace + 1]; return vecSol[coarseSpace + 1];
} }
/// Get the solution vector of the coarse space of a system component.
inline Vec& getVecSolCoarseByComponent(int comp)
{
int vecIndex = componentIthCoarseMap[comp] + 1;
return vecSol[vecIndex];
}
/** \brief /** \brief
* Checks whether a given DOF index in some component is a coarse space DOF. * Checks whether a given DOF index in some component is a coarse space DOF.
* Note (TODO): The specification of both, the component number and FE * Note (TODO): The specification of both, the component number and FE
......
...@@ -530,13 +530,10 @@ namespace AMDiS { ...@@ -530,13 +530,10 @@ namespace AMDiS {
for (DofContainer::iterator it = allBoundaryDofs.begin(); for (DofContainer::iterator it = allBoundaryDofs.begin();
it != allBoundaryDofs.end(); ++it) it != allBoundaryDofs.end(); ++it)
if (meshDistributor->getDofMap()[feSpace].isRankDof(**it)) { if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
MSG("INSERT INTERFACE RANK DOF %d\n", **it);
interfaceDofMap[feSpace].insertRankDof(**it); interfaceDofMap[feSpace].insertRankDof(**it);
} else { else
MSG("INSERT INTERFACE NON-RANK DOF %d\n", **it);
interfaceDofMap[feSpace].insertNonRankDof(**it); interfaceDofMap[feSpace].insertNonRankDof(**it);
}
} }
......
...@@ -136,8 +136,6 @@ namespace AMDiS { ...@@ -136,8 +136,6 @@ namespace AMDiS {
if (!dofMat) if (!dofMat)
continue; continue;
MSG("=============== COMP %d %d ================= \n", rowComponent, colComponent);
ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent]; ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent]; ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
...@@ -163,10 +161,6 @@ namespace AMDiS { ...@@ -163,10 +161,6 @@ namespace AMDiS {
bool isColCoarse = bool isColCoarse =
isCoarseSpace(colComponent, feSpaces[colComponent], col(*icursor)); isCoarseSpace(colComponent, feSpaces[colComponent], col(*icursor));
if (colComponent == 2 && col(*icursor) == 19) {
MSG("HERE WE TEST: %d %d\n", isColCoarse, isRowCoarse);
}
if (isColCoarse) { if (isColCoarse) {
if (isRowCoarse) { if (isRowCoarse) {
cols.push_back(col(*icursor)); cols.push_back(col(*icursor));
...@@ -190,71 +184,56 @@ namespace AMDiS { ...@@ -190,71 +184,56 @@ namespace AMDiS {
// === Set matrix values. === // === Set matrix values. ===
if (isRowCoarse) { if (isRowCoarse) {
MSG("A-GET MAT INDEX %d\n", rowComponent);
int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, *cursor); int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, *cursor);
MSG("B-GET MAT INDEX %d\n", colComponent);
for (unsigned int i = 0; i < cols.size(); i++) for (unsigned int i = 0; i < cols.size(); i++)
cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]); cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
// matCoarseCoarse
//MSG("SET COARSE MAT COMP %d %d WITH %d entry\n", rowComponent, colComponent, cols.size());
MatSetValues(getMatCoarseByComponent(rowComponent, colComponent), MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
1, &rowIndex, cols.size(), 1, &rowIndex, cols.size(),
&(cols[0]), &(values[0]), ADD_VALUES); &(cols[0]), &(values[0]), ADD_VALUES);
if (colsOther.size()) { if (colsOther.size()) {
if (subdomainLevel == 0) { if (subdomainLevel == 0) {
MSG("C0-GET MAT INDEX %d\n", colComponent);
for (unsigned int i = 0; i < colsOther.size(); i++) for (unsigned int i = 0; i < colsOther.size(); i++)
colsOther[i] = colsOther[i] =
interiorMap->getMatIndex(colComponent, colsOther[i]); interiorMap->getMatIndex(colComponent, colsOther[i]);
} else { } else {
MSG("C1-GET MAT INDEX %d\n", colComponent);
for (unsigned int i = 0; i < colsOther.size(); i++) for (unsigned int i = 0; i < colsOther.size(); i++)
colsOther[i] = colsOther[i] =
interiorMap->getMatIndex(colComponent, colsOther[i]) + interiorMap->getMatIndex(colComponent, colsOther[i]) +
rStartInterior; rStartInterior;
} }
// matCoarseInt
// MSG("SET COARSE-INT MAT COMP %d %d WITH %d entry\n", rowComponent, colComponent, cols.size());
MatSetValues(getMatCoarseInteriorByComponent(rowComponent), MatSetValues(getMatCoarseInteriorByComponent(rowComponent),
1, &rowIndex, colsOther.size(), 1, &rowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES); &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
} }
} else { } else {
MSG("D-GET MAT INDEX %d\n", rowComponent);
int localRowIndex = int localRowIndex =
(subdomainLevel == 0 ? (subdomainLevel == 0 ?
interiorMap->getLocalMatIndex(rowComponent, *cursor) : interiorMap->getLocalMatIndex(rowComponent, *cursor) :
interiorMap->getMatIndex(rowComponent, *cursor)); interiorMap->getMatIndex(rowComponent, *cursor));
MSG("E-GET MAT INDEX %d\n", colComponent);
for (unsigned int i = 0; i < cols.size(); i++) { for (unsigned int i = 0; i < cols.size(); i++) {
if (subdomainLevel == 0) if (subdomainLevel == 0)
cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]); cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]);
else else
cols[i] = interiorMap->getMatIndex(colComponent, cols[i]); cols[i] = interiorMap->getMatIndex(colComponent, cols[i]);
} }
MSG("DONE\n");
MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(), MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(),
&(cols[0]), &(values[0]), ADD_VALUES); &(cols[0]), &(values[0]), ADD_VALUES);
if (colsOther.size()) { if (colsOther.size()) {
MSG("F-GET MAT INDEX %d\n", colComponent);
int globalRowIndex = interiorMap->getMatIndex(rowComponent, *cursor); int globalRowIndex = interiorMap->getMatIndex(rowComponent, *cursor);
if (subdomainLevel != 0) if (subdomainLevel != 0)
globalRowIndex += rStartInterior; globalRowIndex += rStartInterior;
MSG("G-GET MAT INDEX %d\n", colComponent);
for (unsigned int i = 0; i < colsOther.size(); i++) for (unsigned int i = 0; i < colsOther.size(); i++)
colsOther[i] = colsOther[i] =
colCoarseSpace->getMatIndex(colComponent, colsOther[i]); colCoarseSpace->getMatIndex(colComponent, colsOther[i]);
// matIntCoarse
// MSG("SET INT-COARSE MAT COMP %d %d WITH %d entry\n", rowComponent, colComponent, cols.size());
MatSetValues(getMatInteriorCoarseByComponent(colComponent), MatSetValues(getMatInteriorCoarseByComponent(colComponent),
1, &globalRowIndex, colsOther.size(), 1, &globalRowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES); &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
...@@ -293,7 +272,8 @@ namespace AMDiS { ...@@ -293,7 +272,8 @@ namespace AMDiS {
// === Transfer values from DOF vector to the PETSc vector. === // === Transfer values from DOF vector to the PETSc vector. ===
if (coarseSpaceMap.size()) { if (coarseSpaceMap.size()) {
for (int i = 0; i < vec->getSize(); i++) for (int i = 0; i < vec->getSize(); i++)
setDofVector(getVecRhsInterior(), getVecRhsCoarse(), vec->getDOFVector(i), i); setDofVector(getVecRhsInterior(),
getVecSolCoarseByComponent(i), vec->getDOFVector(i), i);
} else { } else {
for (int i = 0; i < vec->getSize(); i++) for (int i = 0; i < vec->getSize(); i++)
setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i); setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i);
...@@ -499,7 +479,7 @@ namespace AMDiS { ...@@ -499,7 +479,7 @@ namespace AMDiS {
blockComponents); blockComponents);
int nComponents = static_cast<int>(blockComponents.size()); int nComponents = static_cast<int>(blockComponents.size());
TEST_EXIT(nComponents > 0)("No is block for block %d defined!\n", i); TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
// Check if blocks are continous // Check if blocks are continous
for (int j = 0; j < nComponents; j++) { for (int j = 0; j < nComponents; j++) {
...@@ -709,9 +689,8 @@ namespace AMDiS { ...@@ -709,9 +689,8 @@ namespace AMDiS {
const FiniteElemSpace *feSpace = vec->getFeSpace(); const FiniteElemSpace *feSpace = vec->getFeSpace();
PeriodicMap &perMap = meshDistributor->getPeriodicMap(); PeriodicMap &perMap = meshDistributor->getPeriodicMap();
ParallelDofMapping *rowCoarseSpace = NULL; ParallelDofMapping *rowCoarseSpace =
if (coarseSpaceMap.size()) (coarseSpaceMap.size() ? coarseSpaceMap[nRowVec] : NULL);
rowCoarseSpace = coarseSpaceMap[nRowVec];
// Traverse all used DOFs in the dof vector. // Traverse all used DOFs in the dof vector.
DOFVector<double>::Iterator dofIt(vec, USED_DOFS); DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
......
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