Merge #1411

1411: Improve col integ test, add extruded face to all_spaces r=charleskawczynski a=charleskawczynski

This PR adds `FaceExtrudedFiniteDifferenceSpace` to `all_spaces`, and fixes the test that breaks with it.

This revealed that our definite face column integrals are only 1st order accurate, which is perhaps fine, but used to be unclear from the tests.

Co-authored-by: Charles Kawczynski <kawczynski.charles@gmail.com>

test/Fields/field.jl

@@ -701,12 +701,19 @@ convergence_rate(err, Δh) =
     end
     device = ClimaComms.device()
     context = ClimaComms.SingletonCommsContext(device)
-    for zelem in (2^2, 2^3, 2^4, 2^5)
-        for space in TU.all_spaces(FT; zelem, context)
-            # Filter out spaces without z coordinates:
-            TU.has_z_coordinates(space) || continue
-            # Skip spaces incompatible with Fields.bycolumn:
-            TU.bycolumnable(space) || continue
+    for space_constructor in [
+        TU.ColumnCenterFiniteDifferenceSpace,
+        TU.ColumnFaceFiniteDifferenceSpace,
+        TU.CenterExtrudedFiniteDifferenceSpace,
+        TU.FaceExtrudedFiniteDifferenceSpace,
+    ]
+        device isa ClimaComms.CUDADevice && continue # broken on gpu
+        for zelem in (2^2, 2^3, 2^4, 2^5)
+            space = space_constructor(FT; zelem, context)
+            # # Filter out spaces without z coordinates:
+            # TU.has_z_coordinates(space) || continue
+            # # Skip spaces incompatible with Fields.bycolumn:
+            # TU.bycolumnable(space) || continue
 
             Y = fill((; y = FT(1)), space)
             zcf = Fields.coordinate_field(Y.y).z
@@ -715,7 +722,7 @@ convergence_rate(err, Δh) =
             Δz_1 = CUDA.allowscalar() do
                 parent(Δz_col)[1]
             end
-            key = (space, zelem)
+            key = zelem
             if !haskey(results, key)
                 results[key] =
                     Dict(:maxerr => 0, :Δz_1 => FT(0), :∫y => [], :y => [])
@@ -725,26 +732,41 @@ convergence_rate(err, Δh) =
             y = Y.y
             @. y .= 1 + sin(zcf)
             # Compute max error against analytic solution
-            maxerr = 0
-            Fields.bycolumn(axes(Y.y)) do colidx
-                Operators.column_integral_definite!(∫y[colidx], y[colidx])
-                maxerr = max(
-                    maxerr,
-                    maximum(abs.(parent(∫y[colidx]) .- ∫y_analytic)),
-                )
-                nothing
+            maxerr = FT(0)
+            if space isa Spaces.ExtrudedFiniteDifferenceSpace
+                Fields.bycolumn(axes(Y.y)) do colidx
+                    Operators.column_integral_definite!(∫y[colidx], y[colidx])
+                    maxerr = max(
+                        maxerr,
+                        maximum(abs.(parent(∫y[colidx]) .- ∫y_analytic)),
+                    )
+                    nothing
+                end
+            else
+                Operators.column_integral_definite!(∫y, y)
+                maxerr = max(maxerr, maximum(abs.(parent(∫y) .- ∫y_analytic)))
             end
             results[key][:Δz_1] = Δz_1
             results[key][:maxerr] = maxerr
             push!(results[key][:∫y], ∫y)
             push!(results[key][:y], y)
             nothing
         end
+        maxerrs = map(key -> results[key][:maxerr], collect(keys(results)))
+        Δzs_1 = map(key -> results[key][:Δz_1], collect(keys(results)))
+        cr = convergence_rate(maxerrs, Δzs_1)
+        if nameof(space_constructor) == :ColumnCenterFiniteDifferenceSpace
+            @test 2 < sum(abs.(cr)) / length(cr) < 2.01
+        elseif nameof(space_constructor) == :ColumnFaceFiniteDifferenceSpace
+            @test_broken 2 < sum(abs.(cr)) / length(cr) < 2.01
+        elseif nameof(space_constructor) == :CenterExtrudedFiniteDifferenceSpace
+            @test 2 < sum(abs.(cr)) / length(cr) < 2.01
+        elseif nameof(space_constructor) == :FaceExtrudedFiniteDifferenceSpace
+            @test_broken 2 < sum(abs.(cr)) / length(cr) < 2.01
+        else
+            error("Uncaught case")
+        end
     end
-    maxerrs = map(key -> results[key][:maxerr], collect(keys(results)))
-    Δzs_1 = map(key -> results[key][:Δz_1], collect(keys(results)))
-    cr = convergence_rate(maxerrs, Δzs_1)
-    @test 2 < sum(abs.(cr)) / length(cr) < 2.01
 end
 
 @testset "Allocation tests for integrals" begin

test/Spaces/extruded_cuda.jl

@@ -26,7 +26,7 @@ compare(cpu, gpu, sym) =
     )
     collect(TU.all_spaces(Float64; zelem = 10, context)) # make sure we can construct spaces
     as = collect(TU.all_spaces(Float64; zelem = 10, context))
-    @test length(as) == 7
+    @test length(as) == 8
 end
 
 @testset "copyto! with CuArray-backed extruded spaces" begin

test/TestUtilities/TestUtilities.jl

@@ -155,7 +155,7 @@ function all_spaces(
         ColumnFaceFiniteDifferenceSpace(FT; zelem, context),
         SphereSpectralElementSpace(FT; context),
         CenterExtrudedFiniteDifferenceSpace(FT; zelem, context, helem),
-        # FaceExtrudedFiniteDifferenceSpace(FT; zelem, context, helem), # errors on sum
+        FaceExtrudedFiniteDifferenceSpace(FT; zelem, context, helem),
         # TODO: incorporate this list of spaces somehow:
         #     space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z)
         #     space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad)