Case Study 2: Designing a Chunked, Compressed HDF5 Archive

"Choosing the chunk shape is choosing which reads are fast and which are slow. Decide on purpose." — the operational lesson every large-simulation group learns once

Executive Summary

Case Study 1 audited a mute text pipeline and chose a self-describing target. Now we build the scalable output — and go one level deeper than "just write NetCDF." The simulation here produces a long time series of 2D fields and needs an archive that is compact on disk, fast to read the way the group actually reads it, and organized so a human can navigate a thousand snapshots. That points at HDF5, where we control the hierarchy, the chunk shape, and the compression directly. This case study designs the group layout, builds a per-snapshot writer with a chunking-and-compression property list, then optimizes the two decisions that dominate performance — chunk shape and compression level — against the group's real access pattern, and estimates the payoff. Every design choice is made on purpose, with its tradeoff stated.

Skills applied - Designing an HDF5 group hierarchy for a time series (§25.3) - Building a chunked, deflate-compressed dataset writer with a property list (§25.3) - Choosing a chunk shape from the read/write access pattern — the core optimization (§25.3) - Trading compression level against write speed (§25.3) - Attaching run metadata as attributes for reproducibility (§25.3, §25.5) - Estimating and validating the storage payoff (§25.1, §25.4)

Background

The solver evolves a $1024 \times 1024$ temperature field and saves a snapshot every few steps — call it 1,000 snapshots over a run. The group's analysis habits, gathered by asking (never guess an access pattern), are: usually load one whole snapshot to make a figure; sometimes extract a small region of interest across many snapshots to plot its time history; rarely touch everything at once. Two requirements fall out immediately — whole-snapshot reads must be fast (the common case), and the archive must be small enough to keep online. HDF5 serves both if we choose its parameters deliberately.

Phase 1 — Design the Hierarchy

HDF5's groups let the file document its own organization. We use a shallow, predictable layout — a metadata group and a fields group, one dataset per snapshot under a zero-padded name so they sort:

/                       (root)
  /meta                 (group: run-level attributes)
      nx, ny, dx, dy    (attributes)
      alpha, source     (attributes)
  /fields               (group: the time series)
      /fields/step_000000   (dataset: 1024x1024 double, chunked+compressed)
      /fields/step_000001
      ...
      /fields/step_000999

The design choices, each with a reason: zero-padded names (step_000042) so the thousand datasets sort lexically the way they sort numerically; a separate /meta group so run parameters live in one predictable place rather than being repeated on every dataset; and one dataset per snapshot rather than a single 3D (x, y, time) dataset. The 3D-with-extendible-time design is also valid and more compact for pure whole-array streaming, but it needs extendible dataspaces and hyperslab writes (h5sselect_hyperslab_f, h5dset_extent_f) — more API surface, and the per-snapshot layout is simpler to write, browse, and partially recover if a run aborts. We note the alternative and choose simplicity on purpose.

Phase 2 — Build the Writer

The per-snapshot writer creates a chunked, compressed dataset under /fields. It is the example-03 pattern, generalized to take a step number and build the dataset path. (Requires HDF5 installed; not run here.)

subroutine write_snapshot(file_id, u, step, chunk, level)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use hdf5
  integer(hid_t), intent(in) :: file_id       ! the open /fields group or file
  real(dp),       intent(in) :: u(:,:)
  integer,        intent(in) :: step, chunk(2), level
  integer(hid_t)   :: space_id, dset_id, dcpl
  integer(hsize_t) :: dims(2), chunk_dims(2)
  character(len=16) :: name
  integer :: hdferr

  dims       = [int(size(u,1), hsize_t), int(size(u,2), hsize_t)]
  chunk_dims = [int(chunk(1),  hsize_t), int(chunk(2),  hsize_t)]
  write(name, '(a, i6.6)') 'step_', step        ! step_000042 (internal-file write, Ch.12)

  call h5screate_simple_f(2, dims, space_id, hdferr)
  call h5pcreate_f(H5P_DATASET_CREATE_F, dcpl, hdferr)
  call h5pset_chunk_f(dcpl, 2, chunk_dims, hdferr)     ! chosen in Phase 3
  call h5pset_shuffle_f(dcpl, hdferr)                  ! reorder bytes to help deflate
  call h5pset_deflate_f(dcpl, level, hdferr)           ! chosen in Phase 4
  call h5dcreate_f(file_id, trim(name), H5T_NATIVE_DOUBLE, space_id, &
                   dset_id, hdferr, dcpl)
  call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, u, dims, hdferr)

  call h5pclose_f(dcpl, hdferr)                        ! close every handle
  call h5dclose_f(dset_id, hdferr)
  call h5sclose_f(space_id, hdferr)
