fix: Fracture/3D cell co-location in parallel mesh redistribution by bd713 · Pull Request #4008 · GEOS-DEV/GEOS

bd713 · 2026-03-30T04:05:44Z

Replaces #3985

codecov · 2026-05-05T21:28:48Z

Codecov Report

❌ Patch coverage is 50.19405% with 385 lines in your changes missing coverage. Please review.
⚠️ Please upload report for BASE (develop@506fae7). Learn more about missing BASE report.

Files with missing lines	Patch %	Lines
...nents/mesh/generators/VTKSuperCellPartitioning.cpp	39.38%	257 Missing ⚠️
...rc/coreComponents/mesh/generators/VTKUtilities.cpp	63.30%	102 Missing ⚠️
...reComponents/mesh/generators/ParMETISInterface.cpp	0.00%	18 Missing ⚠️
...coreComponents/mesh/generators/CollocatedNodes.cpp	38.46%	8 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             develop    #4008   +/-   ##
==========================================
  Coverage           ?   64.50%           
==========================================
  Files              ?     1384           
  Lines              ?   123603           
  Branches           ?        0           
==========================================
  Hits               ?    79729           
  Misses             ?    43874           
  Partials           ?        0

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

jafranc

Great ! Thanks, some shortening suggestion here and there

jafranc · 2026-05-06T06:49:33Z

      for( vtkIdType const * q: quad )
      {
-        fracture->InsertNextCell( VTK_QUAD, numPoints, q );
+        fracture->InsertNextCell( VTK_QUAD, 4, q );


isn't it weird to hardcode there and still use numPoints above ?

Made more explicit

jafranc · 2026-05-06T07:30:19Z

+  }
+
+  array1d< real_t > tpwgts( numParts );
+  tpwgts.setValues< serialPolicy >( 1.0f / static_cast< real_t >( numParts ) );


isn't parmetis dividing in equal weights if NULL is pass as tpwgts ? If so, it might avoid a check for zero-div here

Tried this, but it makes the unit tests crash.
But will protect the division.

jafranc · 2026-05-06T07:57:04Z

 CollocatedNodes::CollocatedNodes( string const & faceBlockName,
-                                  vtkSmartPointer< vtkDataSet > faceMesh )
+                                  vtkSmartPointer< vtkDataSet > faceMesh,
+                                  bool performGlobalCheck )
 {
  // The vtk field to the collocated nodes for fractures.
  string const COLLOCATED_NODES = "collocated_nodes";

  vtkIdTypeArray const * collocatedNodes = vtkIdTypeArray::FastDownCast( faceMesh->GetPointData()->GetArray( COLLOCATED_NODES.c_str() ) );

-  // Depending on the parallel split, the vtk face mesh may be empty on a rank.
-  // In that case, vtk will not provide any field for the emtpy mesh.
-  // Therefore, not finding the duplicated nodes field on a rank cannot be interpreted as a globally missing field.
-  // Converting the address into an integer and exchanging it through the MPI ranks let us find out
-  // if the field is globally missing or not.
-  std::uintptr_t const address = MpiWrapper::max( reinterpret_cast< std::uintptr_t >(collocatedNodes) );
-  if( address == 0 )
+  if( performGlobalCheck )
  {
-    GEOS_LOG_RANK_0( "Available point data fields in '" << faceBlockName << "':" );
-    for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )
+    // Depending on the parallel split, the vtk face mesh may be empty on a rank.
+    // In that case, vtk will not provide any field for the empty mesh.
+    // Therefore, not finding the duplicated nodes field on a rank cannot be interpreted as a globally missing field.
+    // Converting the address into an integer and exchanging it through the MPI ranks let us find out
+    // if the field is globally missing or not.
+    std::uintptr_t const address = MpiWrapper::max( reinterpret_cast< std::uintptr_t >(collocatedNodes) );
+    if( address == 0 )
    {
-      GEOS_LOG_RANK_0( " - " << faceMesh->GetPointData()->GetArrayName( i ) << " of type '" << faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() << "'" );
+      GEOS_LOG_RANK_0( GEOS_FMT( "Available point data fields in '{}':", faceBlockName ) );
+      for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )
+      {
+        GEOS_LOG_RANK_0( GEOS_FMT( " - {} of type '{}'",
+                                   faceMesh->GetPointData()->GetArrayName( i ),
+                                   faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() ) );
+      }
+      GEOS_ERROR( GEOS_FMT( "Could not find valid field \"{}\" for fracture \"{}\".",
+                            COLLOCATED_NODES,
+                            faceBlockName ) );
+    }
+  }
+  else
+  {
+    if( !collocatedNodes )
+    {
+      GEOS_LOG_RANK_0( GEOS_FMT( "Available point data fields in '{}':", faceBlockName ) );
+      for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )
+      {
+        GEOS_LOG_RANK_0( GEOS_FMT( " - {} of type '{}'",
+                                   faceMesh->GetPointData()->GetArrayName( i ),
+                                   faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() ) );
+      }
+      GEOS_ERROR( GEOS_FMT( "Could not find valid field \"{}\" for fracture \"{}\" on this rank.",
+                            COLLOCATED_NODES,
+                            faceBlockName ) );


Suggested change

CollocatedNodes::CollocatedNodes( string const & faceBlockName,

vtkSmartPointer< vtkDataSet > faceMesh )

vtkSmartPointer< vtkDataSet > faceMesh,

bool performGlobalCheck )

