Fix downstream ClimaAtmos CI failures

- Compact kernel parameters: adapt every DataLayout view to a
  CompactDeviceView at kernel launch, storing an Int32 offset per
  restricted dimension instead of a SubArray (72-88 vs 128-160 bytes per
  broadcast argument), so that large EDMF broadcasts fit the 4 KiB
  parameter limit of sm_60 GPUs; full arrays pass through unchanged.
- Combine DataScopes with pairwise recursion instead of unrolled
  helpers: broadcasts with many arguments (EDMF tendencies) make
  inference hit its recursion limiter inside UnrolledUtilities' methods,
  which the DataLayouts recursion_relation lift cannot cover, so every
  scope in the launch path widens to Any and dispatches at runtime at
  every slice. Interleaved Larcform1 measurements: median step
  allocations drop from 17118 (10.19 MB) to 3280 (1.18 MB), below
  main's 5298 (1.12 MB), with median step time 9.6 ms -> 8.2 ms.
- QuasiMonotoneLimiter: apply the limit through one scalar view per
  component of ρq, indexing DataLayouts instead of reshaping slab
  parent arrays, since Base.reshape of a device-array SubArray is
  uncompilable in kernels under Julia 1.11.
- Fix the GPU MPI DSS exchange buffer index, which overlapped items
  whenever Nv > 1 and Nf > 1, corrupting every distributed GPU
  weighted_dss! of extruded multi-component fields.
- Restore two behaviors from main that downstream code relies on:
  zero-size fields are hidden from DataLayout propertynames (restart
  comparisons recurse into Tensor bases otherwise), and array2field
  derives its parent shape from any layout (PointSpace included).
- Support the old universal CartesianIndex{5} indexing convention on
  layouts, slabs, and columns, which registered downstream packages
  like ClimaCoreTempestRemap still use.
- Dev-install lib/ClimaCoreMakie in the performance pipeline, since the
  registered version still destructures the old 5-tuple layout size.
- Remove the unused needs_projection methods for AutoBroadcasters.
- unit_loops.jl: seed the RNG, and compare single-threaded reductions
  against a hand-written pairwise reference, since the @simd blocks in
  Base's mapreduce reassociate contiguous values unless bounds checks
  are enabled.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: dennisYatunin

.buildkite/perf/pipeline.yml

@@ -26,6 +26,10 @@ steps:
       # Instantiate and dev install ClimaCore
       - julia --project=$COUPLER_PATH/experiments/AMIP/ -e 'using Pkg; Pkg.instantiate(;verbose=true)'
       - julia --project=$COUPLER_PATH/experiments/AMIP/ -e 'using Pkg; Pkg.develop(path=".")'
+      # The registered ClimaCoreMakie still destructures the old 5-tuple layout
+      # size in plot_triangles, which crashes the AMIP postprocessing plots, so
+      # dev-install the migrated in-repo version alongside ClimaCore.
+      - julia --project=$COUPLER_PATH/experiments/AMIP/ -e 'using Pkg; Pkg.develop(path="lib/ClimaCoreMakie")'
       - julia --project=$COUPLER_PATH/experiments/AMIP/ -e 'using Pkg; Pkg.add("MPI")'
       - julia --project=$COUPLER_PATH/experiments/AMIP/ -e 'using Pkg; Pkg.precompile()'
       - julia --project=$COUPLER_PATH/experiments/AMIP/ -e 'using Pkg; Pkg.status()'

ext/cuda/data_layouts.jl

@@ -6,3 +6,122 @@ import ClimaCore: DataLayouts
 include("scopes.jl")
 include("loops.jl")
 include("data_layouts_threadblock.jl")
