Heat3D — single-purpose drivers¶

MPI drivers for the 3D heat equation \(\partial u/\partial t = D \nabla^2 u\) on a uniform periodic brick \([0,N)^3\) with spacing 1 (CPU, plus HIP twins for compact FD and implicit-Euler spectral). Each method ships as its own executable, and the binaries are arranged as a model hierarchy from “raw pointers in main” to “compose three production stacks”:

Binary	Layer	What it shows
`heat3d_fd_scratch`	Version 0 — from scratch	The bottom of the hierarchy. No `HeatModel`, no `HeatGrads`, no `for_each_` helpers. Three nested `for (k) for (j) for (i)` loops over `[0, n)`, padded `lin = (i+hw)sx + (j+hw)sy + (k+hw)sz` computed by hand, neighbours accessed as `u_ptr[lin ± sx/sy/sz]`, per-step Laplacian in a plain `std::vector<double>` of size `nxnynz` with no halo. The only OpenPFC piece in the hot loop is `pfc::comm::HaloExchange::exchange()`.
`heat3d_fd_manual`	Laboratory-style FD (2nd order, 7-point)	The same physics, lifted one level up. The driver still spells out the stencil and the explicit Euler update, but it now uses `HeatModel::rhs(t, HeatGrads&)` for the per-cell RHS, the `for_each_inner / for_each_border / for_each_owned` lambda iterators for the loops, and `start_halo_exchange / finish_halo_exchange` to overlap halos with the inner-region work.
`heat3d_fd`	Compact FD (orders 2–20)	The production path. Two `pfc::data::Field<double, HostSpace>` fields (`u`, `du`), `pfc::comm::HaloExchange<HostSpace, double> halo(u, decomp, rank, MPI_COMM_WORLD)` and `pfc::gradient::FDGradient<HeatGrads> grad(u, cfg.fd_order)` — both bound to `u`. Each step is `halo.exchange()` / `pfc::data::Field::for_each_owned(du, …) / u += dt * du` (with `du[idx]` writes). For overlap, `halo.start()` / `halo.finish()` are the split API. For an FFT-safe unpadded core + separated face halos, tests use `pfc::sim::stacks::FDCPUStack` (which still bundles those concerns through `pfc::sim::DuField`).
`heat3d_fd_hip`	Compact FD on HIP	Same physics and CLI as `heat3d_fd`, on `FDGPUStack<HIPSpace>`: device `HaloExchange`, `for_each_interior_device` Laplacian, Euler axpy on the GCD. One MPI rank per GPU. Optional `HEAT3D_PROFILE_JSON` writes schema-v4 `wall_step` frames for LUMI-G scaling.
`heat3d_spectral_pointwise`	Point-wise spectral RHS + explicit Euler	The spectral twin of `heat3d_fd`: identical user-facing time loop, but the residual is built by `pfc::field::SpectralGradient<HeatGrads>` (1 fwd + 3 inv FFTs/step) instead of FD stencils. The same `HeatModel::rhs` is applied cell-by-cell through `pfc::sim::DuField`.
`heat3d_spectral`	Implicit Euler in Fourier space	HeFFTe-backed forward + backward FFT per step (`heat3d::SpectralHeatPropagator`).
`heat3d_spectral_hip`	Implicit Euler on HIP	HIP twin of `heat3d_spectral`: `HIPSpectralStack` + device multiply of \(1/(1-\Delta t D k_\mathrm{lap})\), 2 FFTs/step. Same CLI as the CPU spectral driver. `HEAT3D_PROFILE_JSON` / `HEAT3D_WARMUP` match `heat3d_fd_hip`. Built when HIP spectral (rocFFT HeFFTe) is on.
`heat3d_fd_convergence_study`	Verification, not a physics demo	Not part of the model hierarchy above: sweeps `fd_order` x `N` on a single-Fourier-mode problem and reports the observed order of accuracy of `heat3d_fd`’s own `FDGradient` stencil against its design order — see Order of accuracy (FD) below.
`heat3d_spectral_content_study`	Verification + a method recommendation	Also not part of the hierarchy: answers which spatial operator is the cheaper route to a given accuracy, as a function of how much of Nyquist the field’s content reaches — see Spectral against finite difference: where each wins below.

All three CPU FD drivers (scratch, manual, fd) compute the same thing with the same 7-point central stencil. Tests assert l2_scratch == l2_manual == l2_compact to within 1e-7. The higher binaries share HeatModel for physics and heat3d::report for the method / timing / l2_error summary line; heat3d_fd_scratch only uses heat3d::report and inlines its own physics.