{

// The vtk field to the collocated nodes for fractures.

string const COLLOCATED_NODES = "collocated_nodes";

vtkIdTypeArray const * collocatedNodes = vtkIdTypeArray::FastDownCast( faceMesh->GetPointData()->GetArray( COLLOCATED_NODES.c_str() ) );

// Depending on the parallel split, the vtk face mesh may be empty on a rank.

// In that case, vtk will not provide any field for the emtpy mesh.

// Therefore, not finding the duplicated nodes field on a rank cannot be interpreted as a globally missing field.

// Converting the address into an integer and exchanging it through the MPI ranks let us find out

// if the field is globally missing or not.

std::uintptr_t const address = MpiWrapper::max( reinterpret_cast< std::uintptr_t >(collocatedNodes) );

if( address == 0 )

if( performGlobalCheck )

{

GEOS_LOG_RANK_0( "Available point data fields in '" << faceBlockName << "':" );

for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )

// Depending on the parallel split, the vtk face mesh may be empty on a rank.

// In that case, vtk will not provide any field for the empty mesh.

// Therefore, not finding the duplicated nodes field on a rank cannot be interpreted as a globally missing field.

// Converting the address into an integer and exchanging it through the MPI ranks let us find out

// if the field is globally missing or not.

std::uintptr_t const address = MpiWrapper::max( reinterpret_cast< std::uintptr_t >(collocatedNodes) );

if( address == 0 )

{

GEOS_LOG_RANK_0( " - " << faceMesh->GetPointData()->GetArrayName( i ) << " of type '" << faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() << "'" );

GEOS_LOG_RANK_0( GEOS_FMT( "Available point data fields in '{}':", faceBlockName ) );

for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )

{

GEOS_LOG_RANK_0( GEOS_FMT( " - {} of type '{}'",

faceMesh->GetPointData()->GetArrayName( i ),

faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() ) );

}

GEOS_ERROR( GEOS_FMT( "Could not find valid field \"{}\" for fracture \"{}\".",

COLLOCATED_NODES,

faceBlockName ) );

}

}

else