+
+# Kernel parameters are limited to 4 KiB of memory before compute capability
+# 7.0, and a SubArray of a 5-D CuDeviceArray uses 128-160 of those bytes per
+# broadcast argument (64 for the parent array, 48-80 for the index ranges, and
+# 16 for precomputed linear-indexing fields), so broadcasts over a few dozen
+# field views cannot be launched as kernels. Since every extent of the array in
+# a DataLayout is either available from the layout's type or identical to the
+# corresponding extent of the parent array, the index ranges can be replaced
+# with an Int32 offset for every restricted dimension, plus an Int32 extent for
+# every restricted dimension whose extent is not available from the type. The
+# type parameters are the extent of each dimension `E` (with 0 for extents that
+# are only available at runtime), the restricted dimensions `R`, and the tuple
+# lengths `K = length(R)` and `D = count(d -> iszero(E[d]), R)`.
+struct CompactDeviceView{T, N, E, R, K, D, A <: AbstractArray{T, N}} <:
+       AbstractArray{T, N}
+    parent::A
+    offsets::NTuple{K, Int32}
+    dynamic_extents::NTuple{D, Int32}
+end
+
+Base.parent(view::CompactDeviceView) = view.parent
+
+DataLayouts.DataScope(
+    ::Type{<:CompactDeviceView{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, A}},
+) where {A} = DataLayouts.DataScope(A)
+
+@inline function Base.size(
+    view::CompactDeviceView{<:Any, N, E, R},
+) where {N, E, R}
+    return ntuple(Val(N)) do d
+        iszero(E[d]) || return Int(E[d])
+        d in R || return size(parent(view), d)
+        return Int(view.dynamic_extents[count(r -> r <= d && iszero(E[r]), R)])
+    end
+end
+
+# Index into the parent array that corresponds to an index into the view,
+# shifted by the stored offset along each restricted dimension
+@inline parent_index(
+    view::CompactDeviceView{<:Any, N, <:Any, R},
+    index,
+) where {N, R} =
+    ntuple(Val(N)) do d
+        position = findfirst(==(d), R)
+        isnothing(position) ? index[d] : index[d] + Int(view.offsets[position])
+    end
+
+Base.@propagate_inbounds function Base.getindex(
+    view::CompactDeviceView{<:Any, N},
+    index::Vararg{Integer, N},
+) where {N}
+    @boundscheck checkbounds(view, index...)
+    return @inbounds parent(view)[parent_index(view, index)...]
+end
+
+Base.@propagate_inbounds function Base.setindex!(
+    view::CompactDeviceView{<:Any, N},
+    value,
+    index::Vararg{Integer, N},
+) where {N}
+    @boundscheck checkbounds(view, index...)
+    @inbounds parent(view)[parent_index(view, index)...] = value
+    return view
+end
+
+# Index types generated by stable_view for unrestricted and restricted
+# dimensions of the parent array in a DataLayout
+const ViewDimIndex = Union{
+    Base.Slice{<:Base.OneTo{<:Integer}},
+    Base.OneTo{<:Integer},
+    UnitRange{<:Integer},
+}
+
+# Extent of every parent array dimension that is available from a DataLayout's
+# type, with 0 for dimensions whose extents are only available at runtime
+@inline inferred_parent_extents(data, array) = DataLayouts.add_f_dim(
+    map(extent -> something(extent, 0), DataLayouts.inferred_size(data)),
+    DataLayouts.num_basetypes(eltype(array), eltype(data)),
+    Val(DataLayouts.f_dim(data)),
+)
+
+compact_device_view(array, data) = array
+
+@inline function compact_device_view(
+    array::SubArray{<:Any, N, <:CUDA.CuDeviceArray, <:NTuple{N, ViewDimIndex}},
+    data::DataLayouts.DataLayout,
+) where {N}
+    extents = inferred_parent_extents(data, array)
+    length(extents) == N ||
+        throw(DimensionMismatch("DataLayout extents do not match its array"))
+    indices = parentindices(array)
+    restricted =
+        filter(d -> indices[d] isa UnitRange, ntuple(identity, Val(N)))
+    dynamic = filter(d -> iszero(extents[d]), restricted)
+    view = CompactDeviceView{
+        eltype(array),
+        N,
+        extents,
+        restricted,
+        length(restricted),
+        length(dynamic),
+        typeof(parent(array)),
+    }(
+        parent(array),
+        map(d -> Int32(first(indices[d]) - 1), restricted),
+        map(d -> Int32(length(indices[d])), dynamic),
+    )
+    # Validates every extent assumption at launch time, including that any
+    # Base.OneTo indices span their full parent dimensions
+    size(view) == size(array) ||
+        throw(DimensionMismatch("DataLayout extents do not match its array"))
+    return view
+end
+
+Adapt.adapt_structure(to::CUDA.KernelAdaptor, data::DataLayouts.DataLayout) =
+    DataLayouts.rebuild(
+        data,
+        compact_device_view(Adapt.adapt(to, parent(data)), data),
+    )

ext/cuda/topologies_dss.jl

