update some tests

CliMA · Jul 18, 2023 · 237fce1 · 237fce1
1 parent 4e59637
commit 237fce1
Show file tree

Hide file tree

Showing 3 changed files with 54 additions and 53 deletions.
diff --git a/test/Operators/finitedifference/column.jl b/test/Operators/finitedifference/column.jl
@@ -630,7 +630,7 @@ end
     @test conv_curl_sin_f[1] ≤ conv_curl_sin_f[2] ≤ conv_curl_sin_f[3]
 end
 
-@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform) periodic mesh, constant w" begin
+@testset "Upwind3rdOrderC2F + DivergenceF2C on (uniform) periodic mesh, constant w" begin
     FT = Float64
     n_elems_seq = 2 .^ (5, 6, 7, 8)
 
@@ -651,17 +651,15 @@ end
 
         centers = getproperty(Fields.coordinate_field(cs), :z)
 
-        # Upwind3rdOrderBiasedProductC2F Center -> Face operator
+        # Upwind3rdOrderC2F Center -> Face operator
         # Unitary, constant advective velocity
         w = Geometry.WVector.(ones(fs))
         # c = sin(z), scalar field defined at the centers
         c = sin.(centers)
 
-        third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F()
-        third_order_fluxsinᶠ = third_order_fluxᶠ.(w, c)
-
+        U = Operators.Upwind3rdOrderC2F()
         divf2c = Operators.DivergenceF2C()
-        adv_wc = divf2c.(third_order_fluxsinᶠ)
+        adv_wc = divf2c.(w .* U.(w, c))
 
         Δh[k] = cs.face_local_geometry.J[1]
 
@@ -672,15 +670,15 @@ end
     # Check convergence rate
     conv_adv_wc = convergence_rate(err_adv_wc, Δh)
 
-    # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z)
+    # Upwind3rdOrderC2F conv, with f(z) = sin(z)
     @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 5e-4
     @test conv_adv_wc[1] ≈ 3 atol = 0.1
     @test conv_adv_wc[2] ≈ 3 atol = 0.1
     @test conv_adv_wc[3] ≈ 3 atol = 0.1
     @test conv_adv_wc[1] ≤ conv_adv_wc[2] ≤ conv_adv_wc[2]
 end
 
-@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform) periodic mesh, varying sign w" begin
+@testset "Upwind3rdOrderC2F + DivergenceF2C on (uniform) periodic mesh, varying sign w" begin
     FT = Float64
     n_elems_seq = 2 .^ (5, 6, 7, 8)
 
@@ -702,17 +700,15 @@ end
         centers = getproperty(Fields.coordinate_field(cs), :z)
         faces = getproperty(Fields.coordinate_field(fs), :z)
 
-        # Upwind3rdOrderBiasedProductC2F Center -> Face operator
+        # Upwind3rdOrderC2F Center -> Face operator
         # w = cos(z), vertical velocity field defined at the faces
         w = Geometry.WVector.(cos.(faces))
         # c = sin(z), scalar field defined at the centers
         c = sin.(centers)
 
-        third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F()
-        third_order_fluxsinᶠ = third_order_fluxᶠ.(w, c)
-
+        U = Operators.Upwind3rdOrderC2F()
         divf2c = Operators.DivergenceF2C()
-        adv_wc = divf2c.(third_order_fluxsinᶠ)
+        adv_wc = divf2c.(w .* U.(w, c))
 
         Δh[k] = cs.face_local_geometry.J[1]
 
@@ -724,14 +720,14 @@ end
     # Check convergence rate
     conv_adv_wc = convergence_rate(err_adv_wc, Δh)
 
-    # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z), w(z) = cos(z)
+    # Upwind3rdOrderC2F conv, with f(z) = sin(z), w(z) = cos(z)
     @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 4e-3
     @test conv_adv_wc[1] ≈ 2 atol = 0.2
     @test conv_adv_wc[2] ≈ 2 atol = 0.1
     @test conv_adv_wc[3] ≈ 2 atol = 0.1
 end
 
-@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided1stOrder + DivergenceF2C SetValue BCs, constant w" begin
+@testset "Upwind3rdOrderC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided1stOrder + DivergenceF2C SetValue BCs, constant w" begin
     FT = Float64
     n_elems_seq = 2 .^ (4, 6, 8, 10)
     stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0))
@@ -752,15 +748,15 @@ end
 
             centers = getproperty(Fields.coordinate_field(cs), :z)
 
