Skip to content

Commit

Permalink
update some tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@simonbyrne
simonbyrne committed Jul 18, 2023
1 parent 4e59637 commit 6a00f21
Show file tree
Hide file tree
Showing 3 changed files with 55 additions and 54 deletions.
92 changes: 40 additions & 52 deletions test/Operators/finitedifference/column.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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]

Expand All @@ -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)

Expand All @@ -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]

Expand All @@ -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))
Expand All @@ -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(),
)
Expand All @@ -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]

Expand All @@ -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))
Expand All @@ -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(),
)
Expand All @@ -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
Expand All @@ -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)

Expand All @@ -870,21 +866,18 @@ 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ᶠ))
adv_wc = @. divf2c.(first_order_fluxsinᶠ) + corrected_antidiff_flux
corrected_antidiff_flux = @. divf2c(C * w * (U3(w, c) - U1(w, c)))
adv_wc = @. divf2c(w * U1(w, c)) + corrected_antidiff_flux

Δh[k] = cs.face_local_geometry.J[1]

Expand All @@ -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))
Expand All @@ -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(),
)
Expand All @@ -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

Expand All @@ -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))
Expand All @@ -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(),
)
Expand All @@ -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
Expand All @@ -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
Expand Down
15 changes: 13 additions & 2 deletions test/Operators/finitedifference/column_benchmark_kernels.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 2 additions & 0 deletions test/Operators/finitedifference/column_benchmark_utils.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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!)) =
Expand Down

0 comments on commit 6a00f21

Please sign in to comment.