From 1c3dee63f0d2b2514cba99bbad8ece745cd43938 Mon Sep 17 00:00:00 2001
From: Charles Kawczynski <kawczynski.charles@gmail.com>
Date: Thu, 3 Aug 2023 09:31:42 -0700
Subject: [PATCH] Improve col integ test, add extruded face to all_spaces

---
 test/Fields/field.jl                | 59 +++++++++++++++++++----------
 test/TestUtilities/TestUtilities.jl |  2 +-
 2 files changed, 41 insertions(+), 20 deletions(-)

diff --git a/test/Fields/field.jl b/test/Fields/field.jl
index 06e8b62c82..dfa7ea9041 100644
--- a/test/Fields/field.jl
+++ b/test/Fields/field.jl
@@ -701,12 +701,18 @@ 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,
+    ]
+        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 +721,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,14 +731,19 @@ 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
@@ -740,11 +751,21 @@ convergence_rate(err, Δh) =
             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
diff --git a/test/TestUtilities/TestUtilities.jl b/test/TestUtilities/TestUtilities.jl
index 35aac3fded..f9957fcace 100644
--- a/test/TestUtilities/TestUtilities.jl
+++ b/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)