-            # Upwind3rdOrderBiasedProductC2F Center -> Face operator
+            # Upwind3rdOrderC2F Center -> Face operator
             # Unitary, constant advective velocity
             w = Geometry.WVector.(ones(fs))
             # c = sin(z), scalar field defined at the centers
             Δz = FT(2pi / n)
             c = (cos.(centers .- Δz / 2) .- cos.(centers .+ Δz / 2)) ./ Δz
             s = sin.(centers)
 
-            third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F(
+            U = Operators.Upwind3rdOrderC2F(
                 bottom = Operators.OneSided1stOrder(),
                 top = Operators.OneSided1stOrder(),
             )
@@ -774,7 +770,7 @@ end
                 ),
             )
 
-            adv_wc = divf2c.(third_order_fluxᶠ.(w, c))
+            adv_wc = divf2c.(w .* U.(w, c))
 
             Δh[k] = cs.face_local_geometry.J[1]
 
@@ -784,15 +780,15 @@ end
 
         # Check convergence rate
         conv_adv_wc = convergence_rate(err_adv_wc, Δh)
-        # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z), w(z) = 1
+        # Upwind3rdOrderC2F conv, with f(z) = sin(z), w(z) = 1
         @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 0.2006
         @test conv_adv_wc[1] ≈ 0.5 atol = 0.2
         @test conv_adv_wc[2] ≈ 0.5 atol = 0.3
         @test conv_adv_wc[3] ≈ 1.0 atol = 0.55
     end
 end
 
-@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder + DivergenceF2C SetValue BCs, varying sign w" begin
+@testset "Upwind3rdOrderC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder + DivergenceF2C SetValue BCs, varying sign w" begin
     FT = Float64
     n_elems_seq = 2 .^ (4, 6, 8, 10)
     stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0))
@@ -814,13 +810,13 @@ end
             centers = getproperty(Fields.coordinate_field(cs), :z)
             faces = getproperty(Fields.coordinate_field(fs), :z)
 
-            # Upwind3rdOrderBiasedProductC2F Center -> Face operator
+            # Upwind3rdOrderC2F Center -> Face operator
             # w = cos(z), vertical velocity field defined at the faces
             w = Geometry.WVector.(cos.(faces))
             # c = sin(z), scalar field defined at the centers
             c = sin.(centers)
 
-            third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F(
+            U = Operators.Upwind3rdOrderC2F(
                 bottom = Operators.OneSided3rdOrder(),
                 top = Operators.OneSided3rdOrder(),
             )
@@ -829,7 +825,7 @@ end
                 bottom = Operators.SetValue(Geometry.WVector(FT(0.0))),
                 top = Operators.SetValue(Geometry.WVector(FT(0.0))),
             )
-            adv_wc = divf2c.(third_order_fluxᶠ.(w, c))
+            adv_wc = divf2c.(w .* U.(w, c))
 
             Δh[k] = cs.face_local_geometry.J[1]
             # Errors
@@ -840,15 +836,15 @@ end
 
         # Check convergence rate
         conv_adv_wc = convergence_rate(err_adv_wc, Δh)
-        # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z), w(z) = cos(z)
+        # Upwind3rdOrderC2F conv, with f(z) = sin(z), w(z) = cos(z)
         @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 2e-1
         @test conv_adv_wc[1] ≈ 2 atol = 0.1
         @test conv_adv_wc[2] ≈ 2 atol = 0.1
         @test conv_adv_wc[3] ≈ 2 atol = 0.1
     end
 end
 
-@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform) periodic mesh" begin
+@testset "Simple FCT: lin combination of Upwind1stOrderC2F + Upwind3rdOrderC2F on (uniform) periodic mesh" begin
     FT = Float64
     n_elems_seq = 2 .^ (5, 6, 7, 8)
 
@@ -870,20 +866,17 @@ end
         centers = getproperty(Fields.coordinate_field(cs), :z)
         C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding)
 
-        # UpwindBiasedProductC2F & Upwind3rdOrderBiasedProductC2F Center -> Face operator
+        # Upwind1stOrderC2F & Upwind3rdOrderC2F Center -> Face operator
         # Unitary, constant advective velocity
         w = Geometry.WVector.(ones(fs))
         # c = sin(z), scalar field defined at the centers
         c = sin.(centers)
 