end subroutine write_snapshot

Two design decisions are already visible. The chunk shape and compression level are parameters, not hard-coded — because they are exactly what we tune in the next two phases. And h5pset_shuffle_f runs before deflate: shuffle regroups the bytes of the values (all the first bytes, then all the second, …), which makes the near-identical high-order bytes of a smooth field line up and compress far better. The compile command is h5fc -std=f2018 -Wall ... -o writer.

Phase 3 — Optimize the Chunk Shape

This is the decision the epigraph warns about. A chunk is the atomic unit of I/O and compression: to read any element, HDF5 reads and decompresses the whole chunk containing it. So the chunk shape should match how you read.

Chunk shape (of a $1024^2$ field) Fast for Slow for Verdict
$(1024, 1024)$ — one chunk = whole field reading a whole snapshot reading a small region (still fetches all 8 MB) fine if you only read wholes
$(1, 1024)$ — one row per chunk reading a row reading a column or a square region poor for 2D regions
$(128, 128)$ — square tiles whole-field and square regions of interest nothing badly the balanced choice
$(16, 16)$ — tiny tiles pinpoint element access everything else (per-chunk overhead dominates) too small

The group's common case is whole-snapshot reads and its occasional case is a square region across time. Square tiles of $(128, 128)$ serve both: a whole read touches all $8 \times 8 = 64$ chunks in order (essentially as fast as one big chunk), and a $128 \times 128$ region-of-interest touches just one chunk instead of dragging in the entire field. Tiny chunks lose to per-chunk bookkeeping; whole-field chunks make the region case read 64× more data than it needs. We choose $(128, 128)$, and — the point — we can say exactly why.

Design rule. Make the chunk the shape and size of your typical read: big enough that per-chunk overhead is negligible (tens of kB to a few MB), small enough that a partial read does not drag in data you will not use. A $128 \times 128$ double chunk is $128^2 \times 8 = 131{,}072$ bytes — a comfortable 128 kB.

Phase 4 — Optimize the Compression Level

Deflate takes a level from 1 (fast, less shrinkage) to 9 (slow, most shrinkage). The right level depends on whether writing or reading dominates, and smooth scientific fields hit diminishing returns early:

Deflate level Relative write cost Typical ratio on a smooth field Use when
0 (off) 1× 1.0× you are compute-bound and disk is free
1 ~1.2× ~2.5× writing dominates; you want some shrinkage cheaply
4–6 ~2× ~3× the sweet spot for most runs
9 ~5×+ ~3.2× archival, write-once read-many, CPU to spare

The jump from level 1 to 6 buys real ratio; the jump from 6 to 9 buys little ratio for a lot of CPU on smooth data. With shuffle already improving the ratio, level 5 is the deliberate choice: near the best ratio at a modest write cost. (Ratios are Tier-2 illustrative — smooth fields compress well, noisy ones far less; always measure on your data before committing an archive.)

⚡ Optimization, honestly. Compression trades CPU for disk and I/O bandwidth. On a machine where the disk is the bottleneck (common in HPC), compressing can make writing faster overall — fewer bytes to push — even though it costs CPU. On a machine where the CPU is saturated with physics, a high level steals cycles from the simulation. Choose the level for your bottleneck; there is no universal best.

Phase 5 — Estimate and Validate the Payoff

Quantify what the design achieves. A short program computes the storage before and after (pure arithmetic, exact output):

