run-problem

name	run-problem
description	Set up and launch an Aether.jl test problem on CPU or GPU, serial or under MPI — environment recipe, problem-script anatomy, GPU launch and float-type discipline, mpiexec launch pattern, distributed reductions, and verified Sod/KH baselines. Use when running a physics problem, convergence study, GPU production run, or serial-vs-MPI comparison.

Running a problem (CPU and GPU, serial and MPI)

Problem scripts are throwaway drivers — keep them outside the repo tree (scratchpad or another scratch location; AGENTS.md forbids scratch files in the repo). Only problems that belong in the suite go under test/problems/.

Environment

MPI is a weak dependency: using Aether alone never loads the MPI extension. A distributed script needs a scratch environment with both Aether and MPI, which activates MPIExt automatically:

# <scratch>/env/Project.toml
[deps]
Aether = "f2bfbca3-4e38-4f9f-8f4d-cc8c0e3d2b3a"
MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[sources]
Aether = {path = "<path-to-Aether.jl>"}

julia --project=<scratch>/env --startup-file=no -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()'

Serial-only scripts can instead use --project=<path-to-Aether.jl> directly. Precompile once before an mpiexec launch so ranks don't race.

For GPU runs, add CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" to the scratch environment above. A distributed GPU run needs both CUDA and MPI.

Problem script anatomy

test/problems/sod.jl is the template of record. The skeleton:

using Aether
using MPI

MPI.Initialized() || MPI.Init()
communicator = MPI.COMM_WORLD
rank = MPI.Comm_rank(communicator)
nranks = MPI.Comm_size(communicator)

mesh = Mesh(CPU(); size = (N, 1, 1), extent = (1, 1, 1),
            cells_per_block = (N ÷ nblocks, 1, 1),
            communicator = nranks == 1 ? nothing : communicator)
simulation = Simulation(mesh;
                        eos = IdealHydro(1.4, 1e-12, 1e-10),
                        boundary_conditions = BoundaryConditions(ix1 = OutflowBC(),
                                                                 ox1 = OutflowBC()),
                        reconstruction = PLM(),
                        riemann_solver = HLLC(),
                        stepper = RK2(),
                        cfl = 0.8)

# primitives (ρ, v1, v2, v3, e) with e = p / (γ - 1); pointwise at cell centers
set_initial_condition!(simulation) do x1, x2, x3
    ρ, p = x1 < 1/2 ? (1.0, 1.0) : (0.125, 0.1)
    return (ρ, 0, 0, 0, p / (1.4 - 1))
end

run!(simulation; stop_time = 0.2)

Post-processing — loop local blocks, interior indices only:

w = on_architecture(CPU(), simulation.state.w0)          # w[i, j, k, var, m]
geometry = on_architecture(CPU(), mesh.block_geometry)   # local blocks only
indices = mesh.block_indices
(; is, ie, js, je, ks, ke) = indices
for m in 1:mesh.nblocks, i in is:ie
    x = Aether.Meshes.x1v(mesh.coordinates, i, geometry[m], indices)
    ρ = w[i, 1, 1, 1, m]
    p = (γ - 1) * w[i, 1, 1, 5, m]
end

• cfl is the fraction of the stepper's stability limit: 0.8 is fine in 1D; multidimensional runs want roughly 1 / ndimensions.
• nghost(reconstruction) must fit the mesh's ghost width (default 2; WENOZ/PPM need more — pass a wider nghost to Mesh).

Launching

Serial: julia --project=<env> --startup-file=no script.jl.

MPI — never call the artifact mpiexec directly from the shell: it fails with libhwloc.so.15: cannot open shared object file because the JLL library paths only exist inside Julia's environment. Launch through Julia, the same pattern test/runtests.jl uses:

julia --project=<env> --startup-file=no -e '
using MPI
run(`$(MPI.mpiexec()) -n 4 $(Base.julia_cmd()) --startup-file=no
     --project=$(Base.active_project()) script.jl arg1 arg2`)'

Run several resolutions inside one launch (loop in the script) — Julia startup dominates otherwise.

GPU runs

Verified recipe (2048² Kelvin–Helmholtz validation, RTX 4090, 2026-07-11):

import CUDA    ## `import`, not `using`: CUDA also exports `synchronize`

FT = Float64
architecture = GPU(CUDA.CUDABackend())   # guard with CUDA.functional()
mesh = Mesh(architecture, FT; size = (2048, 2048, 1), extent = (1, 1, 1),
            cells_per_block = (512, 512, 1), nghost = 3)
