Skip to content

feat(Fields): add soa_broadcast! fast path for point-wise kernels#2543

Closed
haakon-e wants to merge 1 commit into
mainfrom
he/soa-broadcast
Closed

feat(Fields): add soa_broadcast! fast path for point-wise kernels#2543
haakon-e wants to merge 1 commit into
mainfrom
he/soa-broadcast

Conversation

@haakon-e

@haakon-e haakon-e commented Jul 9, 2026

Copy link
Copy Markdown
Member

Purpose

Add Fields.soa_broadcast!(f, out, args...), a restricted, fast-compiling stand-in for Base.copyto! of a point-wise broadcast.

What

soa_broadcast! applies a point-wise function f to the per-point scalar values of args and writes the result into the pre-existing field out in a single device pass, indexing the underlying arrays linearly rather than through CartesianIndex.

  • Restricted to uniform same-space, unmasked, scalar-field arguments. f returns a scalar (for a scalar out) or a NamedTuple matching eltype(out).
  • CPU method: a serial loop. CUDA method: a plain kernel launch.
  • Marked experimental; not exported.

Why

For many-call-site point-wise functions, for example the quadrature-based 2M+P3 microphysics tendencies in ClimaAtmos, the standard broadcast makes the GPU kernel compile very slowly: minutes for the instantaneous tendency, and more than 30 minutes without completion for the full Rosenbrock 2M+P3 tendency. A compile-time decomposition (see the benchmark comment below) identifies the cause: forcing always_inline on the broadcast kernel balloons the IR and the LLVM compile. soa_broadcast! recovers the fast compile chiefly by not forcing inlining. Two further choices are secondary for compile time and kept for other reasons: linear indexing (instead of CartesianIndex), which does not change compile time but gives a roughly 1.6x faster warm kernel on the heavy tendency and reaches the speed of a hand-written CuArray kernel without leaving the Fields interface; and extracting inputs with map rather than ntuple(_, Val(N)), the safe choice for heterogeneous inputs (the N > 10 boxing fallback cited earlier did not reproduce on the benchmarked Julia/CUDA).

Notes

  • Draft: opened to unblock the downstream ClimaAtmos 2M+P3 EDMF work that consumes this API. The compile-time decomposition is in the comment below.
  • dr/data_style_broadcast_getindex targets a related goal (a specialized Base.getindex for Broadcasted{<:DataStyle}) via the generic broadcast path; this is an orthogonal, restricted fast path and the two do not overlap at the file level.

@haakon-e

haakon-e commented Jul 9, 2026

Copy link
Copy Markdown
Member Author

cc: @imreddyTeja here's the hacky fix that helps my P3 quadrature kernels compiling in reasonable time.

Add `Fields.soa_broadcast!(f, out, args...)`, a restricted point-wise
fill that indexes the underlying arrays linearly rather than through
`CartesianIndex`. This avoids the `CartesianIndex` code generation that
makes many-call-site point-wise kernels (e.g. quadrature integrals)
compile slowly on the GPU. Restricted to uniform same-space, unmasked,
scalar-field arguments; provides CPU and CUDA methods.
@haakon-e haakon-e force-pushed the he/soa-broadcast branch from 27279a3 to 51f7bba Compare July 9, 2026 05:29
Comment thread ext/cuda/soa_broadcast.jl
n > 0 || return nothing
# A plain launch (no forced `always_inline`): forcing inlining of the per-point
# function balloons the IR and the LLVM compile for many-call-site functions.
kernel = CUDA.@cuda launch = false soa_broadcast_kernel!(f, outs, ins, n)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the main improvement in compilation time might come from not having always_inline. I previously found that the Cartesian indexing does slow down compilation, but not severely.

@haakon-e

haakon-e commented Jul 11, 2026

Copy link
Copy Markdown
Member Author

The compile-time decomposition promised above. The standalone benchmark scripts are at /home/haakon/soa-broadcast-bench/ on the clima cluster (the synthetic case needs only ClimaCore + CUDA; the real 2M+P3 case adds CloudMicrophysics); I can push them to a branch here or a gist if that is easier to run.

On a many-field point-wise kernel, the compile-time advantage of soa_broadcast! over the Base broadcast comes almost entirely from one choice: not forcing always_inline on the CUDA kernel. On the production 2M+P3 Rosenbrock tendency the Base broadcast does not compile within 30 minutes, while soa_broadcast! compiles it in 33 s (rosenbrock_manual) or 106 s (rosenbrock_exact). Reverting only the inlining choice reproduces the 30-minute compile. Linear indexing and map-versus-ntuple argument extraction each change compile time by at most a few seconds, and the raw-CuArray kernel ClimaAtmos uses is the same speed as soa_broadcast!. Linear indexing does give a roughly 1.6x faster warm kernel on the heavy tendency.

Compile time, full Rosenbrock 2M+P3 tendency, A100:

variant rosenbrock_manual rosenbrock_exact
soa 33 s 106 s
cartesian (CartesianIndex) 33 s 111 s
ntuple 34 s 107 s
cuarray (raw kernel) 33 s 107 s
base (Base broadcast) > 30 min (did not finish) > 30 min (did not finish)
forceinline > 30 min (did not finish) > 30 min (did not finish)
Method, variant definitions, full tables, per-component decomposition, and caveats