Initial condition matches the diffusion examples: \(u(\mathbf{x},0)=\exp(-|\mathbf{x}|^2/(4D))\) with origin at \((0,0,0)\). The diffusion coefficient \(D\) is hard-pinned to 1.0 via heat3d::kD in include/heat3d/heat_model.hpp so all heat3d binaries share one fixed value (and their L2 outputs stay directly comparable). To experiment with a different coefficient, change the literal there and rebuild.

Problem setup¶

The Problem setup block the report contract (#112) asks every application to carry. There is no science preset here and no JSON: this application is entirely a verification and benchmarking instrument, and no configuration of it is a production science case.

Item	Description
Use case	Method comparison on one fixed PDE: several finite-difference orders and two spectral treatments of the same diffusion problem
Question	For one fixed problem, how do a stencil of order `2m` and a spectral solve compare in accuracy, cost per step, and scaling?
Domain	Uniform periodic brick `[0,N)^3`, `dx = dy = dz = 1`, coordinate origin at index `(0,0,0)`
Grid/time	`N n_steps dt [fd_order]` on the command line; `N >= 8`, `dt > 0`, `fd_order` even in `[2,20]`. No output cadence: the drivers write no snapshots, only an end-of-run summary line
Boundary conditions	Periodic on all three axes. `HeatModel::boundary_value` exists but is empty by default and no driver populates it – there is no Dirichlet or Neumann path here
Initial condition	`u(x,0) = exp(-
Key physical parameters	`D = heat3d::kD = 1.0` exactly, a compile-time constant shared by every binary so their errors are directly comparable. Nondimensional; no SI unit appears anywhere in this app
Observable	`l2_error_vs_R3_analytic_rms` against the closed-form `R^3` Gaussian `(1+t)^(-3/2) exp(-
Model maturity	numerical verification: analytical – every driver is compared against the closed-form Gaussian, and the implicit-Euler symbol, the constant-field fixed point and the reference solution each have their own closed-form checks. Physical completeness: canonical – this is the heat equation, with nothing left out. Calibration: none, by design

Why periodic, and the reference mismatch¶

The domain is periodic but the reference is the infinite-domain Gaussian, which never wraps. For a localized bump the mismatch is small – at D = 1 the profile is already below 1e-7 eight cells from its peak – but it is not zero, and it contributes to the reported error even when the numerics are perfectly consistent. The reported number is therefore an upper bound on discretization error, not a clean measurement of it. That is acceptable because every driver carries the same offset, which is exactly why D is pinned and the grid is identical across binaries.

What is not verified¶

No test performs a grid-refinement study or fits an observed convergence order. The only order-related check asserts that the L2 error at fd_order 8 and 10 is below the second-order value on the same grid – monotonicity, not an observed-order match to the design order.

Build¶

Enabled with OpenPFC_BUILD_APPS=ON (default). Requires HeFFTe (the spectral binaries use the same FFT stack as the library). When CMake finds OpenMP for C++, the CPU binaries link it so the FD paths can use hybrid MPI + OpenMP: the compact driver’s interior loop runs under a single omp parallel for collapse(2) inside pfc::field::for_each_interior; the manual driver opts into the same parallelisation explicitly via pfc::field::for_each_inner_omp and pfc::field::for_each_owned_omp; the from-scratch driver runs serial (it stays single-threaded on purpose so the bare loops do not need OMP-aware indexing). Set OMP_NUM_THREADS to control the thread count for the compact and manual drivers.

Source layout¶

Shared physics, IC, propagator, parser, and reporting headers (live in include/heat3d/):

include/heat3d/heat_model.hpp — the physics, all in one file: the heat3d::HeatGrads per-point grads aggregate (xx, yy, zz) drawn from the OpenPFC catalog so the kernel’s templated evaluators (pfc::field::FDGradient<G>, pfc::field::SpectralGradient<G>) compute exactly those second derivatives and nothing else; the single source-level constexpr heat3d::kD for the diffusion coefficient; the small self-contained heat3d::HeatModel struct (initial-condition lambda + optional boundary-value provider + per-point RHS kD * (g.xx + g.yy + g.zz)). OpenPFC-free (only <cmath> and <functional> are included); this is the only file a physicist edits to define a new heat problem or change D. See docs/extending_openpfc/per_point_grads.md for the per-point grads contract.
include/heat3d/spectral_heat_propagator.hpp — heat-specific implicit-Euler-in-Fourier-space propagator. fill_implicit_euler_symbol builds the 1 / (1 - dt·D·k²) table (Catch2-tested without a GPU). The CPU class exposes step(pfc::data::Field<double, pfc::HostSpace>&) on pfc::fft::IHostFFT.
include/heat3d/spectral_heat_propagator_hip.hpp — HIP twin: uploads that table once, then each step is a device forward FFT, multiply_complex_real_hip_impl, and inverse FFT.
include/heat3d/cli.hpp — RunConfig plus the slim per-binary parsers parse_fd / parse_spectral and their _or_print_usage wrappers. D is not a CLI knob (it lives in heat_model.hpp). Each binary already knows its own discretisation, so the parsers do not consume an argv[1] discriminator. Header-only, MPI-free, OpenPFC-free; trivially unit-testable.
include/heat3d/reporting.hpp — analytic_gaussian (closed-form reference solution on \(\mathbb{R}^3\)), fd_extra_metadata (FD/OpenMP info string), and the rank-0 report template that prints the canonical method / timing / l2_error triplet, shared by all heat3d binaries.
include/heat3d/convergence_study.hpp — the Order of accuracy (FD) study’s shared run_case(fd_order, N): single-Fourier-mode IC, pfc::gradient::FDGradient<HeatGrads> bound the same way heat3d_fd.cpp binds it, and an exact (not time-marched) closed-form evolution to isolate spatial from temporal error. Used by both heat3d_fd_convergence_study and tests/test_heat3d_fd_convergence.cpp.
include/heat3d/spectral_content_study.hpp — Heat3D driver around periodic_spectra.hpp and fd_symbols.hpp: Gaussian RK4 validation of the shipped FD stack, plus occupancy-matched family maps. Used by heat3d_spectral_content_study and tests/test_heat3d_spectral_content.cpp.

Per-binary drivers (live in src/cpu/):

src/cpu/heat3d_fd_scratch.cpp — version 0, the from-scratch driver. Inlines the IC exp(-r²/(4 D)) and the Laplacian; runs three nested for (k) for (j) for (i) loops with manual padded lin = (i+hw)*sx + (j+hw)*sy + (k+hw)*sz and stencil neighbours u_ptr[lin ± sx/sy/sz]; per-step Laplacian held in a std::vector<double> of size nx*ny*nz with no halo (visibly different flat indexing from the padded u). The only OpenPFC pieces in the hot loop are pfc::comm::HaloExchange::exchange(), plus a padded field ctor for storage and domain::create / decomposition::create for geometry. Read this one first.
src/cpu/heat3d_fd.cpp — compact FD driver. After constructing two pfc::data::Field<double, HostSpace> fields (u, du) plus pfc::comm::HaloExchange<HostSpace, double> halo(u, decomp, rank, MPI_COMM_WORLD) and pfc::gradient::FDGradient<HeatGrads> grad(u, cfg.fd_order), the hot loop reads:
```
for (int step = 0; step < cfg.n_steps; ++step) {
  halo.exchange();
  pfc::field::for_each(du, [&](const auto& idx) {
    const auto g = pfc::gradient::evaluate(grad, idx);
    du[idx] = g.xx + g.yy + g.zz;   // D = 1 in this driver
  });
  u += cfg.dt * du;            // explicit Euler, on the page
  t  += cfg.dt;
}
```
Halo, gradient evaluator, and per-cell sweep are three visible primitives composed by the user. halo.exchange() is the blocking one-shot; pfc::gradient::evaluate(grad, idx) keeps the inner loop dimension-agnostic; pfc::data::Field::operator[] is used for element access. The spectral twin still bundles those concerns through pfc::sim::DuField.
src/hip/heat3d_fd_hip.cpp — HIP twin of compact FD. FDGPUStack<HIPSpace> owns the padded device field and HaloExchange<HIPSpace>. Default step on more than one rank (HEAT3D_HALO_OVERLAP=1) calls pfc::comm::GpuHaloOverlap (Waitall): interior on a non-blocking compute stream, Faces pack/MPI on the default stream, then finish() + one-launch boundary shell. A single rank has no MPI halo, so the default falls back to blocking; =0 forces blocking; =2 is Testall (MPI_Testall pump). Built when OpenPFC_ENABLE_HIP is on.
src/hip/heat3d_spectral_hip.cpp — HIP twin of implicit-Euler spectral. HIPSpectralStack + SpectralHeatPropagatorHIP. Built when OpenPFC_ENABLE_HIP_SPECTRAL is on. CTest heat3d-spectral-hip-smoke looks for HEAT3D_SPECTRAL_HIP_CHECKSUM. LUMI submit: docs/lumi_slurm/submit_heat3d_spectral_hip_scaling.sh.

src/cpu/heat3d_fd_manual.cpp — laboratory-style 2nd-order central FD driver. The hot loop reads:

halo.start_halo_exchange(u.data(), u.size());
field::for_each_inner_omp(u, hw, stencil_step);   // overlap
halo.finish_halo_exchange();
field::for_each_border (u, hw, stencil_step);     // touches halos
field::for_each_owned_omp(u, [&](int i, int j, int k) {
  u(i, j, k) += cfg.dt * du(i, j, k);             // explicit Euler
});

with a single stencil_step lambda calling model.rhs(0.0, HeatGrads{xx, yy, zz}). Each stage is wrapped in pfc::runtime::tic(timer, "...") / toc(timer, "...") and print_timing_summary(timer, 0) prints a sorted breakdown on rank 0 at the end.

src/cpu/heat3d_spectral.cpp — implicit-Euler spectral driver. Calls heat3d::SpectralHeatPropagator::step(stack.u()) once per step (forward FFT → diagonal multiply in k-space → inverse FFT). HIP twin: heat3d_spectral_hip.
src/cpu/heat3d_spectral_pointwise.cpp — point-wise spectral RHS: the spectral twin of heat3d_fd. Built on pfc::sim::stacks::SpectralCPUStack + pfc::sim::DuField<HeatGrads, SpectralGradient<HeatGrads>>, so the user-facing time loop reads du.apply(...) / u += dt * du / t += dt; the residual is materialised by pfc::field::SpectralGradient<HeatGrads> (1 forward + 3 inverse FFTs/step) instead of a stencil sweep. The DuField shim hides halo prep + per-cell evaluation here, where the FD twin (heat3d_fd) instead spells those primitives out in main.
src/cpu/heat3d_fd_convergence_study.cpp — the Order of accuracy (FD) sweep: loops heat3d::convergence::run_case(fd_order, N) over every (fd_order, N) pair, prints the design-vs-observed-order table, and writes out/report/data/heat3d_fd_order_convergence.csv.
src/cpu/heat3d_spectral_content_study.cpp — the Spectral against finite difference study: prints the error map, the crossover fractions and the equal-accuracy cost ladder, runs the validation cases, and writes out/report/data/heat3d_spectral_content_{map,crossover,validation}.csv.

Tests:

tests/test_heat3d.cpp — Catch2 unit tests covering: heat3d::kD is pinned to 1.0; HeatModel (default IC at the origin, IC override, the rhs = kD * (xx + yy + zz) formula); the slim per-binary CLI parsers (happy paths + every rejection case); fill_implicit_euler_symbol (DC = 1, first-mode formula, wrapped \(k_y\), size mismatch) and a constant-field fixed point of SpectralHeatPropagator; the analytic reference solution; two single-rank integration tests against a padded field / explicit Euler stepper; two single-rank cases for the manual driver (smoke + L2-vs-analytic, and parity vs the compact FDCPUStack path); and two single-rank cases for the from-scratch driver (smoke + L2-vs-analytic, and parity vs the compact path to within 1e-7). Built into the test_heat3d executable and registered with CTest as heat3d-all-tests whenever OpenPFC_BUILD_TESTS=ON and Catch2 is available (set HEAT3D_ENABLE_TESTS=OFF to skip). HIP smokes: heat3d-fd-hip-smoke (production overlap), heat3d-fd-hip-blocking-smoke, heat3d-spectral-hip-smoke.
tests/test_heat3d_fd_convergence.cpp — regression guard for Order of accuracy (FD): asserts the observed order of accuracy at fd_order 4 and 8 matches the design order within a stated tolerance. Registered as heat3d-fd-convergence; pure FD (no HeFFTe needed).

Where OpenMP runs (FD)¶

In heat3d_fd.cpp the per-step sweep over du is delegated to pfc::field::for_each (include/openpfc/kernel/field/brick_iteration.hpp), which iterates the brick’s owned cells in cache-friendly k-outer / j-middle / i-inner order and yields a pfc::Int3 to the user lambda. Each cell evaluates pfc::gradient::evaluate(grad, idx) against pfc::gradient::FDGradient<HeatGrads> (per-axis stencil application from kernel/field/fd_apply.hpp), writes the per-point residual into du[idx], and the explicit-Euler step u += dt * du then updates the field in place.

In heat3d_fd_manual.cpp the same parallelisation is applied explicitly: pfc::field::for_each_inner_omp parallelises the inner-region stencil over (k, j), pfc::field::for_each_owned_omp parallelises the Euler update; the (cheap) border pass stays serial. Not parallel in either driver: the MPI halo exchange (pfc::comm::SparseExchange on separated-layout stacks elsewhere; pfc::comm::HaloExchange in these heat3d FD binaries).

So htop/top can look single-threaded when wall time is dominated by MPI (many ranks, thin subdomains), or when the launcher has pinned the process to a narrow CPU mask (OpenMP then schedules all threads inside that mask).

Single MPI rank on Linux: heat3d calls pfc::runtime::reset_cpu_affinity_if_single_mpi_rank (which uses sched_setaffinity) after MPI_Init to allow all online logical CPUs for that process, so OpenMP can use the machine without extra mpirun flags. Set OPENPFC_NO_RESET_AFFINITY (any value) to keep the launcher’s mask unchanged (e.g. if your site policy forbids overriding binding).

Several MPI ranks on one node: the reset is not applied (MPI_COMM_WORLD size \(>1\)); use launcher options so each rank gets a disjoint CPU set and set OMP_NUM_THREADS per rank accordingly, for example:

export OMP_NUM_THREADS=4
mpirun --bind-to none -n 4 ./apps/heat3d/heat3d_fd 256 25 1e-6 12
# or: export OMPI_MCA_hwloc_base_binding_policy=none

Rank 0 prints omp_max_threads and omp_get_num_procs() in the summary line. After a single-rank run, omp_get_num_procs() should match the number of CPUs OpenMP may use; if it stays at 1, check binding and whether OPENPFC_NO_RESET_AFFINITY is set.

Usage¶

heat3d_fd_scratch          <N> <n_steps> <dt>
heat3d_fd_manual           <N> <n_steps> <dt>
heat3d_fd                  <N> <n_steps> <dt> <fd_order>
heat3d_fd_hip              <N> <n_steps> <dt> <fd_order>
heat3d_spectral_pointwise  <N> <n_steps> <dt>
heat3d_spectral            <N> <n_steps> <dt>
heat3d_spectral_hip        <N> <n_steps> <dt>
heat3d_fd_convergence_study [output.csv]
heat3d_spectral_content_study [--data-dir DIR] [--no-validate|--held-out-only|--held-out-families]

D is not a CLI knob: it is fixed at heat3d::kD = 1.0 in include/heat3d/heat_model.hpp. Edit the literal there if you want to experiment.
fd_order (compact FD only): even integers 2 through 20. Halo width is fd_order/2 per face. Each local subdomain must be wider than fd_order in every dimension (the driver aborts with MPI_Abort if the interior slab is empty).
The scratch and manual drivers hard-code 2nd-order central (halo width 1) — extend the inline stencil / stencil_step lambda if you want a different order; the goal is to keep the showcase code uncluttered.
Stability (FD): explicit Euler requires a sufficiently small dt (CFL); higher spatial order does not remove this limit.
Spectral: each step applies \(\hat u \leftarrow \hat u / (1 + \Delta t\,D|\mathbf{k}|^2)\) in complex Fourier space (implicit Euler for the linear diffusion operator).

Example¶

From your build directory:

mpirun -n 4 ./apps/heat3d/heat3d_fd_scratch         64 200 0.001
mpirun -n 4 ./apps/heat3d/heat3d_fd_manual          64 200 0.001
mpirun -n 4 ./apps/heat3d/heat3d_fd                 64 200 0.001 4
mpirun -n 4 ./apps/heat3d/heat3d_spectral_pointwise 64 200 0.001
mpirun -n 4 ./apps/heat3d/heat3d_spectral           64 200 0.001
# HIP (one rank per GPU; LUMI-G uses the submit helpers under docs/lumi_slurm/):
#   heat3d_fd_hip       64 200 0.001 2
#   heat3d_spectral_hip 64 200 0.001

All three FD binaries report identical l2_error_vs_R3_analytic_rms (the manual binary additionally prints a per-section breakdown — inner / halo_wait / border / euler — from pfc::runtime::print_timing_summary).

Rank 0 prints timing_s, avg_step_time_s (MPI max across ranks), omp_max_threads when OpenMP is enabled, and an RMS L2 error against the infinite-domain Gaussian reference

\[ u(\mathbf{x},t)=(1+t)^{-3/2}\exp\!\left(-\frac{|\mathbf{x}|^2}{4D(1+t)}\right) \]

The domain is periodic; for a localized Gaussian the mismatch at boundaries contributes to the reported error even when the numerics are consistent.

Benchmarks (FD)¶

All numbers below are avg_step_time_s from rank 0 (same as MPI_MAX across ranks for these runs). Build: Release, GCC 11.2.0, Open MPI 5.0.10, HeFFTe CPU prefix on CMAKE_PREFIX_PATH. Host: login node used for development (timing varies with load; repeat for publication-quality data).

Higher-order FD and MPI (`OMP_NUM_THREADS=1`)¶

Fixed case: N=256, n_steps=25, dt=1e-6 (D = heat3d::kD = 1.0). Higher orders use wider stencils (more flops per step) but still scale down with rank count on this grid.

fd_order	ranks=1	ranks=2	ranks=4	ranks=8
8	0.651413	0.313773	0.156707	0.0852177
10	0.920717	0.460287	0.198479	0.100299
12	0.972677	0.542606	0.264795	0.155862
14	1.03012	0.58818	0.310875	0.151631
16	1.21346	0.617872	0.331	0.164514
18	1.42205	0.682012	0.36659	0.18127
20	1.51928	0.839491	0.474414	0.190416

Reproduce (after module load gcc/11.2.0 openmpi/5.0.10 and configuring with that toolchain):

export OMP_NUM_THREADS=1
mpirun -n 8 --mca btl tcp,self --mca oob tcp ./apps/heat3d/heat3d_fd 256 25 1e-6 20

Sweep example:

export OMP_NUM_THREADS=1
for ranks in 1 2 4 8; do
  for ord in 8 10 12 14 16 18 20; do
    mpirun -n "$ranks" --mca btl tcp,self --mca oob tcp ./apps/heat3d/heat3d_fd 256 25 1e-6 "$ord"
  done
done

OpenMP scaling (`fd_order=12`, one MPI rank)¶

Same grid 256³, n_steps=25, dt=1e-6, mpirun -n 1 (D = heat3d::kD = 1.0).

Correctness: For this case, l2_error_vs_R3_analytic_rms is identical for OMP_NUM_THREADS in {1,2,4,8} (bit-for-bit in repeated runs), so the OpenMP region is not changing the numerics—only scheduling.

Scaling (2026-04-29, shared login node): Real FD run, mpirun -n 1, after the in-process single-rank CPU mask reset on Linux so omp_get_num_procs() matches the node (here 32). l2_error_vs_R3_analytic_rms=9.91404e-21 was identical for \(\{1,2,4,8\}\) threads. Speedup is sub-linear (memory traffic per step + fork/join each step), but 8 threads \(\approx\) 5.3\(\times\) faster than 1 thread on this capture—consistent with a healthy OpenMP Laplacian loop rather than a pinned single-core mask.

OMP_NUM_THREADS	avg_step_time_s
1	2.86396
2	1.53138
4	0.883718
8	0.535633

for t in 1 2 4 8; do
  export OMP_NUM_THREADS=$t
  mpirun -n 1 --mca btl tcp,self --mca oob tcp ./apps/heat3d/heat3d_fd 256 25 1e-6 12
done

Hybrid MPI + OpenMP (FD)¶

Use fewer MPI ranks and OMP_NUM_THREADS matching usable cores so halo exchange stays smaller while threads parallelize over interior \((i_y,i_z)\) lines (the omp parallel for collapse(2) lives in pfc::field::for_each_interior for the compact driver and in pfc::field::for_each_inner_omp / for_each_owned_omp for the manual driver).

export OMP_NUM_THREADS=8
export OMP_PROC_BIND=close   # optional; site-dependent
mpirun --bind-to none -n 4 --mca btl tcp,self ./apps/heat3d/heat3d_fd 256 25 1e-6 12

The spectral binaries are unchanged (MPI over ranks only).

Order of accuracy (FD)¶

heat3d_fd ships central FD stencils at every even order 2–20, but until now nothing measured whether raising the order actually buys accuracy. heat3d_fd_convergence_study does: for each even fd_order in {2, 4, 6, 8, 10, 12} it sweeps N in {16, 24, 32, 48, 64} on a periodic domain of length \(2\pi\) with the single-Fourier-mode initial condition \(u_0(x,y,z) = \cos(3x)\), whose exact PDE solution is \(u(x,t) = e^{-9Dt}\cos(3x)\) — chosen (instead of the Gaussian IC the other heat3d binaries use) because it is periodic to machine precision at every grid resolution, so there is no domain-truncation error competing with the stencil error being measured.

How temporal error was kept out of the fit. Two designs were tried and rejected before the one below, because both reproduced the textbook failure mode this kind of study is warned about — “the curves all flatten at the same floor, and you are measuring the time integrator, not the stencil”:

Explicit Euler with dt = dt_safety * dx^2 / D (the standard parabolic-CFL scaling, re-derived at each case’s own dx). Euler’s O(dt) global temporal error then scales as dx^2 — the same power as the plain 7-point stencil’s own design order — so fd_order 4/6/8/10/12 all flattened onto one common dx^2 floor at fine grids instead of showing their own steeper slopes.
Explicit Euler with one fixed dt for the whole sweep (independent of dx, sized from the finest grid’s stability limit). This removes the floor’s dx^2 growth at coarse grids, but the floor itself (measured: ≈1.8e-6, matching the predicted ~0.5·D²·m⁴·t·dt·e^{-Dm²t} for that dt) was still far above round-off, and reaching a round-off-level floor this way needs a dt some 5 orders of magnitude smaller — i.e. n_steps = t_final/dt about 10⁵× larger, since one dt is shared by every case. That is minutes turning into most of a day for no benefit: a fixed, non-vanishing dt was always going to leave some floor.

What is used instead: no time-stepping loop at all. cos(3x) is an exact eigenfunction of pfc::gradient::FDGradient<HeatGrads>’s stencil under periodic wraparound — any linear, shift-invariant stencil is diagonalized by the discrete Fourier basis (a standard circulant-matrix fact) — so Lap_FD[u0] = -k_eff² · u0 for one k_eff shared by every grid point. The study verifies this numerically rather than assuming it: the measured eigenvalue_residual (RMS of Lap_FD[u0] - eigenvalue·u0, relative to u0) is ~1e-12 to 1e-14 in every case, i.e. round-off, not model error. Because u stays exactly proportional to cos(3x) for all time under the (spatially discretized, not yet time discretized) ODE, that ODE has a closed form, u(x,t) = u0(x)·exp(D·eigenvalue·t), evaluated directly with zero time-discretization error — not “small”, exactly zero, by construction. What is left after comparing against the true PDE solution is purely the stencil’s spatial truncation error. (This mirrors tests/integration/scenarios/time_integration/test_rk3_convergence.cpp in the opposite direction: that test replaces the spatial operator by its exact eigenvalue action to cleanly isolate a temporal order; this study replaces the time evolution by its exact action to cleanly isolate a spatial order.) See apps/heat3d/include/heat3d/convergence_study.hpp for the full derivation and code.

A cross-check against real time-stepping. Before discovering that a fixed dt cannot cheaply reach round-off, design 2 above was run to completion. Away from its own floor its numbers agree closely with the exact-evolution numbers below — e.g. fd_order 4, N=32: 2.6146e-4 (Euler, dt≈1.93e-6, 25939 steps) vs 2.6321e-4 (exact evolution) — a 0.7% difference consistent with Euler’s own residual O(dt) error at that dt, confirming the exact-evolution shortcut is measuring the same thing real time-stepping would, just without paying for the steps.

Measured table (LUMI, CPU, Release; out/report/data/heat3d_fd_order_convergence.csv has the full per-N breakdown):

fd_order (design)	L2 error, N=16	L2 error, N=64	observed order (range across N)
2	2.297e-02	1.465e-03	1.97 – 2.00
4	3.854e-03	1.683e-05	3.84 – 3.98
6	7.817e-04	2.333e-07	5.71 – 5.96
8	1.745e-04	3.576e-09	7.57 – 7.95
10	4.132e-05	5.820e-11	9.43 – 9.94
12	1.019e-05	1.013e-12	11.28 – 11.84

Every design order is reproduced within a few tenths of a unit, tightening towards the design value as N grows (finite-N sub-leading truncation terms, not a defect). fd_order 12 at N=64 sits on the double-precision round-off floor: its L2 error (1.013e-12) is within an order of magnitude of the measured eigenvalue_residual there (8.95e-13), so the last pairwise ratio for fd_order 12 is measuring floating-point noise as much as the stencil — expected and reported rather than hidden, per the usual caveat for high-order stencils on fine grids. fd_order 10 stays clean (well above its own floor) across the whole N range tested. test_heat3d_fd_convergence.cpp pins fd_order 4 (measured avg ≈3.92) and fd_order 8 (measured avg ≈7.80) as a regression guard, using all five grid sizes since neither gets near the round-off floor in this range.

Figure: out/report/figures/heat3d_fd_order_convergence.svg (log-log L2 error vs dx, one line per order, dotted reference slopes for orders 2, 8, 12, and the round-off floor marked). Regenerate after re-running the study with:

./apps/heat3d/heat3d_fd_convergence_study out/report/data/heat3d_fd_order_convergence.csv
python3 out/report/figures/make_figures.py   # needs matplotlib; see that file's docstring

A boundary-shell bug found while building this study (worth knowing if you use the FDCPUStack + pfc::sim::steppers::create pattern elsewhere): the first version of this study reused the same construction as test_heat3d.cpp’s Gaussian heat3d_gaussian_l2_rms helper — FDCPUStack + pfc::field::create<G>(stack.u(), order) + pfc::sim::steppers::create. FDCPUStack::u() is an unpadded Field (storage_halo() == 0); binding FDGradient to it via the raw-pointer constructor restricts the evaluator to [fd_order/2, N-fd_order/2) per axis, and pfc::sim::for_each_interior (what the steppers::create factory uses) leaves everything outside that interior untouched in du — a boundary shell of width fd_order/2 per face is silently never updated by the Euler step, for every N and every fd_order. FDCPUStack::du<G>() covers that shell separately through its face-halo buffers, but nothing does if the evaluator + stepper are built and driven directly, which is exactly what the existing Gaussian helper (and this study’s first draft) does. A Gaussian IC decays to ~0 near the domain edge, so the frozen shell is invisible in its L2 error; a periodic single-mode IC has full amplitude at the edge, so the same construction produced a flat ~3–6% error with no grid-resolution dependence at all — immediately visible, and the tell that something structural (not a tuning problem) was wrong. The fix used here is to match heat3d_fd.cpp’s own construction instead (padded Field + HaloExchange + FDGradient bound to the padded field, whose stencil covers the whole owned domain) — see apps/heat3d/include/heat3d/convergence_study.hpp. Anyone else composing FDCPUStack with pfc::sim::steppers::create directly (rather than via FDCPUStack::du<G>()) should check for the same gap.

Spectral against finite difference: where each wins¶

The order sweep above is the single-mode accuracy measurement, and one smooth mode is the most favourable case a high-order stencil can be handed. It says nothing about a field carrying real content near the grid scale, so on its own it cannot decide between the FD and spectral paths. heat3d_spectral_content_study closes that, and the derivation lives in include/heat3d/spectral_content_study.hpp.

The idea. The heat equation is linear on a periodic box, so every operator is diagonal in the Fourier basis and each mode decays independently: exactly as exp(-D|k|^2 t), and under FD as exp(D*lambda_p(k)*t) with the stencil’s own symbol. The spectral operator uses -|k|^2 itself, so its spatial error is identically zero. By Parseval the L2 error of any initial field is therefore a closed-form sum over that field’s spectrum — no simulation required to draw the map.

The parameter. Occupancy f is the fraction of Nyquist at which a family’s amplitude prescription reaches 1e-3 of its peak (a top-hat uses the same f as its support cutoff). The original map used a periodic Gaussian, which is one family in periodic_spectra.hpp. The Gaussian map depends on f alone, not on the grid — pinned by tests/test_heat3d_spectral_content.cpp and tests/unit/kernel/field/test_periodic_spectra.cpp. Other families can share f while placing different energy near high-error stencil modes; heat3d_spectral_content_study writes that comparison without a new timing campaign.

The result. Under an N^3 cost model at a fixed step count, FD order p is the cheaper route to accuracy eps iff the coarsest grid it can use still keeps the content above cbrt(cost_fd / cost_spectral) of Nyquist. With the measured per-step costs (out/report/data/heat3d_method_cost.csv) that threshold is 0.32 (FD-2) to 0.49 (FD-12): a 32x cheaper step buys only 32^(1/3) = 3.2x in grid spacing. The crossover in accuracy is at L2 ~ 6e-7, owned by FD-12; below 1e-2 FD-2 costs a twentieth of the spectral path. The full tables are in the applications catalog (Heat3D comparison; formerly out/report/16_scalability.qmd).

Validated, not asserted. Ten points of the map were re-measured by running the shipped FD stack (padded Field + HaloExchange + FDGradient<HeatGrads>) under RK4 against the exact solution of the same sampled field. Prediction and measurement agree to between eight and ten significant figures.

One implementation note worth stealing. The map needs the stencil’s dispersion defect theta^2 + lambda_p*h^2, which for order 12 at theta = 0.1 is 1.2e-19 while the two terms being subtracted are 1e-2: computed directly it is pure round-off. The way out is that the order-2M central second-difference stencil is exactly the series (2 asin(delta/2))^2 = sum_m a_m delta^(2m) truncated at m = M, with delta^2 = 2 - 2cos(theta) and a_m = 2/(m^2 * C(2m,m)). The defect is then the tail of that series — all-positive terms, no cancellation at all. The test pins the identity against the shipped EvenCentralD2 tables, so a future edit to fd_stencils.hpp cannot break it silently.

Reproduce (single rank, login node, no allocation needed):

./apps/heat3d/heat3d_spectral_content_study --no-validate
./apps/heat3d/heat3d_spectral_content_study --families-only
./apps/heat3d/heat3d_spectral_content_study --validate-families --no-validate --no-families
# Held-out family RK4 (compute node; cosine-sum ICs). The driver writes
# the CSV then exits nonzero if max |ratio-1| is not below the frozen 1e-6:
#   sbatch apps/heat3d/slurm/heldout_families.sbatch
python3 out/report/figures/make_figures.py            # needs matplotlib

Figure: out/report/figures/heat3d_spectral_content.svg.