{

if( !collocatedNodes )

{

GEOS_LOG_RANK_0( GEOS_FMT( "Available point data fields in '{}':", faceBlockName ) );

for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )

{

GEOS_LOG_RANK_0( GEOS_FMT( " - {} of type '{}'",

faceMesh->GetPointData()->GetArrayName( i ),

faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() ) );

}

GEOS_ERROR( GEOS_FMT( "Could not find valid field \"{}\" for fracture \"{}\" on this rank.",

COLLOCATED_NODES,

faceBlockName ) );

CollocatedNodes::CollocatedNodes( string const & faceBlockName,

vtkSmartPointer< vtkDataSet > faceMesh,

bool performGlobalCheck )

{

string const COLLOCATED_NODES = "collocated_nodes";

vtkIdTypeArray const * collocatedNodes = vtkIdTypeArray::FastDownCast(

faceMesh->GetPointData()->GetArray( COLLOCATED_NODES.c_str() ) );

// Determine whether the field is missing:

// - if globalCheck: reduce across ranks so an empty-mesh detection is any ranks

bool const isMissing = performGlobalCheck

? ( MpiWrapper::max( reinterpret_cast< std::uintptr_t >( collocatedNodes ) ) == 0 )

: ( collocatedNodes == nullptr );

if( isMissing )

{

GEOS_LOG_RANK_0( GEOS_FMT( "Available point data fields in '{}':", faceBlockName ) );

for( int i = 0; i < faceMesh->GetPointData()->GetNumberOfArrays(); ++i )

{

GEOS_LOG_RANK_0( GEOS_FMT( " - {} of type '{}'",

faceMesh->GetPointData()->GetArrayName( i ),

faceMesh->GetPointData()->GetArray( i )->GetDataTypeAsString() ) );

}

GEOS_ERROR( GEOS_FMT( "Could not find valid field \"{}\" for fracture \"{}\"{}.",

COLLOCATED_NODES,

faceBlockName,

performGlobalCheck ? "" : " on this rank" ) );

}

else

init( collocatedNodes );

}

jafranc · 2026-05-06T08:12:13Z

+  // Options: [use_defaults, log_level, seed, coupling]
+  // PARMETIS_PSR_UNCOUPLED = 0 (default - uses PartitionSmallGraph)
+  // PARMETIS_PSR_COUPLED = 1 (forces distributed algorithm)
+  idx_t options[4] = { 1, 0, 2022, 0 };


Suggested change

idx_t options[4] = { 1, 0, 2022, 0 };

idx_t options[4] = { 1, 0, 2022, PARMETIS_PSR_UNCOUPLED };

jafranc · 2026-05-06T09:08:49Z


  // GetCellAtId() is conditionally thread-safe, use POLICY argument
-  forAll< POLICY >( num3dCells, [&cells, &globalPointId, elemToNodes = elemToNodes.toView()] ( localIndex const cellIdx )
+  forAll< POLICY >( numCells, [&cells, &globalPointId, elemToNodes = elemToNodes.toView()] ( localIndex const cellIdx )


why aren't we using GetCellPoints() here ?
https://vtk.org/doc/nightly/html/classvtkDataSet.html#ab4d677c257a58e4eb4a80757f9b371ea

This might the need for passing cells through the interface and the IsShareable logic

yes, that part was refactored.
@jafranc, could you review the changes and confirm that they correctly addressed your comments.

jafranc · 2026-05-06T12:24:31Z

+  int myEdgeCount = static_cast< int >( myEdgesSrc.size() );
+
+  // Gather counts from all ranks
+  array1d< int > allEdgeCounts( numRanks );
+  MpiWrapper::allgather( &myEdgeCount, 1, allEdgeCounts.data(), 1, comm );
+
+  // Compute displacements for symmetrization
+  array1d< int > edgeDispls( numRanks + 1 );
+  edgeDispls[0] = 0;
+  for( int r = 0; r < numRanks; ++r )
+  {
+    edgeDispls[r + 1] = edgeDispls[r] + allEdgeCounts[r];
+  }


Suggested change

int myEdgeCount = static_cast< int >( myEdgesSrc.size() );

// Gather counts from all ranks

array1d< int > allEdgeCounts( numRanks );

MpiWrapper::allgather( &myEdgeCount, 1, allEdgeCounts.data(), 1, comm );

// Compute displacements for symmetrization

array1d< int > edgeDispls( numRanks + 1 );

edgeDispls[0] = 0;

for( int r = 0; r < numRanks; ++r )

{

edgeDispls[r + 1] = edgeDispls[r] + allEdgeCounts[r];

}

int myEdgeCount = static_cast< int >( myEdgesSrc.size() );

// Gather counts from all ranks

array1d< int > allEdgeCounts( numRanks );

MpiWrapper::allgather( &myEdgeCount, 1, allEdgeCounts.data(), 1, comm );

// Compute displacements for symmetrization

array1d< int > edgeDispls( numRanks + 1 );

edgeDispls[0] = 0;

std::partial_sum(allEdgeCounts.begin(), allEdgeCounts.end(),

edgeDispls.begin() + 1);

jafranc · 2026-05-06T12:47:51Z

+  // Check 6: Graph symmetry for local edges
+  // -----------------------------------------------------------------------
+  {
+    std::unordered_set< uint64_t > localOutgoingNeighbors;


can we refactor this duplication into a single locally declared function ?

Duplication removed

jafranc · 2026-05-06T12:49:35Z

+// SECTION 5: PARTITION UNPACKING (after ParMETIS)
+// =============================================================================================
+array1d< int64_t >
+unpackSuperCellPartitioning(


as packing/unpacking is rather specific in GEOS context, unless related, I would suggest renaming

renamed expandSuperCellPartitioningToCells

jafranc · 2026-05-06T12:52:18Z

+  // -----------------------------------------------------------------------
+  // Step 1: Build 3D cell connectivity graph via fractures
+  // -----------------------------------------------------------------------
+  std::map< vtkIdType, std::set< vtkIdType > > fractureGraph;


isn't the use of stdMap rather than std::map a directive ?

yes, changed all of them in favor of StdContainerWrappers

jafranc · 2026-05-06T12:59:51Z

+  vtkIdType const numCells = cells3D->GetNumberOfCells();
+  array1d< int64_t > cellPartitioning( numCells );
+
+  for( vtkIdType i = 0; i < numCells; ++i )
+  {
+    vtkIdType superCellId = superCellIdArray->GetValue( i );
+    auto it = superCellIdToLocalIdx.find( superCellId );
+
+    GEOS_ERROR_IF( it == superCellIdToLocalIdx.end(),
+                   GEOS_FMT( "Rank {}: Cell {} has unknown super-cell ID {}",
+                             rank, i, superCellId ) );
+
+    cellPartitioning[i] = superPartitioning[it->second];
+  }
+
+  // -----------------------------------------------------------------------
+  // Step 3: Validate that super-cells weren't split
+  // -----------------------------------------------------------------------
+  stdMap< vtkIdType, std::set< int64_t > > superCellToRanks;
+
+  for( vtkIdType i = 0; i < numCells; ++i )
+  {
+    vtkIdType const scId = superCellIdArray->GetValue( i );
+    superCellToRanks.get_inserted( scId ).insert( cellPartitioning[i] );
+  }
+
+  vtkIdType numSplitSuperCells = 0;
+  for( auto const & [superCellId, ranks] : superCellToRanks )
+  {
+    if( ranks.size() > 1 )
+    {
+      ++numSplitSuperCells;
+    }
+  }
+
+  vtkIdType totalSplitSuperCells = MpiWrapper::sum( numSplitSuperCells, comm );
+
+  GEOS_ERROR_IF( totalSplitSuperCells > 0,
+                 GEOS_FMT( "Partitioning failed: {:L} super-cells were split across ranks!",
+                           totalSplitSuperCells ) );
+
+  return cellPartitioning;


Can't we group these do-check architecture into a single loop leveraging set's lookup insertion idiom (like above) and disjunction on the second value of the tuple ?

DENEL Bertrand added 4 commits March 28, 2026 22:46

Draft

3669b67

working version

435e91c

Clean

e97bf72

Clean

0feb386

bd713 self-assigned this Mar 30, 2026

bd713 added the type: bug Something isn't working label Mar 30, 2026

bd713 requested review from MelReyCG, OmarDuran, corbett5, dkachuma, rrsettgast and wrtobin as code owners March 30, 2026 04:05

bd713 added ci: run CUDA builds Allows to triggers (costly) CUDA jobs ci: run integrated tests Allows to run the integrated tests in GEOS CI ci: run code coverage enables running of the code coverage CI jobs labels Mar 30, 2026

DENEL Bertrand added 4 commits March 30, 2026 13:18

Fixed map

397bfa7

Clean

32b26b5

gpu

3bef343

Clean Omar

768299d

bd713 mentioned this pull request Mar 31, 2026

fix: Fracture/3D cell co-location in parallel mesh redistribution #3985

Closed

dkachuma approved these changes Apr 21, 2026

View reviewed changes

Merge branch 'develop' into fix/bd713/fractureDistribution

4f04ed0

OmarDuran requested a review from jafranc as a code owner May 4, 2026 19:03

Merge branch 'develop' into fix/bd713/fractureDistribution

536ae04

jafranc approved these changes May 6, 2026

View reviewed changes

DENEL Bertrand added 3 commits May 8, 2026 19:42

JF Review

e9b6c01

Missing protection

5e0009c

stdContainer 1

bb5b023

	idx_t options[4] = { 1, 0, 2022, 0 };
	idx_t options[4] = { 1, 0, 2022, PARMETIS_PSR_UNCOUPLED };

Conversation

bd713 commented Mar 30, 2026

Uh oh!

codecov Bot commented May 5, 2026

Codecov Report

Uh oh!

jafranc left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants