diff --git a/AMDiS/src/parallel/ParallelDofMapping.cc b/AMDiS/src/parallel/ParallelDofMapping.cc index e62df99a1b64bd47bda94336f88135b4fc2f621d..fc651b0fcc16f689cd0daf29eaf5efdc738910d2 100644 --- a/AMDiS/src/parallel/ParallelDofMapping.cc +++ b/AMDiS/src/parallel/ParallelDofMapping.cc @@ -252,7 +252,8 @@ namespace AMDiS { rStartDofs = computeStartDofs(); // And finally, compute the matrix indices. - computeMatIndex(); + if (needMatIndex) + computeMatIndex(needMatIndexFromGlobal); } @@ -271,10 +272,10 @@ namespace AMDiS { } - void ParallelDofMapping::computeMatIndex() + void ParallelDofMapping::computeMatIndex(bool globalIndex) { FUNCNAME("ParallelDofMapping::computeMatIndex()"); - + dofToMatIndex.clear(); // The offset is always added to the local matrix index. The offset for the @@ -293,7 +294,10 @@ namespace AMDiS { for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) { if (data[feSpaces[i]].isRankDof(it->first)) { int globalMatIndex = it->second.local + offset; - dofToMatIndex.add(i, it->first, globalMatIndex); + if (globalIndex) + dofToMatIndex.add(i, it->second.global, globalMatIndex); + else + dofToMatIndex.add(i, it->first, globalMatIndex); } } @@ -318,7 +322,10 @@ namespace AMDiS { for (; !it.endDofIter(); it.nextDof()) if (dofMap.count(it.getDofIndex())) - sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex())); + if (globalIndex) + sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global)); + else + sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex())); stdMpi.send(it.getRank(), sendGlobalDofs); } @@ -336,7 +343,10 @@ namespace AMDiS { for (; !it.endDofIter(); it.nextDof()) { if (dofMap.count(it.getDofIndex())) { DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++]; - dofToMatIndex.add(i, it.getDofIndex(), d); + if (globalIndex) + dofToMatIndex.add(i, dofMap[it.getDofIndex()].global, d); + else + dofToMatIndex.add(i, it.getDofIndex(), d); } } } diff --git a/AMDiS/src/parallel/ParallelDofMapping.h b/AMDiS/src/parallel/ParallelDofMapping.h index 40d603fed034c9c90151591eb75a6716231afb25..2a1ef9e0f7d4152482d6255b813b1fed6d879ea2 100644 --- a/AMDiS/src/parallel/ParallelDofMapping.h +++ b/AMDiS/src/parallel/ParallelDofMapping.h @@ -260,6 +260,8 @@ namespace AMDiS { sendDofs(NULL), recvDofs(NULL), hasNonLocalDofs(false), + needMatIndex(false), + needMatIndexFromGlobal(false), feSpaces(0), nRankDofs(-1), nLocalDofs(-1), @@ -374,6 +376,15 @@ namespace AMDiS { ERROR_EXIT("MUST BE IMPLEMENTED!\n"); } + void setComputeMatIndex(bool b, bool global = false) + { + needMatIndex = b; + needMatIndexFromGlobal = global; + } + + /// Compute local and global matrix indices. + void computeMatIndex(bool globalIndex); + protected: /// Insert a new FE space DOF mapping for a given FE space. void addFeSpace(const FiniteElemSpace* feSpace); @@ -390,9 +401,6 @@ namespace AMDiS { /// Compute \ref rStartDofs. int computeStartDofs(); - /// Compute local and global matrix indices. - void computeMatIndex(); - private: /// MPI communicator object; MPI::Intracomm mpiComm; @@ -405,6 +413,16 @@ namespace AMDiS { /// are also owned by this rank. bool hasNonLocalDofs; + /// If true, matrix indeces for the stored DOFs are computed, see + /// \ref computeMatIndex. + bool needMatIndex; + + /// If matrix indices should be computed, this variable defines if the + /// mapping from DOF indices to matrix row indices is defined on local + /// or global DOF indices. If true, the mapping is to specify and to use + /// on global ones, otherwise on local DOF indices. + bool needMatIndexFromGlobal; + /// Maps from FE space pointers to DOF mappings. map data; diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc index a184448ace289541534d83567ea56e3ce6d14e8e..8affcba7bb80d40e56edcc50ccf8777769935483 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.cc +++ b/AMDiS/src/parallel/PetscSolverFeti.cc @@ -244,15 +244,25 @@ namespace AMDiS { primalDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(), true, true); + primalDofMap.setComputeMatIndex(true); + dualDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(), false, false); + dualDofMap.setComputeMatIndex(true); + lagrangeMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(), true, true); + lagrangeMap.setComputeMatIndex(true); + localDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(), false, false); - if (fetiPreconditioner == FETI_DIRICHLET) + localDofMap.setComputeMatIndex(true); + + if (fetiPreconditioner == FETI_DIRICHLET) { interiorDofMap.init(mpiComm, feSpaces, meshDistributor->getFeSpaces(), false, false); + interiorDofMap.setComputeMatIndex(true); + } } diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc index 3c43669edfae09554350f6b2494145f0d95f0151..c2a53b8d9ace865f474a5586f8e85ced02601863 100644 --- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc +++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc @@ -26,37 +26,20 @@ namespace AMDiS { TEST_EXIT_DBG(mat)("No DOF matrix defined!\n"); double wtime = MPI::Wtime(); - vector feSpaces = getFeSpaces(mat); - dofMap->update(feSpaces); - - int nRankRows = dofMap->getRankDofs(); - int nOverallRows = dofMap->getOverallDofs(); - - // === Create PETSc vector (solution and a temporary vector). === - - VecCreateMPI(mpiComm, nRankRows, nOverallRows, &petscSolVec); - VecCreateMPI(mpiComm, nRankRows, nOverallRows, &petscTmpVec); - - int testddd = 1; - Parameters::get("block size", testddd); - if (testddd > 1) { - VecSetBlockSize(petscSolVec, testddd); - VecSetBlockSize(petscTmpVec, testddd); - } + + // === Check if mesh was changed and, in this case, recompute matrix === + // === nnz structure and matrix indices. === int recvAllValues = 0; int sendValue = static_cast(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz); mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM); - recvAllValues = 1; - if (!d_nnz || recvAllValues != 0 || alwaysCreateNnzStructure) { - - // Global DOF mapping must only be recomputed, if the NNZ structure has - // changed (thus, also the mesh has changed). - createGlobalDofMapping(feSpaces); + vector feSpaces = getFeSpaces(mat); + dofMap->setComputeMatIndex(true, true); + dofMap->update(feSpaces); if (d_nnz) { delete [] d_nnz; @@ -70,6 +53,22 @@ namespace AMDiS { } + // === Create PETSc vector (solution and a temporary vector). === + + int nRankRows = dofMap->getRankDofs(); + int nOverallRows = dofMap->getOverallDofs(); + + VecCreateMPI(mpiComm, nRankRows, nOverallRows, &petscSolVec); + VecCreateMPI(mpiComm, nRankRows, nOverallRows, &petscTmpVec); + + int testddd = 1; + Parameters::get("block size", testddd); + if (testddd > 1) { + VecSetBlockSize(petscSolVec, testddd); + VecSetBlockSize(petscTmpVec, testddd); + } + + // === Create PETSc matrix with the computed nnz data structure. === MatCreateMPIAIJ(mpiComm, nRankRows, nRankRows, @@ -134,7 +133,6 @@ namespace AMDiS { TEST_EXIT_DBG(vec)("No DOF vector defined!\n"); TEST_EXIT_DBG(dofMap)("No parallel DOF map defined!\n"); - vector feSpaces = getFeSpaces(vec); int nRankRows = dofMap->getRankDofs(); int nOverallRows = dofMap->getOverallDofs(); @@ -279,7 +277,7 @@ namespace AMDiS { // === Row DOF index is not periodic. === // Get PETSc's mat row index. - int rowIndex = dofToMatIndex.get(nRowMat, globalRowDof); + int rowIndex = dofMap->getMatIndex(nRowMat, globalRowDof); cols.clear(); values.clear(); @@ -292,7 +290,7 @@ namespace AMDiS { // Test if the current col dof is a periodic dof. bool periodicCol = perMap.isPeriodic(colFe, globalColDof); // Get PETSc's mat col index. - int colIndex = dofToMatIndex.get(nColMat, globalColDof); + int colIndex = dofMap->getMatIndex(nColMat, globalColDof); // Ignore all zero entries, expect it is a diagonal entry. if (value(*icursor) == 0.0 && rowIndex != colIndex) @@ -341,7 +339,7 @@ namespace AMDiS { } for (unsigned int i = 0; i < newCols.size(); i++) { - cols.push_back(dofToMatIndex.get(nColMat, newCols[i])); + cols.push_back(dofMap->getMatIndex(nColMat, newCols[i])); values.push_back(scaledValue); } } @@ -428,8 +426,8 @@ namespace AMDiS { // === Translate the matrix entries to PETSc's matrix. for (unsigned int i = 0; i < entry.size(); i++) { - int rowIdx = dofToMatIndex.get(nRowMat, entry[i].first); - int colIdx = dofToMatIndex.get(nColMat, entry[i].second); + int rowIdx = dofMap->getMatIndex(nRowMat, entry[i].first); + int colIdx = dofMap->getMatIndex(nColMat, entry[i].second); colsMap[rowIdx].push_back(colIdx); valsMap[rowIdx].push_back(scaledValue); @@ -474,7 +472,7 @@ namespace AMDiS { DegreeOfFreedom globalRowDof = (*dofMap)[feSpace][dofIt.getDOFIndex()].global; // Get PETSc's mat index of the row DOF. - int index = dofToMatIndex.get(nRowVec, globalRowDof); + int index = dofMap->getMatIndex(nRowVec, globalRowDof); if (perMap.isPeriodic(feSpace, globalRowDof)) { std::set& perAsc = perMap.getAssociations(feSpace, globalRowDof); @@ -484,7 +482,7 @@ namespace AMDiS { for (std::set::iterator perIt = perAsc.begin(); perIt != perAsc.end(); ++perIt) { int mappedDof = perMap.map(feSpace, *perIt, globalRowDof); - int mappedIndex = dofToMatIndex.get(nRowVec, mappedDof); + int mappedIndex = dofMap->getMatIndex(nRowVec, mappedDof); VecSetValues(petscVec, 1, &mappedIndex, &value, ADD_VALUES); } } else { @@ -577,7 +575,7 @@ namespace AMDiS { int globalRowDof = (*dofMap)[feSpaces[i]][*cursor].global; // The corresponding global matrix row index of the current row DOF. - int petscRowIdx = dofToMatIndex.get(i, globalRowDof); + int petscRowIdx = dofMap->getMatIndex(i, globalRowDof); if ((*dofMap)[feSpaces[i]].isRankDof(*cursor)) { // === The current row DOF is a rank DOF, so create the === @@ -601,7 +599,7 @@ namespace AMDiS { for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { int globalColDof = (*dofMap)[feSpaces[j]][col(*icursor)].global; - int petscColIdx = dofToMatIndex.get(j, globalColDof); + int petscColIdx = dofMap->getMatIndex(j, globalColDof); if (value(*icursor) != 0.0 || petscRowIdx == petscColIdx) { // The row DOF is a rank DOF, if also the column is a rank DOF, @@ -629,7 +627,7 @@ namespace AMDiS { if (value(*icursor) != 0.0) { int globalColDof = (*dofMap)[feSpaces[j]][col(*icursor)].global; - int petscColIdx = dofToMatIndex.get(j, globalColDof); + int petscColIdx = dofMap->getMatIndex(j, globalColDof); sendMatrixEntry[sendToRank]. push_back(make_pair(petscRowIdx, petscColIdx)); @@ -686,99 +684,4 @@ namespace AMDiS { d_nnz[i] = std::min(d_nnz[i], nRankRows); } - - void PetscSolverGlobalMatrix::createGlobalDofMapping(vector &feSpaces) - { - FUNCNAME("PetscSolverGlobalMatrix::createGlobalDofMapping()"); - - int offset = dofMap->getStartDofs(); - Mesh *mesh = feSpaces[0]->getMesh(); - - dofToMatIndex.clear(); - - for (unsigned int i = 0; i < feSpaces.size(); i++) { - - // === Create indices for all DOFs in rank' domain. === - std::set rankDofSet; - mesh->getAllDofs(feSpaces[i], rankDofSet); - for (std::set::iterator it = rankDofSet.begin(); - it != rankDofSet.end(); ++it) - if ((*dofMap)[feSpaces[i]].isRankDof(**it)) { - int globalIndex = (*dofMap)[feSpaces[i]][**it].global; - - int globalMatIndex = - globalIndex - (*dofMap)[feSpaces[i]].rStartDofs + offset; - - dofToMatIndex.add(i, globalIndex, globalMatIndex); - } - - // === Communicate interior boundary DOFs between domains. === - - StdMpi > stdMpi(mpiComm); - - for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpaces[i]); - !it.end(); it.nextRank()) { - vector sendGlobalDofs; - - for (; !it.endDofIter(); it.nextDof()) { - int globalIndex = (*dofMap)[feSpaces[i]][it.getDofIndex()].global; - int globalMatIndex = dofToMatIndex.get(i, globalIndex); - sendGlobalDofs.push_back(globalMatIndex); - } - - stdMpi.send(it.getRank(), sendGlobalDofs); - } - - for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpaces[i]); - !it.end(); it.nextRank()) - stdMpi.recv(it.getRank()); - - stdMpi.startCommunication(); - - for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpaces[i]); - !it.end(); it.nextRank()) - for (; !it.endDofIter(); it.nextDof()) { - int globalIndex = (*dofMap)[feSpaces[i]][it.getDofIndex()].global; - int globalMatIndex = - stdMpi.getRecvData(it.getRank())[it.getDofCounter()]; - - dofToMatIndex.add(i, globalIndex, globalMatIndex); - } - - - // === Communicate DOFs on periodic boundaries. === - - stdMpi.clear(); - - for (DofComm::Iterator it(meshDistributor->getPeriodicDofs(), feSpaces[i]); - !it.end(); it.nextRank()) { - vector sendGlobalDofs; - - for (; !it.endDofIter(); it.nextDof()) { - int ind0 = (*dofMap)[feSpaces[i]][it.getDofIndex()].global; - int ind1 = dofToMatIndex.get(i, ind0); - - sendGlobalDofs.push_back(ind0); - sendGlobalDofs.push_back(ind1); - } - - stdMpi.send(it.getRank(), sendGlobalDofs); - stdMpi.recv(it.getRank()); - } - - stdMpi.startCommunication(); - - for (DofComm::Iterator it(meshDistributor->getPeriodicDofs(), feSpaces[i]); - !it.end(); it.nextRank()) - for (; !it.endDofIter(); it.nextDof()) { - int ind0 = stdMpi.getRecvData(it.getRank())[it.getDofCounter() * 2]; - int ind1 = stdMpi.getRecvData(it.getRank())[it.getDofCounter() * 2 + 1]; - - dofToMatIndex.add(i, ind0, ind1); - } - - // === Update offset. === - offset += (*dofMap)[feSpaces[i]].nRankDofs; - } - } } diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h index 87e3913b2466744598eacf1aea5e428e0359baa6..b95f68dd02966052c041fd7f14234c1a76aa9f6b 100644 --- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h +++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h @@ -68,8 +68,6 @@ namespace AMDiS { void setDofVector(Vec& petscVec, DOFVector* vec, int nRowVec, bool rankOnly = false); - void createGlobalDofMapping(vector &feSpaces); - protected: /// Arrays definig the non zero pattern of Petsc's matrix. int *d_nnz, *o_nnz; @@ -91,10 +89,6 @@ namespace AMDiS { /// operators using DOFVectors from old timestep containing many zeros due to /// some phase fields. bool alwaysCreateNnzStructure; - - /// Mapping from global DOF indices to global matrix indices under - /// consideration of possibly multiple components. - DofToMatIndex dofToMatIndex; };