@@ -252,7 +252,7 @@ function Topologies.fill_send_buffer!(
             # Avoid indexing with a colon, which would allocate in the kernel.
             (h, p) = (send_buf_idx[send_index, 1], send_buf_idx[send_index, 2])
             item = perimeter_data[v, p, 1, h]
-            buffer_index = v + (send_index - 1) * Nv * Nf
+            buffer_index = (v - 1) * Nf + 1 + (send_index - 1) * Nv * Nf
             DataLayouts.set_struct!(send_data, item, buffer_index, Val(1))
         end
         nothing
@@ -279,7 +279,7 @@ function Topologies.load_from_recv_buffer!(
             (v, recv_index) = CartesianIndices((Nv, nrecv))[gidx].I
             # Avoid indexing with a colon, which would allocate in the kernel.
             (h, p) = (recv_buf_idx[recv_index, 1], recv_buf_idx[recv_index, 2])
-            buffer_index = v + (recv_index - 1) * Nv * Nf
+            buffer_index = (v - 1) * Nf + 1 + (recv_index - 1) * Nv * Nf
             item_view = DataLayouts.view_struct(recv_data, T, buffer_index, Val(1))
             parent_view = parent(view(perimeter_data, v, p, 1, h))
             for f in 1:Nf

src/DataLayouts/DataLayouts.jl

@@ -274,7 +274,13 @@ function replace_basetype(data::DataLayout, ::Type{B}) where {B}
     return similar_layout(data, T, B)
 end
 
-@inline Base.propertynames(data::DataLayout) = fieldnames(eltype(data))
+# Hide zero-size fields (e.g. the singleton bases of Tensors) from property
+# iteration, so that generic code which recursively walks propertynames only
+# encounters properties that contain data. Zero-size fields are still
+# accessible through getproperty, which returns a view with an empty F axis.
+@inline function Base.propertynames(::DataLayout{T}) where {T}
+    filter(name -> sizeof(fieldtype(T, name)) > 0, fieldnames(T))
+end
 
 # Wrap the field index in a Val as soon as it is available, resolving field
 # views through specialization rather than constant propagation. Making the Val
@@ -532,10 +538,11 @@ include("deprecated.jl")
 
 # Drop the recursion limits of this module's Core.kwcall methods and recursive
 # DataScope functions, so that kwarg functions like fill! and column_reduce! can
-# be composed and slice_subscope can repeatedly partition a scope. The default
-# limit makes the compiler widen argument types, leading to dynamic dispatch.
+# be composed, multiple scopes can be combined, and is_subscope/slice_subscope
+# can repeatedly partition a scope. The default limit makes the compiler widen
+# argument types, leading to dynamic dispatch and runtime allocations.
 @static if hasfield(Method, :recursion_relation)
-    for f in (Core.kwcall, slice_subscope, is_subscope), method in methods(f)
+    for f in (Core.kwcall, DataScope, is_subscope, slice_subscope), method in methods(f)
         method.module === (@__MODULE__) || continue
         method.recursion_relation = Returns(true)
     end

src/DataLayouts/deprecated.jl

@@ -31,3 +31,26 @@ const IJHF{S, Nij} = VIJHF{S, 1, Nij, Nij, nothing}
 # These were exported on main, and packages like ClimaCoreTempestRemap access
 # them through `using ClimaCore.DataLayouts`.
 export AbstractData, IJFH, IJHF
+
+# Old: getindex/setindex! with a 5-D "universal index" CartesianIndex(i, j, f,
+# v, h), whose f component was ignored because the entire value was read or
+# written. New: Cartesian indices have one coordinate per array dimension, in
+# the order (v, i, j, h). Packages like ClimaCoreTempestRemap still index full
+# layouts, slabs, and columns with universal indices (e.g.,
+# `slab(data, h)[CartesianIndex(i, j, 1, 1, 1)]`), so the methods below
+# translate 5-D indices from the old convention to the new one. They override
+# the generic AbstractArray interpretation of a `CartesianIndex{5}` as a 4-D
+# index with a trailing singleton coordinate, so new code must not add trailing
+# coordinates when indexing into 4-dimensional layouts. Remove these methods
+# together with the aliases above.
+@inline universal_index_shim(I::CartesianIndex{5}) =
+    CartesianIndex(I[4], I[1], I[2], I[5])
+@propagate_inbounds Base.getindex(data::VIJHWithF, I::CartesianIndex{5}) =
+    getindex(data, universal_index_shim(I))
+@propagate_inbounds Base.setindex!(data::VIJHWithF, value, I::CartesianIndex{5}) =
+    setindex!(data, value, universal_index_shim(I))
+
+# Old: every component of a universal index was ignored for 0-dimensional data.
+@propagate_inbounds Base.getindex(data::DataF, I::CartesianIndex{5}) = data[]
+@propagate_inbounds Base.setindex!(data::DataF, value, I::CartesianIndex{5}) =
+    setindex!(data, value)

src/DataLayouts/scopes.jl

@@ -46,15 +46,15 @@ DataScope(::Type{<:Base.ReshapedArray{<:Any, <:Any, A}}) where {A} = DataScope(A
 # Infer parent types of other AbstractArrays (constant-folding not guaranteed).
 DataScope(::Type{A}) where {A <: AbstractArray} = DataScope(return_type(parent, Tuple{A}))
 
-DataScope(arg1, arg2, args...) =
-    unrolled_reduce(unrolled_map(DataScope, (arg1, arg2, args...))) do scope1, scope2
-        is_subscope(scope1, scope2) ? scope1 :
-        is_subscope(scope2, scope1) ? scope2 :
-        throw(ArgumentError(non_overlapping_scopes_string(scope1, scope2)))
-    end
+DataScope(scope1::DataScope, scope2::DataScope) =
+    is_subscope(scope1, scope2) ? scope1 :
+    is_subscope(scope2, scope1) ? scope2 :
+    throw(ArgumentError(non_overlapping_scopes_string(scope1, scope2)))
 @generated non_overlapping_scopes_string(::S1, ::S2) where {S1, S2} =
     "$S1 and $S2 do not overlap, so they cannot be put in the same DataScope"
 
+DataScope(arg1, arg2, args...) = DataScope(DataScope(arg1), DataScope(arg2, args...))
+
 """
     partition(scope)
 
@@ -194,11 +194,11 @@ struct ThisThreadPool <: DataScope end
 # _sym_to_tpid (0 = :interactive, 1 = :default); fall back to the public API if the
 # internals change.
 @static if isdefined(Threads, :_nthreads_in_pool)
-    @inline default_pool_size() = Int(Threads._nthreads_in_pool(Int8(1)))
-    @inline interactive_pool_size() = Int(Threads._nthreads_in_pool(Int8(0)))
+    default_pool_size() = Int(Threads._nthreads_in_pool(Int8(1)))
+    interactive_pool_size() = Int(Threads._nthreads_in_pool(Int8(0)))
 else
-    @inline default_pool_size() = Threads.threadpoolsize(:default)
-    @inline interactive_pool_size() = Threads.threadpoolsize(:interactive)
+    default_pool_size() = Threads.threadpoolsize(:default)
+    interactive_pool_size() = Threads.threadpoolsize(:interactive)
 end
 
 # Threads.threading_run compiles faster than an equivalent static Threads.@threads loop
@@ -219,8 +219,8 @@ end
 
 # Task-local storage marks ClimaCore-launched threads, distinguishing them from external
 # threaded loops; storage is nothing until first set, so storage-less threads are external.
-@inline running_in_threaded_loop() = !iszero(ccall(:jl_in_threaded_region, Cint, ()))
-@inline function running_in_external_threaded_loop()
+running_in_threaded_loop() = !iszero(ccall(:jl_in_threaded_region, Cint, ()))
+function running_in_external_threaded_loop()
     running_in_threaded_loop() || return false
     storage = current_task().storage
     return isnothing(storage) ||

src/Fields/Fields.jl

@@ -584,10 +584,7 @@ The struct type of the resulting `Field` is set to the array's element type.
 """
 function array2field(array, space)
     data = Spaces.local_geometry_data(space)
-    (; Nv, Ni, Nj, Nh, F) = DataLayouts.shape_params(data)
-    Nh_dynamic = isnothing(Nh) ? DataLayouts.nelems(data) : Nh
-    array_size =
-        DataLayouts.add_f_dim((Nv, Ni, Nj, Nh_dynamic), 1, Val(F))
+    array_size = DataLayouts.add_f_dim(size(data), 1, Val(DataLayouts.f_dim(data)))
     parent_array = reshape(array, array_size)
     return Field(DataLayouts.rebuild(data, parent_array, eltype(array)), space)
 end

src/Fields/field_iterator.jl

@@ -38,8 +38,6 @@ function flattened_property_chains!(prop_chains, f::Field, pc = ())
     else
         for pn in propertynames(f)
             p = getproperty(f, pn)
-            # Skip properties of singleton types, which do not contain data
-            sizeof(eltype(p)) == 0 && continue
             flattened_property_chains!(prop_chains, p, (pc..., pn))
         end
     end

src/Geometry/auto_broadcaster_methods.jl

@@ -20,16 +20,6 @@ mul_with_projection(x::AutoBroadcaster, y, lg) =
 mul_with_projection(x, y::AutoBroadcaster, lg) =
     nested_broadcast(y -> mul_with_projection(x, y, lg), y)
 
-needs_projection(
-    ::Type{X},
-    ::Type{Y},
-) where {X <: AutoBroadcaster, Y <: AutoBroadcaster} =
-    needs_projection(eltype(X), eltype(Y))
-needs_projection(::Type{X}, ::Type{Y}) where {X <: AutoBroadcaster, Y} =
-    needs_projection(eltype(X), Y)
-needs_projection(::Type{X}, ::Type{Y}) where {X, Y <: AutoBroadcaster} =
-    needs_projection(X, eltype(Y))
-
 mul_return_type(
     ::Type{X},
     ::Type{Y},

src/Limiters/Limiters.jl

@@ -1,6 +1,7 @@
 module Limiters
 
 import ..DataLayouts, ..Topologies, ..Spaces, ..Fields
+using UnrolledUtilities
 import ..DebugOnly: call_post_op_callback, post_op_callback
 import ClimaCore: slab