program archive_payoff
  use, intrinsic :: iso_fortran_env, only: dp => real64, int64
  implicit none
  integer(int64), parameter :: nx = 1024, ny = 1024, n_snaps = 1000
  real(dp),       parameter :: ratio = 3.0_dp
  integer(int64) :: raw_bytes
  real(dp) :: raw_gb, comp_gb

  raw_bytes = nx * ny * 8_int64 * n_snaps
  raw_gb    = real(raw_bytes, dp) / 1.0e9_dp
  comp_gb   = raw_gb / ratio
  print '(a, i0)',   'bytes/snapshot     : ', nx * ny * 8_int64
  print '(a, f0.1)', 'raw archive  (GB)  : ', raw_gb
  print '(a, f0.1)', 'gzip5 archive (GB) : ', comp_gb
end program archive_payoff

bytes/snapshot     : 8388608
raw archive  (GB)  : 8.4
gzip5 archive (GB) : 2.8

Hand-check: $1024^2 \times 8 = 8{,}388{,}608$ bytes per snapshot; $\times 1000 = 8.388\ldots \times 10^9$ bytes $= 8.4$ GB raw; at 3× that is $2.8$ GB on disk. The chunking-and-compression design turns an 8.4 GB archive into a 2.8 GB one that also supports fast region-of-interest reads — the storage win and the access win in a single decision.

Finally, a reader that lists the archive's contents validates the structure (requires HDF5; not run). Using h5gn_members_f to count the datasets under /fields:

integer :: nmembers, hdferr
call h5gn_members_f(file_id, '/fields', nmembers, hdferr)
print '(a, i0)', 'snapshots in archive: ', nmembers

snapshots in archive: 1000

Sanity check. A correct build reports 1,000 members under /fields, each an $8 \times 8$-chunk, shuffle+deflate-5 dataset of shape $(1024, 1024)$ that h5dump -H confirms — and a total file well under the 3 GB the estimate predicts. If the count is short, a run aborted; if a dataset is uncompressed, the property list was not passed to h5dcreate_f.

Discussion Questions

  1. The group reads whole snapshots often and square regions occasionally. Justify $(128, 128)$ chunks over $(1024, 1024)$ using both access patterns. What would you choose if they only ever read wholes?
  2. Why does h5pset_shuffle_f improve the ratio on real(dp) fields specifically, and why must it come before deflate in the filter pipeline?
  3. On an I/O-bound cluster, compression can make writing faster despite costing CPU. Explain the mechanism, and name the machine where the opposite is true.
  4. We chose one-dataset-per-snapshot over a single 3D dataset with an extendible time axis. Give one advantage of each design, and the API cost that pushed us toward the simpler one.

Your Turn: Extensions

  • Option A (build). Implement the /meta group: create it with h5gcreate_f, then attach nx, ny, dx, dy, and alpha as attributes (h5acreate_f/h5awrite_f). Why store run parameters once in /meta rather than as attributes on every snapshot dataset?
  • Option B (optimize). Add a chunk and level sweep: write the same field with chunk shapes $\{(1024,1024), (128,128), (16,16)\}$ and levels $\{0, 5, 9\}$, and predict (from the tables) which combination minimizes file size and which minimizes write time. Which single combination would you ship?
  • Option C (design, ⭐⭐⭐). Sketch the 3D-extendible-dataset alternative: a single temperature(x, y, time) with time unlimited, written a snapshot at a time via a hyperslab. List the extra calls it needs versus Phase 2, and state one workload where it clearly beats one-dataset-per-step.

Key Takeaways

  • The access pattern is the design input. Chunk shape decides which reads are cheap. Match the chunk to the typical read — square tiles for mixed whole/region access — and say why. Guessing the pattern is how archives end up slow.
  • Compression is a tunable trade, not a switch. Level 5 with shuffle is the smooth-field sweet spot; level 9 buys little extra ratio for a lot of CPU, and level 0 is right only when disk is free. Choose for your bottleneck, and measure on your data — the ratios here are illustrative.
  • Shuffle before deflate. Byte-reordering aligns the near-identical high bytes of smooth floating-point fields, and it must precede the compressor in the pipeline to help it.
  • Hierarchy makes a big archive navigable. A /meta group and zero-padded /fields/step_NNNNNN datasets let both a human and a program find their way through a thousand snapshots — the "filesystem in a file" earning its keep.
  • Estimate the payoff in numbers. 8.4 GB → 2.8 GB and fast region reads, from two deliberate parameters. Design decisions you can defend with arithmetic are the ones that survive review.