sandbox/acastillo/output_fields/readme

Output Fields

This module provides routines for field output and specialized quantity computations for post-processing simulation results.

Quantity Computation Routines

reference_height()

Location: available_potential.h

Description: Routines useful to evaluate the reference height required for the available potential energy in stratified flows.

Example: Examples are available in test_available_potential.c and heating.c.

probability_distribution_1D() and probability_distribution_2D()

Location: histograms1D.h and histograms2D.h

Description: Histogram, PDF and (in 1D) CDF of a scalar field, volume- or f-weighted, via GSL. available_potential.h and available_potential_vof.h both build on the 1D version. The 2D version stores no CDF, a joint histogram having no single natural cumulative distribution.

Example: Examples are available in tests_histograms/.

strain_rate_sq()

Location: strain_rate.h

Description: Squared magnitude of the strain-rate tensor, S^2 = S_{ij}S_{ij}, from centred differences on the cell-centred velocity, mainly for dissipation diagnostics (\varepsilon = 2\mu S_{ij}S_{ij}). Metric factors are not included, so this is Cartesian-only; the trace is not removed, consistent with \nabla\cdot\mathbf{u} = 0.

Example: Examples are available in test_strain_rate_affine.c (exactness on an affine field) and test_strain_rate_smooth.c (second-order convergence).

count_phase()

Location: count_phase.h

Description: Tags the connected regions of one phase of a VOF field f and writes, per region, its volume, centre of mass, volume-weighted average velocity and interfacial area (perimeter in 2D). above selects droplets (f > threshold) or bubbles (f < threshold). Works in 2D and 3D; uses the globals f and u.

Example: Examples are available in test_count_phase.c and test_count_phase_bubbles.c.

spectrum_scalar_stack() and spectrum_vector_stack()

Location: spectra/spectra.h

Description: Shell-averaged horizontal spectra of a list of fields, on one or more planes of constant z, written as ASCII or appended to a single HDF5 file along its time axis. A plane is sampled onto a regular lattice with spectra_sample.h, transformed with GSL and binned in |k| by spectra_shell.h; bin 0 holds the squared mean, so \sum_k E(k) is the mean square of the plane exactly. Requires GSL, and HDF5 for the second writer. 3D only, and for a bi-periodic domain: the planes cover one period, so no window and no detrending are needed.

Sample at the finest grid size, m = 1 << maxlevel, and keep the planes inside the refined region. Any coarser lattice lands on cell faces rather than centres, and on an adaptive grid the resulting displacement varies across the plane, so it does not cancel in |\hat c|^2.

Example: Examples are available in tests_spectra/, covering the sampling (test_spectra_sample.c), the transform and binning (test_spectra_modes.c), the sampling rule above (test_spectra_amr.c) and the two writers (test_spectra_ascii.c and test_spectra_hdf5.c).

profile_foreach_region(), profile_product_foreach_region(), profile_dissipation_foreach_region(), profile_dissipation_scalar_foreach_region()

Location: profiles/profiles.h

Description: Extracts profiles (mean, product, velocity-gradient/strain-rate dissipation) from specified regions of the computational domain. Adapted from Antoon’s profile utilities.

Example: Examples are available in tests_profiles/.

File Format Output Routines

output_vtu() and output_vtu_box()

Location: vtu/output_vtu.h and vtu/output_vtu_box.h

Description: Compatible with the VTK XML unstructured grid format which can be read using ParaView. These functions write scalar and vector fields in binary format for efficiency. The unstructured grid format is required for quadtree and octree meshes. In MPI mode, each process writes its own file, linked together using .pvtu master files.

Key Features:

• Binary data output for performance
• MPI support with .pvtu master files
  
• Box region selection with output_vtu_box()

Examples: Examples are available in the tests_outputs/ directory:

• test_output_vtu.c - Basic VTU output
• test_output_vtu_box.c - Region-specific output

output_xdmf()

Location: xdmf/output_xdmf.h

Description: Compatible with the XDMF Model and Format which can be read using ParaView or Visit. Data is split into Light data (XML metadata) and Heavy data (HDF5 binary). This used parallel I/O and is compatible with multiple visualization tools including ParaView, Visit, Matlab, and Tecplot.

Installation Requirements:

# Serial HDF5
sudo apt install libhdf5-dev hdf5-tools
export HDF5_INCDIR=/usr/include/hdf5/serial
export HDF5_LIBDIR=/usr/lib/x86_64-linux-gnu/hdf5/serial

# Parallel HDF5  
sudo apt install libhdf5-mpi-dev hdf5-tools
export HDF5_INCDIR=/usr/include/hdf5/openmpi
export HDF5_LIBDIR=/usr/lib/x86_64-linux-gnu/hdf5/openmpi

which have to be linked, for instance by adding

CFLAGS+= -disable-dimensions
HDFLIBS=-I$(HDF5_INCDIR) -L$(HDF5_LIBDIR) -lhdf5 -lhdf5_hl
LIBS+=-I$(BASILISK) $(HDFLIBS) 

