Skip to content

Faster and parallel reading; velocity-aware GMSH preprocessing - #95

Merged
davschneller merged 19 commits into
masterfrom
davschneller/modernize
Sep 22, 2026
Merged

davschneller merged 19 commits into
masterfrom
davschneller/modernize

Conversation

@davschneller

@davschneller davschneller commented Sep 22, 2026 •

Copy link
Copy Markdown
Contributor

Also add a lot of tests.

The new CTest suite runs pumgen with 1 to 4 ranks on small gmsh-generated
fixtures and compares the written datasets with h5diff against reference PUML
files. fixtures/generate.py computes the references from the gmsh model
directly, i.e. independently of PUMGen. The element and vertex counts of the
layered fixture do not divide evenly among 3 or 4 ranks. CI now runs the tests
for the Release and the Debug (AddressSanitizer) build.

AI-generated. Model: Claude Opus 5.5
Vertices and elements were distributed in ceil-sized chunks, while
SerialMeshFile allocates the chunks given by getChunksize. Whenever the counts
do not divide evenly (e.g. 815206 elements on 3 ranks), ranks received more
data than their buffers hold; pumgen then aborted or hung. Group numbers and
boundary conditions were routed to their elements with the same mismatch, and
readBoundaries only read the first nBoundaries/P faces in many cases, so with
2 ranks half of the boundary conditions were silently dropped.

All four readers now use getChunksize/getChunksum, a new getChunkOwner maps
elements to ranks, and groups and boundaries are read completely and ended by
an explicit terminator message. The serial reader no longer reads past a
group section that ends exactly at a chunk boundary or past the last section.

AI-generated. Model: Claude Opus 5.5
--filter-enable failed for every mesh: the chunk shape of the one-dimensional
datasets (group, boundary) was set with rank 2 and a zero second dimension,
which HDF5 rejects. The chunk rows were also derived from the local size,
although the collective dataset creation needs the same shape on all ranks.
The chunk shape now uses the dataset rank and the global size.

With more ranks than cells, the empty ranks dereferenced the end iterator of
their empty insphere list. They now contribute +inf to the minimum and take
part in the collective writes with an empty selection.

AI-generated. Model: Claude Opus 5.5
The write chunk size was read from --filter-chunksize, so --chunksize had no
effect and e.g. --filter-chunksize 4096 shrank the write buffer to 4 KiB.

--order took an optional argument, so "-o 2" left the order empty (read as 0,
reported as "Order too high: 0") and consumed the input file name as a
positional argument; the order is now a required argument and has to be at
least 1. --compactify-datatypes is a flag, but required an argument and thus
swallowed the next command line argument.

AI-generated. Model: Claude Opus 5.5
Writing with compact types aborted at the geometry: the datatype was closed
whenever the option was set, including the immutable predefined float type.
Besides, the reduced integer size was derived from the number of cells rather
than from the stored values. The bit-packed boundary codes then did not fit
and HDF5 clipped them (306 of 1597 cells in the layered test mesh), and vertex
indices could exceed the reduced range as well.

The size is now the smallest one that holds the global value range of the
dataset, including a sign bit for signed types; data with negative values is
not reduced. Only copied types are closed.

AI-generated. Model: Claude Opus 5.5
Parse errors were only appended to a message and parsing continued with the
value 0, e.g. a coordinate "abc" was silently read as 0 and pumgen succeeded;
an erroneous node tag became 0 and wrapped around when shifted to a 0-based
index. Errors now throw, and parseFile returns false with the first error
message.

Further messages and checks:
- binary MSH files are rejected with a hint instead of "Expected 0"
- MSH 4.0 is rejected, since the parser implements the 4.1 layout
- unknown element types are reported instead of indexing past the node count
  table (which aborted with std::out_of_range)
- missing $Entities/$Nodes/$Elements sections are reported
The node buffer is sized by the largest element (up to 1000 nodes) instead of
the number of element types (136), which overflowed the stack in the msh2
parser, and msh2 elements without tags get tag 0 instead of the tag of the
previous element.

AI-generated. Model: Claude Opus 5.5
Only rank 0 parses the GMSH file, but every rank asked its own (empty)
builder whether the mesh has a vertex identification. Rank 0 therefore
scattered and wrote the identification alone, and pumgen hung in the
collective HDF5 write. The flag is now broadcast.

Adds a periodic fixture (a box, periodic in x) whose reference
identification is computed from gmsh's periodic node pairs; the tests
compare the identify dataset as well.

AI-generated. Model: Claude Opus 5.5
The msh4 parser used tag - 1 as the vertex index and only checked it with an
assert, so any file whose tags do not form 1..N wrote out of bounds, e.g.
dense tags starting at 1001. Element and periodic node references were not
checked in either parser, and duplicate tags left vertices undefined.

Now the node tags have to be contiguous; in msh4 files they may start at any
positive tag. Duplicate tags and references to unknown nodes are reported,
and $Nodes has to precede $Elements and $Periodic. Gmsh writes contiguous
tags also for models with volumes outside of any physical group, so
non-contiguous tags are rejected instead of remapped.

Parametric node coordinates (Mesh.SaveParametric = 1) are skipped instead of
being read as the next node, and facet node indices use std::size_t like the
cell node indices.

AI-generated. Model: Claude Opus 5.5
The msh4 parser picked the physical tag of an element from the surface map
only for element type 2 (linear triangles) and from the volume map for all
other types. The boundary conditions of higher-order meshes were thus looked
up with surface entity tags in the volume map (208 of 320 cells wrong in the
order-2 test mesh). The map is now chosen by the entity dimension of the
element block, and lookups no longer insert entries.

Elements of surfaces or volumes without a physical group still get 0, but
pumgen now warns about them. Partitioned MSH files, whose physical groups
live in $PartitionedEntities, are rejected instead of being read with all
groups and boundary conditions 0.

The minimum insphere assumed 4 nodes per cell and reported 0 for order 2; it
now uses the four vertices of each cell.

AI-generated. Model: Claude Opus 5.5
The lexer read the file character by character through std::istream::get and
converted numbers with strtol/strtod; on a mesh with 3.6M cells (170 MB MSH
4.1, 196 MB MSH 2.2) parsing took 2.18 s and 2.82 s. The new MshInput reads
the file in 16 MiB blocks and parses the numbers in place with std::from_chars
(with a strtod fallback for standard libraries without floating-point
from_chars), which brings parsing down to 0.41 s and 0.54 s (5.3x); the
results are identical.

The parsers become part of src/meshreader and read the values they expect
instead of generic tokens. Unneeded sections are skipped as a whole. Errors
are still reported with line and column; binary MSH 2 files get a message of
their own. A new unit test parses the fixtures with a buffer of 256 bytes, so
that tokens and section markers straddle the buffer boundaries.

AI-generated. Model: Claude Opus 5.5
Binary files are parsed as the ASCII ones, but with the values read as raw
4-byte ints, doubles and data-size integers; node tags, coordinates and
element nodes are read in bulk. Other byte orders (detected by the marker in
$MeshFormat) and a data size of 4 are supported as well. On the mesh with
3.6M cells, parsing the binary file takes 0.18 s instead of 0.43 s for the
ASCII file (and 2.18 s for the ASCII file before the buffered reader).

Errors in binary data are reported with their byte offset. The binary
fixtures are written by gmsh; fixtures/generate.py derives big-endian and
4-byte-size variants and first checks that its converter reproduces gmsh's
file byte for byte. The references of the binary tests are compared with a
tolerance, since they are built from the ASCII files (%.16g).

AI-generated. Model: Claude Opus 5.5
In one of about 90 runs, a test expecting an error did not see pumgen's error
message: mpiexec may drop the output of a process which ends with MPI_Abort.
The tests expecting an error on one rank now start pumgen as MPI singleton,
which also shortens them from about 2.3 s to 0.35 s each.

AI-generated. Model: Claude Opus 5.5
The matching kept a list of adjacent surface facets for every vertex and
intersected these lists for each face of each cell. It now marks the vertices
of the surface mesh and looks up only the faces of cells with at least three
such vertices, in the surface faces sorted by their vertices. On the mesh with
3.6M cells (151016 boundary faces), 161886 faces remain to be looked up; the
matching takes 0.105 s instead of 0.21 s, and the peak memory of reading the
mesh drops from 230 MB to 215 MB. The results are identical.

Adds a test for a surface face given twice, which remains an error.