simulation = Simulation(mesh; eos = IdealHydro(FT(5//3), FT(1e-12), FT(1e-10)),
                        reconstruction = PPM(), riemann_solver = HLLC(),
                        stepper = RK3(FT), cfl = FT(0.4))

• Float-type discipline: the constructor validates that the EOS and the stepper carry the mesh's FT — build them as IdealHydro(FT(γ), FT(...), ...) and RK3(FT), and pass cfl = FT(...). IC functions may return plain Float64 literals (converted on assignment to the host staging array).
• Nothing else changes relative to CPU: set_initial_condition! fills on the host with one transfer, and the post-processing pattern above (host copy via on_architecture(CPU(), ...)) is architecture-agnostic. Blocks all launch in a single kernel sweep, so a handful of large blocks (e.g. 16 × 512² in 2D, 128³ in 3D) is fine.
• Time series: run! lands on stop_time exactly, so sample diagnostics by looping run!(simulation; stop_time = target) over targets. A full-array host copy per sample costs ~1–2 s at 2048² — negligible at Δt = 0.02 cadence.
• VRAM: HydroState holds 6 padded 5-variable arrays; ≈ 6 × 5 × (n1+2ng)(n2+2ng)(n3+2ng) × sizeof(FT) — ~1 GB at 2048² Float64.
• gpu_aware_mpi = true hands device buffers straight to MPI (untested here; the default host-staged path is what MPI runs get).

Measured throughput, RTX 4090 (Mzone-cycles/s; memory-bandwidth-bound — Float32 runs exactly 2× Float64, so FP64 is fine even on consumer GPUs):

configuration	Float64	Float32
2D PPM + HLLC + RK3 (2048²/4096²)	93–94	193
3D PPM + HLLC + RK3 (256³)	68	—
3D PLM + HLLC + RK2 (256³)	189	401

Wall-time rule of thumb: cycles = stop_time / (cfl_fraction × limit × dx / max wavespeed); 2048² KH to t = 2 (26.6k cycles) ran in 20.5 min Float64. Doubling 2D resolution costs 8× (4× zones, 2× cycles). Reference points: 4096² KH to t = 1.5 is ~1.9 h Float64 / ~57 min Float32. First kernel compilation adds ~1–2 min to any fresh Julia session; warm up with a few stop_cycle steps before timing anything.

Distributed rules

• cells_per_block must divide size per dimension, and total blocks ≥ nranks (every rank needs at least one block). Blocks are dealt to ranks in contiguous Z-order slices.
• mesh.nblocks and mesh.block_geometry are rank-local; never index by global block id.
• Norms: accumulate locally (divide by the global cell count), then MPI.Allreduce(partial, +, communicator).
• Assembling a solution: MPI.gather(rows, communicator) (lowercase, serializing) → rank 0 gets a vector of per-rank vectors; reduce(vcat, ...) and sort by coordinate. Only rank 0 prints or writes files.

Verification anchors (Sod, PLM + HLLC + RK2, CFL 0.8, t = 0.2)

• Results are bitwise invariant to block decomposition and rank count (verified serial 1-block ≡ serial 4-block ≡ 4-rank at N = 1024): ghost exchange copies exact values and the dt min-reduction is exact. cmp on full-precision solution dumps is a legitimate MPI regression check.
• L1 errors vs the exact Riemann solution at N = 1024: ρ 1.145301e-3, v 1.665293e-3, p 6.773966e-4 (560 cycles).
• L1 convergence orders over N = 512 → 8192: ρ ≈ 0.81 (contact-limited, expect ⅔–1), v ≈ 1.06 and p ≈ 1.02 (shock-limited, expect ≈ 1) — the correct pattern for a second-order TVD scheme on a discontinuous solution.

Verification anchors (Kelvin–Helmholtz, PPM + HLLC + RK3, CFL 0.4)

McNally, Lyra & Passy (2012, ApJS 201, 18) smooth double shear layer; test/problems/kelvin_helmholtz.jl is the in-suite template of record and has the exact ICs and the mode-amplitude diagnostic. The paper's 4096² Pencil Code reference data (M(t) and max ½ρvy², Δt = 0.02, GCI uncertainty ~2e-6): https://www.colinmcnally.ca/khcomp/khmode_rev0.txt and khener_rev1.txt.

• Diagnostic trap: M is normalized by the sum of the exponential weights e^{-4πΔ} with Δ mirrored about the mid-plane (both interfaces), so M(0) = 0.01 exactly — normalizing by cell count, or weighting only the y = 1/4 interface, silently breaks the reference comparison.
• Max |M/Mref − 1| over t ∈ [0, 1.5], Float64: 2048² → 0.01%, 512² → 0.14%, 256² → 0.59%. Transient dip minimum at t = 0.18 (M ≈ 8.1439e-3) matches to 5 digits; peak ½ρvy² at t = 1.5 within −0.12%.
• Effective growth rate d ln M/dt over t ∈ [0.6, 1.0] ≈ 2.996 (reference: 2.9966); the pure linear eigenvalue for this profile at k = 4π is σ = 2.834 (Fourier-collocation eigensolver of the linearized compressible Euler equations) — the fitted slope sits above σ because the mode is an overstability (Re ω ≠ 0) and beats during growth. Don't fit t < 0.4.