diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 7f6d9754336db6bd0a4901f2069cd532227d17f4..6104d8ee28a9fe9ca219d8b8c9c9e38dec0ef5d1 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -238,18 +238,18 @@ namespace AMDiS {
if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- if (dofMap->isRankDof(rowIndices[i]))
- applyDBCs.insert(static_cast<int>(row));
+ // if (dofMap->isRankDof(rowIndices[i])) {
+ applyDBCs.insert(row);
+ // dirichletDofs.push_back(row);
+ // }
#else
- applyDBCs.insert(static_cast<int>(row));
+ applyDBCs.insert(row);
#endif
}
} else {
for (int j = 0; j < nCol; j++) {
DegreeOfFreedom col = colIndices[j];
- if (fabs(elMat[i][j]) > 1e-10) {
- ins[row][col] += elMat[i][j];
- }
+ ins[row][col] += elMat[i][j];
}
}
}
@@ -493,11 +493,16 @@ namespace AMDiS {
{
FUNCNAME("DOFMatrix::removeRowsWithDBC()");
+ // Do the following only in sequential code. In parallel mode, the specific
+ // solver method must care about dirichlet boundary conditions.
+
+#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
inserter_type &ins = *inserter;
for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
ins[*it][*it] = 1.0;
rows->clear();
+#endif
}
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index 5105828e941cfe9c67adec713fb4fe00925ebb22..84410b014e334742cb75e8936ba98dc78b8e047f 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -455,27 +455,25 @@ namespace AMDiS {
/// Maps local col indices of an element to global matrix indices.
vector<DegreeOfFreedom> colIndices;
- /* \brief
- * A set of row indices. When assembling the DOF matrix, all rows, that
- * correspond to a dof at a dirichlet boundary, are ignored and the row is
- * left blank. After assembling, the diagonal element of the matrix must be
- * set to 1. The indices of all rows, where this should be done, are stored
- * in this set.
- */
+ /// A set of row indices. When assembling the DOF matrix, all rows, that
+ /// correspond to a dof at a dirichlet boundary, are ignored and the row is
+ /// left blank. After assembling, the diagonal element of the matrix must be
+ /// set to 1. The indices of all rows, where this should be done, are stored
+ /// in this set.
std::set<int> applyDBCs;
- /* \brief
- * Number of non zero entries per row (average). For instationary problems this
- * information may be used in the next timestep to accelerate insertion of
- * elemnts during assembling.
- */
- int nnzPerRow;
+ /// Number of non zero entries per row (average). For instationary problems this
+ /// information may be used in the next timestep to accelerate insertion of
+ /// elemnts during assembling.
+ int nnzPerRow;
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
/// Is used in parallel computations to check whether specific DOFs in the
/// row FE spaces are owned by the rank or not. This is used to ensure that
/// Dirichlet BC is handled correctly in parallel computations.
FeSpaceDofMap *dofMap;
+
+ // std::set<DegreeOfFreedom> dirichletDofs;
#endif
/// Inserter object: implemented as pointer, allocated and deallocated as needed
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 36ac7118450e2caff0f6fe525265fccc09080b12..e05c35a45c5eca0519368ae1a97164cbdba7ded8 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -219,10 +219,10 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
- #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
double localResult = result;
MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
- #endif
+#endif
return result;
}
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 41afc730400dd5447dff825fc3b5a2d7aeccbfc0..d334ab56116745e8bfe735d28d6b2e49166a61bb 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -250,6 +250,17 @@ namespace AMDiS {
return dofMap->isRankDof(dof);
}
+
+ inline void setDirichletDofValue(DegreeOfFreedom dof,
+ T value)
+ {
+ dirichletDofValues[dof] = value;
+ }
+
+ map<DegreeOfFreedom, T>& getDirichletValues()
+ {
+ return dirichletDofValues;
+ }
#endif
protected:
@@ -282,10 +293,12 @@ namespace AMDiS {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
FeSpaceDofMap *dofMap;
-#endif
+ map<DegreeOfFreedom, T> dirichletDofValues;
+#endif
};
+
/// Specifies which operation should be done after coarsening
typedef enum{
NO_OPERATION = 0,
diff --git a/AMDiS/src/DirichletBC.cc b/AMDiS/src/DirichletBC.cc
index a739647d71268304ad26490c7ac3ebdb016efb67..a5f5ec79e2fe4d5d572d067317e4e341f0b7bb22 100644
--- a/AMDiS/src/DirichletBC.cc
+++ b/AMDiS/src/DirichletBC.cc
@@ -65,15 +65,22 @@ namespace AMDiS {
for (int i = 0; i < nBasFcts; i++) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- if (vector->isRankDof(dofIndices[i]))
+ // if (vector->isRankDof(dofIndices[i]))
#endif
if (localBound[i] == boundaryType) {
+ double value = 0.0;
if (f) {
elInfo->coordToWorld(*(basFcts->getCoords(i)), worldCoords);
- (*vector)[dofIndices[i]] = (*f)(worldCoords);
+ value = (*f)(worldCoords);
}
if (dofVec)
- (*vector)[dofIndices[i]] = (*dofVec)[dofIndices[i]];
+ value = (*dofVec)[dofIndices[i]];
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ vector->setDirichletDofValue(dofIndices[i], value);
+#else
+ (*vector)[dofIndices[i]] = value;
+#endif
}
}
}
diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc
index 07633a87c61a2fa880d83590c57f91245c7afeaa..f212abf6e18d54b7b80db9e04e34d4b705325418 100644
--- a/AMDiS/src/Element.cc
+++ b/AMDiS/src/Element.cc
@@ -76,13 +76,30 @@ namespace AMDiS {
}
- void Element::createNewDofPtrs()
+ void Element::createNewDofPtrs(bool setDofs)
{
FUNCNAME("Element::setDofPtrs()");
TEST_EXIT_DBG(mesh)("no mesh!\n");
dof = mesh->createDofPtrs();
+
+ if (setDofs) {
+ for (int i = 0; i < mesh->getGeo(VERTEX); i++)
+ dof[i] = mesh->getDof(VERTEX);
+
+ if (mesh->getNumberOfDofs(EDGE))
+ for (int i = 0; i < mesh->getGeo(EDGE); i++)
+ dof[mesh->getNode(EDGE) + i] = mesh->getDof(EDGE);
+
+ if (mesh->getNumberOfDofs(FACE))
+ for (int i = 0; i < mesh->getGeo(FACE); i++)
+ dof[mesh->getNode(FACE) + i] = mesh->getDof(FACE);
+
+ if (mesh->getNumberOfDofs(CENTER))
+ for (int i = 0; i < mesh->getGeo(CENTER); i++)
+ dof[mesh->getNode(CENTER) + i] = mesh->getDof(CENTER);
+ }
}
diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h
index a13d3c95f8df720337391eb342caac3cc6de0ecf..607eff4cf82684a5e1c4cff1ce002bac158df516 100644
--- a/AMDiS/src/Element.h
+++ b/AMDiS/src/Element.h
@@ -538,7 +538,7 @@ namespace AMDiS {
void deserialize(std::istream &in);
/// Sets Element's \ref dof pointer.
- void createNewDofPtrs();
+ void createNewDofPtrs(bool setDofs = false);
protected:
/// Sets Element's \ref index. Used by friend class Mesh.
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index 1bf6e97404efc6bdb9e13f9315fd22c5611dacc8..fa9d755d8b73108c70a99d54d184d0527de8aa96 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -803,7 +803,7 @@ namespace AMDiS {
DimVec<int> nDof;
/// Number of nodes on a single element where DOFs are located. Needed for
- /// the (de-) allocation of the dof-vector on the element (\ref Element::dof).
+ /// the (de-) allocation of the DOF-vector on the element (\ref Element::dof).
/// Here "node" is equivalent to the number of basis functions on the element.
int nNodeEl;
diff --git a/AMDiS/src/parallel/InteriorBoundary.cc b/AMDiS/src/parallel/InteriorBoundary.cc
index 5451e74aafede2c0e849eca1b1fd0c7fd430e467..f907b32cc437db8502169e101126cb3e620e5880 100644
--- a/AMDiS/src/parallel/InteriorBoundary.cc
+++ b/AMDiS/src/parallel/InteriorBoundary.cc
@@ -12,6 +12,7 @@
#include "parallel/InteriorBoundary.h"
#include "parallel/ElementObjectDatabase.h"
+#include "parallel/MeshDistributor.h"
#include "FiniteElemSpace.h"
#include "BasisFunction.h"
#include "Serializer.h"
@@ -177,8 +178,18 @@ namespace AMDiS {
bound.type = it->second;
- AtomicBoundary& b = getNewPeriodic(otherElementRank);
- b = bound;
+ if (MeshDistributor::sebastianMode) {
+ if (bound.rankObj.elIndex == 3 &&
+ bound.rankObj.ithObj == 1 ||
+ bound.rankObj.elIndex == 78 &&
+ bound.rankObj.ithObj == 2) {
+ AtomicBoundary& b = getNewPeriodic(otherElementRank);
+ b = bound;
+ }
+ } else {
+ AtomicBoundary& b = getNewPeriodic(otherElementRank);
+ b = bound;
+ }
}
}
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 644c52ed4d8429eb76cf2ab410fde711224c824a..aae499703c286b5ccc0f045a121f08226027aa67 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -68,6 +68,7 @@ namespace AMDiS {
return (*dof1 < *dof2);
}
+ bool MeshDistributor::sebastianMode = false;
MeshDistributor::MeshDistributor()
: problemStat(0),
@@ -1322,10 +1323,6 @@ namespace AMDiS {
// === Add new macro elements to mesh. ===
- TEST_EXIT_DBG(feSpaces.size() == 1)("Not yet implemented!\n");
- TEST_EXIT_DBG(feSpaces[0]->getBasisFcts()->getDegree() == 1)
- ("Not yet implemented!\n");
-
for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
it != newMacroEl.end(); ++it) {
MacroElement *mel = *it;
@@ -1336,10 +1333,7 @@ namespace AMDiS {
mel->setNeighbour(i, NULL);
// Create new DOFs for the macro element.
- mel->getElement()->createNewDofPtrs();
-
- for (int i = 0; i < mesh->getGeo(VERTEX); i++)
- mel->getElement()->setDof(i, mesh->getDof(VERTEX));
+ mel->getElement()->createNewDofPtrs(true);
// Push the macro element to all macro elements in mesh.
mesh->getMacroElements().push_back(mel);
diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h
index dff4bfb2bea2bfab44914fa1ba154092db922513..75c412dd9f5b476dbcf058dfcedf0c667e4f13b8 100644
--- a/AMDiS/src/parallel/MeshDistributor.h
+++ b/AMDiS/src/parallel/MeshDistributor.h
@@ -550,6 +550,8 @@ namespace AMDiS {
MeshLevelData levelData;
public:
+ static bool sebastianMode;
+
/// The boundary DOFs are sorted by subobject entities, i.e., first all
/// face DOFs, edge DOFs and to the last vertex DOFs will be set to
/// communication structure vectors, \ref sendDofs and \ref recvDofs.
diff --git a/AMDiS/src/parallel/ParallelDofMapping.cc b/AMDiS/src/parallel/ParallelDofMapping.cc
index bc306d6ccdd9a49347791e294dc337dad4a432ea..21d16e390a17527d36d8eed370ff47b1ad862df1 100644
--- a/AMDiS/src/parallel/ParallelDofMapping.cc
+++ b/AMDiS/src/parallel/ParallelDofMapping.cc
@@ -382,25 +382,74 @@ namespace AMDiS {
stdMpi.startCommunication();
- {
- for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpaces[i]);
- !it.end(); it.nextRank()) {
- int rank = it.getRank();
- if (meshLevel > 0)
- rank = levelData->mapRank(rank, 0, meshLevel);
-
- int counter = 0;
- for (; !it.endDofIter(); it.nextDof()) {
- if (dofMap.count(it.getDofIndex())) {
- DegreeOfFreedom d = stdMpi.getRecvData(rank)[counter++];
- if (globalIndex)
- dofToMatIndex.add(i, dofMap[it.getDofIndex()].global, d);
- else
- dofToMatIndex.add(i, it.getDofIndex(), d);
- }
+ for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpaces[i]);
+ !it.end(); it.nextRank()) {
+ int rank = it.getRank();
+ if (meshLevel > 0)
+ rank = levelData->mapRank(rank, 0, meshLevel);
+
+ int counter = 0;
+ for (; !it.endDofIter(); it.nextDof()) {
+ if (dofMap.count(it.getDofIndex())) {
+ DegreeOfFreedom d = stdMpi.getRecvData(rank)[counter++];
+ if (globalIndex)
+ dofToMatIndex.add(i, dofMap[it.getDofIndex()].global, d);
+ else
+ dofToMatIndex.add(i, it.getDofIndex(), d);
}
}
}
+
+
+ // === Communicate DOFs on periodic boundaries. ===
+
+ stdMpi.clear();
+
+ for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, feSpaces[i]);
+ !it.end(); it.nextRank()) {
+ vector<DegreeOfFreedom> sendGlobalDofs;
+
+ for (; !it.endDofIter(); it.nextDof())
+ if (dofMap.count(it.getDofIndex())) {
+ if (globalIndex) {
+ sendGlobalDofs.push_back(dofMap[it.getDofIndex()].global);
+ sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global));
+ } else {
+ sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex()));
+ }
+ }
+
+ int rank = it.getRank();
+ if (meshLevel > 0)
+ rank = levelData->mapRank(rank, 0, meshLevel);
+
+ stdMpi.send(rank, sendGlobalDofs);
+ stdMpi.recv(rank);
+ }
+
+ stdMpi.startCommunication();
+
+ for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, feSpaces[i]);
+ !it.end(); it.nextRank()) {
+
+ int rank = it.getRank();
+ if (meshLevel > 0)
+ rank = levelData->mapRank(rank, 0, meshLevel);
+
+ int counter = 0;
+
+ for (; !it.endDofIter(); it.nextDof())
+ if (dofMap.count(it.getDofIndex())) {
+ if (globalIndex) {
+ dofToMatIndex.add(i,
+ stdMpi.getRecvData(it.getRank())[counter++],
+ stdMpi.getRecvData(it.getRank())[counter++]);
+ } else {
+ dofToMatIndex.add(i, it.getDofIndex(),
+ stdMpi.getRecvData(it.getRank())[counter++]);
+ }
+ }
+ }
}
}
diff --git a/AMDiS/src/parallel/PetscSolver.h b/AMDiS/src/parallel/PetscSolver.h
index 022a50b39fc77fd765ce4dd9a2e21ec91d29d2fc..231a9203046edd6e2191ce0e80950cbaaa1aa540 100644
--- a/AMDiS/src/parallel/PetscSolver.h
+++ b/AMDiS/src/parallel/PetscSolver.h
@@ -125,8 +125,9 @@ namespace AMDiS {
DegreeOfFreedom dof)
{
FUNCNAME("PetscSolver::isCoarseSpace()");
- TEST_EXIT_DBG(coarseSpaceMap)
- ("Subdomain solver does not contain a coarse space!\n");
+
+ if (coarseSpaceMap == NULL)
+ return false;
return (*coarseSpaceMap)[feSpace].isSet(dof);
}
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 9d4b86ec332d669e487fb05a97491783ed8bd273..c7de58a86fd0e3281bb60b05036c2860b696cdd6 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -22,6 +22,12 @@ namespace AMDiS {
using namespace std;
+ double FetiTimings::fetiSolve = 0.0;
+ double FetiTimings::fetiSolve01 = 0.0;
+ double FetiTimings::fetiSolve02 = 0.0;
+ double FetiTimings::fetiPreconditioner = 0.0;
+
+
// y = mat * x
int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
{
@@ -44,25 +50,44 @@ namespace AMDiS {
// y = mat * x
int petscMultMatFeti(Mat mat, Vec x, Vec y)
{
+ FUNCNAME("petscMultMatFeti()");
+
// F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
// => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L)
+ double wtime = MPI::Wtime();
+
void *ctx;
MatShellGetContext(mat, &ctx);
FetiData* data = static_cast<FetiData*>(ctx);
MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
+
+ double wtime01 = MPI::Wtime();
data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
+
+ FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);
+
MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
+
+ wtime01 = MPI::Wtime();
KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
+ FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);
+
MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
+
+ wtime01 = MPI::Wtime();
data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
+ FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);
+
MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
+ FetiTimings::fetiSolve += (MPI::Wtime() - wtime);
+
return 0;
}
@@ -70,6 +95,8 @@ namespace AMDiS {
// y = PC * x
PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
{
+ double wtime = MPI::Wtime();
+
// Get data for the preconditioner
void *ctx;
PCShellGetContext(pc, &ctx);
@@ -125,6 +152,8 @@ namespace AMDiS {
// Multiply with scaled Lagrange constraint matrix.
MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);
+ FetiTimings::fetiPreconditioner += (MPI::Wtime() - wtime);
+
return 0;
}
@@ -1452,6 +1481,7 @@ namespace AMDiS {
MSG("FETI-DP timing 09: %.5f seconds (create rhs vector)\n",
MPI::Wtime() - wtime);
wtime = MPI::Wtime();
+ FetiTimings::reset();
}
@@ -1463,6 +1493,11 @@ namespace AMDiS {
MSG("FETI-DP timing 10: %.5f seconds (application of FETI-DP operator)\n",
MPI::Wtime() - wtime);
wtime = MPI::Wtime();
+
+ MSG("FETI-DP timing 10a: %.5f [ %.5f %.5f ] seconds (FETI-DP KSP Solve)\n",
+ FetiTimings::fetiSolve, FetiTimings::fetiSolve01, FetiTimings::fetiSolve02);
+ MSG("FETI-DP timing 10b: %.5f seconds (FETI-DP KSP Solve)\n",
+ FetiTimings::fetiPreconditioner);
}
// === Solve for u_primals. ===
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index 0a3cf2771876b14582ee641425306b37f8ee5551..cbe71f8713cfd5f2b24c4e38d5e05fce938b060b 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -269,6 +269,27 @@ namespace AMDiS {
bool printTimings;
};
+
+ class FetiTimings {
+ private:
+ FetiTimings() {}
+
+ public:
+ static void reset()
+ {
+ fetiSolve = 0.0;
+ fetiSolve01 = 0.0;
+ fetiSolve02 = 0.0;
+ fetiPreconditioner = 0.0;
+ }
+
+ public:
+ static double fetiSolve;
+ static double fetiSolve01;
+ static double fetiSolve02;
+
+ static double fetiPreconditioner;
+ };
}
#endif
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 78146dc3aac96ad5fc38fabd951584c1a1825a81..d7b3d3ac81603e6d190afc8a7ecd9564e6cfe893 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -102,7 +102,14 @@ namespace AMDiS {
MatAssemblyBegin(matIntInt, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(matIntInt, MAT_FINAL_ASSEMBLY);
+
+ // === Remove Dirichlet BC DOFs. ===
+
+ removeDirichletBcDofs(mat);
+
+
// === Init PETSc solver. ===
+
KSPCreate(mpiCommGlobal, &kspInterior);
KSPGetPC(kspInterior, &pcInterior);
KSPSetOperators(kspInterior, matIntInt, matIntInt, SAME_NONZERO_PATTERN);
@@ -296,6 +303,10 @@ namespace AMDiS {
}
+ // === Remove Dirichlet BC DOFs. ===
+
+ removeDirichletBcDofs(mat);
+
// === Create solver for the non primal (thus local) variables. ===
KSPCreate(mpiCommLocal, &kspInterior);
@@ -342,6 +353,10 @@ namespace AMDiS {
}
+ // === Remove Dirichlet BC DOFs. ===
+ removeDirichletBcDofs(vec);
+
+
VecAssemblyBegin(rhsInterior);
VecAssemblyEnd(rhsInterior);
@@ -350,6 +365,9 @@ namespace AMDiS {
VecAssemblyEnd(rhsCoarseSpace);
}
+
+ // === Remove null space, if requested. ===
+
if (removeRhsNullSpace) {
TEST_EXIT_DBG(coarseSpaceMap == NULL)("Not supported!\n");
@@ -409,6 +427,9 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");
+ double wtime = MPI::Wtime();
+ double t0 = 0.0, t1 = 0.0;
+
Vec tmp;
if (mpiCommLocal.Get_size() == 1)
VecCreateSeq(mpiCommLocal, interiorMap->getRankDofs(), &tmp);
@@ -428,7 +449,13 @@ namespace AMDiS {
VecRestoreArray(rhs, &rhsValues);
VecRestoreArray(tmp, &tmpValues);
+ t0 = MPI::Wtime() - wtime;
+
+ wtime = MPI::Wtime();
KSPSolve(kspInterior, tmp, tmp);
+ t1 = MPI::Wtime() - wtime;
+
+ wtime = MPI::Wtime();
VecGetArray(tmp, &tmpValues);
VecGetArray(sol, &rhsValues);
@@ -440,6 +467,9 @@ namespace AMDiS {
VecRestoreArray(tmp, &tmpValues);
VecDestroy(&tmp);
+ t0 += MPI::Wtime() - wtime;
+
+ MSG("TIMEING: %.5f %.5f\n", t0, t1);
}
@@ -712,7 +742,7 @@ namespace AMDiS {
if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
continue;
- if (coarseSpaceMap && isCoarseSpace(feSpace, dofIt.getDOFIndex())) {
+ if (isCoarseSpace(feSpace, dofIt.getDOFIndex())) {
TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");
int index = coarseSpaceMap->getMatIndex(nRowVec, dofIt.getDOFIndex());
@@ -752,6 +782,85 @@ namespace AMDiS {
}
+ void PetscSolverGlobalMatrix::removeDirichletBcDofs(Matrix<DOFMatrix*> *mat)
+ {
+ FUNCNAME("PetscSolverGlobalMatrix::removeDirichletBcDofs()");
+
+ vector<int> dofsInterior, dofsCoarse;
+
+ int nComponents = mat->getNumRows();
+ for (int i = 0; i < nComponents; i++) {
+ if ((*mat)[i][i]) {
+ const FiniteElemSpace *feSpace = (*mat)[i][i]->getRowFeSpace();
+
+ std::set<DegreeOfFreedom> &dirichletDofs = *((*mat)[i][i]->getApplyDBCs());
+
+ MSG("DIRICHLET DOFS: %d %d -> %d\n", i, i, dirichletDofs.size());
+
+ for (std::set<DegreeOfFreedom>::iterator it = dirichletDofs.begin();
+ it != dirichletDofs.end(); ++it) {
+ if (isCoarseSpace(feSpace, *it)) {
+ if ((*coarseSpaceMap)[feSpace].isRankDof(*it))
+ dofsCoarse.push_back(coarseSpaceMap->getMatIndex(i, *it));
+ } else {
+ if ((*interiorMap)[feSpace].isRankDof(*it))
+ dofsInterior.push_back(interiorMap->getMatIndex(i, *it));
+ }
+ }
+ } else {
+ MSG("NO MAT DIAG in %d\n", i);
+ }
+ }
+
+ MatZeroRows(matIntInt, dofsInterior.size(), &(dofsInterior[0]), 1.0,
+ PETSC_NULL, PETSC_NULL);
+
+ if (coarseSpaceMap != NULL)
+ MatZeroRows(matCoarseCoarse, dofsCoarse.size(), &(dofsCoarse[0]), 1.0,
+ PETSC_NULL, PETSC_NULL);
+ }
+
+
+ void PetscSolverGlobalMatrix::removeDirichletBcDofs(SystemVector *vec)
+ {
+ FUNCNAME("PetscSolverGlobalMatrix::removeDirichletBcDofs()");
+
+ int nComponents = vec->getSize();
+ for (int i = 0; i < nComponents; i++) {
+ const FiniteElemSpace *feSpace = vec->getDOFVector(i)->getFeSpace();
+
+ map<DegreeOfFreedom, double> &dirichletValues =
+ vec->getDOFVector(i)->getDirichletValues();
+
+ MSG("DIRICHLET DOFS: %d -> %d\n", i, dirichletValues.size());
+
+ MSG("MAT IS GLOBAL: %d\n", interiorMap->isMatIndexFromGlobal());
+// if (interiorMap->isMatIndexFromGlobal())
+// index =
+// interiorMap->getMatIndex(nRowVec, globalRowDof) + rStartInterior;
+// else
+// index =
+// interiorMap->getMatIndex(nRowVec, dofIt.getDOFIndex()) + rStartInterior;
+
+
+ for (map<DegreeOfFreedom, double>::iterator it = dirichletValues.begin();
+ it != dirichletValues.end(); ++it) {
+ if (isCoarseSpace(feSpace, it->first)) {
+ if ((*coarseSpaceMap)[feSpace].isRankDof(it->first))
+ VecSetValue(rhsCoarseSpace, coarseSpaceMap->getMatIndex(i, it->first),
+ it->second, INSERT_VALUES);
+ } else {
+ if ((*interiorMap)[feSpace].isRankDof(it->first)) {
+ MSG("REMOVE: %d %d %d\n", i, it->first, interiorMap->getMatIndex(i, it->first));
+ VecSetValue(rhsInterior, interiorMap->getMatIndex(i, it->first),
+ it->second, INSERT_VALUES);
+ }
+ }
+ }
+ }
+ }
+
+
void PetscSolverGlobalMatrix::createPetscNnzStructure(Matrix<DOFMatrix*> *mat)
{
FUNCNAME("PetscSolverGlobalMatrix::createPetscNnzStructure()");
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
index 94d2117b46bb9026a775ae2256275be80a3a9edf..ac2631300c9e61b104431099f53ed6448d8b2048 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
@@ -75,6 +75,10 @@ namespace AMDiS {
DOFVector<double>* vec,
int nRowVec, bool rankOnly = false);
+ void removeDirichletBcDofs(Matrix<DOFMatrix*> *mat);
+
+ void removeDirichletBcDofs(SystemVector *vec);
+
inline void setDofVector(Vec vecInterior,
DOFVector<double>* vec,
int nRowVec, bool rankOnly = false)
@@ -89,12 +93,10 @@ namespace AMDiS {
Vec petscSolVec;
- /** \brief
- * Stores the mesh change index of the mesh the nnz structure was created for.
- * Therefore, if the mesh change index is higher than this value, we have to create
- * a new nnz structure for PETSc matrices, because the mesh has been changed and
- * therefore also the assembled matrix structure.
- */
+ /// Stores the mesh change index of the mesh the nnz structure was created for.
+ /// Therefore, if the mesh change index is higher than this value, we have to create
+ /// a new nnz structure for PETSc matrices, because the mesh has been changed and
+ /// therefore also the assembled matrix structure.
int lastMeshNnz;
bool zeroStartVector;