to your Makefile.

HDF5 Setup for Basilisk:

Currently, HDF5 can’t pass through qcc without adding

echo "@include <hdf5.h>" > $BASILISK/ast/std/hdf5.h
echo "@include <hdf5_hl.h>" > $BASILISK/ast/std/hdf5_hl.h
echo "typedef hid_t, hsize_t, herr_t, H5L_info_t;" >> $BASILISK/ast/defaults.h

and disabling dimensional analysis. This is more an HDF5 problem than a Basilisk problem, but it’s still worth using.

Examples: Examples are available in tests_outputs/test_output_xmf*.c.

output_vtkhdf() and output_vtkhdf_box()

Location: vtkhdf/output_vtkhdf.h and vtkhdf/output_vtkhdf_box.h

Description: Compatible with the modern VTK-HDF format which can be read using ParaView or Visit. This format provides robust parallel capabilities. VTK-HDF files use a structured HDF5 layout with VTKHDF groups containing datasets for points, connectivity, and field data.

Examples: Examples are available in tests_outputs/test_output_vtkhdf*.c.

Running the tests

tests_profiles/, tests_quantities/ and tests_spectra/ each carry a run_tests.py driver:

load_gcc
cd tests_quantities && ./run_tests.py all

It is only a driver: each test is still run by make <test>.tst, which stays the single source of truth for how a test is compiled, executed and compared against its reference. A test may be run against several grids, the grid name selecting the dimension too (multigrid3D/octree are 3D); each combination has its own frozen <test>.<grid>.ref, since a coarser grid legitimately resolves some quantities less well. Note that multigrid3D decomposes the domain into octants and so requires 8^i MPI ranks.

The tests themselves compute their norms in C, using foreach reductions so that the numbers do not depend on the rank count, and write them to a small table. The thresholds, convergence orders and structural checks live in a Python checker beside the sources, whose report is what the .ref freezes. A test whose dependencies are missing is skipped by the driver rather than allowed to build into a different program: tests_spectra/ skips its HDF5 test when hdf5.h is not on the include path.

Note: Routines in the deprecated/ folder may be outdated and are provided for compatibility purposes only. ## Writing large snapshots quickly: stripe the output directory

dump() on a Lustre filesystem is limited by the stripe count of the directory the file is created in, and on TGCC/irene scratch that default is 1 — every dump lands on a single OST no matter how many ranks write it. Measured on sirene1, level 11 multigrid3D, 1024 ranks, 192 GiB per dump:

4.3x, for one line in the job script, and nothing to compile:

mkdir -p $OUTDIR
lfs setstripe -c 8 $OUTDIR      # BEFORE anything is written into it

Lustre fixes a file’s layout when the file is created, and the layout is inherited from the parent directory. Striping a directory that already holds files does not restripe them, so this has to run before the first dump — in the job script, next to the mkdir.

-c -1 asks for every OST in the directory’s pool, which is not necessarily every OST in the filesystem: sirene1 has 40, but pool t0_D0da23 exposes 8, so -c 8 and -c -1 are the same request there. Check with

lfs getstripe -d $OUTDIR        # what a new file will inherit
lfs df . | grep -c OST          # how many OSTs the pool offers

The same applies to any large output — VTU, XDMF, HDF5 — not just to dump().

Two cautions. Striping widely is right for a few large files and wrong for many small ones, where it multiplies metadata cost; the rule of thumb is one stripe per gigabyte or so, capped by the pool. And these numbers were single unrepeated samples taken while the filesystem was otherwise idle: the 8 vs -1 difference is within their noise, so read the table as 1 is slow, more than 1 is much better, not as a ranking of 8 against -1.

Three things learned getting there

This recipe is what survived an attempt to make dump() faster by rewriting it — a packed, chunked, collective MPI-IO writer that reached the point of being byte-identical to dump() at 1024 ranks and was then 30x slower (1462 s against 49 s for a 192 GiB dump). It was abandoned. What is worth carrying forward:

Measure the machine before believing a diagnosis. The rewrite was justified by a production figure of 71 MiB/s. Stock dump() on the target filesystem measures 1000–4300 MiB/s; the 71 was striping and contention, not the access pattern. Reproducing that one number first — an afternoon, no code — would have prevented the whole project.

A rank’s cells are contiguous in Z-order. On trees as much as on multigrids, z_indexing() gives each rank one unbroken run, so stock dump() does one seek per rank, not one per cell. Any optimisation premised on scattered per-cell seeks is aimed at something that does not happen.

How to time a writer on a shared filesystem. One job at a time — two jobs overlapping made each the other’s noise and produced a result that was later retracted. fsync() inside the timed region, or the clock stops when the bytes reach the page cache and you measure memcpy. Alternate which writer runs first. Repeat. Applying all four reversed the sign of every earlier number.

The code, its tests, the full measurement log and the plan are archived outside this repository, in basilisk_faramix2/diagnostics/dump/.