Speed up the MPAS-to-XDMF converter on large meshes

[xylar]

The problem

mpas_to_xdmf writes one 2D field per index of every extra dimension, so a 3D field such as temperature(Time, nCells, nVertLevels) becomes 80 separate fields on an 80-layer mesh. The converter was reading those fields one at a time.

That is the wrong way round for how MPAS files are stored. A single vertical level is strided across the whole variable on disk, so reading one level means reading the entire variable. Reading 80 levels one at a time therefore reads the variable 80 times.

Measured on a 4,016,915-cell, 80-layer global ocean initial condition:

	time
read one level of temperature	1.61 s
read all 80 levels together	3.85 s

So roughly an 80x read amplification, of which ~24x was pure waste. In a real polaris run this made the export of one init.nc take 62.7 minutes (50:02 for the cell fields, 12:44 for the edge fields). The HDF5 writing was never the problem — that side was measured at over 900 MB/s.

The change

Read as many indices as possible in a single pass over each variable and slice them apart in memory, instead of going back to disk for each one.

Writing one 3D variable on the mesh above:

199 s → 6.5 s (30x)

How much is read at a time is capped by a new max_read_bytes argument (--max-read-gb on the command line, 2 GB by default) so that memory use stays bounded. Variables bigger than the cap are read in as few passes as it allows, so the behavior degrades gradually rather than falling off a cliff. Peak memory is roughly two to three times the cap, because masking the fill values makes a copy.

The main structural change is that unwrapping the extra dimensions moved out of _process_extra_dims and into the writer. Previously the dataset was exploded into one variable per index before anything was written, which threw away the knowledge of which fields came from the same variable — exactly the knowledge needed to read them together. Now variables keep their extra dimensions until they are written.

Three smaller improvements came along with it:

• A new float32 option (-f/--float32) writes floating-point fields in single precision. This halves both the HDF5 file size and the time ParaView needs to load it. The mesh geometry deliberately stays in double precision so that cell shapes are unaffected, and integer fields are left alone.
• Vertex-centered fields map onto the kites of the dual mesh. That mapping was being applied as a lazy fancy-index read from netCDF, which is slow; it is now applied after the data has been read.
• The connectivity array is written as 32-bit integers when the indices fit, halving the largest static array in each file.

An alternative was tried and rejected: reading contiguous blocks of cells and writing them into pre-created HDF5 datasets with hyperslabs. That gives the fastest possible reads, but the resulting scattered writes cost far more than they save (19.8 s versus 6.5 s for the same variable), so the simpler whole-variable approach won.

Compatibility

Output is unchanged. The HDF5 datasets produced before and after this change were verified byte-for-byte identical — same names, dtypes, shapes and checksums — across cell, edge and vertex files, with multiple time steps, two extra dimensions at once, and non-contiguous index lists. The only difference is the order in which attributes are listed in the .xdmf file; the attribute-to-dataset mapping for each time step is identical, and ParaView looks attributes up by name.

All new options default to the old behavior, so existing callers and scripts need no changes.

Testing

• Existing tests pass, plus five new ones covering unwrapped field names and values, vertex-to-kite mapping, max_read_bytes not affecting results, float32 casting, and the new _process_extra_dims behavior. The three tests that target new behavior were confirmed to fail on the pre-change code.
• Spot checks against a direct read of the real 39 GB init.nc match exactly.

Note for polaris

The VizInitStep._export_xdmf step added in E3SM-Project/polaris#612 exports every level of every 3D field, which is ~40 GB for the mesh above. Once this lands in an MPAS-Tools release, passing float32=True would halve that, and restricting extra_dims to the levels actually of interest would cut it much further.