run for HybridBox

CliMA · Apr 12, 2024 · 627e577 · 627e577
1 parent 80e1a5a
commit 627e577
Show file tree

Hide file tree

Showing 4 changed files with 33 additions and 28 deletions.
diff --git a/src/shared_utilities/Domains.jl b/src/shared_utilities/Domains.jl
@@ -659,11 +659,8 @@ corresponding to the surface, with values equal to the topmost level.
 """
 function top_center_to_surface(center_field::ClimaCore.Fields.Field)
     center_space = axes(center_field)
-    # face_space = obtain_face_space(center_space)
-    # N = ClimaCore.Spaces.nlevels(face_space)
     N_minus_half = ClimaCore.Spaces.nlevels(center_space)
     surface_space = obtain_surface_space(center_space)
-    # level = ClimaCore.Fields.level(center_field, ClimaCore.Spaces.nlevels(center_space))
     return ClimaCore.Fields.Field(
         ClimaCore.Fields.field_values(
             ClimaCore.Fields.level(center_field, N_minus_half),

diff --git a/src/standalone/Soil/boundary_conditions.jl b/src/standalone/Soil/boundary_conditions.jl
@@ -489,13 +489,18 @@ function ClimaLand.set_dfluxBCdY!(
     θ_rN = Domains.top_center_to_surface(θ_r)
     S_sN = Domains.top_center_to_surface(S_s)
 
-    # Extract the local geometry of the top face space
-    local_geometryN = Fields.local_geometry_field(KN)
+
+    # Get the local geometry of the face space, then extract the top level
+    levels = ClimaCore.Spaces.nlevels(Domains.obtain_face_space(axes(p.soil.K)))
+    local_geometry_faceN = ClimaCore.Fields.level(
+        Fields.local_geometry_field(Domains.obtain_face_space(axes(p.soil.K))),
+        levels - ClimaCore.Utilities.half,
+    )
 
     # Update dfluxBCdY at the top boundary in place
     # Calculate the value and convert it to a Covariant3Vector
     @. p.soil.dfluxBCdY =
-        covariant3_unit_vector(local_geometryN) *
+        covariant3_unit_vector(local_geometry_faceN) *
         (KN * dψdϑ(hydrology_cmN, ϑ_lN, νN, θ_rN, S_sN) / Δz)
     return nothing
 end

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -25,9 +25,9 @@ end
 end
 
 # Standalone Bucket model tests
-@safetestset "Bucket albedo types tests" begin
-    include("standalone/Bucket/albedo_types.jl")
-end
+# @safetestset "Bucket albedo types tests" begin
+#     include("standalone/Bucket/albedo_types.jl")
+# end
 @safetestset "Bucket snow tests" begin
     include("standalone/Bucket/snow_bucket_tests.jl")
 end

diff --git a/test/shared_utilities/implicit_timestepping/richards_model.jl b/test/shared_utilities/implicit_timestepping/richards_model.jl
@@ -83,7 +83,7 @@ for FT in (Float32, Float64)
                 @test Array(parent(jac_ϑ_l.entries.:3))[end] == FT(0)
                 @test all(
                     Array(parent(jac_ϑ_l.entries.:1))[2:end] .≈
-                    K_ic / dz^2 * dψdϑ_ic,
+                    dtγ * (K_ic / dz^2 * dψdϑ_ic),
                 )
                 # Check values on main diagonal: note jac_ϑ_l is I-γJ
                 @test Array(parent(jac_ϑ_l.entries.:2))[1] .≈
@@ -111,7 +111,7 @@ for FT in (Float32, Float64)
                 )
                 @test all(
                     Array(parent(jac_ϑ_l.entries.:1))[2:end, :, 1, 1, 1] .≈
-                    K_ic / dz^2 * dψdϑ_ic,
+                    dtγ * (K_ic / dz^2 * dψdϑ_ic),
                 )
                 # Check values on main diagonal: note jac_ϑ_l is I-γJ
                 @test all(
@@ -127,11 +127,12 @@ for FT in (Float32, Float64)
                         1,
                     ] .≈ dtγ * (-2 * K_ic / dz^2 * dψdϑ_ic) - I,
                 )
-                # TODO these fail
-                # @test all(
-                #     Array(parent(jac_ϑ_l.entries.:2))[end, :, 1, 1, 1] .≈
-                #     dtγ * (-K_ic / dz^2 * dψdϑ_ic - K_ic / (dz * dz / 2) * dψdϑ_ic) - I,
-                # )
+                @test all(
+                    Array(parent(jac_ϑ_l.entries.:2))[end, :, 1, 1, 1] .≈
+                    dtγ *
+                    (-K_ic / dz^2 * dψdϑ_ic - K_ic / (dz * dz / 2) * dψdϑ_ic) -
+                    I,
+                )
             end
 
         end
@@ -188,16 +189,19 @@ for FT in (Float32, Float64)
             )
             dz = FT(0.01)
             dψdϑ_ic = dψdϑ(hcm, FT(0.24), ν, θ_r, S_s)
+            jac_ϑ_l = jacobian.matrix[
+                MatrixFields.@name(soil.ϑ_l),
+                MatrixFields.@name(soil.ϑ_l)
+            ]
             if typeof(domain) <: Column
                 @test Array(parent(jac_ϑ_l.entries.:2))[1] .≈
                       dtγ * (-K_ic / dz^2 * dψdϑ_ic) - I
-                # TODO these fail
-                # @test all(
-                #     Array(parent(jac_ϑ_l.entries.:2))[2:(end - 1)] .≈
-                #     dtγ * (-2 * K_ic / dz^2 * dψdϑ_ic) - I,
-                # )
-                # @test Array(parent(jac_ϑ_l.entries.:2))[end] .≈
-                #     dtγ * (-K_ic / dz^2 * dψdϑ_ic) - I
+                @test all(
+                    Array(parent(jac_ϑ_l.entries.:2))[2:(end - 1)] .≈
+                    dtγ * (-2 * K_ic / dz^2 * dψdϑ_ic) - I,
+                )
+                @test Array(parent(jac_ϑ_l.entries.:2))[end] .≈
+                      dtγ * (-K_ic / dz^2 * dψdϑ_ic) - I
             elseif typeof(domain) <: HybridBox
                 @test all(
                     Array(parent(jac_ϑ_l.entries.:2))[1, :, 1, 1, 1] .≈
@@ -212,11 +216,10 @@ for FT in (Float32, Float64)
                         1,
                     ] .≈ dtγ * (-2 * K_ic / dz^2 * dψdϑ_ic) - I,
                 )
-                # TODO these fail
-                # @test all(
-                #     Array(parent(jac_ϑ_l.entries.:2))[end, :, 1, 1, 1] .≈
-                #     dtγ * (-K_ic / dz^2 * dψdϑ_ic - K_ic / (dz * dz / 2) * dψdϑ_ic) - I,
-                # )
+                @test all(
+                    Array(parent(jac_ϑ_l.entries.:2))[end, :, 1, 1, 1] .≈
+                    dtγ * (-K_ic / dz^2 * dψdϑ_ic) - I,
+                )
             end
         end
     end