Why is mesh initialization still very slow after converting my Cubit mesh to HDF5?

PyLith
1

Hi there,

I am running a large 3D model with PyLith 5.0.1. Because reading a Cubit Exodus mesh is serial, I converted the mesh to PETSc HDF5 so that PyLith could read it in parallel. I expected the HDF5 workflow to make initialization faster, but the run now remains at Inserting cohesive cells for several hours.

I would like to know whether my workflow and configuration are correct, and why converting the mesh to HDF5 does not improve this part of the initialization.

My mesh has approximately 16.7 million tetrahedra and 2.9 million vertices. The fault is a Cubit side set named FAULT containing 54,754 triangular faces. FAULT_EDGE is a node set containing 635 vertices.

First, I converted the Exodus mesh using:

pylith_convertmesh \
  --mesh_initializer.phases.read_mesh.reader=pylith.meshio.MeshIOCubit \
  --mesh_initializer.phases.read_mesh.reader.filename=sumatra_v5.exo \
  --mesh_initializer.phases.write_mesh.writer.filename=sumatra_v5.h5

I then configured PyLith to read the HDF5 mesh in parallel:

[pylithapp.problem.mesh_initializer]
phases = [read_mesh, distribute_mesh, insert_interfaces]

phases.read_mesh = pylith.initializers.MeshReader
phases.distribute_mesh = pylith.initializers.MeshDistributor
phases.insert_interfaces = pylith.initializers.MeshInsertInterfaces

[pylithapp.problem.mesh_initializer.phases.read_mesh]
reader = pylith.meshio.MeshIOPetsc
reader.filename = sumatra_v5.h5
reader.coordsys.space_dim = 3

I did not include reorder_mesh because pylith_convertmesh had already reordered the mesh. I did not include refine_mesh because I want to use the original mesh resolution.

The fault is configured as:

[pylithapp.problem]
interfaces = [fault]

[pylithapp.problem.interfaces.fault]
label = FAULT
label_value = 1
edge = FAULT_EDGE
edge_value = 1
ref_dir_1 = [-0.1, -0.1, 0.9]

I run PyLith with:

pylith sumatra_v5.cfg --nodes=24

The HDF5 mesh appears to be read and distributed successfully. PyLith then prints:

-- Inserting cohesive cells.

and remains there for more than three hours. All 24 MPI processes continue to use one CPU core each. The workstation has 251 GiB of RAM, about 216 GiB remains available, and there is no swapping, so it does not appear to be waiting for memory or disk I/O.

I also attached gdb to one of the running MPI processes. It was inside:

DMPlexTransformCheckImpingingPoint_Internal
DMPlexTransformCheckImpingingStratum_Internal
DMPlexTransformCheck_Cohesive
DMPlexTransformCheck

My questions are:

  1. Is the Exodus-to-HDF5 conversion command above correct?
  2. Are [read_mesh, distribute_mesh, insert_interfaces] the correct phases for a converted HDF5 mesh when I do not want refinement?
  3. Does HDF5 only improve mesh I/O, while insert_interfaces still performs a separate operation on the complete mesh?
  4. Is more than three hours in DMPlexTransformCheck_Cohesive() expected for a mesh of this size?
  5. Is there a recommended way to reduce this initialization time, especially when running many Green’s-function simulations with the same mesh and fault geometry?

I have checked the FAULT side set for duplicate faces, invalid element-side entries, disconnected components, and nonmanifold edges, and did not find any problems. I can provide the complete configuration files, parameter JSON, log, gdb backtrace, and mesh metadata if needed.

Thank you.

2

Thank you for the detailed description. Optimizing the new cohesive cell insertion implementation is high on our to-do list.

  1. The Exodus-to-HDF5 conversion command looks correct.
  2. When reading in parallel, you are correct in that you should skip reordering the cells. The phases you have are similar to the predefined pylith.initializers.Parallel component.
  3. Reading the mesh using the PETSc HDF5 format only changes how the mesh is read. When you distribute the mesh before inserting cohesive cells, the performance for inserting cohesive cells should be nearly identical.
  4. We are aware that the new algorithm for inserting cohesive cells is not as fast as the previous algorithm, which did not work in parallel. We have not yet had a chance to do benchmark scaling tests with large meshes.
  5. Running on a cluster with multiple compute nodes will better distribute the work and would likely be faster than running on a single workstation with the same number of CPUs. We have verified that the new algorithm running in parallel does scale well for small- to medium-sized meshes compared to running it in serial.
3

Thank you, this confirms that my mesh conversion and initialization phases are set up correctly.

I collected GDB backtraces from all 24 MPI ranks after the program had remained at “Inserting cohesive cells” for about 9 hours. The results suggest that the main bottleneck may be in the consistency check after the cohesive transform, rather than mesh reading.

Four ranks were still executing:

PetscFindInt
ISLocate_General
DMLabelGetValue
DMPlexTransformGetSourcePoint
DMPlexTransformCheckImpingingPoint_Internal
DMPlexTransformCheck_Cohesive
DMPlexTransformCheck

One backtrace showed PetscFindInt() operating on an index set containing 121,283 entries.

The other 20 ranks had reached DMPlexDistributeOwnership() inside Distributor::distributeOverlap() and were waiting in MPI_Allreduce. This looks like a large load imbalance during or immediately after DMPlexTransformCheck().

The mesh contains 16,712,488 tetrahedral cells and 2,894,795 vertices. The fault side set contains 54,754 triangular faces. I am using 24 MPI ranks on one workstation.

The PETSc revision reported by the PyLith binary is:

v3.25.2-216-g1c60be80d13

Is DMPlexTransformCheck() required for every production run, or is there an option to disable this check after the mesh has already been validated?

4

Skipping DMPlexTransformCheck() is one option. From a developer’s perspective, fixing the performance issue may be nearly as easy.

5

Matt fixed the bottleneck in DMPlexTransformCheck(). It is now over 1000x faster.

We are preparing to release v5.0.2 in the next week or so. We will incorporate this fix.

1 Like
6

That sounds like a worthwhile patch! Isn’t it funny how sometimes someone just needs to point you to a function to take a look and fix an algorithmic bottleneck you didn’t know about!

Cheers
WB

1 Like
7

This is excellent news! Thank you for the update, and thank Matt for fixing the problem, I’m looking forward to using v5.0.2.