Trilinos
Problem with sumIntoGlobalValues when assembling a Tepetra::FECrsMatrix
#13943
-
|
I'm transitioning some Epetra code to Tpetra. During the global assembly element loop, ghost DOFs are not summed into the matrix. I followed the example in I set up a small example with 3 elements, with rank 1 owning dof 4 which is also part of element 2 owned by rank 0 (it is a ghost in rank 0). During the the assembly, the arguments of But the final matrix, after It did not sum up the contribution from rank 0? I checked the graph, created with, graph_ =
Teuchos::rcp(new Tpetra::FECrsGraph<LO, GO, NO>(rowMap_, owned_plus_shared_map, maxNumEntriesPerRow));The row map contains only owned dofs (domain and range are the same). The column map has owned plus shared: Any ideas on how to debug this? @trilinos/tpetra |
Beta Was this translation helpful? Give feedback.
Answered by
nasseralkmim
Apr 7, 2025
Replies: 2 comments 4 replies
-
Can you show the full code for this small example? |
Beta Was this translation helpful? Give feedback.
1 reply
-
|
For this example, with each node with 4 DOFS for 3 fields. I assume that node 2 is a ghost for rank 0, since it owns element 2 but the node is in rank 1. The owned map is created here: std::vector<int> myGlobalElements;
for(const auto& entry : DofToAdjacentDofs) {
myGlobalElements.push_back(entry.first);
}
Teuchos::ArrayView<const GO> myGlobalElementsView(myGlobalElements);
const Tpetra::global_size_t numGlobalElements = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
const GO indexBase = 0;
rowMap_ = Teuchos::rcp(new Tpetra::Map<LO, GO, NO>(numGlobalElements, myGlobalElementsView, indexBase, comm));Ghost dofs here, std::set<GO> myGhostGlobalElementsSet; // Use set for automatic uniqueness
int nodekey = -1;
Node* theNode = NULL;
theGrid.ResetNodeLoopToGhosted(); // Start iterating only ghost nodes
while (theGrid.NextNode(nodekey, theNode, false) !=
0) // include_ghosted=false is fine here as we started with ghost iterator
{
EntityDofs ghostDofEntity = theNumberer.GetNodeDofs(theNode->GetNodeId());
std::vector<int> ghost_available_dofs;
std::vector<int> ghost_gids_int; // Use int first as GetDofs returns int
ghostDofEntity.GetAvailableDofTypes(ghost_available_dofs);
ghostDofEntity.GetDofs(ghost_available_dofs, ghost_gids_int);
for (int gid_int : ghost_gids_int) {
myGhostGlobalElementsSet.insert(static_cast<GO>(gid_int));
}
}This gives as output: I assume this is correct. I combine the owned plus shared here: // Combine owned and ghost DOF GIDs into a map
std::vector<GO> myAllGlobalElements;
myAllGlobalElements.reserve(myGlobalElements.size() + myGhostGlobalElementsSet.size());
myAllGlobalElements.insert(myAllGlobalElements.end(), myGlobalElements.begin(), myGlobalElements.end());
myAllGlobalElements.insert(myAllGlobalElements.end(), myGhostGlobalElementsSet.begin(), myGhostGlobalElementsSet.end());
Teuchos::ArrayView<const GO> myAllGlobalElementsView(myAllGlobalElements);
auto owned_plus_shared_map = Teuchos::rcp(new Tpetra::Map<LO, GO, NO>(numGlobalElements, myAllGlobalElementsView, indexBase, comm));I find the max number of entries per row with this: GO maxNumEntriesPerRow = 0;
for (const auto& entry : DofToAdjacentDofs) {
// entry.first is the global row number (dof global index)
// entry.second is the set of adjacent dofs (column indices)
const auto& adjacentDofs = entry.second;
if (adjacentDofs.size() + 1 > maxNumEntriesPerRow) {
maxNumEntriesPerRow = adjacentDofs.size();
}
}Which results in 24 (2 elements sharing one node, results in 16 dofs from one element, 8 from the other). Then I create the graph: graph_ =
Teuchos::rcp(new Tpetra::FECrsGraph<LO, GO, NO>(rowMap_, owned_plus_shared_map, maxNumEntriesPerRow));
// Populate the graph_
graph_->beginAssembly();
for(const auto& entry : DofToAdjacentDofs) {
GO globalRow = entry.first;
const std::set<int>& adjacentDofs = entry.second;
// Create a view of the column indices for Tpetra
// Include the diagonal entry itself
std::vector<GO> indicesVec;
indicesVec.reserve(adjacentDofs.size() + 1);
indicesVec.push_back(globalRow); // Add diagonal
for (int col : adjacentDofs) {
indicesVec.push_back(static_cast<GO>(col));
}
Teuchos::ArrayView<const GO> indicesView(indicesVec);
graph_->insertGlobalIndices(globalRow, indicesView);
}
graph_->endAssembly();Which results in this row and column maps. Domain and range are the same as row map. The matrix is created with the graph ( A_ = Teuchos::rcp(new CrsMatrix(graph_));
A_->setAllToScalar(0.0);The assembly is done like this: A_->beginAssembly();
// Loop over OWNED elements (rank 0: 2 and rank 1: 0 and 1)
while(theGrid.NextElement(elekey,theEle,include_ghosted)!=0) {
// dofN: number of dofs
// dofIi: converts local to global id
// elemat: element matrix
for(int i=0;i<dofN;++i)
{
GO global_row = dofIi[i];
for(int j=0;j<dofN;++j)
{
GO global_col = dofIi[j];
SC value = elemat[i][j];
A_->sumIntoGlobalValues(global_row, 1, &value, &global_col);
}
}
}
A_->endAssembly();I output the values in this loop and I confirmed that they are correct. I can give more details if you think it is needed. |
Beta Was this translation helpful? Give feedback.
-
Unless you've used a I'm also not sure what to make of this: std::vector<int> myGlobalElements;
for(const auto& entry : DofToAdjacentDofs) {
myGlobalElements.push_back(entry.first);
}Since I don't know what I would strongly recommend creating a self-contained minimal reproducible example. Use what would normally be an absurdly simplistic representation of element connectivity (say, a vector of vectors of |
Beta Was this translation helpful? Give feedback.
3 replies
-
|
@jjramsey Thanks for taking a look into it. I set up a minimal example with the same data, and it works as expected. Now I'm more puzzled. Thanks! |
Beta Was this translation helpful? Give feedback.
-
|
I found one difference between the MRE and my app. Whereas in the MRE: Do you think that could be the reason, @jjramsey? |
Beta Was this translation helpful? Give feedback.
-
|
I found the culprit. |
Beta Was this translation helpful? Give feedback.
Answer selected by
nasseralkmim
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I found the culprit.
I was creating the graph based on owned dofs only.
The solution was to loop over owned elements and include ghost dofs in the graph.