Method. Compile time is the first kernel launch in a fresh Julia process (CUDA compiles lazily on first launch), with a warm launch subtracted to remove device runtime. Each configuration runs in its own fresh process, so the CUDA compile cache and Julia codegen are cold. Cheap variants are repeated and the median reported with the min-max spread; the intractable base and forceinline variants run once under a 30-minute timeout. All compiling variants produce identical output checksums. A100-SXM4-80GB, Julia 1.12.5, CUDA.jl 5.11.3, ClimaCore he/soa-broadcast; single CenterFiniteDifferenceSpace column, 1000 levels.

Variants. Each reverts exactly one choice of soa_broadcast!, so a variant delta isolates that component.

variant indexing inlining arg extraction path
base CartesianIndex forced Base broadcast ClimaCore Base.copyto!
soa linear not forced map Fields.soa_broadcast!
cartesian CartesianIndex not forced map custom kernel (isolates indexing)
forceinline linear forced map custom kernel (isolates inlining)
ntuple linear not forced ntuple(_, Val(N)) custom kernel (isolates arg extraction)
cuarray linear not forced map kernel over raw CuArrays

Full Rosenbrock tendency. BMT.bulk_microphysics_tendencies over nsub linearized-implicit substeps; 12 field inputs, 9-field NamedTuple output; the mode, dt, and nsub are passed as isbits values (see the first caveat).

rosenbrock_manual (analytic 8x8 Jacobian):

variant compile (s) spread (s) warm (ms)
soa 33.1 33.1-33.5 177
cuarray 33.1 33.1-34.3 177
cartesian 32.5 32.5-32.8 177
ntuple 33.9 33.9-36.8 177
base > 30 min (did not finish) - -
forceinline > 30 min (did not finish) - -

rosenbrock_exact (ForwardDiff Jacobian):

variant compile (s) spread (s) warm (ms)
soa 105.6 105.6-105.9 426
cuarray 106.8 106.8-111.0 425
cartesian 110.7 108.3-110.7 702
ntuple 107.3 107.3-107.9 699
base > 30 min (did not finish) - -
forceinline > 30 min (did not finish) - -

On the exact kernel the warm runtime separates: cartesian and ntuple run about 1.65x slower than soa (702 and 699 ms versus 426 ms), so linear indexing and map help runtime while leaving compile time unchanged. The manual kernel shows no such difference (all 177 ms). This is a shared-node measurement and is indicative.

Instantaneous tendency (lighter reference). The instantaneous 2M+P3 tendency the Rosenbrock step wraps compiles for every variant, so it gives the full decomposition the intractable variants cannot.

variant compile (s) warm (ms)
soa 17.3 86.9
cuarray 18.8 86.8
cartesian 18.4 86.6
ntuple 19.4 87.0
base 304.6 53.4
forceinline 293.3 61.6

Per-component, relative to soa: inlining (forceinline - soa) +276 s, 96% of the base - soa gap; the N-dimensional data-layout indexing plus Base argument handling (base - forceinline) +11 s, 4%; indexing flip (cartesian - soa) +1 s; arg extraction (ntuple - soa) +2 s; raw-CuArray kernel (cuarray - soa) +2 s.

Synthetic case (ClimaCore only). A tunable point-wise function of N inputs and arithmetic depth D reproduces the result without CloudMicrophysics. With 8-field NamedTuple output at N=16, base 19.6 s versus soa 6.1 s, and inlining accounts for 95% of the difference (forceinline 19.0 s). With scalar output all variants compile within about a second of each other, so the effect is triggered by multi-field output.

Caveats.

  • The full Rosenbrock 2M+P3 tendency compiles allocation-free (rosenbrock_manual and rosenbrock_exact; zero allocation intrinsics in the device LLVM). An earlier run reporting InvalidIRError (gpu_gc_pool_alloc) on every path was an artifact of the benchmark setup: the per-point closure captured the mode, dt, and nsub as non-const globals, which forces GPU dynamic dispatch. Passing them as isbits values, as ClimaAtmos does with the mode in a concrete type parameter and n_substeps in a concrete field, compiles cleanly. Reproducers must do the same.
  • cartesian reverts indexing on the flat parent arrays; base additionally includes ClimaCore's full N-dimensional data-layout indexing, which is the residual base - forceinline difference (11 s on the instantaneous tendency).
  • The ntuple(_, Val(N)) boxing regression cited in the PR (runtime-loop fallback for N > 10 triggering gpu_gc_pool_alloc) did not reproduce on this Julia/CUDA: ntuple(_, Val(16)) compiled cleanly for both homogeneous and heterogeneous input tuples. For same-space fields with identical parent array types the fallback stays type-stable, so map versus ntuple is not a measurable compile-time factor. map remains the safer choice for heterogeneous inputs.

@haakon-e

Copy link
Copy Markdown
Member Author

I'll close this since this isn't actually being merged.

@haakon-e haakon-e closed this Jul 13, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants