Fix comments and docstrings that no longer match the code

Comments only; no behaviour change.

Stale claims corrected:

- ext/cuda/operators_finite_difference.jl: the launch-config comment said the
  kernel "uses block and grid indices instead of computing cartesian indices
  from a linear index" (it does the opposite now), that the block size is
  `(n_face_levels, Ni, 1)` (it is `(n_face_levels, threads_dim_y, 1)`), and that
  the path is used only when a block holds 32-256 threads (that guard is gone).
  Rewrote to describe the current launch configuration, noted that masked spaces
  are supported too, explained the `threads_dim_y` heuristic, and dropped the
  dead commented-out `32 <= n_face_levels * Ni <= 256` line.

- Both files referred to a `us::UniversalSize` argument that does not exist. The
  kernel reads the horizontal extents off the output layout's type parameters
  via `vijh_params`; the "compile-time constant divisor" point still holds, so
  it is now stated in terms of what the code actually does.

- ext/cuda/operators_fd_eager.jl: "A `SetBoundaryOperator` only ever outputs to
  faces" is false since `return_space(::SetBoundaryOperator, space) = space` --
  the conversion also puts one on a center output (DivergenceF2C +
  SetDivergence) and on a face input (GradientF2C + SetValue). The reason for
  gating `idx` on the compile-time staggering is unchanged and kept.

- src/MatrixFields/operator_matrices.jl: "the boundary row is empty (`rzero`)".
  `rzero` returns a zero-filled row of the same bandwidth, not an empty row;
  it is the multiply's band clipping that makes the wider row harmless.

- src/Operators/finitedifference.jl: six comments described inlined stencils as
  belonging to "(deleted)" operators. `LeftBiasedC2F`/`RightBiasedC2F` are still
  exported operators -- only their `stencil_interior` methods went away -- so
  that wording was misleading; the 3rd-order variants really were removed, which
  is now said explicitly.

Outdated signatures and documentation:

- All twelve docstring signatures in operators_fd_eager.jl were missing the
  `hidx` argument added to `calc_level_val`/`get_op_row`/`project_row2_for_mul`,
  or the `mask` argument added to `eager_copyto_stencil_kernel!`.

- `SetBoundaryOperator`'s docstring listed only `SetValue`, though it accepts
  `SetGradient`, `SetCurl` and `SetDivergence` and projects the first two onto
  the `Covariant3`/`Contravariant123` axes. All four are now documented with
  their projection semantics, along with the space-preserving behaviour and the
  fact that a side without a condition is left untouched. The four `@ref`
  targets are all listed in docs/src/operators.md, so the docs build resolves
  them.

Verified: docstring parses and renders, and CPU MatrixFields/operator_matrices.jl
(275+41) and finitedifference/unit_column.jl still pass.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01BLMKMmk7JezBvagh3AybPQ

Author: imreddyTeja

ext/cuda/operators_fd_eager.jl

@@ -104,10 +104,10 @@ compile-time constant.
     space.staggering isa Spaces.CellCenter && !Topologies.isperiodic(space)
 
 """
-    eager_copyto_stencil_kernel!(out, bc::BC, space)
+    eager_copyto_stencil_kernel!(out, bc::BC, mask, space)
 
 CUDA kernel to compute the value of a `Broadcasted` or `StencilBroadcasted` at a single index.
-This calls `calc_level_val(bc, space)`, which  computes the value of the broadcasted
+This calls `calc_level_val(bc, hidx, space)`, which computes the value of the broadcasted
 expression at the given index, and then copies the result into `out`.
 """
 Base.@propagate_inbounds function eager_copyto_stencil_kernel!(
@@ -119,8 +119,10 @@ Base.@propagate_inbounds function eager_copyto_stencil_kernel!(
     v = threadIdx().x
     col_idx = threadIdx().y + (blockIdx().x - 1) * blockDim().y
     (i, j, h) = if mask isa NoMask
-        # `Nij` comes from the type of `us`, so it is a compile-time constant and
-        # the `CartesianIndices` decomposition below uses a fixed-divisor `divrem`.
+        # `Ni` and `Nj` are read off the output layout's type parameters (see
+        # `vijh_params`), so they are compile-time constants and the `CartesianIndices`
+        # decomposition below is a fixed-divisor `divrem`. Only `Nh` is a runtime value,
+        # and being the last extent it is never divided by.
         size_params = ClimaCore.DataLayouts.vijh_params(ClimaCore.Fields.field_values(out))
         Nj = size_params.Nj
         Ni = size_params.Ni
@@ -149,7 +151,7 @@ end
 # All the functions below this line should not be used outside of this file
 
 """
-    calc_level_val(bc, space)
+    calc_level_val(bc, hidx, space)
 
 Call `calc_level_val` on all the arguments of `bc`, and then apply the function `bc.f` to the results.
 """
@@ -166,7 +168,7 @@ Base.@propagate_inbounds function calc_level_val(
 end
 
 """
-    reconstruct_space_and_call_calc_level_val(arg, space)
+    reconstruct_space_and_call_calc_level_val(arg, (hidx, space))
 
 If `arg` is a `Broadcasted`, `StencilBroadcasted`, or `Field`,
 reconstruct the space for the argument and call `calc_level_val` on it. This allows
@@ -189,7 +191,7 @@ Base.@propagate_inbounds reconstruct_space_and_call_calc_level_val(
 ) where {A, S} = @inbounds @inline calc_level_val(arg, space_idx_tpl[1], space_idx_tpl[2])
 
 """
-    calc_level_val(val::T, space)
+    calc_level_val(val::T, hidx, space)
 
 If `val` is not a `Broadcasted`, `StencilBroadcasted`, or `Field`, just return `val`.
 If it is a `Ref`, return `val[]`. If it is a one element tuple, return the element.
@@ -200,7 +202,7 @@ Base.@propagate_inbounds calc_level_val(val::T, hidx, space) where {V, T <: Tupl
 Base.@propagate_inbounds calc_level_val(arg::S, hidx, space) where {S} = arg
 
 """
-    calc_level_val(bc::StencilBroadcasted{<:Any, <: MultiplyColumnwiseBandMatrixField}, space)
+    calc_level_val(bc::StencilBroadcasted{<:Any, <: MultiplyColumnwiseBandMatrixField}, hidx, space)
 
 Call `calc_level_val` on both args of `bc`, place the result of the second arg into shared memory,
 and then perform the multiplication.
@@ -284,12 +286,12 @@ Base.@propagate_inbounds function calc_level_val(
 end
 
 """
-    calc_level_val(bc::StencilBroadcasted{<:Any, <: SetBoundaryOperator}, space)
+    calc_level_val(bc::StencilBroadcasted{<:Any, <: SetBoundaryOperator}, hidx, space)
 
-A `SetBoundaryOperator` only modifies the two boundary faces, and is the identity
-in the interior. At the boundaries we dispatch to `stencil_left_boundary` /
-`stencil_right_boundary` (which extract and project the boundary value), and in the
-interior we reuse the eagerly-computed value of the argument.
+A `SetBoundaryOperator` only modifies the two boundary levels of the space it is applied
+to, and is the identity in the interior. At the boundaries we dispatch to
+`stencil_left_boundary` / `stencil_right_boundary` (which extract and project the
+boundary value), and in the interior we reuse the eagerly-computed value of the argument.
 """
 Base.@propagate_inbounds function calc_level_val(
     bc::BC,
@@ -303,13 +305,14 @@ Base.@propagate_inbounds function calc_level_val(
     op = bc.op
     v = threadIdx().x
     val_no_bcs = @inline @inbounds calc_level_val(bc.args[1i32], hidx, space)
-    # A `SetBoundaryOperator` only ever outputs to faces. Gating the boundary logic on
-    # the (compile-time) staggering type ensures that when this method is compiled for a
-    # center-output space, the whole face-boundary branch -- including
-    # `should_call_left_boundary`, whose `idx < left_interior_idx` comparison would
-    # otherwise mix a `PlusHalf` face index with an integer center index and pull in
-    # non-GPU-compatible error-formatting code -- is dropped as dead code.
-
+    # A `SetBoundaryOperator` is space-preserving (`return_space(op, space) = space`), so
+    # this method is compiled for both staggerings: the automatic conversion puts one on a
+    # face output (InterpolateC2F + SetValue), on a center output (DivergenceF2C +
+    # SetDivergence), and on a face input (GradientF2C + SetValue). Deriving `idx` from
+    # the compile-time staggering type keeps the two apart, so the `PlusHalf` face index
+    # only reaches `should_call_*_boundary` when compiling for a face space and its
+    # `idx < left_interior_idx` comparison never mixes a `PlusHalf` with an integer center
+    # index -- which would pull in non-GPU-compatible error-formatting code.
     idx = space.staggering isa Spaces.CellFace ? (v - half) : v
     if Operators.should_call_left_boundary(idx, space, op, bc.args...)
         lbw = Operators.left_boundary_window(space)
@@ -338,7 +341,7 @@ Base.@propagate_inbounds function calc_level_val(
 end
 
 """
-    calc_level_val(bc::StencilBroadcasted{<:Any, <: LinVanLeerC2F}, space)
+    calc_level_val(bc::StencilBroadcasted{<:Any, <: LinVanLeerC2F}, hidx, space)
 
 Special case of `calc_level_val` for `LinVanLeerC2F`s, which makes the
 top and bottom face values not use the fallback `Operators.getidx`, since that
@@ -358,7 +361,7 @@ Base.@propagate_inbounds function calc_level_val(
 end
 
 """
-    calc_level_val(bc::StencilBroadcasted, space)
+    calc_level_val(bc::StencilBroadcasted, hidx, space)
 
 Fallback case of `calc_level_val` that calls `Operators.getidx`. This is used for
 affine BCs or values that won't fit in shmmem.
@@ -378,7 +381,7 @@ Base.@propagate_inbounds function calc_level_val(
 end
 
 """
-    calc_level_val(f::Field, space)
+    calc_level_val(arg::Field, hidx, space)
 
 Returns the value of the field `f` at the thread's index.
 When the staggering of `space` is `CellCenter`, the thread with `v == CUDA.blockDim().x` returns `new(eltype(f))`
@@ -400,7 +403,7 @@ Base.@propagate_inbounds function calc_level_val(
 end
 
 """
-    calc_level_val(bc::StencilBroadcasted{<:Any, <: FDOperatorMatrix}, space)
+    calc_level_val(bc::StencilBroadcasted{<:Any, <: FDOperatorMatrix}, hidx, space)
 
 Return the correct row of the operator matrix for the current thread
 """
@@ -420,7 +423,7 @@ Base.@propagate_inbounds function calc_level_val(
 end
 
 """
-    get_op_row(op, args, space)
+    get_op_row(op, args, hidx, space)
 
 Get the correct row of the operator matrix for the current thread, taking into account boundary conditions.
 """
@@ -476,7 +479,7 @@ end
 
 
 """
-    project_row2_for_mul
+    project_row2_for_mul(mat1_row, mat2_row, hidx, space)
 
 Projects `mat2_row` onto the correct axis for multiplication with `mat1_row` if necessary, and returns the projected row.
 """

ext/cuda/operators_finite_difference.jl

@@ -77,28 +77,29 @@ function Base.copyto!(
     else
         bc′ = disable_shmem_style(bc)
         (_, Ni, Nj, Nh) = size(out_fv)
-        # This uses block and grid indices instead of computing cartesian indices from a
-        # linear index. The launch configuration is optimized for common use case of 64 face
-        # levels and Ni = Nj = 4. Both periodic and non-periodic vertical topologies are
-        # supported (`has_padding_thread` accounts for the extra face level of non-periodic
-        # spaces); high-resolution columns fall through to `copyto_stencil_kernel!` below.
-        # `eager_copyto_stencil_kernel!` requires a  block size of (n_face_levels, Ni, 1)
-        # this block config is better for VIJFH. It is only used when the total number of
-        # threads in a block is between 32 and 256 to avoid underutilization of the GPU and
-        # errors due to too many registers used when the block size is too large.
+        # `eager_copyto_stencil_kernel!` requires one x-thread per face level, so a block
+        # is (n_face_levels, columns_per_block, 1) and the grid indexes the remaining
+        # columns. Each thread derives a linear column index from its y-thread and block
+        # index and decomposes it into `(i, j, h)`; this layout suits VIJFH, where the
+        # vertical axis is contiguous. Both periodic and non-periodic vertical topologies
+        # are supported (`has_padding_thread` accounts for the extra face level of
+        # non-periodic spaces), as are masked spaces. High-resolution columns (more face
+        # levels than fit in a block) fall through to `copyto_stencil_kernel!` below.
         # TODO: auto reduce max reg usage when needed because of high res columns
         # Size the dynamic shared memory to fit the largest single expression result in
         # the broadcasted tree (see `max_eager_shmem_per_thread`). If even that does not
         # fit in the device's per-block shared memory, there is no way to eagerly evaluate
         # the expression, so error out instead of silently falling back.
         eager_shmem_per_thread = max_eager_shmem_per_thread(bc′)
         if !high_resolution
-            #    32 <= n_face_levels * Ni <= 256
             # mask.N holds the active column count in a one-element device array;
             # reading it on the host needs @allowscalar.
             n_columns =
                 mask isa NoMask ? Ni * Nj * Nh :
                 CUDA.@allowscalar(mask.N[1])
+            # One column per block keeps register pressure low, which matters more than
+            # occupancy until there are enough columns to saturate the device; past that,
+            # pack as many columns into each block as 256 threads allow.
             # 108 is the number of SMs in an A100. TODO: get this value from CUDA.jl to better optimize for different GPUs
             threads_dim_y = n_columns > 256 * 108 ? div(256, n_face_levels) : 1
             block_dim_x = div(n_columns, threads_dim_y, RoundUp)
@@ -112,11 +113,10 @@ function Base.copyto!(
                  Split the expression into smaller sub-expressions so that each \
                  intermediate matrix/vector result is smaller.",
             )
-            # `us` (a `UniversalSize`) encodes `Nij` in its type, so the kernel
-            # decomposes the linear column index into `(i, j, h)` using a
-            # `CartesianIndices` whose horizontal extents are compile-time
-            # constants. This keeps the per-thread `divrem` a cheap fixed-divisor
-            # operation
+            # `axes(out)` is passed so the kernel can recover the output layout's
+            # horizontal extents from its type parameters (see `vijh_params`). The
+            # `CartesianIndices` the kernel builds from them therefore divides by
+            # compile-time constants, keeping the per-thread `divrem` cheap.
             args = (
                 strip_space(out, space),
                 strip_space(bc′, space),

src/MatrixFields/operator_matrices.jl

@@ -550,8 +550,10 @@ Operators.stencil_right_boundary(
 # An operator matrix can only capture the *linear* part of the operator; the
 # constant contributed by the boundary value is zeroed out (see `has_affine_bc`).
 # For every operator except GradientF2C/DivergenceF2C, a value-fixing condition
-# prescribes the *output* at the boundary as a pure constant, so its linear part is
-# zero and the boundary row is empty (`rzero`).
+# prescribes the *output* at the boundary as a pure constant, so its linear part is zero
+# and the boundary row is all zeros (`rzero` of the row type, which keeps the row's
+# bandwidth and zeroes its entries; the multiply clips the out-of-range band entries at
+# the column ends, so the row need not be narrowed).
 const ValueFixingBoundaryCondition = Union{
     Operators.SetValue,
     Operators.SetGradient,

src/Operators/finitedifference.jl

@@ -1007,9 +1007,11 @@ Base.@propagate_inbounds function stencil_left_boundary(
         getidx(space, velocity, idx, hidx),
         Geometry.LocalGeometry(space, idx, hidx),
     )
-    # `a⁻` is the (deleted) `LeftBiasedC2F` interior stencil, inlined here.
+    # `a⁻` is `LeftBiasedC2F`'s interior stencil, inlined here (the operator remains,
+    # but its `stencil_interior` method is now derived from its operator matrix).
     a⁻ = getidx(space, arg, idx - half, hidx)
-    # `a⁺` is the (deleted) `RightBiased3rdOrderC2F` interior stencil, inlined here.
+    # `a⁺` is the removed `RightBiased3rdOrderC2F` operator's interior stencil, inlined
+    # here because this is now its only remaining use.
     a⁺ =
         (
             4 * getidx(space, arg, idx - half, hidx) +
@@ -1034,14 +1036,16 @@ Base.@propagate_inbounds function stencil_right_boundary(
         getidx(space, velocity, idx, hidx),
         Geometry.LocalGeometry(space, idx, hidx),
     )
-    # `a⁻` is the (deleted) `LeftBiased3rdOrderC2F` interior stencil, inlined here.
+    # `a⁻` is the removed `LeftBiased3rdOrderC2F` operator's interior stencil, inlined
+    # here because this is now its only remaining use.
     a⁻ =
         (
             -2 * getidx(space, arg, idx - 1 - half, hidx) +
             10 * getidx(space, arg, idx - half, hidx) +
             4 * getidx(space, arg, idx + half, hidx)
         ) / 12
-    # `a⁺` is the (deleted) `RightBiasedC2F` interior stencil, inlined here.
+    # `a⁺` is `RightBiasedC2F`'s interior stencil, inlined here (the operator remains,
+    # but its `stencil_interior` method is now derived from its operator matrix).
     a⁺ = getidx(space, arg, idx + half, hidx)
     return Geometry.Contravariant3Vector(upwind_biased_product(v, a⁻, a⁺))
 
@@ -1063,7 +1067,8 @@ Base.@propagate_inbounds function stencil_left_boundary(
         getidx(space, velocity, idx, hidx),
         Geometry.LocalGeometry(space, idx, hidx),
     )
-    # `a` is the (deleted) `RightBiased3rdOrderC2F` interior stencil, inlined here.
+    # `a` is the removed `RightBiased3rdOrderC2F` operator's interior stencil, inlined
+    # here because this is now its only remaining use.
     a =
         (
             4 * getidx(space, arg, idx - half, hidx) +
@@ -1090,7 +1095,8 @@ Base.@propagate_inbounds function stencil_right_boundary(
         getidx(space, velocity, idx, hidx),
         Geometry.LocalGeometry(space, idx, hidx),
     )
-    # `a` is the (deleted) `LeftBiased3rdOrderC2F` interior stencil, inlined here.
+    # `a` is the removed `LeftBiased3rdOrderC2F` operator's interior stencil, inlined
+    # here because this is now its only remaining use.
     a =
         (
             -2 * getidx(space, arg, idx - 1 - half, hidx) +
@@ -1632,9 +1638,23 @@ abstract type BoundaryOperator <: FiniteDifferenceOperator end
 """
     SetBoundaryOperator(;boundaries...)
 
-This operator only modifies the values at the boundary faces, or the center cells adjacent to the boundary faces:
+This operator is the identity in the interior, and replaces the value at each boundary
+for which a condition is given. It preserves the space of its argument, so it modifies
+the boundary faces of a face field or the boundary center cells of a center field. A side
+with no condition is left untouched.
+
+The following boundary conditions are supported:
 
  - [`SetValue(val)`](@ref): set the value to be `val` on the boundary.
+ - [`SetGradient(val)`](@ref): set the value to be `val` on the boundary, projected onto
+   the `Covariant3` axis.
+ - [`SetCurl(val)`](@ref): set the value to be `val` on the boundary, projected onto the
+   `Contravariant123` axis.
+ - [`SetDivergence(val)`](@ref): set the value to be `val` on the boundary.
+
+The projecting conditions exist so that this operator can reapply the boundary conditions
+of the operator it was derived from when a broadcast is rewritten as an operator matrix
+multiply; see `MatrixFields.modifies_output`.
 """
 struct SetBoundaryOperator{BCS} <: BoundaryOperator
     bcs::BCS