CliMA · dennisYatunin · May 1, 2026 · Jan 10, 2026
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
@@ -57,13 +57,6 @@ steps:
         agents:
           slurm_gpus: 1
 
-  - group: "Unit: RecursiveApply"
-    steps:
-
-      - label: "Unit: RecursiveApply"
-        key: unit_recursive_apply
-        command: "julia --color=yes --check-bounds=yes --project=.buildkite test/RecursiveApply/unit_recursive_apply.jl"
-
   - group: "Unit: CUDA utils"
     steps:
 
@@ -93,6 +86,10 @@ steps:
         key: unit_plushalf
         command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Utilities/unit_plushalf.jl"
 
+      - label: "Unit: auto_broadcaster"
+        key: unit_auto_broadcaster
+        command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Utilities/unit_auto_broadcaster.jl"
+
   - group: "Unit: DataLayouts"
     steps:
 
@@ -227,9 +224,9 @@ steps:
         key: unit_axistensors
         command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Geometry/axistensors.jl"
 
-      - label: "Unit: rmul_with_projection"
-        key: unit_rmul_with_projection
-        command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Geometry/rmul_with_projection.jl"
+      - label: "Unit: mul_with_projection"
+        key: unit_mul_with_projection
+        command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Geometry/mul_with_projection.jl"
 
   - group: "Unit: Meshes"
     steps:

diff --git a/benchmarks/3d/se_kernels.jl b/benchmarks/3d/se_kernels.jl
@@ -13,9 +13,7 @@ import ClimaCore:
     Spaces,
     Quadratures,
     Topologies,
-    DataLayouts,
-    RecursiveApply
-
+    DataLayouts
 const C1 = ClimaCore.Geometry.Covariant1Vector
 const C2 = ClimaCore.Geometry.Covariant2Vector
 const C3 = ClimaCore.Geometry.Covariant3Vector
@@ -25,8 +23,6 @@ const CT123 = Geometry.Contravariant123Vector
 const ᶜinterp = Operators.InterpolateF2C()
 const ᶠinterp = Operators.InterpolateC2F()
 
-const ⊞ = RecursiveApply.radd
-
 init_uθ(ϕ, z, R) = 1.0 / R
 init_vθ(ϕ, z, R) = 1.0 / R
 init_w(ϕ, z) = 1.0

diff --git a/benchmarks/bickleyjet/bickleyjet_dg.jl b/benchmarks/bickleyjet/bickleyjet_dg.jl
@@ -7,9 +7,6 @@ import ClimaCore.Operators
 using ClimaCore.Geometry
 import ClimaCore.Geometry: Abstract2DPoint
 
-using ClimaCore.RecursiveApply
-
-
 const parameters = (
     ϵ = 0.1,  # perturbation size for initial condition
     l = 0.5, # Gaussian width
@@ -60,7 +57,7 @@ roe_average(ρ⁻, ρ⁺, var⁻, var⁺) =
     (sqrt(ρ⁻) * var⁻ + sqrt(ρ⁺) * var⁺) / (sqrt(ρ⁻) + sqrt(ρ⁺))
 
 function roeflux(n, (y⁻, parameters⁻), (y⁺, parameters⁺))
-    Favg = RecursiveApply.rdiv(flux(y⁻, parameters⁻) ⊞ flux(y⁺, parameters⁺), 2)
+    Favg = (flux(y⁻, parameters⁻) + flux(y⁺, parameters⁺)) / 2
 
     λ = sqrt(parameters⁻.g)
 
@@ -115,7 +112,7 @@ function roeflux(n, (y⁻, parameters⁻), (y⁺, parameters⁺))
     fluxᵀn_ρθ = ((w1 + w2) * θ + w5) * 0.5
 
     Δf = (ρ = -fluxᵀn_ρ, ρu = -fluxᵀn_ρu, ρθ = -fluxᵀn_ρθ)
-    RecursiveApply.rmap(f -> f' * n, Favg) ⊞ Δf
+    return Favg' * n + Δf
 end
 
 function volume!(dydt, y, (parameters,), t)

diff --git a/benchmarks/bickleyjet/core_vs_ref.jl b/benchmarks/bickleyjet/core_vs_ref.jl
@@ -30,8 +30,7 @@ for Nq in Nqs
     volume!(dydt, y0, (parameters,), 0.0)
     # TODO: move this to volume!
     dydt_data = Fields.field_values(dydt)
-    dydt_data .=
-        RecursiveApply.rdiv.(dydt_data, Spaces.local_geometry_data(space).WJ)
+    dydt_data ./= Spaces.local_geometry_data(space).WJ
 
     # setup reference
     X = coordinates(Val(Nq), n1, n2)
@@ -84,8 +83,7 @@ for Nq in Nqs
     add_face!(dydt, y0, (parameters,), 0.0)
     # TODO: move this to volume!
     dydt_data = Fields.field_values(dydt)
-    dydt_data .=
-        RecursiveApply.rdiv.(dydt_data, Spaces.local_geometry_data(space).WJ)
+    dydt_data ./= Spaces.local_geometry_data(space).WJ
 
     fill!(dydt_ref, 0.0)
     add_face_ref!(dydt_ref, y0_ref, (n1, n2, parameters, Val(Nq)), 0.0)

diff --git a/docs/make.jl b/docs/make.jl
@@ -100,7 +100,6 @@ withenv("GKSwstype" => "nul") do
                 "Limiters" => "APIs/limiters_api.md",
                 "InputOutput" => "APIs/input_output_api.md",
                 "Remapping" => "APIs/remapping_api.md",
-                "RecursiveApply" => "APIs/recursive_apply_api.md",
                 "Devices" => "APIs/devices_api.md",
                 "DebugOnly" => "APIs/debug_only_api.md",
             ],

diff --git a/docs/src/APIs/geometry_api.md b/docs/src/APIs/geometry_api.md
@@ -21,6 +21,7 @@ Geometry.LocalGeometry
 
 ```@docs
 Geometry.Δz_metric_component
+Geometry.:⊗
 ```
 
 ## Coordinates

diff --git a/docs/src/APIs/recursive_apply_api.md b/docs/src/APIs/recursive_apply_api.md
diff --git a/docs/src/APIs/utilities_api.md b/docs/src/APIs/utilities_api.md
@@ -5,11 +5,30 @@ CurrentModule = ClimaCore
 ```
 
 ```@docs
-Utilities.PlusHalf
-Utilities.half
+Utilities.unionall_type
 Utilities.replace_type_parameter
 Utilities.fieldtype_vals
 Utilities.new
+Utilities.unsafe_eltype
+Utilities.safe_eltype
+```
+
+## Utilities.PlusHalf
+
+```@docs
+Utilities.PlusHalf
+Utilities.half
+```
+
+## Utilities.AutoBroadcaster
+
+```@docs
+Utilities.AutoBroadcaster
+Utilities.is_auto_broadcastable
+Utilities.add_auto_broadcasters
+Utilities.drop_auto_broadcasters
+Utilities.auto_broadcasted
+Utilities.nested_broadcast
 ```
 
 ## Utilities.Cache

diff --git a/docs/src/geometry.md b/docs/src/geometry.md
@@ -6,7 +6,5 @@ CurrentModule = ClimaCore.Geometry
 
 ```@docs
 mul_with_projection
-rmul_with_projection
 mul_return_type
-rmul_return_type
 ```
diff --git a/docs/src/matrix_fields.md b/docs/src/matrix_fields.md
@@ -210,15 +210,11 @@ J = MatrixFields.FieldMatrix((@name(f), @name(g))=> ∂f_∂g)
 
 ## Optimizations
 
-Each entry of a `FieldMatrix` can be a `ColumnwiseBandMatrixField`, a `DiagonalMatrixRow`, or an
-`UniformScaling`.
-
-A `ColumnwiseBandMatrixField` is a `Field` with a `BandMatrixRow` at each point. It is intended
-to represent a collection of banded matrices, where there is one band matrix for each column
-of the space the `Field` is on. Beyond only storing the diagonals of the band matrix, an `entry`
-can be optimized to use less memory. Each optimized representation can be indexed equivalently to
-non optimized representations, and used in addition, subtraction, matrix-vector multiplication,
-Matrix-matrix multiplication, `RecursiveApply`, and `FieldMatrixSolver`.
+Each entry of a `FieldMatrix` can be a `ColumnwiseBandMatrixField`, a `DiagonalMatrixRow`, or a
+`UniformScaling`. A `ColumnwiseBandMatrixField` is a `Field` with a `BandMatrixRow` at each point.
+It represents a collection of banded matrices, with each column of the `Field` corresponding to a
+specific matrix. If all columns correspond a constant multiple of the identity matrix, the `Field`
+may be replaced with a `ScalingFieldMatrixEntry` (i.e., a `DiagonalMatrixRow` or `UniformScaling`).
 
 For the following sections, `space` is a column space with $N_v$ levels. A column space is
 used for simplicity in this example, but the optimizations work with any space with columns.

diff --git a/docs/src/operators.md b/docs/src/operators.md
@@ -112,8 +112,6 @@ column_accumulate!
 ## Internal APIs
 
 ```@docs
-getidx_return_type
-stencil_return_type
 return_eltype
 return_space
 stencil_interior_width

diff --git a/examples/bickleyjet/bickleyjet_dg.jl b/examples/bickleyjet/bickleyjet_dg.jl
@@ -7,12 +7,10 @@ import ClimaCore:
     Geometry,
     Meshes,
     Operators,
-    RecursiveApply,
     Spaces,
     Quadratures,
     Topologies
 import ClimaCore.Geometry: ⊗
-import ClimaCore.RecursiveApply: ⊞, rdiv, rmap
 
 using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33
 
@@ -106,8 +104,6 @@ roe_average(ρ⁻, ρ⁺, var⁻, var⁺) =
     (sqrt(ρ⁻) * var⁻ + sqrt(ρ⁺) * var⁺) / (sqrt(ρ⁻) + sqrt(ρ⁺))
 
 function roeflux(n, (y⁻, parameters⁻), (y⁺, parameters⁺))
-    Favg = rdiv(flux(y⁻, parameters⁻) ⊞ flux(y⁺, parameters⁺), 2)
-
     λ = sqrt(parameters⁻.g)
 
     ρ⁻, ρu⁻, ρθ⁻ = y⁻.ρ, y⁻.ρu, y⁻.ρθ
@@ -159,9 +155,11 @@ function roeflux(n, (y⁻, parameters⁻), (y⁺, parameters⁺))
         (w1 * (u - c * n) + w2 * (u + c * n) + w3 * u + w4 * (Δu - Δuₙ * n)) *
         0.5
     fluxᵀn_ρθ = ((w1 + w2) * θ + w5) * 0.5
-
     Δf = (ρ = -fluxᵀn_ρ, ρu = -fluxᵀn_ρu, ρθ = -fluxᵀn_ρθ)
-    rmap(f -> f' * n, Favg) ⊞ Δf
+
+    return map(flux(y⁻, parameters⁻), flux(y⁺, parameters⁺), Δf) do F⁻, F⁺, Δf
+        ((F⁻ + F⁺) / 2)' * n + Δf
+    end
 end
 
 
@@ -205,9 +203,8 @@ function rhs!(dydt, y, (parameters, numflux), t)
     end
 
     # 6. Solve for final result
-    dydt_data = Fields.field_values(dydt)
-    dydt_data .=
-        RecursiveApply.rdiv.(dydt_data, Spaces.local_geometry_data(space).WJ)
+    dydt_data =
+        Fields.field_values(dydt) ./ Spaces.local_geometry_data(space).WJ
     M = Quadratures.cutoff_filter_matrix(
         Float64,
         Spaces.quadrature_style(space),