From 237fce14e771cc011ce9caf02781643d603236a6 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Tue, 18 Jul 2023 10:24:12 -0700 Subject: [PATCH] update some tests --- test/Operators/finitedifference/column.jl | 90 ++++++++----------- .../column_benchmark_kernels.jl | 15 +++- .../column_benchmark_utils.jl | 2 + 3 files changed, 54 insertions(+), 53 deletions(-) diff --git a/test/Operators/finitedifference/column.jl b/test/Operators/finitedifference/column.jl index 3f23a42e8a..bfa7523b97 100644 --- 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,7 +670,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-4 @test conv_adv_wc[1] ≈ 3 atol = 0.1 @test conv_adv_wc[2] ≈ 3 atol = 0.1 @@ -680,7 +678,7 @@ end @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,7 +748,7 @@ 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 @@ -760,7 +756,7 @@ end 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,7 +780,7 @@ 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 @@ -792,7 +788,7 @@ end 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,7 +836,7 @@ 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 @@ -848,7 +844,7 @@ end 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,7 +888,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-4 @test conv_adv_wc[1] ≈ 3 atol = 0.1 @test conv_adv_wc[2] ≈ 3 atol = 0.1 @@ -903,7 +896,7 @@ end @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,7 +918,7 @@ 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 @@ -933,11 +926,11 @@ end 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,7 +956,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] ≤ 0.2006 @test conv_adv_wc[1] ≈ 0.5 atol = 0.2 @test conv_adv_wc[2] ≈ 0.5 atol = 0.3 @@ -973,7 +964,7 @@ end 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 index d320a52a85..a5587e2da0 100644 --- 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 index 7e1bd62606..9c9fc14527 100644 --- 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!)) =