-        first_order_fluxᶠ = Operators.UpwindBiasedProductC2F()
-        third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F()
-        first_order_fluxsinᶠ = first_order_fluxᶠ.(w, c)
-        third_order_fluxsinᶠ = third_order_fluxᶠ.(w, c)
+        U1 = Operators.Upwind1stOrderC2F()
+        U3 = Operators.Upwind3rdOrderC2F()
 
         divf2c = Operators.DivergenceF2C()
-        corrected_antidiff_flux =
-            @. divf2c(C * (third_order_fluxsinᶠ - first_order_fluxsinᶠ))
+        corrected_antidiff_flux = @. divf2c(C * w * (U1(w, c) - U3(w, c)))
         adv_wc = @. divf2c.(first_order_fluxsinᶠ) + corrected_antidiff_flux
 
         Δh[k] = cs.face_local_geometry.J[1]
@@ -895,15 +888,15 @@ end
     # Check convergence rate
     conv_adv_wc = convergence_rate(err_adv_wc, Δh)
 
-    # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z)
+    # Upwind3rdOrderC2F conv, with f(z) = sin(z)
     @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 5e-4
     @test conv_adv_wc[1] ≈ 3 atol = 0.1
     @test conv_adv_wc[2] ≈ 3 atol = 0.1
     @test conv_adv_wc[3] ≈ 3 atol = 0.1
     @test conv_adv_wc[1] ≤ conv_adv_wc[2] ≤ conv_adv_wc[2]
 end
 
-@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with OneSided1stOrder BCs" begin
+@testset "Simple FCT: lin combination of Upwind1stOrderC2F + Upwind3rdOrderC2F on (uniform and stretched) non-periodic mesh, with OneSided1stOrder BCs" begin
     FT = Float64
     n_elems_seq = 2 .^ (4, 6, 8, 10)
     stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0))
@@ -925,19 +918,19 @@ end
             centers = getproperty(Fields.coordinate_field(cs), :z)
             C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding)
 
-            # UpwindBiasedProductC2F & Upwind3rdOrderBiasedProductC2F Center -> Face operator
+            # Upwind1stOrderC2F & Upwind3rdOrderC2F Center -> Face operator
             # Unitary, constant advective velocity
             w = Geometry.WVector.(ones(fs))
             # c = sin(z), scalar field defined at the centers
             Δz = FT(2pi / n)
             c = (cos.(centers .- Δz / 2) .- cos.(centers .+ Δz / 2)) ./ Δz
             s = sin.(centers)
 
-            first_order_fluxᶠ = Operators.UpwindBiasedProductC2F(
+            U1 = Operators.Upwind1stOrderC2F(
                 bottom = Operators.Extrapolate(),
                 top = Operators.Extrapolate(),
             )
-            third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F(
+            U3 = Operators.Upwind3rdOrderC2F(
                 bottom = Operators.OneSided1stOrder(),
                 top = Operators.OneSided1stOrder(),
             )
@@ -951,9 +944,7 @@ end
                 ),
             )
 
-            corrected_antidiff_flux = @. divf2c(
-                C * (third_order_fluxᶠ(w, c) - first_order_fluxᶠ(w, c)),
-            )
+            corrected_antidiff_flux = @. divf2c(C * w * (U3(w, c) - U1(w, c)))
             adv_wc =
                 @. divf2c.(first_order_fluxᶠ(w, c)) + corrected_antidiff_flux
 
@@ -965,15 +956,15 @@ end
 
         # Check convergence rate
         conv_adv_wc = convergence_rate(err_adv_wc, Δh)
-        # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z)
+        # Upwind3rdOrderC2F conv, with f(z) = sin(z)
         @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 0.2006
         @test conv_adv_wc[1] ≈ 0.5 atol = 0.2
         @test conv_adv_wc[2] ≈ 0.5 atol = 0.3
         @test conv_adv_wc[3] ≈ 1.0 atol = 0.55
     end
 end
 
-@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder BCs" begin
+@testset "Simple FCT: lin combination of Upwind1stOrderC2F + Upwind3rdOrderC2F on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder BCs" begin
     FT = Float64
     n_elems_seq = 2 .^ (4, 6, 8, 10)
     stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0))
@@ -995,17 +986,17 @@ end
             centers = getproperty(Fields.coordinate_field(cs), :z)
             C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding)
 
-            # UpwindBiasedProductC2F & Upwind3rdOrderBiasedProductC2F Center -> Face operator
+            # Upwind1stOrderC2F & Upwind3rdOrderC2F Center -> Face operator
             # Unitary, constant advective velocity
             w = Geometry.WVector.(ones(fs))
             # c = sin(z), scalar field defined at the centers
             c = sin.(centers)
 
