[MPI] Bug in triangulation preconditioning for explicit triangulations

[eve-le-guillou]

Hi all,

A bug has been brought to my attention for Explicit Triangulations. I've reproduced here using the pipeline below. Though here it uses the DiscreteGradient filter, the problem seems to be located in the preconditioning of the triangulation, rather than the filter.

To reproduce:

  import os
  from paraview.simple import *
  from paraview import vtk
  from paraview import options as pv_options
  # import ttk
  
  from paraview import vtk
  ctrl  = vtk.vtkMultiProcessController.GetGlobalController()
  rank  = ctrl.GetLocalProcessId()
  size  = ctrl.GetNumberOfProcesses()
  
  print(f"Running with {size} MPI ranks")
  
  pv_options.streaming = 1
  
  # Rank detection
  rank = int(os.environ.get("PMI_RANK") or os.environ.get("OMPI_COMM_WORLD_RANK") or -1)
  
  # Load and process mesh
  mesh = Wavelet(registrationName='Wavelet1')
  
  mesh.UpdatePipeline()
  
  tetrahedralize1 = Tetrahedralize(registrationName='Tetrahedralize1', Input=mesh)
  tetrahedralize1.UpdatePipeline()
  
  print(f"[Rank {rank}] Input mesh loaded. Applying TTK filters...")
  
  # create a new 'Redistribute DataSet'
  redistributeDataSet1 = RedistributeDataSet(registrationName='RedistributeDataSet1', Input=tetrahedralize1)
  redistributeDataSet1.UpdatePipeline()
  
  # Ensure consistent ordering for parallel runs
  arrayPreconditioning = TTKArrayPreconditioning(Input=redistributeDataSet1)
  arrayPreconditioning.UpdatePipeline()
  
  # Discrete gradient computation
  discreteGradient = TTKDiscreteGradient(Input=arrayPreconditioning)
  discreteGradient.UpdatePipeline()

To compute this pipeline use:

OMP_NUM_THREADS=1 mpirun -n 4 pvbatch pipeline.py

Expected behavior
The pipeline should generate an Unstructured Grid, distribute it across all processes and compute the Discrete Gradient. Instead, it will either hang indefinitely or crash at the preconditioning stage of the Discrete Gradient. The bug happens regardless of the number of threads.

System:

• OS: Linux
• TTK latest
• MPICH/4.3.0 or OpenMPI/4.0.3 (both fail)
• ParaView 5.13

[julien-tierny]

@pierre-guillou as far as I remember you authored the MPI part for the explicit triangulation, am I correct?
Do you know where this could come from?

[pierre-guillou]

I can reproduce this bug with TTK's up-to-date dev branch (and ParaView 5.12).
Two processes seem to be enough to trigger it.

gdb shows something akin to an infinite loop in ExplicitTriangulation::preconditionDistributedTriangles around

ttk/core/base/explicitTriangulation/ExplicitTriangulation.cpp

Line 1373 in 7bea8bd

const auto lcid{this->cellGidToLid_[j]};

for rank 0. Other ranks are waiting further down in the same method at

ttk/core/base/explicitTriangulation/ExplicitTriangulation.cpp

Line 1398 in 7bea8bd

this->exchangeDistributedInternal(

.

The infinite loop is caused by the ranges.begin & ranges.end values being both at SIZE_MAX (or -1L), plus integer overflow.

These values for the range object are computed above, in ExplicitTriangulation::preconditionDistributedCellRanges at

ttk/core/base/explicitTriangulation/ExplicitTriangulation.cpp

Line 958 in 7bea8bd

this->localCellRanges_.emplace_back(

.

Running this method in gdb shows discrepancies between the values of the cellGid_ and the cellRankArray_ arrays: for a local cell id i, we may have cellGid_[i] == -1L while cellRankArray_[i] == 0L. As a consequence, rank 0 sees all local cells as its own and fails at exchanging data on the sub-domain boundaries.

This problem only seem to affect rank 0. Other ranks seem fine.

IMHO this may stem from the computation of cellGid_ and cellRankArray_.
Do you confirm my findings?

[julien-tierny]

Thanks a lot @pierre-guillou for this quick and precise feedback!
Do you know where this inconsistency on rank 0 can come from?

[pierre-guillou]

No idea for now.

Some other interesting findings: I tried to change the default value of cellRankArray_ from 0 to -1 in

ttk/core/base/explicitTriangulation/ExplicitTriangulation.cpp

Line 148 in 7bea8bd

this->cellRankArray_.resize(this->cellNumber_, 0);

to reduce the effect this default value could have, but now this seems to mess with every ranks.

I think with the 0 default value, other ranks may wrongly assign non-owned cells to rank 0. Switching to -1 leaves cells with cellGid -1 and cellRank -1 in every rank. vertexRankArray_ might also be affected.

I think we need to make sure that cellGid_ and cellRankArray_ are well-defined for every local cell of every rank. I don't know the meaning of the -1 values in cellGid_ (apart from the default value set in

ttk/core/vtk/ttkAlgorithm/ttkAlgorithm.cpp

Line 647 in 7bea8bd

vtkCellIdentifiers->Fill(-1);

).