[certik]

Summary

One comprehensive performance/discussion draft containing all current 2D advection-diffusion improvements, kept in focused commits so their CPU/GPU trade-offs can be discussed together. No platform-dependent layout extension is implemented, and no claim of universal performance optimality is made.

Depends on #70. Its two unchanged correctness commits are included while that PR remains unmerged. This PR does not modify #70.

Scope

• Coefficient folding: replace the two negate/reverse passes in the three 1D operator helpers with one equivalent array assignment.
• Resident Fortran RK4 workspace: retain fields and flux storage across stages/timesteps, with two kernels per RK4 stage. Preserve the order-4 mimetic gradient/interpolation/divergence closures, RK4 accumulation and fixed boundaries. The original formulation remains available with --solver operators; order 4 defaults to workspace and order 2 uses operators.
• Benchmark controls: configurable workload, the substantial 2048x2048/1000-step Gaussian case, warmup/state restoration, separate setup/solve/I/O reporting, optional GPU telemetry, raw input/output and normal CSV output.
• Native reference: the FP32 C CPU oracle and Metal implementation, bounded two-deep scheduling and signed-index comparison. Native sources remain under benchmarks/, outside FPM's automatic example discovery.
• Reproducible tooling and validation: a Formal-owned Pixi manifest/lock with FPM 0.13.0, explicit compiler/SDK selection, guarded GPU launches, full-field/evolution/boundary comparisons, portable callbacks/API exports, and checked size arithmetic, finite timesteps and strict metadata/tolerances.

No performance component is dropped. No compiler source, CI version removal, test skip or speculative layout-switching feature is added.

Measurements

Clean merged-main GPU result

Apple M4 Max, actual Gaussian input, 2048x2048 cells, 1000 steps, three interleaved pairs. Formal and native Metal use identical raw input and FP32 timestep; native geometry is direct 2D indexing with 128x1 threadgroups.

Implementation	Median solve wall	Median completed GPU time
Formal + LFortran workspace	2.167132 s	2.158323 s
Native C + Metal	2.188441 s	2.187208 s

The compiler is an unpatched merged-main checkout at lfortran/lfortran@2247d35082ea190408486b93c12ec2ead4771735, built with LLVM 11.1.0. Maximum full-field difference was 5.96e-8, relative L2 approximately 5.98e-9; fields evolved and fixed boundaries were preserved.

Solve timings exclude compilation, initialization, warmup/JIT and field I/O. Native dispatch/synchronization counts in the detailed report are schedule-derived, not observed counters. GPU command-buffer time is not hardware occupancy, and bandwidth estimates are models rather than DRAM-counter measurements.

Controlled CPU comparisons

Each compiler used three interleaved comparisons with identical inputs and flags within its workload. The original/folded operator variants differ only in the three coefficient-helper implementations.

Compiler and workload	Original operators	Folded-coefficient operators	Workspace speedup vs folded operators
GFortran 16.2, 51x51/2000 steps	0.769787 s	0.769291 s	2.9067x
LFortran CPU --fast, 51x51/1500 steps	1.865753 s	1.866891 s	7.1055x

The coefficient-only changes are about 0.06% in opposite directions, within the observed sample spread. These measurements do not establish a CPU regression, speedup or equivalence across architectures/workloads. The workloads differ between compilers, so the table is not an absolute GFortran-versus-LFortran comparison.

The earlier GPU coefficient-only result, 9.603519 -> 4.057951 s (2.3666x) for 51x51/500 steps with byte-identical fields, remains explicitly historical: it used an older optimized compiler branch, not this clean-main measurement.

Verification and reproduction

The measurements above were captured at 70e73ee; subsequent changes harden host input validation and comparison tools rather than changing the numerical kernels. Kernel equivalence was checked, and four real small GPU executions were repeated at the current ef4fb6e head, with CPU-oracle, finite-field, evolution and fixed-boundary checks. Gaussian maximum difference was 2.33e-10; the nonzero-boundary case agreed exactly.

Additional checks include allocation-free overflow/cap probes, non-finite timestep rejection, strict metadata/tolerance negative cases in normal Python and python -O, and LFortran's 32/32 FPM tests in normal and fast modes. The independent code/measurement review has no remaining blocker or rework finding.

For setup, exact compiler flags, medians/ranges, input provenance and paired commands, see:

• Benchmark README
• Performance and reproduction details

The managed tooling selects the real Apple SDK compiler after Pixi activation; it does not patch LFortran or install a fake clang wrapper. Default Pixi Clang 11 still has a current-Foundation SDK incompatibility; the separate compiler-driver work is tracked in lfortran/lfortran#12628.

Known timing follow-up: raw dumps remain outside final_io_seconds. That counter is not total final-output time, and the reported solve timings are not whole-process wall time. This limitation is documented; the timing code is unchanged.

CI limitations

Current CI at ef4fb6e passes 11/17 compiler jobs, including current LFortran and the other Flang configurations:

• LFortran 0.64-0.66: duplicate LLVM types for the example's numerical_arguments_t.
• Flang 21: compiler crash while lowering the LOCAL clause.
• macOS LFortran 0.67: Homebrew/Z3 loader failure at lfortran --version, before Formal compilation.

These remain documented compiler/packaging limitations; versions and checks have not been removed to make CI green. Pre-existing GFortran 3D failures are also kept separate from this 2D work. The PR remains a discussion draft.

Performance portability discussion

1. The coefficient change is not a storage-layout switch. Ap = -A(m:1:-1,n:1:-1) preserves coefficient values, shape and Fortran storage order. Compiler-dependent traversal and temporary generation can still affect performance.
2. Separate the effects. Workspace reuse, fusion/kernel granularity, allocation and transfer costs differ from coefficient folding. Two velocity scratch arrays are still allocated/copied once per advance call inside the solve; this is not allocation-free execution.
3. Compiler mapping versus physical storage. The merged work in gpu: coalesce contiguous stores and add optional Metal profiling lfortran/lfortran#15506 maps adjacent GPU threads to the contiguous subscript without changing the Fortran array layout. No such compiler implementation is part of this Formal PR.

Which changes benefit other CPUs and workloads? Can the compiler own more temporary elimination, legal loop/thread mapping, fusion or layout transformations, avoiding separate source strategies? If annotations or a language extension become necessary, which measured trade-off and guarantees should motivate them?

Any automatic physical-layout transformation must preserve observable Fortran element order, storage association and ABI/interoperability requirements, including SEQUENCE, BIND(C), EQUIVALENCE, TRANSFER and C pointers. Those are open discussion questions, not a proposed extension hidden in this implementation.