AI-generated. Model: Claude Opus 5.5
Rank 0 kept the complete parsed mesh until pumgen ended, i.e. also while
computing the inspheres and writing. Each part (vertices, elements, groups,
boundary conditions, identification) is now released right after it has been
scattered, and the surface mesh right after the boundary conditions have been
matched. On the mesh with 3.6M cells and one rank, the memory after reading
drops from 392 MB to 192 MB and the peak of the whole run from 534 MB to
388 MB; the peak is now reached while the elements are scattered.

The scatter no longer dereferences the data of empty vectors on the other
ranks.

AI-generated. Model: Claude Opus 5.5
The settings of velocity-aware meshing (refinement cuboids, frequency, easi
file, elements per wavelength) and their XML parsing move from MeshAttributes
to src/sizing, as does the mesh size computation from EasiMeshSize, which now
evaluates the easi model for many points at once. EasiMeshSize keeps only the
SimModSuite part, i.e. finding the region of a point, and passes single points
to the shared computation; the settings and the size formula stay the same.

easi becomes a dependency of its own (option EASI, implied by SIMMETRIX), so
that tools without SimModSuite can use the same settings and sizes.

AI-generated. Model: Claude Opus 5.5
CellVertices fetches the coordinates of the vertices of the local cells which
are owned by other ranks once and provides the vertices per cell, so that
further per-cell computations can share them with the insphere calculation.

AI-generated. Model: Claude Opus 5.5
pumgen-sizefield reads the VelocityAwareMeshing element of a mesh attributes
file and writes, for each refinement cuboid, a binary grid for gmsh's
Structured field with the mesh size of the cuboid's group (given by
bypassFindRegionAndUseGroup, which is required here). The easi model is
evaluated for all grid nodes at once, distributed over the ranks in x slabs;
each rank writes its slab with MPI-IO.

Each node takes the minimum of its 3x3x3 neighbourhood, so that gmsh's
trilinear interpolation never exceeds the smallest size sampled at the
corners of a grid cell; the grid covers the cuboid with a margin of two
spacings, so that all grid cells which intersect the cuboid have defined
corners. A .geo file restricts each field to the physical volume of its group
and uses their minimum as background field; it refers to the grids relative
to its own directory.

On the layered test model (20 x 20 x 10 km, 2 km sediment), gmsh's HXT
meshes with these fields have the same edge length distribution as meshes
from the exact per-point size callback (816945 vs. 815206 cells).

AI-generated. Model: Claude Opus 5.5
With --velocity-check <mesh attributes>, pumgen compares the edge lengths of
the cells in the refinement regions of the VelocityAwareMeshing settings with
the mesh size at their barycentre, for the material of the cell's group. It
logs the distribution of the mean and of the longest edge length over the mesh
size (quantiles from a histogram with 0.23 % wide bins, the maximum, and the
number of cells above 1). This works for meshes from any source, so that
meshes of different generators can be compared by the same measure.

For the layered test model, meshed by gmsh with a size callback, the check
reproduces the independently computed quantiles (mean edge: p50 1.37, p90
1.57, p99 1.71, max 2.14).

AI-generated. Model: Claude Opus 5.5
For binary MSH 4.1 files, rank 0 now only reads the structure of the file
(entities, block headers, periodic nodes) and skips the node and element
data. With the resulting offsets, every rank reads a share of the nodes,
cells and boundary faces. The nodes are sent to the ranks owning them in the
chunk distribution; the cells are read directly by their owners, in file
order as before. The boundary faces go to the ranks given by a hash of their
vertices, their vertices are marked in a bit set shared by all ranks, and
only cell faces whose vertices all lie on the surface mesh are looked up.
Periodic nodes are joined into classes on rank 0 and sent to the owners.

On the mesh with 3.6M cells, the output is identical to the one of the
serial reader for 1, 2, 4 and 8 ranks. The highest peak memory of a rank
drops from 255 MB to 161 MB with 4 ranks and from 235 MB to 129 MB with 8
ranks; with one rank, time and memory stay about the same (1.07 s vs.
1.11 s, 404 MB vs. 378 MB).

--gmsh-reader=serial keeps reading on rank 0; ASCII files are always read
there.

AI-generated. Model: Claude Opus 5.5
@davschneller
davschneller merged commit 1f1a918 into master Sep 22, 2026
3 checks passed
@davschneller
davschneller deleted the davschneller/modernize branch September 22, 2026 20:50
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant