Fix unit_column.jl GPU failures

The reference values in the "Replacements for removed operators and
boundary conditions" testset are built with scalar indexing, so extract
the data they are built from (and the results they are compared
against) on the CPU.

Also apply the biased operators to `ᶠlg.J` rather than to `ᶠlg`: the GPU
compiles the never-taken `NullBoundaryCondition` branch of the stencil,
which multiplies by `NaN` and so cannot handle a `LocalGeometry`
element type.

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

Author: imreddyTeja

test/Operators/finitedifference/unit_column.jl

@@ -368,13 +368,17 @@ end
     wᶜ = Geometry.WVector.(1 .+ sin.(Fields.coordinate_field(cs).z))
     wᶠ = Geometry.WVector.(1 .+ sin.(Fields.coordinate_field(fs).z))
 
+    # the reference values below are built with scalar indexing, so the data
+    # they are built from has to be moved to the CPU first
+    cpu_parent(field) = Array(parent(field))
+
     # the quantities that the removed stencils were written in terms of
-    tᶜ = parent(θᶜ)[:]
-    tᶠ = parent(θᶠ)[:]
-    w³ᶜ = parent(Geometry.contravariant3.(wᶜ, ᶜlg))[:]
-    w³ᶠ = parent(Geometry.contravariant3.(wᶠ, ᶠlg))[:]
-    Jᶜ = parent(ᶜlg.J)[:]
-    Jᶠ = parent(ᶠlg.J)[:]
+    tᶜ = cpu_parent(θᶜ)[:]
+    tᶠ = cpu_parent(θᶠ)[:]
+    w³ᶜ = cpu_parent(Geometry.contravariant3.(wᶜ, ᶜlg))[:]
+    w³ᶠ = cpu_parent(Geometry.contravariant3.(wᶠ, ᶠlg))[:]
+    Jᶜ = cpu_parent(ᶜlg.J)[:]
+    Jᶠ = cpu_parent(ᶠlg.J)[:]
     θ₀ = FT(0.5) # boundary value
     interpf2c = Operators.InterpolateF2C()
 
@@ -395,7 +399,7 @@ end
                 Geometry.Covariant3Vector.(2 .* (θ₀ .- θ_top)),
             ),
         )
-        @test parent(gradc2f.(θᶜ))[:] ≈ ref
+        @test cpu_parent(gradc2f.(θᶜ))[:] ≈ ref
     end
 
     @testset "DivergenceC2F, SetValue(0)" begin
@@ -407,24 +411,26 @@ end
         ]
         # `LeftBiasedF2C(x)[i] = x[i-half]`, so its first level is the bottom
         # face; `RightBiasedF2C(x)[i] = x[i+half]`, so its last level is the top.
-        lg_bot_face = Fields.level(Operators.LeftBiasedF2C().(ᶠlg), 1)
-        lg_top_face = Fields.level(Operators.RightBiasedF2C().(ᶠlg), n)
+        # These are applied to `J` rather than to the local geometry itself, so
+        # that the stencil only ever operates on scalars.
+        J_bot_face = Fields.level(Operators.LeftBiasedF2C().(ᶠlg.J), 1)
+        J_top_face = Fields.level(Operators.RightBiasedF2C().(ᶠlg.J), n)
         set_bcs = Operators.SetBoundaryOperator(
             bottom = Operators.SetValue(
                 Geometry.Jcontravariant3.(
                     Fields.level(wᶜ, 1),
                     Fields.level(ᶜlg, 1),
-                ) .* (2 .* inv.(lg_bot_face.J)),
+                ) .* (2 .* inv.(J_bot_face)),
             ),
             top = Operators.SetValue(
                 Geometry.Jcontravariant3.(
                     Fields.level(wᶜ, n),
                     Fields.level(ᶜlg, n),
-                ) .* (-2 .* inv.(lg_top_face.J)),
+                ) .* (-2 .* inv.(J_top_face)),
             ),
         )
         divc2f = Operators.DivergenceC2F()
-        @test parent(@. set_bcs(divc2f(wᶜ)))[:] ≈ ref
+        @test cpu_parent(@. set_bcs(divc2f(wᶜ)))[:] ≈ ref
     end
 
     @testset "UpwindBiasedProductC2F, SetValue" begin
@@ -437,8 +443,8 @@ end
             Operators.upwind_biased_product(w³ᶠ[i], tᶜ[i - 1], tᶜ[i])
             for i in 1:(n + 1)
         ]
-        v_bot = parent(Geometry.contravariant3.(wᶠ, ᶠlg))[1]
-        v_top = parent(Geometry.contravariant3.(wᶠ, ᶠlg))[n + 1]
+        v_bot = w³ᶠ[1]
+        v_top = w³ᶠ[n + 1]
         set_bcs = Operators.SetBoundaryOperator(
             bottom = Operators.SetValue(
                 Geometry.Contravariant3Vector(
@@ -452,7 +458,7 @@ end
             ),
         )
         upwind = Operators.UpwindBiasedProductC2F()
-        @test parent(@. set_bcs(upwind(wᶠ, θᶜ)))[:] ≈ ref
+        @test cpu_parent(@. set_bcs(upwind(wᶠ, θᶜ)))[:] ≈ ref
     end
 
     @testset "AdvectionC2C, SetValue" begin
@@ -474,7 +480,7 @@ end
         new = @. interpf2c(
             Geometry.dot(Geometry.Contravariant3Vector(wᶠ), gradc2f(θᶜ)),
         )
-        @test parent(new)[:] ≈ ref
+        @test cpu_parent(new)[:] ≈ ref
     end
 
     @testset "AdvectionF2F" begin
@@ -486,7 +492,7 @@ end
             Geometry.Contravariant3Vector(wᶠ),
             interpc2f(gradf2c(θᶠ)),
         )
-        @test parent(new)[2:n] ≈ ref
+        @test cpu_parent(new)[2:n] ≈ ref
     end
 
     @testset "FluxCorrectionC2C, Extrapolate" begin
@@ -512,6 +518,6 @@ end
         )
         # `GradientF2C` returns a `Covariant3Vector`, whose component is the
         # difference of the fluxes across the cell
-        @test parent(gradf2c.(flux))[:] ≈ ref
+        @test cpu_parent(gradf2c.(flux))[:] ≈ ref
     end
 end