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).

strain_and_vorticity_sq()

Location: strain_and_vorticity.h

Description: Squared magnitudes of both the strain-rate and the rotation-rate tensors, S^2 = S_{ij}S_{ij} and \Omega^2 = \Omega_{ij}\Omega_{ij}, from one set of centred differences instead of differencing the velocity twice. Together they give the Q-criterion for vortex identification, Q = \frac{1}{2}(\Omega^2 - S^2). The S^2 is the same quantity as strain_rate_sq()’s, agreeing to roundoff but not bit-for-bit — the off-diagonal terms are summed in a different order — and strain_rate_sq() stays the right call when \Omega^2 is not wanted. Same caveats as that function: Cartesian only, and neither tensor is trace-removed.

Example: Examples are available in test_strain_and_vorticity_affine.c (exactness of both invariants on an affine field, plus agreement with strain_rate_sq()) and test_strain_and_vorticity_smooth.c (second-order convergence of \Omega^2).

interfacial_power_curvature(), interfacial_power_acceleration()

Location: interfacial_power.h

Description: The rate at which the interfacial force does work on the flow, \Psi_\sigma = \int \mathbf{u}\cdot\phi\nabla f\,dV with \phi = \sigma\kappa. Since the transport term integrates to zero on a bi-periodic or wall-bounded box, in the continuum \Psi_\sigma = -d(\sigma A)/dt exactly, so the gap between the two is discretisation error and nothing else — which the interfacial area alone cannot tell you, being a state variable that cannot separate an interface leak from an advection leak.

Two routes compute it. interfacial_power_curvature() rebuilds \phi = \sigma\kappa and is the one to use by default: independent of where it is called from, usable under REDUCED, and it returns the per-cell density in its output field so a caller wanting a histogram needs no second curvature() pass. interfacial_power_acceleration() reads back the acceleration the solver actually applied and is the ground truth, but must be called from a projection event, needs the uniform body acceleration passed in, and cannot be used under REDUCED. Include after tension.h. On trees the two differ wherever the interface sits on a refinement boundary.

It is not a more accurate route to the interfacial area, and should not be used as one. On a static drop, both against their closed forms, with an exact (Vofi) fraction:

PLIC reconstructs the geometry within each cut cell, whereas \Psi_\sigma is a surface integral evaluated as a smeared volume sum weighted by \sigma\kappa, and the height-function curvature degrades first. The residual \Psi_\sigma + \sigma\,dA/dt is therefore not symmetric discretisation error: at these resolutions it is dominated by \Psi_\sigma. The one regime where the area misleads instead is a long run whose net area change is small, where VOF advection drift fails to cancel in \Delta A.

Example: Examples are available in test_interfacial_power_static.c (second-order convergence to the closed form on a static drop, and the null of a uniform velocity) and test_interfacial_power_routes.c (the two routes agreeing to roundoff on an oscillating drop), with test_interfacial_power_amr.c pinning the tree behaviour from both sides: they still agree while a refined band contains the interface, and separate by ~1e-3 once a level jump runs through it.

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/.

profile_scalar_slab(), profile_product_slab(), profile_dissipation_slab()

Location: profiles/profiles_slab.h

Description: The same three profiles as profiles/profiles.h, in the same file format, computed in one traversal of the grid instead of one per slab. Those routines call foreach_region() for each of the n slabs, each carrying its own reduction, so a profile costs n MPI collectives and the price scales with resolution in the profile direction rather than with the work per slab; here every cell is binned into its slab by coordinate in a single foreach() and the profile is reduced once.

Not a general drop-in. The two agree only where every cell of the sampled region is at the same resolution and hmin, hmax and n are chosen so that a slab is exactly one row of those cells — foreach_region() point-samples at cell centres without interpolation, so under that condition it visits the cells this traversal visits. Over a region spanning several levels, or with slabs thicker than a cell, the two weight cells differently.

Example: test_profiles_slab.c runs both implementations on one field under the required geometry and compares the files column by column (they agree to ~1e-18); test_profiles_slab_bias.c adds a refined band so the region spans two resolutions and checks that they then diverge (~1e-4), so the precondition is pinned from both sides.

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/.