-            first_order_fluxᶠ = Operators.UpwindBiasedProductC2F(
+            U1 = Operators.Upwind1stOrderC2F(
                 bottom = Operators.Extrapolate(),
                 top = Operators.Extrapolate(),
             )
-            third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F(
+            U3 = Operators.Upwind3rdOrderC2F(
                 bottom = Operators.OneSided3rdOrder(),
                 top = Operators.OneSided3rdOrder(),
             )
@@ -1014,11 +1005,8 @@ end
                 bottom = Operators.SetValue(Geometry.WVector(FT(0.0))),
                 top = Operators.SetValue(Geometry.WVector(FT(0.0))),
             )
-            corrected_antidiff_flux = @. divf2c(
-                C * (third_order_fluxᶠ(w, c) - first_order_fluxᶠ(w, c)),
-            )
-            adv_wc =
-                @. divf2c.(first_order_fluxᶠ(w, c)) + corrected_antidiff_flux
+            corrected_antidiff_flux = @. divf2c(C * w * (U3(w, c) - U1(w, c)))
+            adv_wc = @. divf2c.(U1(w, c)) + corrected_antidiff_flux
 
             Δh[k] = cs.face_local_geometry.J[1]
             # Errors
@@ -1028,7 +1016,7 @@ end
 
         # Check convergence rate
         conv_adv_wc = convergence_rate(err_adv_wc, Δh)
-        # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z)
+        # Upwind3rdOrderC2F conv, with f(z) = sin(z)
         @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 5e-1
         @test conv_adv_wc[1] ≈ 2.5 atol = 0.1
         @test conv_adv_wc[2] ≈ 2.5 atol = 0.1

diff --git a/test/Operators/finitedifference/column_benchmark_kernels.jl b/test/Operators/finitedifference/column_benchmark_kernels.jl
@@ -99,13 +99,24 @@ function op_divgrad_uₕ!(c, f, bcs)
     @. c.uₕ2 = div(f.y * grad(c.uₕ))
     return nothing
 end
-function op_divUpwind3rdOrderBiasedProductC2F!(c, f, bcs)
+function op_divUpwind3rdOrderBiasedProductC2F_old!(c, f, bcs)
     upwind = Operators.Upwind3rdOrderBiasedProductC2F(bcs.inner)
     divf2c = Operators.DivergenceF2C(bcs.outer)
     @. c.y = divf2c(upwind(f.w, c.x))
     return nothing
 end
-
+function op_divUpwind1sOrderBiasedProductC2F!(c, f, bcs)
+    upwind = Operators.Upwind1stOrderC2F(bcs.inner)
+    divf2c = Operators.DivergenceF2C(bcs.outer)
+    @. c.y = divf2c(f.w * upwind(f.w, c.x))
+    return nothing
+end
+function op_divUpwind3rdOrderBiasedProductC2F!(c, f, bcs)
+    upwind = Operators.Upwind3rdOrderC2F(bcs.inner)
+    divf2c = Operators.DivergenceF2C(bcs.outer)
+    @. c.y = divf2c(f.w * upwind(f.w, c.x))
+    return nothing
+end
 function op_broadcast_example0!(c, f, bcs)
     Fields.bycolumn(axes(f.ᶠu³)) do colidx
         CT3 = Geometry.Contravariant3Vector

diff --git a/test/Operators/finitedifference/column_benchmark_utils.jl b/test/Operators/finitedifference/column_benchmark_utils.jl
@@ -298,6 +298,8 @@ bcs_tested(c, ::typeof(op_UpwindBiasedProductC2F!)) = (set_value_bcs(c), extrapo
 bcs_tested(c, ::typeof(op_Upwind3rdOrderBiasedProductC2F!)) = (set_upwind_biased_3_bcs(c), extrapolate_bcs(c))
 
 # Composed operators (bcs handled case-by-case)
+bcs_tested(c, ::typeof(op_divUpwind1stOrderBiasedProductC2F!)) =
+    ((; inner = set_upwind_biased_3_bcs(c), outer = set_value_contra3_bcs(c)), )
 bcs_tested(c, ::typeof(op_divUpwind3rdOrderBiasedProductC2F!)) =
     ((; inner = set_upwind_biased_3_bcs(c), outer = set_value_contra3_bcs(c)), )
 bcs_tested(c, ::typeof(op_divgrad_CC!)) =