diff --git a/src/Models/HydrostaticFreeSurfaceModels/set_hydrostatic_free_surface_model.jl b/src/Models/HydrostaticFreeSurfaceModels/set_hydrostatic_free_surface_model.jl index e256f653dc..877dbc678e 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/set_hydrostatic_free_surface_model.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/set_hydrostatic_free_surface_model.jl @@ -62,7 +62,7 @@ model.velocities.u throw(ArgumentError("name $fldname not found in model.velocities, model.tracers, or model.free_surface")) end - @apply_regionally set!(ϕ, value) + set!(ϕ, value) end # initialize!(model) diff --git a/src/MultiRegion/multi_region_boundary_conditions.jl b/src/MultiRegion/multi_region_boundary_conditions.jl index b0011be7a6..acb6a032b7 100644 --- a/src/MultiRegion/multi_region_boundary_conditions.jl +++ b/src/MultiRegion/multi_region_boundary_conditions.jl @@ -13,7 +13,7 @@ using Oceananigans.BoundaryConditions: import Oceananigans.BoundaryConditions: fill_halo_regions!, fill_halo_event! -@inline bc_str(::MultiRegionObject) = "MultiRegion Boundary Conditions" +@inline bc_str(::MultiRegionObject) = "MultiRegionObject" @inline function fill_halo_regions!(fields::NamedTuple, grid::ConformalCubedSphereGridOfSomeKind, args...; kwargs...) u = haskey(fields, :u) ? fields.u : nothing diff --git a/src/MultiRegion/multi_region_field.jl b/src/MultiRegion/multi_region_field.jl index 09040ec9f5..3772876bb6 100644 --- a/src/MultiRegion/multi_region_field.jl +++ b/src/MultiRegion/multi_region_field.jl @@ -1,6 +1,6 @@ using Oceananigans.AbstractOperations: AbstractOperation, compute_computed_field! using Oceananigans.BoundaryConditions: default_auxiliary_bc -using Oceananigans.Fields: FunctionField, data_summary, AbstractField, instantiated_location +using Oceananigans.Fields: FunctionField, data_summary, AbstractField, instantiated_location, set_to_function! using Oceananigans.Operators: assumed_field_location using Oceananigans.OutputWriters: output_indices @@ -9,8 +9,17 @@ using Base: @propagate_inbounds import Oceananigans.BoundaryConditions: regularize_field_boundary_conditions import Oceananigans.Diagnostics: hasnan import Oceananigans.DistributedComputations: reconstruct_global_field, CommunicationBuffers -import Oceananigans.Fields: set!, compute!, compute_at!, interior, communication_buffers, - validate_field_data, validate_boundary_conditions, validate_indices +import Oceananigans.Fields: set_to_array!, + set_to_field!, + compute!, + compute_at!, + interior, + communication_buffers, + validate_field_data, + validate_boundary_conditions, + validate_indices, + set! + import Oceananigans.Grids: xnodes, ynodes import Base: fill!, axes @@ -111,15 +120,20 @@ function reconstruct_global_indices(indices, p::YPartition, N) return (idx_x, idx_y, idx_z) end -## Functions applied regionally -set!(mrf::MultiRegionField, v) = apply_regionally!(set!, mrf, v) -fill!(mrf::MultiRegionField, v) = apply_regionally!(fill!, mrf, v) +# Set fields regionally +set_to_array!(mrf::MultiRegionField, a) = apply_regionally!(set_to_array!, mrf, a) +set_to_field!(mrf::MultiRegionField, v) = apply_regionally!(set_to_field!, mrf, v) -set!(mrf::MultiRegionField, a::Number) = apply_regionally!(set!, mrf, a) +# Exporting some set! methods +# set! a function field can be safely done regionally +set!(mrf::MultiRegionField, f::FunctionField) = apply_regionally!(set!, mrf, f) +set!(mrf::MultiRegionField, f::Function) = set_to_function!(mrf, f) +set!(mrf::MultiRegionField, mro::MultiRegionObject) = apply_regionally!(set!, mrf, mro) + +# Fill fields regionally +fill!(mrf::MultiRegionField, v) = apply_regionally!(fill!, mrf, v) fill!(mrf::MultiRegionField, a::Number) = apply_regionally!(fill!, mrf, a) -set!(mrf::MultiRegionField, f::Function) = apply_regionally!(set!, mrf, f) -set!(u::MultiRegionField, v::MultiRegionField) = apply_regionally!(set!, u, v) compute!(mrf::GriddedMultiRegionField, time=nothing) = apply_regionally!(compute!, mrf, time) # Disambiguation (same as computed_field.jl:64) @@ -216,3 +230,23 @@ ynodes(ψ::AbstractField{<:Any, <:Any, <:Any, <:OrthogonalSphericalShellGrid}) = # Convenience @propagate_inbounds Base.getindex(mrf::MultiRegionField, r::Int) = getregion(mrf, r) @propagate_inbounds Base.lastindex(mrf::MultiRegionField) = lastindex(mrf.grid) + +import Base: == + +function ==(a::MultiRegionField, b::MultiRegionField) + if regions(a) == regions(b) + return all(a[r] == b[r] for r in regions(a)) + else + return false + end +end + +import Base: isapprox + +function isapprox(a::MultiRegionField, b::MultiRegionField; kw...) + if regions(a) == regions(b) + return all(isapprox(a[r], b[r]; kw...) for r in regions(a)) + else + return false + end +end diff --git a/src/OrthogonalSphericalShellGrids/conformal_cubed_sphere_panel.jl b/src/OrthogonalSphericalShellGrids/conformal_cubed_sphere_panel.jl index 35b2f76985..bb7b07ff68 100644 --- a/src/OrthogonalSphericalShellGrids/conformal_cubed_sphere_panel.jl +++ b/src/OrthogonalSphericalShellGrids/conformal_cubed_sphere_panel.jl @@ -34,7 +34,8 @@ function Adapt.adapt_structure(to, conformal_mapping::CubedSphereConformalMappin adapt(to, conformal_mapping.ξᶠᵃᵃ), adapt(to, conformal_mapping.ηᵃᶠᵃ), adapt(to, conformal_mapping.ξᶜᵃᵃ), - adapt(to, conformal_mapping.ηᵃᶜᵃ)) + adapt(to, conformal_mapping.ηᵃᶜᵃ) + ) end const ConformalCubedSpherePanelGrid{FT, TX, TY, TZ, CZ, CC, FC, CF, FF, Arch} = diff --git a/src/Utils/multi_region_transformation.jl b/src/Utils/multi_region_transformation.jl index d6912624c0..c2c192fe1a 100644 --- a/src/Utils/multi_region_transformation.jl +++ b/src/Utils/multi_region_transformation.jl @@ -14,6 +14,12 @@ struct MultiRegionObject{R} regional_objects :: R end +function Base.summary(mo::MultiRegionObject) + obj = summary(first(mo.regional_objects)) + N = length(mo.regional_objects) + return "MultiRegionObject with $N × $obj" +end + ##### ##### Convenience structs ##### diff --git a/test/runtests.jl b/test/runtests.jl index 6c7c2dd450..7873dc8185 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -232,7 +232,6 @@ CUDA.allowscalar() do end end - # Tests for Enzyme extension if group == :enzyme || group == :all @testset "Enzyme extension tests" begin diff --git a/test/test_cubed_sphere_circulation.jl b/test/test_cubed_sphere_circulation.jl deleted file mode 100644 index faeb9c210c..0000000000 --- a/test/test_cubed_sphere_circulation.jl +++ /dev/null @@ -1,173 +0,0 @@ -using Oceananigans.Operators: Γᶠᶠᶜ, ζ₃ᶠᶠᶜ - -function diagnose_velocities_from_streamfunction(ψ, arch, grid) - ψᶠᶠᶜ = Field{Face, Face, Center}(grid) - uᶠᶜᶜ = Field{Face, Center, Center}(grid) - vᶜᶠᶜ = Field{Center, Face, Center}(grid) - - for (f, grid_face) in enumerate(grid.faces) - Nx, Ny, Nz = size(grid_face) - - ψᶠᶠᶜ_face = ψᶠᶠᶜ.data.faces[f] - uᶠᶜᶜ_face = uᶠᶜᶜ.data.faces[f] - vᶜᶠᶜ_face = vᶜᶠᶜ.data.faces[f] - - for i in 1:Nx+1, j in 1:Ny+1 - ψᶠᶠᶜ_face[i, j, 1] = ψ(grid_face.λᶠᶠᵃ[i, j], grid_face.φᶠᶠᵃ[i, j]) - end - - for i in 1:Nx+1, j in 1:Ny - uᶠᶜᶜ_face[i, j, 1] = (ψᶠᶠᶜ_face[i, j, 1] - ψᶠᶠᶜ_face[i, j+1, 1]) / grid.faces[f].Δyᶠᶜᵃ[i, j] - end - - for i in 1:Nx, j in 1:Ny+1 - vᶜᶠᶜ_face[i, j, 1] = (ψᶠᶠᶜ_face[i+1, j, 1] - ψᶠᶠᶜ_face[i, j, 1]) / grid.faces[f].Δxᶜᶠᵃ[i, j] - end - end - - return uᶠᶜᶜ, vᶜᶠᶜ, ψᶠᶠᶜ -end - -function set_velocities_from_streamfunction!(u, v, ψ, arch, grid) - uψ, vψ, ψ₀ = diagnose_velocities_from_streamfunction(ψ, arch, grid) - set!(u, uψ) - set!(v, vψ) - return nothing -end - -for arch in archs - - # These tests cause an undefined `Bound Access Error` on GPU's CI with the new CUDA version. - # The error is not reproducible neither on Tartarus nor on Sverdrup. - # These are excised for the moment (PR #2253) as Cubed sphere will be reworked - if !(arch isa GPU) - @testset "Cubed sphere circulation [$(typeof(arch))]" begin - @info " Testing cubed sphere circulation [$(typeof(arch))]..." - - grid = ConformalCubedSphereGrid(cs32_filepath, arch, Nz=1, z=(-1, 0)) - - FT = eltype(grid) - Nx, Ny, Nz, Nf = size(grid) - R = grid.faces[1].radius - - u_field = XFaceField(grid) - v_field = YFaceField(grid) - - ψ(λ, φ) = R * sind(φ) - @allowscalar set_velocities_from_streamfunction!(u_field, v_field, ψ, arch, grid) - - fill_horizontal_velocity_halos!(u_field, v_field, arch) - - grid_faces = [grid.faces[f] for f in 1:Nf] - u_faces = [get_face(u_field, f) for f in 1:Nf] - v_faces = [get_face(v_field, f) for f in 1:Nf] - - circulation(i, j, f) = Γᶠᶠᶜ(i, j, 1, grid_faces[f], u_faces[f], v_faces[f]) - vorticity(i, j, f) = ζ₃ᶠᶠᶜ(i, j, 1, grid_faces[f], u_faces[f], v_faces[f]) - - CUDA.allowscalar(true) - - @testset "Face 1-2 boundary" begin - # Quick test at a single point. - Γ1 = circulation(33, 16, 1) - Γ2 = circulation(1, 16, 2) - @test Γ1 == Γ2 - - # Check the entire boundary. - Γ1 = [circulation(33, j, 1) for j in 1:Ny+1] - Γ2 = [circulation(1, j, 2) for j in 1:Ny+1] - - # There should be no truncation errors here. - @test Γ1 == Γ2 - end - - @testset "Face 1-3 boundary" begin - Γ1 = circulation(16, 33, 1) - Γ3 = circulation(1, 16, 3) - - # We expect a truncation error of around 3ϵ here. - ϵ = eps(maximum(abs, [Γ1, Γ3])) - @test isapprox(Γ1, Γ3, atol=3ϵ) - - # ffc locations are shifted by one along this boundary (and go in reverse). - Γ1 = [circulation(i, 33, 1) for i in 1:Nx] - Γ3 = [circulation(1, j, 3) for j in Ny+1:-1:2] - - # We expect truncation errors at this boundary. - ϵ = eps(maximum(abs, [Γ1; Γ3])) - @test all(isapprox.(Γ1, Γ3, atol=4ϵ)) - end - - @testset "Face 1-5 boundary" begin - Γ1 = [circulation(1, j, 1) for j in 1:Ny] - Γ5 = [circulation(i, 33, 5) for i in Nx+1:-1:2] - - # Hmmm there's more truncation error here. - ϵ = eps(maximum(abs, [Γ1; Γ5])) - @test all(isapprox.(Γ1, Γ5, atol=32ϵ)) - end - - @testset "Face 1-6 boundary" begin - Γ1 = [circulation(i, 1, 1) for i in 1:Nx+1] - Γ6 = [circulation(i, 33, 6) for i in 1:Nx+1] - @test Γ1 == Γ6 - end - - @testset "Face 1-2-3 corner" begin - Γ1 = circulation(33, 33, 1) - Γ2 = circulation(1, 33, 2) - Γ3 = circulation(1, 1, 3) - @test Γ1 == Γ2 == Γ3 - - ζ1 = vorticity(33, 33, 1) - ζ2 = vorticity(1, 33, 2) - ζ3 = vorticity(1, 1, 3) - @test ζ1 == ζ2 == ζ3 - end - - @testset "Face 1-3-5 corner" begin - Γ1, Γ3, Γ5 = Γ = [circulation(1, 33, f) for f in (1, 3, 5)] - ϵ = eps(maximum(abs, Γ)) - @test isapprox(Γ1, Γ3, atol=32ϵ) - @test isapprox(Γ1, Γ5, atol=32ϵ) - @test isapprox(Γ3, Γ5, atol=32ϵ) - end - - @testset "Face 1-5-6 corner" begin - Γ1 = circulation(1, 1, 1) - Γ5 = circulation(33, 33, 5) - Γ6 = circulation(1, 33, 6) - @test Γ1 == Γ5 == Γ6 - end - - @testset "Face 1-2-6 corner" begin - Γ1 = circulation(33, 1, 1) - Γ2 = circulation(1, 1, 2) - Γ6 = circulation(33, 33, 6) - @test Γ1 == Γ2 == Γ6 - end - - @testset "Face 2-3-4 corner" begin - Γ2 = circulation(33, 33, 2) - Γ3 = circulation(33, 1, 3) - Γ4 = circulation(1, 1, 4) - Γ = [Γ2, Γ3, Γ4] - - ϵ = eps(maximum(abs, Γ)) - @test isapprox(Γ2, Γ3, atol=32ϵ) - @test isapprox(Γ2, Γ4, atol=32ϵ) - @test isapprox(Γ3, Γ4, atol=32ϵ) - end - - @testset "Face 2-4-6 corner" begin - Γ2, Γ4, Γ6 = Γ = [circulation(33, 1, f) for f in (2, 4, 6)] - ϵ = eps(maximum(abs, Γ)) - @test isapprox(Γ2, Γ4, atol=32ϵ) - @test isapprox(Γ2, Γ6, atol=32ϵ) - @test isapprox(Γ4, Γ6, atol=32ϵ) - end - - CUDA.allowscalar(false) - end - end -end diff --git a/test/test_cubed_sphere_halo_exchange.jl b/test/test_cubed_sphere_halo_exchange.jl deleted file mode 100644 index ff05acc7d8..0000000000 --- a/test/test_cubed_sphere_halo_exchange.jl +++ /dev/null @@ -1,676 +0,0 @@ -using Oceananigans.CubedSpheres -using Oceananigans.BoundaryConditions: fill_halo_regions! -using Oceananigans.CubedSpheres: west_halo, east_halo, south_halo, north_halo, fill_horizontal_velocity_halos!, get_face - -# Opposite of the `Base.digits` function -# Source: https://stackoverflow.com/a/55529778 -function undigits(d; base=10) - (s, b) = promote(zero(eltype(d)), base) - mult = one(s) - for val in d - s += val * mult - mult *= b - end - return s -end - -cs32_filepath = datadep"cubed_sphere_32_grid/cubed_sphere_32_grid.jld2" - -for arch in archs - - # These tests cause an undefined `Bound Access Error` on GPU's CI with the new CUDA version. - # The error is not reproducible neither on Tartarus nor on Sverdrup. - # These are excised for the moment (PR #2253) as Cubed sphere will be reworked - if !(arch isa GPU) - @testset "Cubed sphere tracer halo exchange [$(typeof(arch))]" begin - @info " Testing cubed sphere tracer halo exchange [$(typeof(arch))]..." - - grid = ConformalCubedSphereGrid(cs32_filepath, arch, Nz=1, z=(-1, 0)) - field = CenterField(grid) - - ## We will fill each grid point with a 5-digit integer "fiijj" where - ## the f digit is the face number, the ii digits are the i index, and - ## the jj digits are the j index. We then check that the halo exchange - ## happened correctly. - - face_digit(n) = digits(abs(Int(n)))[5] - i_digits(n) = digits(abs(Int(n)))[3:4] |> undigits - j_digits(n) = digits(abs(Int(n)))[1:2] |> undigits - - for (face_number, field_face) in enumerate(faces(field)) - for j in 1:field_face.grid.Ny, i in 1:field_face.grid.Nx, - @allowscalar field_face[i, j, 1] = parse(Int, @sprintf("%d%02d%02d", face_number, i, j)) - end - end - - fill_halo_regions!(field) - - allowscalar(true) - - @testset "Source and destination faces are correct" begin - for (face_number, field_face) in enumerate(faces(field)) - west_halo_vals = west_halo(field_face, include_corners=false) |> Array - east_halo_vals = east_halo(field_face, include_corners=false) |> Array - south_halo_vals = south_halo(field_face, include_corners=false) |> Array - north_halo_vals = north_halo(field_face, include_corners=false) |> Array - - @test all(face_digit.(west_halo_vals) .== grid.face_connectivity[face_number].west.face) - @test all(face_digit.(east_halo_vals) .== grid.face_connectivity[face_number].east.face) - @test all(face_digit.(south_halo_vals) .== grid.face_connectivity[face_number].south.face) - @test all(face_digit.(north_halo_vals) .== grid.face_connectivity[face_number].north.face) - end - end - - @testset "1W halo <- 5N boundary halo exchange" begin - # Grid point [i, j] = [0, 1] in 1W halo should be from [i, j] = [32, 32] in 5N boundary. - west_halo_south_value = get_face(field, 1)[0, 1, 1] - @test face_digit(west_halo_south_value) == 5 - @test i_digits(west_halo_south_value) == 32 - @test j_digits(west_halo_south_value) == 32 - - # Grid point [i, j] = [0, 32] in 1W halo should be from [i, j] = [1, 32] in 5N boundary. - west_halo_north_value = get_face(field, 1)[0, 32, 1] - @test face_digit(west_halo_north_value) == 5 - @test i_digits(west_halo_north_value) == 1 - @test j_digits(west_halo_north_value) == 32 - - west_halo_values = west_halo(get_face(field, 1), include_corners=false)[:] |> Array - @test all(face_digit.(west_halo_values) .== 5) - @test all(i_digits.(west_halo_values) .== reverse(1:32)) - @test all(j_digits.(west_halo_values) .== 32) - end - - @testset "1E halo <- 2W boundary halo exchange" begin - # Grid point [i, j] = [33, 1] in 1E halo should be from [i, j] = [1, 1] in 2W boundary. - east_halo_south_value = get_face(field, 1)[33, 1, 1] - @test face_digit(east_halo_south_value) == 2 - @test i_digits(east_halo_south_value) == 1 - @test j_digits(east_halo_south_value) == 1 - - # Grid point [i, j] = [33, 32] in 1E halo should be from [i, j] = [1, 32] in 2W boundary. - east_halo_north_value = get_face(field, 1)[33, 32, 1] - @test face_digit(east_halo_north_value) == 2 - @test i_digits(east_halo_north_value) == 1 - @test j_digits(east_halo_north_value) == 32 - - east_halo_values = east_halo(get_face(field, 1), include_corners=false)[:] |> Array - @test all(face_digit.(east_halo_values) .== 2) - @test all(i_digits.(east_halo_values) .== 1) - @test all(j_digits.(east_halo_values) .== 1:32) - end - - @testset "1S halo <- 6N boundary halo exchange" begin - # Grid point [i, j] = [1, 0] in 1S halo should be from [i, j] = [1, 32] in 6N boundary. - south_halo_west_value = get_face(field, 1)[1, 0, 1] - @test face_digit(south_halo_west_value) == 6 - @test i_digits(south_halo_west_value) == 1 - @test j_digits(south_halo_west_value) == 32 - - # Grid point [i, j] = [32, 0] in 1S halo should be from [i, j] = [32, 32] in 6N boundary. - south_halo_east_value = get_face(field, 1)[32, 0, 1] - @test face_digit(south_halo_east_value) == 6 - @test i_digits(south_halo_east_value) == 32 - @test j_digits(south_halo_east_value) == 32 - - south_halo_values = south_halo(get_face(field, 1), include_corners=false)[:] |> Array - @test all(face_digit.(south_halo_values) .== 6) - @test all(i_digits.(south_halo_values) .== 1:32) - @test all(j_digits.(south_halo_values) .== 32) - end - - @testset "1N halo <- 3W boundary halo exchange" begin - # Grid point [i, j] = [1, 33] in 1N halo should be from [i, j] = [1, 32] in 3W boundary. - north_halo_west_value = get_face(field, 1)[1, 33, 1] - @test face_digit(north_halo_west_value) == 3 - @test i_digits(north_halo_west_value) == 1 - @test j_digits(north_halo_west_value) == 32 - - # Grid point [i, j] = [32, 33] in 1N halo should be from [i, j] = [1, 1] in 3W boundary. - north_halo_east_value = get_face(field, 1)[32, 33, 1] - @test face_digit(north_halo_east_value) == 3 - @test i_digits(north_halo_east_value) == 1 - @test j_digits(north_halo_east_value) == 1 - - north_halo_values = north_halo(get_face(field, 1), include_corners=false)[:] |> Array - @test all(face_digit.(north_halo_values) .== 3) - @test all(i_digits.(north_halo_values) .== 1) - @test all(j_digits.(north_halo_values) .== reverse(1:32)) - end - - @testset "2W halo <- 1E boundary halo exchange" begin - # Grid point [i, j] = [0, 1] in 2W halo should be from [i, j] = [32, 1] in 1E boundary. - west_halo_south_value = get_face(field, 2)[0, 1, 1] - @test face_digit(west_halo_south_value) == 1 - @test i_digits(west_halo_south_value) == 32 - @test j_digits(west_halo_south_value) == 1 - - # Grid point [i, j] = [0, 32] in 2W halo should be from [i, j] = [32, 32] in 1E boundary. - west_halo_north_value = get_face(field, 2)[0, 32, 1] - @test face_digit(west_halo_north_value) == 1 - @test i_digits(west_halo_north_value) == 32 - @test j_digits(west_halo_north_value) == 32 - - west_halo_values = west_halo(get_face(field, 2), include_corners=false)[:] |> Array - @test all(face_digit.(west_halo_values) .== 1) - @test all(i_digits.(west_halo_values) .== 32) - @test all(j_digits.(west_halo_values) .== 1:32) - end - - @testset "2E halo <- 4S boundary halo exchange" begin - # Grid point [i, j] = [33, 1] in 2E halo should be from [i, j] = [32, 1] in 4S boundary. - east_halo_south_value = get_face(field, 2)[33, 1, 1] - @test face_digit(east_halo_south_value) == 4 - @test i_digits(east_halo_south_value) == 32 - @test j_digits(east_halo_south_value) == 1 - - # Grid point [i, j] = [33, 32] in 2E halo should be from [i, j] = [1, 1] in 4S boundary. - east_halo_north_value = get_face(field, 2)[33, 32, 1] - @test face_digit(east_halo_north_value) == 4 - @test i_digits(east_halo_north_value) == 1 - @test j_digits(east_halo_north_value) == 1 - - east_halo_values = east_halo(get_face(field, 2), include_corners=false)[:] |> Array - @test all(face_digit.(east_halo_values) .== 4) - @test all(i_digits.(east_halo_values) .== reverse(1:32)) - @test all(j_digits.(east_halo_values) .== 1) - end - - @testset "2S halo <- 6E boundary halo exchange" begin - # Grid point [i, j] = [1, 0] in 2S halo should be from [i, j] = [32, 32] in 6E boundary. - south_halo_west_value = get_face(field, 2)[1, 0, 1] - @test face_digit(south_halo_west_value) == 6 - @test i_digits(south_halo_west_value) == 32 - @test j_digits(south_halo_west_value) == 32 - - # Grid point [i, j] = [32, 0] in 2S halo should be from [i, j] = [32, 1] in 6E boundary. - south_halo_east_value = get_face(field, 2)[32, 0, 1] - @test face_digit(south_halo_east_value) == 6 - @test i_digits(south_halo_east_value) == 32 - @test j_digits(south_halo_east_value) == 1 - - south_halo_values = south_halo(get_face(field, 2), include_corners=false)[:] |> Array - @test all(face_digit.(south_halo_values) .== 6) - @test all(i_digits.(south_halo_values) .== 32) - @test all(j_digits.(south_halo_values) .== reverse(1:32)) - end - - @testset "2N halo <- 3S boundary halo exchange" begin - # Grid point [i, j] = [1, 33] in 2N halo should be from [i, j] = [1, 1] in 3S boundary. - north_halo_west_value = get_face(field, 2)[1, 33, 1] - @test face_digit(north_halo_west_value) == 3 - @test i_digits(north_halo_west_value) == 1 - @test j_digits(north_halo_west_value) == 1 - - # Grid point [i, j] = [32, 33] in 2N halo should be from [i, j] = [32, 1] in 3S boundary. - north_halo_east_value = get_face(field, 2)[32, 33, 1] - @test face_digit(north_halo_east_value) == 3 - @test i_digits(north_halo_east_value) == 32 - @test j_digits(north_halo_east_value) == 1 - - north_halo_values = north_halo(get_face(field, 2), include_corners=false)[:] |> Array - @test all(face_digit.(north_halo_values) .== 3) - @test all(i_digits.(north_halo_values) .== 1:32) - @test all(j_digits.(north_halo_values) .== 1) - end - - @testset "3W halo <- 1N boundary halo exchange" begin - # Grid point [i, j] = [0, 1] in 3W halo should be from [i, j] = [32, 32] in 1N boundary. - west_halo_south_value = get_face(field, 3)[0, 1, 1] - @test face_digit(west_halo_south_value) == 1 - @test i_digits(west_halo_south_value) == 32 - @test j_digits(west_halo_south_value) == 32 - - # Grid point [i, j] = [0, 32] in 3W halo should be from [i, j] = [1, 32] in 1N boundary. - west_halo_north_value = get_face(field, 3)[0, 32, 1] - @test face_digit(west_halo_north_value) == 1 - @test i_digits(west_halo_north_value) == 1 - @test j_digits(west_halo_north_value) == 32 - - west_halo_values = west_halo(get_face(field, 3), include_corners=false)[:] |> Array - @test all(face_digit.(west_halo_values) .== 1) - @test all(i_digits.(west_halo_values) .== reverse(1:32)) - @test all(j_digits.(west_halo_values) .== 32) - end - - @testset "3E halo <- 4W boundary halo exchange" begin - # Grid point [i, j] = [33, 1] in 3E halo should be from [i, j] = [1, 1] in 4W boundary. - east_halo_south_value = get_face(field, 3)[33, 1, 1] - @test face_digit(east_halo_south_value) == 4 - @test i_digits(east_halo_south_value) == 1 - @test j_digits(east_halo_south_value) == 1 - - # Grid point [i, j] = [33, 32] in 3E halo should be from [i, j] = [1, 32] in 4W boundary. - east_halo_north_value = get_face(field, 3)[33, 32, 1] - @test face_digit(east_halo_north_value) == 4 - @test i_digits(east_halo_north_value) == 1 - @test j_digits(east_halo_north_value) == 32 - - east_halo_values = east_halo(get_face(field, 3), include_corners=false)[:] |> Array - @test all(face_digit.(east_halo_values) .== 4) - @test all(i_digits.(east_halo_values) .== 1) - @test all(j_digits.(east_halo_values) .== 1:32) - end - - @testset "3S halo <- 2N boundary halo exchange" begin - # Grid point [i, j] = [1, 0] in 3S halo should be from [i, j] = [1, 32] in 2N boundary. - south_halo_west_value = get_face(field, 3)[1, 0, 1] - @test face_digit(south_halo_west_value) == 2 - @test i_digits(south_halo_west_value) == 1 - @test j_digits(south_halo_west_value) == 32 - - # Grid point [i, j] = [32, 0] in 3S halo should be from [i, j] = [32, 32] in 2N boundary. - south_halo_east_value = get_face(field, 3)[32, 0, 1] - @test face_digit(south_halo_east_value) == 2 - @test i_digits(south_halo_east_value) == 32 - @test j_digits(south_halo_east_value) == 32 - - south_halo_values = south_halo(get_face(field, 3), include_corners=false)[:] |> Array - @test all(face_digit.(south_halo_values) .== 2) - @test all(i_digits.(south_halo_values) .== 1:32) - @test all(j_digits.(south_halo_values) .== 32) - end - - @testset "3N halo <- 5W boundary halo exchange" begin - # Grid point [i, j] = [1, 33] in 3N halo should be from [i, j] = [1, 32] in 5W boundary. - north_halo_west_value = get_face(field, 3)[1, 33, 1] - @test face_digit(north_halo_west_value) == 5 - @test i_digits(north_halo_west_value) == 1 - @test j_digits(north_halo_west_value) == 32 - - # Grid point [i, j] = [32, 33] in 3N halo should be from [i, j] = [1, 1] in 5W boundary. - north_halo_east_value = get_face(field, 3)[32, 33, 1] - @test face_digit(north_halo_east_value) == 5 - @test i_digits(north_halo_east_value) == 1 - @test j_digits(north_halo_east_value) == 1 - - north_halo_values = north_halo(get_face(field, 3), include_corners=false)[:] |> Array - @test all(face_digit.(north_halo_values) .== 5) - @test all(i_digits.(north_halo_values) .== 1) - @test all(j_digits.(north_halo_values) .== reverse(1:32)) - end - - @testset "4W halo <- 3E boundary halo exchange" begin - # Grid point [i, j] = [0, 1] in 4W halo should be from [i, j] = [32, 1] in 3E boundary. - west_halo_south_value = get_face(field, 4)[0, 1, 1] - @test face_digit(west_halo_south_value) == 3 - @test i_digits(west_halo_south_value) == 32 - @test j_digits(west_halo_south_value) == 1 - - # Grid point [i, j] = [0, 32] in 4W halo should be from [i, j] = [32, 32] in 1N boundary. - west_halo_north_value = get_face(field, 4)[0, 32, 1] - @test face_digit(west_halo_north_value) == 3 - @test i_digits(west_halo_north_value) == 32 - @test j_digits(west_halo_north_value) == 32 - - west_halo_values = west_halo(get_face(field, 4), include_corners=false)[:] |> Array - @test all(face_digit.(west_halo_values) .== 3) - @test all(i_digits.(west_halo_values) .== 32) - @test all(j_digits.(west_halo_values) .== 1:32) - end - - @testset "4E halo <- 6S boundary halo exchange" begin - # Grid point [i, j] = [33, 1] in 4E halo should be from [i, j] = [32, 1] in 6S boundary. - east_halo_south_value = get_face(field, 4)[33, 1, 1] - @test face_digit(east_halo_south_value) == 6 - @test i_digits(east_halo_south_value) == 32 - @test j_digits(east_halo_south_value) == 1 - - # Grid point [i, j] = [33, 32] in 4E halo should be from [i, j] = [1, 1] in 6S boundary. - east_halo_north_value = get_face(field, 4)[33, 32, 1] - @test face_digit(east_halo_north_value) == 6 - @test i_digits(east_halo_north_value) == 1 - @test j_digits(east_halo_north_value) == 1 - - east_halo_values = east_halo(get_face(field, 4), include_corners=false)[:] |> Array - @test all(face_digit.(east_halo_values) .== 6) - @test all(i_digits.(east_halo_values) .== reverse(1:32)) - @test all(j_digits.(east_halo_values) .== 1) - end - - @testset "4S halo <- 2E boundary halo exchange" begin - # Grid point [i, j] = [1, 0] in 4S halo should be from [i, j] = [32, 32] in 2E boundary. - south_halo_west_value = get_face(field, 4)[1, 0, 1] - @test face_digit(south_halo_west_value) == 2 - @test i_digits(south_halo_west_value) == 32 - @test j_digits(south_halo_west_value) == 32 - - # Grid point [i, j] = [32, 0] in 4S halo should be from [i, j] = [32, 1] in 2E boundary. - south_halo_east_value = get_face(field, 4)[32, 0, 1] - @test face_digit(south_halo_east_value) == 2 - @test i_digits(south_halo_east_value) == 32 - @test j_digits(south_halo_east_value) == 1 - - south_halo_values = south_halo(get_face(field, 4), include_corners=false)[:] |> Array - @test all(face_digit.(south_halo_values) .== 2) - @test all(i_digits.(south_halo_values) .== 32) - @test all(j_digits.(south_halo_values) .== reverse(1:32)) - end - - @testset "4N halo <- 5S boundary halo exchange" begin - # Grid point [i, j] = [1, 33] in 4N halo should be from [i, j] = [1, 1] in 5S boundary. - north_halo_west_value = get_face(field, 4)[1, 33, 1] - @test face_digit(north_halo_west_value) == 5 - @test i_digits(north_halo_west_value) == 1 - @test j_digits(north_halo_west_value) == 1 - - # Grid point [i, j] = [32, 33] in 4N halo should be from [i, j] = [32, 1] in 5S boundary. - north_halo_east_value = get_face(field, 4)[32, 33, 1] - @test face_digit(north_halo_east_value) == 5 - @test i_digits(north_halo_east_value) == 32 - @test j_digits(north_halo_east_value) == 1 - - north_halo_values = north_halo(get_face(field, 4), include_corners=false)[:] |> Array - @test all(face_digit.(north_halo_values) .== 5) - @test all(i_digits.(north_halo_values) .== 1:32) - @test all(j_digits.(north_halo_values) .== 1) - end - - ## TODO: Test faces 5 and 6 (or generalize the test to all faces)! - - CUDA.allowscalar(false) - - end - end -end - -for arch in archs - - # These tests cause an undefined `Bound Access Error` on GPU's CI with the new CUDA version. - # The error is not reproducible neither on Tartarus nor on Sverdrup. - # These are excised for the moment (PR #2253) as Cubed sphere will be reworked - if !(arch isa GPU) - @testset "Cubed sphere velocity halo exchange [$(typeof(arch))]" begin - @info " Testing cubed sphere velocity halo exchange [$(typeof(arch))]..." - - grid = ConformalCubedSphereGrid(cs32_filepath, arch, Nz=1, z=(-1, 0)) - - u_field = XFaceField(grid) - v_field = YFaceField(grid) - - ## We will fill each grid point with a 6-digit integer "ufiijj" where - ## the u digit is 1 for u and 2 for v, the f digit is the face number, - ## the ii digits are the i index, and the jj digits are the j index. - ## We then check that the halo exchange happened correctly (including - ## any sign changes). - - U_DIGIT = 1 - V_DIGIT = 2 - - uv_digit(n) = digits(abs(Int(n)))[6] - face_digit(n) = digits(abs(Int(n)))[5] - i_digits(n) = digits(abs(Int(n)))[3:4] |> undigits - j_digits(n) = digits(abs(Int(n)))[1:2] |> undigits - - for (face_number, u_field_face) in enumerate(faces(u_field)) - for j in 1:u_field_face.grid.Ny, i in 1:u_field_face.grid.Nx+1 - @allowscalar u_field_face[i, j, 1] = parse(Int, @sprintf("%d%d%02d%02d", U_DIGIT, face_number, i, j)) - end - end - - for (face_number, v_field_face) in enumerate(faces(v_field)) - for j in 1:v_field_face.grid.Ny+1, i in 1:v_field_face.grid.Nx - @allowscalar v_field_face[i, j, 1] = parse(Int, @sprintf("%d%d%02d%02d", V_DIGIT, face_number, i, j)) - end - end - - fill_horizontal_velocity_halos!(u_field, v_field, arch) - - CUDA.allowscalar(true) - - # @testset "Source and destination faces are correct" begin - # for field in (u_field, v_field), (face_number, field_face) in enumerate(faces(field)) - # west_halo_vals = west_halo(field_face, include_corners=false) - # east_halo_vals = east_halo(field_face, include_corners=false) - # south_halo_vals = south_halo(field_face, include_corners=false) - # north_halo_vals = north_halo(field_face, include_corners=false) - - # @test all(face_digit.(west_halo_vals) .== grid.face_connectivity[face_number].west.face) - # @test all(face_digit.(east_halo_vals) .== grid.face_connectivity[face_number].east.face) - # @test all(face_digit.(south_halo_vals) .== grid.face_connectivity[face_number].south.face) - # @test all(face_digit.(north_halo_vals) .== grid.face_connectivity[face_number].north.face) - # end - # end - - @testset "(1W halo <- 5N boundary) horizontal velocity halo exchange" begin - # Grid point u[i, j] = u[0, 1] in 1W halo should be from +v[i, j] = +v[32, 32] in 5N boundary. - u_west_halo_south_value = get_face(u_field, 1)[0, 1, 1] - @test uv_digit(u_west_halo_south_value) == V_DIGIT - @test sign(u_west_halo_south_value) == +1 - @test face_digit(u_west_halo_south_value) == 5 - @test i_digits(u_west_halo_south_value) == 32 - @test j_digits(u_west_halo_south_value) == 32 - - # Grid point u[i, j] = u[0, 32] in 1W halo should be from +v[i, j] = +v[1, 32] in 5N boundary. - u_west_halo_north_value = get_face(u_field, 1)[0, 32, 1] - @test uv_digit(u_west_halo_north_value) == V_DIGIT - @test sign(u_west_halo_north_value) == +1 - @test face_digit(u_west_halo_north_value) == 5 - @test i_digits(u_west_halo_north_value) == 1 - @test j_digits(u_west_halo_north_value) == 32 - - u_west_halo_values = west_halo(get_face(u_field, 1), include_corners=false)[:] |> Array - @test all(uv_digit.(u_west_halo_values) .== V_DIGIT) - @test all(sign.(u_west_halo_values) .== +1) - @test all(face_digit.(u_west_halo_values) .== 5) - @test all(i_digits.(u_west_halo_values) .== reverse(1:32)) - @test all(j_digits.(u_west_halo_values) .== 32) - - # Note: The 1W v halo <- 5N -u boundary exchange involves a sign change and an index shift. - - # Grid point v[i, j] = v[0, 1] in 1W halo should be from -u[i, j] = -u[1, 32] in 6NE corner! - v_west_halo_south_value = get_face(v_field, 1)[0, 1, 1] - @test uv_digit(v_west_halo_south_value) == U_DIGIT - @test sign(v_west_halo_south_value) == -1 - @test face_digit(v_west_halo_south_value) == 6 - @test i_digits(v_west_halo_south_value) == 1 - @test j_digits(v_west_halo_south_value) == 32 - - # Grid point v[i, j] = u[0, 33] in 1W halo should be from -u[i, j] = -u[2, 32] in 5N boundary. - v_west_halo_north_value = get_face(v_field, 1)[0, 32, 1] - @test uv_digit(v_west_halo_north_value) == U_DIGIT - @test sign(v_west_halo_north_value) == -1 - @test face_digit(v_west_halo_north_value) == 5 - @test i_digits(v_west_halo_north_value) == 2 - @test j_digits(v_west_halo_north_value) == 32 - - v_west_halo_values = west_halo(get_face(v_field, 1), include_corners=false)[:] |> Array - @test all(uv_digit.(v_west_halo_values) .== U_DIGIT) - @test all(sign.(v_west_halo_values) .== -1) - # @test all(face_digit.(v_west_halo_values) .== 5) - # @test all(i_digits.(v_west_halo_values) .== reverse(2:33)) - @test all(j_digits.(v_west_halo_values) .== 32) - end - - @testset "(1E halo <- 2W boundary) horizontal velocity halo exchange" begin - # Grid point u[i, j] = u[33, 1] in 1E halo should be from +u[i, j] = +u[1, 1] in 2W boundary. - u_east_halo_south_value = get_face(u_field, 1)[33, 1, 1] - @test uv_digit(u_east_halo_south_value) == U_DIGIT - @test sign(u_east_halo_south_value) == +1 - @test face_digit(u_east_halo_south_value) == 2 - @test i_digits(u_east_halo_south_value) == 1 - @test j_digits(u_east_halo_south_value) == 1 - - # Grid point u[i, j] = u[33, 32] in 1E halo should be from +u[i, j] = +u[1, 32] in 2W boundary. - u_east_halo_north_value = get_face(u_field, 1)[33, 32, 1] - @test uv_digit(u_east_halo_north_value) == U_DIGIT - @test sign(u_east_halo_north_value) == +1 - @test face_digit(u_east_halo_north_value) == 2 - @test i_digits(u_east_halo_north_value) == 1 - @test j_digits(u_east_halo_north_value) == 32 - - u_east_halo_values = east_halo(get_face(u_field, 1), include_corners=false)[:] |> Array - @test all(uv_digit.(u_east_halo_values) .== U_DIGIT) - @test all(sign.(u_east_halo_values) .== +1) - @test all(face_digit.(u_east_halo_values) .== 2) - @test all(i_digits.(u_east_halo_values) .== 1) - @test all(j_digits.(u_east_halo_values) .== 1:32) - - # Grid point v[i, j] = v[33, 1] in 1E halo should be from +v[i, j] = +v[1, 1] in 2W boundary. - v_east_halo_south_value = get_face(v_field, 1)[33, 1, 1] - @test uv_digit(v_east_halo_south_value) == V_DIGIT - @test sign(v_east_halo_south_value) == +1 - @test face_digit(v_east_halo_south_value) == 2 - @test i_digits(v_east_halo_south_value) == 1 - @test j_digits(v_east_halo_south_value) == 1 - - # # Grid point v[i, j] = v[33, 32] in 1E halo should be from +v[i, j] = +v[1, 32] in 2W boundary. - v_east_halo_north_value = get_face(v_field, 1)[33, 32, 1] - @test uv_digit(v_east_halo_north_value) == V_DIGIT - @test sign(v_east_halo_north_value) == +1 - @test face_digit(v_east_halo_north_value) == 2 - @test i_digits(v_east_halo_north_value) == 1 - @test j_digits(v_east_halo_north_value) == 32 - - v_east_halo_values = east_halo(get_face(v_field, 1), include_corners=false)[:] |> Array - @test all(uv_digit.(v_east_halo_values) .== V_DIGIT) - @test all(sign.(v_east_halo_values) .== +1) - @test all(face_digit.(v_east_halo_values) .== 2) - @test all(i_digits.(v_east_halo_values) .== 1) - @test all(j_digits.(v_east_halo_values) .== 1:32) - end - - @testset "(1S halo <- 6N boundary) horizontal velocity halo exchange" begin - # Grid point u[i, j] = u[1, 0] in 1S halo should be from +u[i, j] = +u[1, 32] in 6N boundary. - u_south_halo_west_value = get_face(u_field, 1)[1, 0, 1] - @test uv_digit(u_south_halo_west_value) == U_DIGIT - @test sign(u_south_halo_west_value) == +1 - @test face_digit(u_south_halo_west_value) == 6 - @test i_digits(u_south_halo_west_value) == 1 - @test j_digits(u_south_halo_west_value) == 32 - - # Grid point u[i, j] = u[32, 0] in 1S halo should be from +u[i, j] = +u[32, 32] in 6N boundary. - u_south_halo_east_value = get_face(u_field, 1)[32, 0, 1] - @test uv_digit(u_south_halo_east_value) == U_DIGIT - @test sign(u_south_halo_east_value) == +1 - @test face_digit(u_south_halo_east_value) == 6 - @test i_digits(u_south_halo_east_value) == 32 - @test j_digits(u_south_halo_east_value) == 32 - - u_south_halo_values = south_halo(get_face(u_field, 1), include_corners=false)[:] |> Array - @test all(uv_digit.(u_south_halo_values) .== U_DIGIT) - @test all(sign.(u_south_halo_values) .== +1) - @test all(face_digit.(u_south_halo_values) .== 6) - @test all(i_digits.(u_south_halo_values) .== 1:32) - @test all(j_digits.(u_south_halo_values) .== 32) - - # Grid point v[i, j] = v[1, 0] in 1S halo should be from +v[i, j] = +v[1, 32] in 6N boundary. - v_south_halo_west_value = get_face(v_field, 1)[1, 0, 1] - @test uv_digit(v_south_halo_west_value) == V_DIGIT - @test sign(v_south_halo_west_value) == +1 - @test face_digit(v_south_halo_west_value) == 6 - @test i_digits(v_south_halo_west_value) == 1 - @test j_digits(v_south_halo_west_value) == 32 - - # # Grid point v[i, j] = v[32, 0] in 1S halo should be from +v[i, j] = +v[32, 32] in 6N boundary. - v_south_halo_east_value = get_face(v_field, 1)[32, 0, 1] - @test uv_digit(v_south_halo_east_value) == V_DIGIT - @test sign(v_south_halo_east_value) == +1 - @test face_digit(v_south_halo_east_value) == 6 - @test i_digits(v_south_halo_east_value) == 32 - @test j_digits(v_south_halo_east_value) == 32 - - v_south_halo_values = south_halo(get_face(v_field, 1), include_corners=false)[:] |> Array - @test all(uv_digit.(v_south_halo_values) .== V_DIGIT) - @test all(sign.(v_south_halo_values) .== +1) - @test all(face_digit.(v_south_halo_values) .== 6) - @test all(i_digits.(v_south_halo_values) .== 1:32) - @test all(j_digits.(v_south_halo_values) .== 32) - end - - @testset "(1N halo <- 3W boundary) horizontal velocity halo exchange" begin - - # Note: The 1N u halo <- 3W -v boundary exchange involves a sign change and an index shift. - - # Grid point u[i, j] = u[1, 33] in 1N halo should be from -v[i, j] = -u[1, 32] in 5NE corner! - u_nouth_halo_west_value = get_face(u_field, 1)[1, 33, 1] - @test uv_digit(u_nouth_halo_west_value) == U_DIGIT - @test sign(u_nouth_halo_west_value) == -1 - @test face_digit(u_nouth_halo_west_value) == 5 - @test i_digits(u_nouth_halo_west_value) == 1 - @test j_digits(u_nouth_halo_west_value) == 32 - - # Grid point u[i, j] = u[32, 33] in 1N halo should be from -v[i, j] = -v[1, 2] in 3W boundary. - u_north_halo_east_value = get_face(u_field, 1)[32, 33, 1] - @test uv_digit(u_north_halo_east_value) == V_DIGIT - @test sign(u_north_halo_east_value) == -1 - @test face_digit(u_north_halo_east_value) == 3 - @test i_digits(u_north_halo_east_value) == 1 - @test j_digits(u_north_halo_east_value) == 2 - - u_north_halo_values = north_halo(get_face(u_field, 1), include_corners=false)[:] |> Array - # @test all(uv_digit.(u_north_halo_values) .== V_DIGIT) - @test all(sign.(u_north_halo_values) .== -1) - # @test all(face_digit.(u_north_halo_values) .== 3) - @test all(i_digits.(u_north_halo_values) .== 1) - # @test all(j_digits.(u_north_halo_values) .== reverse(2:33)) - - # Grid point v[i, j] = v[1, 33] in 1N halo should be from +u[i, j] = +u[1, 32] in 3W boundary. - v_north_halo_west_value = get_face(v_field, 1)[1, 33, 1] - @test uv_digit(v_north_halo_west_value) == U_DIGIT - @test sign(v_north_halo_west_value) == +1 - @test face_digit(v_north_halo_west_value) == 3 - @test i_digits(v_north_halo_west_value) == 1 - @test j_digits(v_north_halo_west_value) == 32 - - # # Grid point v[i, j] = v[32, 33] in 1N halo should be from +u[i, j] = +u[1, 1] in 3W boundary. - v_north_halo_east_value = get_face(v_field, 1)[32, 33, 1] - @test uv_digit(v_north_halo_east_value) == U_DIGIT - @test sign(v_north_halo_east_value) == +1 - @test face_digit(v_north_halo_east_value) == 3 - @test i_digits(v_north_halo_east_value) == 1 - @test j_digits(v_north_halo_east_value) == 1 - - v_north_halo_values = north_halo(get_face(v_field, 1), include_corners=false)[:] |> Array - @test all(uv_digit.(v_north_halo_values) .== U_DIGIT) - @test all(sign.(v_north_halo_values) .== +1) - @test all(face_digit.(v_north_halo_values) .== 3) - @test all(i_digits.(v_north_halo_values) .== 1) - @test all(j_digits.(v_north_halo_values) .== reverse(1:32)) - end - - ## TODO: Test the other faces, especially face 2? I hope we can generalize this... - - @testset "Velocities at the corners (MITgcm regression)" begin - # Face 1 - @test get_face(u_field, 1)[1, 33, 1] == -150132 - @test get_face(u_field, 1)[33, 0, 1] == 220101 - - @test get_face(v_field, 1)[0, 1, 1] == -160132 - @test get_face(v_field, 1)[0, 33, 1] == -150132 - @test get_face(v_field, 1)[33, 33, 1] == 230101 - - # Face 2 - @test get_face(u_field, 2)[1, 0, 1] == -213201 - @test get_face(u_field, 2)[33, 33, 1] == 140101 - - @test get_face(v_field, 2)[33, 1, 1] == -263201 - @test get_face(v_field, 2)[0, 33, 1] == 130101 - @test get_face(v_field, 2)[33, 33, 1] == -140101 - - # Face 3 - @test get_face(u_field, 3)[1, 33, 1] == -110132 - - @test get_face(v_field, 3)[0, 33, 1] == -110132 - @test get_face(v_field, 3)[33, 33, 1] == 250101 - - # Face 4 - @test get_face(u_field, 4)[33, 33, 1] == 160101 - - @test get_face(v_field, 4)[33, 1, 1] == -223201 - @test get_face(v_field, 4)[0, 33, 1] == 150101 - @test get_face(v_field, 4)[33, 33, 1] == -160101 - - # Face 5 - @test get_face(v_field, 5)[0, 33, 1] == -130132 - @test get_face(v_field, 5)[33, 33, 1] == 210101 - - # Face 6 - @test get_face(v_field, 6)[0, 33, 1] == 110101 - @test get_face(v_field, 6)[33, 33, 1] == -120101 - end - - CUDA.allowscalar(false) - - end - end -end diff --git a/test/test_multi_region_cubed_sphere.jl b/test/test_multi_region_cubed_sphere.jl index 2c8999045a..a7c87f5432 100644 --- a/test/test_multi_region_cubed_sphere.jl +++ b/test/test_multi_region_cubed_sphere.jl @@ -247,603 +247,693 @@ panel_sizes = ((8, 8, 1), (9, 9, 2)) end end -@testset "Immersed conformal cubed sphere construction" begin - for FT in float_types, arch in archs, non_uniform_conformal_mapping in (false, true) - Nx, Ny, Nz = 9, 9, 9 - cm = non_uniform_conformal_mapping ? "non-uniform conformal mapping" : "uniform conformal mapping" - - @info " Testing immersed conformal cubed sphere grid [$FT, $(typeof(arch)), $cm]..." - - underlying_grid = ConformalCubedSphereGrid(arch, FT; - panel_size = (Nx, Ny, Nz), z = (-1, 0), radius = 1, - non_uniform_conformal_mapping) - @inline bottom(x, y) = ifelse(abs(y) < 30, - 2, 0) - immersed_grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom); active_cells_map = true) - - # Test that the grid is constructed correctly. - for panel in 1:6 - grid = getregion(immersed_grid, panel) - - if panel == 3 || panel == 6 # North and South panels should be completely immersed. - @test isempty(grid.interior_active_cells) - else # Other panels should have some active cells. - @test !isempty(grid.interior_active_cells) +@testset "Immersed cubed sphere construction" begin + for FT in float_types + for arch in archs + for non_uniform_conformal_mapping in (false, true) + A = typeof(arch) + @info " Testing immersed cubed sphere grid [$FT, $A]..." + + Nx, Ny, Nz = 9, 9, 9 + underlying_grid = ConformalCubedSphereGrid(arch, FT; + panel_size = (Nx, Ny, Nz), z = (-1, 0), radius = 1, + non_uniform_conformal_mapping) + + @inline bottom(x, y) = ifelse(abs(y) < 30, - 2, 0) + immersed_grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom); active_cells_map = true) + + # Test that the grid is constructed correctly. + for panel in 1:6 + grid = getregion(immersed_grid, panel) + + if panel == 3 || panel == 6 # North and South panels should be completely immersed. + @test isempty(grid.interior_active_cells) + else # Other panels should have some active cells. + @test !isempty(grid.interior_active_cells) + end + end end end end end @testset "Testing conformal cubed sphere fill halos for tracers" begin - for FT in float_types, arch in archs, non_uniform_conformal_mapping in (false, true) - cm = non_uniform_conformal_mapping ? "non-uniform conformal mapping" : "uniform conformal mapping" - @info " Testing fill halos for tracers [$FT, $(typeof(arch)), $cm]..." - - Nx, Ny, Nz = 9, 9, 1 - - underlying_grid = ConformalCubedSphereGrid(arch, FT; - panel_size = (Nx, Ny, Nz), z = (0, 1), radius = 1, - horizontal_direction_halo = 3, non_uniform_conformal_mapping) - @inline bottom(x, y) = ifelse(abs(y) < 30, - 2, 0) - immersed_grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom); active_cells_map = true) - - grids = (underlying_grid, immersed_grid) - - for grid in grids - c = CenterField(grid) - - region = Iterate(1:6) - @apply_regionally data = create_c_test_data(grid, region) - set!(c, data) - fill_halo_regions!(c) - - Hx, Hy, Hz = halo_size(c.grid) - - west_indices = 1:Hx, 1:Ny - south_indices = 1:Nx, 1:Hy - east_indices = Nx-Hx+1:Nx, 1:Ny - north_indices = 1:Nx, Ny-Hy+1:Ny - - # Confirm that the tracer halos were filled according to connectivity described at ConformalCubedSphereGrid docstring. - @allowscalar begin - @test get_halo_data(getregion(c, 1), West()) == reverse(create_c_test_data(grid, 5)[north_indices...], dims=1)' - @test get_halo_data(getregion(c, 1), East()) == create_c_test_data(grid, 2)[west_indices...] - @test get_halo_data(getregion(c, 1), South()) == create_c_test_data(grid, 6)[north_indices...] - @test get_halo_data(getregion(c, 1), North()) == reverse(create_c_test_data(grid, 3)[west_indices...], dims=2)' - - @test get_halo_data(getregion(c, 2), West()) == create_c_test_data(grid, 1)[east_indices...] - @test get_halo_data(getregion(c, 2), East()) == reverse(create_c_test_data(grid, 4)[south_indices...], dims=1)' - @test get_halo_data(getregion(c, 2), South()) == reverse(create_c_test_data(grid, 6)[east_indices...], dims=2)' - @test get_halo_data(getregion(c, 2), North()) == create_c_test_data(grid, 3)[south_indices...] - - @test get_halo_data(getregion(c, 3), West()) == reverse(create_c_test_data(grid, 1)[north_indices...], dims=1)' - @test get_halo_data(getregion(c, 3), East()) == create_c_test_data(grid, 4)[west_indices...] - @test get_halo_data(getregion(c, 3), South()) == create_c_test_data(grid, 2)[north_indices...] - @test get_halo_data(getregion(c, 3), North()) == reverse(create_c_test_data(grid, 5)[west_indices...], dims=2)' - - @test get_halo_data(getregion(c, 4), West()) == create_c_test_data(grid, 3)[east_indices...] - @test get_halo_data(getregion(c, 4), East()) == reverse(create_c_test_data(grid, 6)[south_indices...], dims=1)' - @test get_halo_data(getregion(c, 4), South()) == reverse(create_c_test_data(grid, 2)[east_indices...], dims=2)' - @test get_halo_data(getregion(c, 4), North()) == create_c_test_data(grid, 5)[south_indices...] - - @test get_halo_data(getregion(c, 5), West()) == reverse(create_c_test_data(grid, 3)[north_indices...], dims=1)' - @test get_halo_data(getregion(c, 5), East()) == create_c_test_data(grid, 6)[west_indices...] - @test get_halo_data(getregion(c, 5), South()) == create_c_test_data(grid, 4)[north_indices...] - @test get_halo_data(getregion(c, 5), North()) == reverse(create_c_test_data(grid, 1)[west_indices...], dims=2)' - - @test get_halo_data(getregion(c, 6), West()) == create_c_test_data(grid, 5)[east_indices...] - @test get_halo_data(getregion(c, 6), East()) == reverse(create_c_test_data(grid, 2)[south_indices...], dims=1)' - @test get_halo_data(getregion(c, 6), South()) == reverse(create_c_test_data(grid, 4)[east_indices...], dims=2)' - @test get_halo_data(getregion(c, 6), North()) == create_c_test_data(grid, 1)[south_indices...] - end # CUDA.@allowscalar + for FT in float_types + for arch in archs + for non_uniform_conformal_mapping in (false, true) + A = typeof(arch) + @info " Testing fill halos for tracers [$FT, $A]..." + + Nx, Ny, Nz = 9, 9, 1 + underlying_grid = ConformalCubedSphereGrid(arch, FT; + panel_size = (Nx, Ny, Nz), z = (0, 1), radius = 1, + horizontal_direction_halo = 3, non_uniform_conformal_mapping) + + @inline bottom(x, y) = ifelse(abs(y) < 30, - 2, 0) + immersed_grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom); active_cells_map = true) + + grids = (underlying_grid, immersed_grid) + + for grid in grids + c = CenterField(grid) + region = Iterate(1:6) + @apply_regionally data = create_c_test_data(grid, region) + set!(c, data) + fill_halo_regions!(c) + + Hx, Hy, Hz = halo_size(c.grid) + + west_indices = 1:Hx, 1:Ny + south_indices = 1:Nx, 1:Hy + east_indices = Nx-Hx+1:Nx, 1:Ny + north_indices = 1:Nx, Ny-Hy+1:Ny + + # Confirm that the tracer halos were filled according to connectivity described at ConformalCubedSphereGrid docstring. + @allowscalar begin + @test get_halo_data(getregion(c, 1), West()) == reverse(create_c_test_data(grid, 5)[north_indices...], dims=1)' + @test get_halo_data(getregion(c, 1), East()) == create_c_test_data(grid, 2)[west_indices...] + @test get_halo_data(getregion(c, 1), South()) == create_c_test_data(grid, 6)[north_indices...] + @test get_halo_data(getregion(c, 1), North()) == reverse(create_c_test_data(grid, 3)[west_indices...], dims=2)' + + @test get_halo_data(getregion(c, 2), West()) == create_c_test_data(grid, 1)[east_indices...] + @test get_halo_data(getregion(c, 2), East()) == reverse(create_c_test_data(grid, 4)[south_indices...], dims=1)' + @test get_halo_data(getregion(c, 2), South()) == reverse(create_c_test_data(grid, 6)[east_indices...], dims=2)' + @test get_halo_data(getregion(c, 2), North()) == create_c_test_data(grid, 3)[south_indices...] + + @test get_halo_data(getregion(c, 3), West()) == reverse(create_c_test_data(grid, 1)[north_indices...], dims=1)' + @test get_halo_data(getregion(c, 3), East()) == create_c_test_data(grid, 4)[west_indices...] + @test get_halo_data(getregion(c, 3), South()) == create_c_test_data(grid, 2)[north_indices...] + @test get_halo_data(getregion(c, 3), North()) == reverse(create_c_test_data(grid, 5)[west_indices...], dims=2)' + + @test get_halo_data(getregion(c, 4), West()) == create_c_test_data(grid, 3)[east_indices...] + @test get_halo_data(getregion(c, 4), East()) == reverse(create_c_test_data(grid, 6)[south_indices...], dims=1)' + @test get_halo_data(getregion(c, 4), South()) == reverse(create_c_test_data(grid, 2)[east_indices...], dims=2)' + @test get_halo_data(getregion(c, 4), North()) == create_c_test_data(grid, 5)[south_indices...] + + @test get_halo_data(getregion(c, 5), West()) == reverse(create_c_test_data(grid, 3)[north_indices...], dims=1)' + @test get_halo_data(getregion(c, 5), East()) == create_c_test_data(grid, 6)[west_indices...] + @test get_halo_data(getregion(c, 5), South()) == create_c_test_data(grid, 4)[north_indices...] + @test get_halo_data(getregion(c, 5), North()) == reverse(create_c_test_data(grid, 1)[west_indices...], dims=2)' + + @test get_halo_data(getregion(c, 6), West()) == create_c_test_data(grid, 5)[east_indices...] + @test get_halo_data(getregion(c, 6), East()) == reverse(create_c_test_data(grid, 2)[south_indices...], dims=1)' + @test get_halo_data(getregion(c, 6), South()) == reverse(create_c_test_data(grid, 4)[east_indices...], dims=2)' + @test get_halo_data(getregion(c, 6), North()) == create_c_test_data(grid, 1)[south_indices...] + end # CUDA.@allowscalar + end + end end end end @testset "Testing conformal cubed sphere fill halos for horizontal velocities" begin - for FT in float_types, arch in archs, non_uniform_conformal_mapping in (false, true) - cm = non_uniform_conformal_mapping ? "non-uniform conformal mapping" : "uniform conformal mapping" - @info " Testing fill halos for horizontal velocities [$FT, $(typeof(arch)), $cm]..." - - Nx, Ny, Nz = 9, 9, 1 - - underlying_grid = ConformalCubedSphereGrid(arch, FT; - panel_size = (Nx, Ny, Nz), z = (0, 1), radius = 1, - horizontal_direction_halo = 3, non_uniform_conformal_mapping) - @inline bottom(x, y) = ifelse(abs(y) < 30, - 2, 0) - immersed_grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom); active_cells_map = true) - - grids = (underlying_grid, immersed_grid) - - for grid in grids - u = XFaceField(grid) - v = YFaceField(grid) - - region = Iterate(1:6) - @apply_regionally u_data = create_u_test_data(grid, region) - @apply_regionally v_data = create_v_test_data(grid, region) - set!(u, u_data) - set!(v, v_data) - - fill_halo_regions!((u, v); signed = true) - - Hx, Hy, Hz = halo_size(u.grid) - - south_indices = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=nothing, index=:all) - east_indices = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=nothing, index=:all) - north_indices = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=nothing, index=:all) - west_indices = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=nothing, index=:all) - - south_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:endpoint, index=:first) - south_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:endpoint, index=:last) - east_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:endpoint, index=:first) - east_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:endpoint, index=:last) - north_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:endpoint, index=:first) - north_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:endpoint, index=:last) - west_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:endpoint, index=:first) - west_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:endpoint, index=:last) - - south_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:subset, index=:first) - south_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:subset, index=:last) - east_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:subset, index=:first) - east_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:subset, index=:last) - north_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:subset, index=:first) - north_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:subset, index=:last) - west_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:subset, index=:first) - west_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:subset, index=:last) - - # Confirm that the zonal velocity halos were filled according to connectivity described at ConformalCubedSphereGrid docstring. - @allowscalar begin - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(u, 1), West()) == reverse(create_v_test_data(grid, 5)[north_indices...], dims=1)' - @test get_halo_data(getregion(u, 1), East()) == create_u_test_data(grid, 2)[west_indices...] - @test get_halo_data(getregion(u, 1), South()) == create_u_test_data(grid, 6)[north_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(u, 1), North(); - operation=:subset, - index=:first) == - reverse(create_v_test_data(grid, 3)[west_indices_subset_skip_first_index...], dims=2)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(u, 1), North(); - operation=:endpoint, - index=:first) == - reverse(create_u_test_data(grid, 5)[north_indices_first...]) - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(u, 2), West()) == create_u_test_data(grid, 1)[east_indices...] - @test get_halo_data(getregion(u, 2), East()) == reverse(create_v_test_data(grid, 4)[south_indices...], dims=1)' - @test get_halo_data(getregion(u, 2), North()) == create_u_test_data(grid, 3)[south_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(u, 2), South(); - operation=:subset, - index=:first) == - reverse(create_v_test_data(grid, 6)[east_indices_subset_skip_first_index...], dims=2)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(u, 2), South(); - operation=:endpoint, - index=:first) == - create_v_test_data(grid, 1)[east_indices_first...] - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(u, 3), West()) == reverse(create_v_test_data(grid, 1)[north_indices...], dims=1)' - @test get_halo_data(getregion(u, 3), East()) == create_u_test_data(grid, 4)[west_indices...] - @test get_halo_data(getregion(u, 3), South()) == create_u_test_data(grid, 2)[north_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(u, 3), North(); - operation=:subset, - index=:first) == - reverse(create_v_test_data(grid, 5)[west_indices_subset_skip_first_index...], dims=2)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(u, 3), North(); - operation=:endpoint, - index=:first) == - reverse(create_u_test_data(grid, 1)[north_indices_first...]) - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(u, 4), West()) == create_u_test_data(grid, 3)[east_indices...] - @test get_halo_data(getregion(u, 4), East()) == reverse(create_v_test_data(grid, 6)[south_indices...], dims=1)' - @test get_halo_data(getregion(u, 4), North()) == create_u_test_data(grid, 5)[south_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(u, 4), South(); - operation=:subset, - index=:first) == - reverse(create_v_test_data(grid, 2)[east_indices_subset_skip_first_index...], dims=2)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(u, 4), South(); - operation=:endpoint, - index=:first) == - create_v_test_data(grid, 3)[east_indices_first...] - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(u, 5), West()) == reverse(create_v_test_data(grid, 3)[north_indices...], dims=1)' - @test get_halo_data(getregion(u, 5), East()) == create_u_test_data(grid, 6)[west_indices...] - @test get_halo_data(getregion(u, 5), South()) == create_u_test_data(grid, 4)[north_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(u, 5), North(); - operation=:subset, - index=:first) == - reverse(create_v_test_data(grid, 1)[west_indices_subset_skip_first_index...], dims=2)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(u, 5), North(); - operation=:endpoint, - index=:first) == - reverse(create_u_test_data(grid, 3)[north_indices_first...]) - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(u, 6), West()) == create_u_test_data(grid, 5)[east_indices...] - @test get_halo_data(getregion(u, 6), East()) == reverse(create_v_test_data(grid, 2)[south_indices...], dims=1)' - @test get_halo_data(getregion(u, 6), North()) == create_u_test_data(grid, 1)[south_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(u, 6), South(); - operation=:subset, - index=:first) == - reverse(create_v_test_data(grid, 4)[east_indices_subset_skip_first_index...], dims=2)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(u, 6), South(); - operation=:endpoint, - index=:first) == - create_v_test_data(grid, 5)[east_indices_first...] - end # CUDA.@allowscalar - - # Confirm that the meridional velocity halos were filled according to connectivity described at - # ConformalCubedSphereGrid docstring. - @allowscalar begin - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(v, 1), East()) == create_v_test_data(grid, 2)[west_indices...] - @test get_halo_data(getregion(v, 1), South()) == create_v_test_data(grid, 6)[north_indices...] - @test get_halo_data(getregion(v, 1), North()) == reverse(create_u_test_data(grid, 3)[west_indices...], dims=2)' - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(v, 1), West(); - operation=:subset, - index=:first) == - reverse(create_u_test_data(grid, 5)[north_indices_subset_skip_first_index...], dims=1)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(v, 1), West(); - operation=:endpoint, - index=:first) == - create_u_test_data(grid, 6)[north_indices_first...] - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(v, 2), West()) == create_v_test_data(grid, 1)[east_indices...] - @test get_halo_data(getregion(v, 2), South()) == reverse(create_u_test_data(grid, 6)[east_indices...], dims=2)' - @test get_halo_data(getregion(v, 2), North()) == create_v_test_data(grid, 3)[south_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(v, 2), East(); - operation=:subset, - index=:first) == - reverse(create_u_test_data(grid, 4)[south_indices_subset_skip_first_index...], dims=1)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(v, 2), East(); - operation=:endpoint, - index=:first) == - reverse(create_v_test_data(grid, 6)[east_indices_first...]) - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(v, 3), East()) == create_v_test_data(grid, 4)[west_indices...] - @test get_halo_data(getregion(v, 3), South()) == create_v_test_data(grid, 2)[north_indices...] - @test get_halo_data(getregion(v, 3), North()) == reverse(create_u_test_data(grid, 5)[west_indices...], dims=2)' - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(v, 3), West(); - operation=:subset, - index=:first) == - reverse(create_u_test_data(grid, 1)[north_indices_subset_skip_first_index...], dims=1)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(v, 3), West(); - operation=:endpoint, - index=:first) == - create_u_test_data(grid, 2)[north_indices_first...] - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(v, 4), West()) == create_v_test_data(grid, 3)[east_indices...] - @test get_halo_data(getregion(v, 4), South()) == reverse(create_u_test_data(grid, 2)[east_indices...], dims=2)' - @test get_halo_data(getregion(v, 4), North()) == create_v_test_data(grid, 5)[south_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(v, 4), East(); - operation=:subset, - index=:first) == - reverse(create_u_test_data(grid, 6)[south_indices_subset_skip_first_index...], dims=1)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(v, 4), East(); - operation=:endpoint, - index=:first) == - reverse(create_v_test_data(grid, 2)[east_indices_first...]) - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(v, 5), East()) == create_v_test_data(grid, 6)[west_indices...] - @test get_halo_data(getregion(v, 5), South()) == create_v_test_data(grid, 4)[north_indices...] - @test get_halo_data(getregion(v, 5), North()) == reverse(create_u_test_data(grid, 1)[west_indices...], dims=2)' - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(v, 5), West(); - operation=:subset, - index=:first) == - reverse(create_u_test_data(grid, 3)[north_indices_subset_skip_first_index...], dims=1)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(v, 5), West(); - operation=:endpoint, - index=:first) == - create_u_test_data(grid, 4)[north_indices_first...] - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(v, 6), West()) == create_v_test_data(grid, 5)[east_indices...] - @test get_halo_data(getregion(v, 6), South()) == reverse(create_u_test_data(grid, 4)[east_indices...], dims=2)' - @test get_halo_data(getregion(v, 6), North()) == create_v_test_data(grid, 1)[south_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(v, 6), East(); - operation=:subset, - index=:first) == - reverse(create_u_test_data(grid, 2)[south_indices_subset_skip_first_index...], dims=1)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(v, 6), East(); - operation=:endpoint, - index=:first) == - reverse(create_v_test_data(grid, 4)[east_indices_first...]) - end # CUDA.@allowscalar + for FT in float_types + for arch in archs + for non_uniform_conformal_mapping in (false, true) + A = typeof(arch) + @info " Testing fill halos for horizontal velocities [$FT, $A]..." + + Nx, Ny, Nz = 9, 9, 1 + underlying_grid = ConformalCubedSphereGrid(arch, FT; + panel_size = (Nx, Ny, Nz), z = (0, 1), radius = 1, + horizontal_direction_halo = 3, non_uniform_conformal_mapping) + + @inline bottom(x, y) = ifelse(abs(y) < 30, - 2, 0) + immersed_grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom); active_cells_map = true) + + grids = (underlying_grid, immersed_grid) + + for grid in grids + u = XFaceField(grid) + v = YFaceField(grid) + + region = Iterate(1:6) + @apply_regionally u_data = create_u_test_data(grid, region) + @apply_regionally v_data = create_v_test_data(grid, region) + set!(u, u_data) + set!(v, v_data) + + fill_halo_regions!((u, v); signed = true) + + Hx, Hy, Hz = halo_size(u.grid) + + south_indices = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=nothing, index=:all) + east_indices = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=nothing, index=:all) + north_indices = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=nothing, index=:all) + west_indices = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=nothing, index=:all) + + south_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:endpoint, index=:first) + south_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:endpoint, index=:last) + east_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:endpoint, index=:first) + east_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:endpoint, index=:last) + north_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:endpoint, index=:first) + north_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:endpoint, index=:last) + west_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:endpoint, index=:first) + west_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:endpoint, index=:last) + + south_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:subset, index=:first) + south_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:subset, index=:last) + east_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:subset, index=:first) + east_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:subset, index=:last) + north_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:subset, index=:first) + north_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:subset, index=:last) + west_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:subset, index=:first) + west_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:subset, index=:last) + + # Confirm that the zonal velocity halos were filled according to connectivity described at ConformalCubedSphereGrid docstring. + @allowscalar begin + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(u, 1), West()) == reverse(create_v_test_data(grid, 5)[north_indices...], dims=1)' + @test get_halo_data(getregion(u, 1), East()) == create_u_test_data(grid, 2)[west_indices...] + @test get_halo_data(getregion(u, 1), South()) == create_u_test_data(grid, 6)[north_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(u, 1), North(); + operation=:subset, + index=:first) == - reverse(create_v_test_data(grid, 3)[west_indices_subset_skip_first_index...], dims=2)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(u, 1), North(); + operation=:endpoint, + index=:first) == - reverse(create_u_test_data(grid, 5)[north_indices_first...]) + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(u, 2), West()) == create_u_test_data(grid, 1)[east_indices...] + @test get_halo_data(getregion(u, 2), East()) == reverse(create_v_test_data(grid, 4)[south_indices...], dims=1)' + @test get_halo_data(getregion(u, 2), North()) == create_u_test_data(grid, 3)[south_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(u, 2), South(); + operation=:subset, + index=:first) == - reverse(create_v_test_data(grid, 6)[east_indices_subset_skip_first_index...], dims=2)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(u, 2), South(); + operation=:endpoint, + index=:first) == - create_v_test_data(grid, 1)[east_indices_first...] + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(u, 3), West()) == reverse(create_v_test_data(grid, 1)[north_indices...], dims=1)' + @test get_halo_data(getregion(u, 3), East()) == create_u_test_data(grid, 4)[west_indices...] + @test get_halo_data(getregion(u, 3), South()) == create_u_test_data(grid, 2)[north_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(u, 3), North(); + operation=:subset, + index=:first) == - reverse(create_v_test_data(grid, 5)[west_indices_subset_skip_first_index...], dims=2)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(u, 3), North(); + operation=:endpoint, + index=:first) == - reverse(create_u_test_data(grid, 1)[north_indices_first...]) + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(u, 4), West()) == create_u_test_data(grid, 3)[east_indices...] + @test get_halo_data(getregion(u, 4), East()) == reverse(create_v_test_data(grid, 6)[south_indices...], dims=1)' + @test get_halo_data(getregion(u, 4), North()) == create_u_test_data(grid, 5)[south_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(u, 4), South(); + operation=:subset, + index=:first) == - reverse(create_v_test_data(grid, 2)[east_indices_subset_skip_first_index...], dims=2)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(u, 4), South(); + operation=:endpoint, + index=:first) == - create_v_test_data(grid, 3)[east_indices_first...] + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(u, 5), West()) == reverse(create_v_test_data(grid, 3)[north_indices...], dims=1)' + @test get_halo_data(getregion(u, 5), East()) == create_u_test_data(grid, 6)[west_indices...] + @test get_halo_data(getregion(u, 5), South()) == create_u_test_data(grid, 4)[north_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(u, 5), North(); + operation=:subset, + index=:first) == - reverse(create_v_test_data(grid, 1)[west_indices_subset_skip_first_index...], dims=2)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(u, 5), North(); + operation=:endpoint, + index=:first) == - reverse(create_u_test_data(grid, 3)[north_indices_first...]) + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(u, 6), West()) == create_u_test_data(grid, 5)[east_indices...] + @test get_halo_data(getregion(u, 6), East()) == reverse(create_v_test_data(grid, 2)[south_indices...], dims=1)' + @test get_halo_data(getregion(u, 6), North()) == create_u_test_data(grid, 1)[south_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(u, 6), South(); + operation=:subset, + index=:first) == - reverse(create_v_test_data(grid, 4)[east_indices_subset_skip_first_index...], dims=2)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(u, 6), South(); + operation=:endpoint, + index=:first) == - create_v_test_data(grid, 5)[east_indices_first...] + end # CUDA.@allowscalar + + # Confirm that the meridional velocity halos were filled according to connectivity described at + # ConformalCubedSphereGrid docstring. + @allowscalar begin + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(v, 1), East()) == create_v_test_data(grid, 2)[west_indices...] + @test get_halo_data(getregion(v, 1), South()) == create_v_test_data(grid, 6)[north_indices...] + @test get_halo_data(getregion(v, 1), North()) == reverse(create_u_test_data(grid, 3)[west_indices...], dims=2)' + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(v, 1), West(); + operation=:subset, + index=:first) == - reverse(create_u_test_data(grid, 5)[north_indices_subset_skip_first_index...], dims=1)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(v, 1), West(); + operation=:endpoint, + index=:first) == - create_u_test_data(grid, 6)[north_indices_first...] + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(v, 2), West()) == create_v_test_data(grid, 1)[east_indices...] + @test get_halo_data(getregion(v, 2), South()) == reverse(create_u_test_data(grid, 6)[east_indices...], dims=2)' + @test get_halo_data(getregion(v, 2), North()) == create_v_test_data(grid, 3)[south_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(v, 2), East(); + operation=:subset, + index=:first) == - reverse(create_u_test_data(grid, 4)[south_indices_subset_skip_first_index...], dims=1)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(v, 2), East(); + operation=:endpoint, + index=:first) == - reverse(create_v_test_data(grid, 6)[east_indices_first...]) + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(v, 3), East()) == create_v_test_data(grid, 4)[west_indices...] + @test get_halo_data(getregion(v, 3), South()) == create_v_test_data(grid, 2)[north_indices...] + @test get_halo_data(getregion(v, 3), North()) == reverse(create_u_test_data(grid, 5)[west_indices...], dims=2)' + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(v, 3), West(); + operation=:subset, + index=:first) == - reverse(create_u_test_data(grid, 1)[north_indices_subset_skip_first_index...], dims=1)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(v, 3), West(); + operation=:endpoint, + index=:first) == - create_u_test_data(grid, 2)[north_indices_first...] + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(v, 4), West()) == create_v_test_data(grid, 3)[east_indices...] + @test get_halo_data(getregion(v, 4), South()) == reverse(create_u_test_data(grid, 2)[east_indices...], dims=2)' + @test get_halo_data(getregion(v, 4), North()) == create_v_test_data(grid, 5)[south_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(v, 4), East(); + operation=:subset, + index=:first) == - reverse(create_u_test_data(grid, 6)[south_indices_subset_skip_first_index...], dims=1)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(v, 4), East(); + operation=:endpoint, + index=:first) == - reverse(create_v_test_data(grid, 2)[east_indices_first...]) + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(v, 5), East()) == create_v_test_data(grid, 6)[west_indices...] + @test get_halo_data(getregion(v, 5), South()) == create_v_test_data(grid, 4)[north_indices...] + @test get_halo_data(getregion(v, 5), North()) == reverse(create_u_test_data(grid, 1)[west_indices...], dims=2)' + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(v, 5), West(); + operation=:subset, + index=:first) == - reverse(create_u_test_data(grid, 3)[north_indices_subset_skip_first_index...], dims=1)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(v, 5), West(); + operation=:endpoint, + index=:first) == - create_u_test_data(grid, 4)[north_indices_first...] + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(v, 6), West()) == create_v_test_data(grid, 5)[east_indices...] + @test get_halo_data(getregion(v, 6), South()) == reverse(create_u_test_data(grid, 4)[east_indices...], dims=2)' + @test get_halo_data(getregion(v, 6), North()) == create_v_test_data(grid, 1)[south_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(v, 6), East(); + operation=:subset, + index=:first) == - reverse(create_u_test_data(grid, 2)[south_indices_subset_skip_first_index...], dims=1)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(v, 6), East(); + operation=:endpoint, + index=:first) == - reverse(create_v_test_data(grid, 4)[east_indices_first...]) + end # CUDA.@allowscalar + end + end end end end @testset "Testing conformal cubed sphere fill halos for Face-Face-Any field" begin - for FT in float_types, arch in archs, non_uniform_conformal_mapping in (false, true) - cm = non_uniform_conformal_mapping ? "non-uniform conformal mapping" : "uniform conformal mapping" - @info " Testing fill halos for streamfunction [$FT, $(typeof(arch)), $cm]..." - - Nx, Ny, Nz = 9, 9, 1 - - underlying_grid = ConformalCubedSphereGrid(arch, FT; - panel_size = (Nx, Ny, Nz), z = (0, 1), radius = 1, - horizontal_direction_halo = 3, non_uniform_conformal_mapping) - @inline bottom(x, y) = ifelse(abs(y) < 30, - 2, 0) - immersed_grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom); active_cells_map = true) - - grids = (underlying_grid, immersed_grid) - - for grid in grids - ψ = Field{Face, Face, Center}(grid) - - region = Iterate(1:6) - @apply_regionally data = create_ψ_test_data(grid, region) - set!(ψ, data) - - fill_halo_regions!(ψ) - - Hx, Hy, Hz = halo_size(ψ.grid) - - south_indices = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=nothing, index=:all) - east_indices = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=nothing, index=:all) - north_indices = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=nothing, index=:all) - west_indices = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=nothing, index=:all) - - south_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:endpoint, index=:first) - south_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:endpoint, index=:last) - east_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:endpoint, index=:first) - east_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:endpoint, index=:last) - north_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:endpoint, index=:first) - north_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:endpoint, index=:last) - west_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:endpoint, index=:first) - west_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:endpoint, index=:last) - - south_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:subset, index=:first) - south_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:subset, index=:last) - east_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:subset, index=:first) - east_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:subset, index=:last) - north_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:subset, index=:first) - north_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:subset, index=:last) - west_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:subset, index=:first) - west_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:subset, index=:last) - - # Confirm that the tracer halos were filled according to connectivity described at ConformalCubedSphereGrid docstring. - @allowscalar begin - # Panel 1 - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(ψ, 1), East()) == create_ψ_test_data(grid, 2)[west_indices...] - @test get_halo_data(getregion(ψ, 1), South()) == create_ψ_test_data(grid, 6)[north_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 1), North(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 3)[west_indices_subset_skip_first_index...], dims=2)' - # Currently we do not have any test for the point of intersection of the northwest (halo) corners of panels 1, 3, and 5. - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 1), West(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 5)[north_indices_subset_skip_first_index...], dims=1)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(ψ, 1), West(); - operation=:endpoint, - index=:first) == create_ψ_test_data(grid, 6)[north_indices_first...] - - # Panel 2 - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(ψ, 2), West()) == create_ψ_test_data(grid, 1)[east_indices...] - @test get_halo_data(getregion(ψ, 2), North()) == create_ψ_test_data(grid, 3)[south_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 2), East(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 4)[south_indices_subset_skip_first_index...], dims=1)' - # Currently we do not have any test for the point of intersection of the southeast (halo) corners of panels 2, 4, and 6. - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 2), South(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 6)[east_indices_subset_skip_first_index...], dims=2)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(ψ, 2), South(); - operation=:endpoint, - index=:first) == create_ψ_test_data(grid, 1)[east_indices_first...] - - # Panel 3 - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(ψ, 3), East()) == create_ψ_test_data(grid, 4)[west_indices...] - @test get_halo_data(getregion(ψ, 3), South()) == create_ψ_test_data(grid, 2)[north_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 3), West(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 1)[north_indices_subset_skip_first_index...], dims=1)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(ψ, 3), West(); - operation=:endpoint, - index=:first) == create_ψ_test_data(grid, 2)[north_indices_first...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 3), North(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 5)[west_indices_subset_skip_first_index...], dims=2)' - # Currently we do not have any test for the point of intersection of the northwest (halo) corners of panels 1, 3, and 5. - - # Panel 4 - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(ψ, 4), West()) == create_ψ_test_data(grid, 3)[east_indices...] - @test get_halo_data(getregion(ψ, 4), North()) == create_ψ_test_data(grid, 5)[south_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 4), East(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 6)[south_indices_subset_skip_first_index...], dims=1)' - # Currently we do not have any test for the point of intersection of the southeast (halo) corners of panels 2, 4, and 6. - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 4), South(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 2)[east_indices_subset_skip_first_index...], dims=2)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(ψ, 4), South(); - operation=:endpoint, - index=:first) == create_ψ_test_data(grid, 3)[east_indices_first...] - - # Panel 5 - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(ψ, 5), East()) == create_ψ_test_data(grid, 6)[west_indices...] - @test get_halo_data(getregion(ψ, 5), South()) == create_ψ_test_data(grid, 4)[north_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 5), West(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 3)[north_indices_subset_skip_first_index...], dims=1)' - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(ψ, 5), West(); - operation=:endpoint, - index=:first) == create_ψ_test_data(grid, 4)[north_indices_first...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 5), North(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 1)[west_indices_subset_skip_first_index...], dims=2)' - # Currently we do not have any test for the point of intersection of the northwest (halo) corners of panels 1, 3, and 5. - - # Panel 6 - - # Trivial halo checks with no off-set in index - @test get_halo_data(getregion(ψ, 6), West()) == create_ψ_test_data(grid, 5)[east_indices...] - @test get_halo_data(getregion(ψ, 6), North()) == create_ψ_test_data(grid, 1)[south_indices...] - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 6), East(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 2)[south_indices_subset_skip_first_index...], dims=1)' - # Currently we do not have any test for the point of intersection of the southeast (halo) corners of panels 2, 4, and 6. - - # Non-trivial halo checks with off-set in index - @test get_halo_data(getregion(ψ, 6), South(); - operation=:subset, - index=:first) == reverse(create_ψ_test_data(grid, 4)[east_indices_subset_skip_first_index...], dims=2)' - - # The index appearing on the LHS above is the index to be skipped. - @test get_halo_data(getregion(ψ, 6), South(); - operation=:endpoint, - index=:first) == create_ψ_test_data(grid, 5)[east_indices_first...] - end # CUDA.@allowscalar + for FT in float_types + for arch in archs + for non_uniform_conformal_mapping in (false, true) + A = typeof(arch) + @info " Testing fill halos for streamfunction [$FT, $A]..." + + Nx, Ny, Nz = 9, 9, 1 + underlying_grid = ConformalCubedSphereGrid(arch, FT; + panel_size = (Nx, Ny, Nz), z = (0, 1), radius = 1, + horizontal_direction_halo = 3, non_uniform_conformal_mapping) + + @inline bottom(x, y) = ifelse(abs(y) < 30, - 2, 0) + immersed_grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom); active_cells_map = true) + + grids = (underlying_grid, immersed_grid) + + for grid in grids + ψ = Field{Face, Face, Center}(grid) + + region = Iterate(1:6) + @apply_regionally data = create_ψ_test_data(grid, region) + set!(ψ, data) + + fill_halo_regions!(ψ) + + Hx, Hy, Hz = halo_size(ψ.grid) + + south_indices = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=nothing, index=:all) + east_indices = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=nothing, index=:all) + north_indices = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=nothing, index=:all) + west_indices = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=nothing, index=:all) + + south_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:endpoint, index=:first) + south_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:endpoint, index=:last) + east_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:endpoint, index=:first) + east_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:endpoint, index=:last) + north_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:endpoint, index=:first) + north_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:endpoint, index=:last) + west_indices_first = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:endpoint, index=:first) + west_indices_last = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:endpoint, index=:last) + + south_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:subset, index=:first) + south_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, South(); operation=:subset, index=:last) + east_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:subset, index=:first) + east_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, East(); operation=:subset, index=:last) + north_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:subset, index=:first) + north_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, North(); operation=:subset, index=:last) + west_indices_subset_skip_first_index = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:subset, index=:first) + west_indices_subset_skip_last_index = get_boundary_indices(Nx, Ny, Hx, Hy, West(); operation=:subset, index=:last) + + # Confirm that the tracer halos were filled according to connectivity described at ConformalCubedSphereGrid docstring. + @allowscalar begin + # Panel 1 + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(ψ, 1), East()) == create_ψ_test_data(grid, 2)[west_indices...] + @test get_halo_data(getregion(ψ, 1), South()) == create_ψ_test_data(grid, 6)[north_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 1), North(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 3)[west_indices_subset_skip_first_index...], dims=2)' + # Currently we do not have any test for the point of intersection of the northwest (halo) corners of panels 1, 3, and 5. + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 1), West(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 5)[north_indices_subset_skip_first_index...], dims=1)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(ψ, 1), West(); + operation=:endpoint, + index=:first) == create_ψ_test_data(grid, 6)[north_indices_first...] + + # Panel 2 + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(ψ, 2), West()) == create_ψ_test_data(grid, 1)[east_indices...] + @test get_halo_data(getregion(ψ, 2), North()) == create_ψ_test_data(grid, 3)[south_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 2), East(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 4)[south_indices_subset_skip_first_index...], dims=1)' + # Currently we do not have any test for the point of intersection of the southeast (halo) corners of panels 2, 4, and 6. + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 2), South(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 6)[east_indices_subset_skip_first_index...], dims=2)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(ψ, 2), South(); + operation=:endpoint, + index=:first) == create_ψ_test_data(grid, 1)[east_indices_first...] + + # Panel 3 + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(ψ, 3), East()) == create_ψ_test_data(grid, 4)[west_indices...] + @test get_halo_data(getregion(ψ, 3), South()) == create_ψ_test_data(grid, 2)[north_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 3), West(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 1)[north_indices_subset_skip_first_index...], dims=1)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(ψ, 3), West(); + operation=:endpoint, + index=:first) == create_ψ_test_data(grid, 2)[north_indices_first...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 3), North(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 5)[west_indices_subset_skip_first_index...], dims=2)' + # Currently we do not have any test for the point of intersection of the northwest (halo) corners of panels 1, 3, and 5. + + # Panel 4 + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(ψ, 4), West()) == create_ψ_test_data(grid, 3)[east_indices...] + @test get_halo_data(getregion(ψ, 4), North()) == create_ψ_test_data(grid, 5)[south_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 4), East(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 6)[south_indices_subset_skip_first_index...], dims=1)' + # Currently we do not have any test for the point of intersection of the southeast (halo) corners of panels 2, 4, and 6. + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 4), South(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 2)[east_indices_subset_skip_first_index...], dims=2)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(ψ, 4), South(); + operation=:endpoint, + index=:first) == create_ψ_test_data(grid, 3)[east_indices_first...] + + # Panel 5 + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(ψ, 5), East()) == create_ψ_test_data(grid, 6)[west_indices...] + @test get_halo_data(getregion(ψ, 5), South()) == create_ψ_test_data(grid, 4)[north_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 5), West(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 3)[north_indices_subset_skip_first_index...], dims=1)' + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(ψ, 5), West(); + operation=:endpoint, + index=:first) == create_ψ_test_data(grid, 4)[north_indices_first...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 5), North(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 1)[west_indices_subset_skip_first_index...], dims=2)' + # Currently we do not have any test for the point of intersection of the northwest (halo) corners of panels 1, 3, and 5. + + # Panel 6 + + # Trivial halo checks with no off-set in index + @test get_halo_data(getregion(ψ, 6), West()) == create_ψ_test_data(grid, 5)[east_indices...] + @test get_halo_data(getregion(ψ, 6), North()) == create_ψ_test_data(grid, 1)[south_indices...] + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 6), East(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 2)[south_indices_subset_skip_first_index...], dims=1)' + # Currently we do not have any test for the point of intersection of the southeast (halo) corners of panels 2, 4, and 6. + + # Non-trivial halo checks with off-set in index + @test get_halo_data(getregion(ψ, 6), South(); + operation=:subset, + index=:first) == reverse(create_ψ_test_data(grid, 4)[east_indices_subset_skip_first_index...], dims=2)' + + # The index appearing on the LHS above is the index to be skipped. + @test get_halo_data(getregion(ψ, 6), South(); + operation=:endpoint, + index=:first) == create_ψ_test_data(grid, 5)[east_indices_first...] + end # CUDA.@allowscalar + end + end end end end @testset "Testing simulation on conformal and immersed conformal cubed sphere grids" begin - for non_uniform_conformal_mapping in (false, true) - cm = non_uniform_conformal_mapping ? "non-uniform conformal mapping" : "uniform conformal mapping" - cm_suffix = non_uniform_conformal_mapping ? "NUCM" : "UCM" - for FT in float_types, arch in archs - Nx, Ny, Nz = 18, 18, 9 - - underlying_grid = ConformalCubedSphereGrid(arch, FT; - panel_size = (Nx, Ny, Nz), z = (0, 1), radius = 1, - horizontal_direction_halo = 6, non_uniform_conformal_mapping) - @inline bottom(x, y) = ifelse(abs(y) < 30, - 2, 0) - immersed_grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom); - active_cells_map = true) - - grids = (underlying_grid, immersed_grid) - - for grid in grids - if grid == underlying_grid - @info " Testing simulation on conformal cubed sphere grid [$FT, $(typeof(arch)), $cm]..." - grid_suffix = "UG" - else - @info " Testing simulation on immersed boundary conformal cubed sphere grid [$FT, $(typeof(arch)), $cm]..." - grid_suffix = "IG" - end - - model = HydrostaticFreeSurfaceModel(; grid, - momentum_advection = WENOVectorInvariant(FT; order=5), - tracer_advection = WENO(FT; order=5), - free_surface = SplitExplicitFreeSurface(grid; substeps=12), - coriolis = HydrostaticSphericalCoriolis(FT), - tracers = :b, - buoyancy = BuoyancyTracer()) - - simulation = Simulation(model, Δt=1minute, stop_time=10minutes) - - save_fields_interval = 2minute - checkpointer_interval = 4minutes + for FT in float_types + for arch in archs + for non_uniform_conformal_mapping in (false, true) + + Nx, Ny, Nz = 9, 9, 9 + underlying_grid = ConformalCubedSphereGrid(arch, FT; + panel_size = (Nx, Ny, Nz), z = (0, 1), radius = 1, + horizontal_direction_halo = 6, non_uniform_conformal_mapping) + + @inline bottom(x, y) = ifelse(abs(y) < 30, - 2, 0) + immersed_grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom); + active_cells_map = true) + + grids = (underlying_grid, immersed_grid) + + for grid in grids + if grid == underlying_grid + grid_msg = "" + grid_suffix = "UG" + else + grid_msg = "immersed boundary " + grid_suffix = "IG" + end - filename_checkpointer = - "cubed_sphere_checkpointer_$(FT)_$(typeof(arch))_" * cm_suffix * "_" * grid_suffix - filename_output_writer = "cubed_sphere_output_$(FT)_$(typeof(arch))_" * cm_suffix * "_" * grid_suffix + A = typeof(arch) + cm = non_uniform_conformal_mapping ? "non-uniform conformal mapping" : "uniform conformal mapping" + cm_suffix = non_uniform_conformal_mapping ? "NUCM" : "UCM" + + @info " Testing simulation on $grid_msg conformal cubed sphere grid [$FT, $A, $cm]..." + + model = HydrostaticFreeSurfaceModel(; grid, + momentum_advection = WENOVectorInvariant(FT; order=5), + tracer_advection = WENO(FT; order=5), + free_surface = SplitExplicitFreeSurface(grid; substeps=12), + coriolis = HydrostaticSphericalCoriolis(FT), + tracers = :b, + buoyancy = BuoyancyTracer()) + + ϵ(x, y, z) = rand() + N² = (1 / 3600)^2 + bᵢ(x, y, z) = N² * z + N² * 1e-6 * ϵ(x, y, z) + uᵢ(x, y, z) = 1e-9 * ϵ(x, y, z) + set!(model, b=bᵢ, u=uᵢ, v=uᵢ) + + Δt = 1e-3 + stop_time = 10Δt + simulation = Simulation(model; Δt, stop_time) + + filename_checkpointer = + "cubed_sphere_checkpointer_$(FT)_$(A)_" * cm_suffix * "_" * grid_suffix + filename_output_writer = "cubed_sphere_output_$(FT)_$(A)_" * cm_suffix * "_" * grid_suffix + + # If previous run produced these files, remove them now to ensure a clean test. + for f in readdir(".") + if f == filename_output_writer * ".jld2" || occursin(r"^" * filename_checkpointer * r"_.*\.jld2$", f) + rm(f; force=true) + end + end - # If previous run produced these files, remove them now to ensure a clean test. - for f in readdir(".") - if f == filename_output_writer * ".jld2" || occursin(r"^" * filename_checkpointer * r"_.*\.jld2$", f) - rm(f; force=true) + simulation.output_writers[:checkpointer] = Checkpointer(model, + schedule = IterationInterval(4), + prefix = filename_checkpointer, + overwrite_existing = true) + + outputs = fields(model) + simulation.output_writers[:fields] = JLD2Writer(model, outputs; + schedule = IterationInterval(2), + filename = filename_output_writer, + array_type = Array{FT}, + verbose = false, + overwrite_existing = true) + + u2, v2, w2, b2 = nothing, nothing, nothing, nothing + function save_at_2(sim) + arch = sim.model.grid.architecture + if iteration(sim) == 2 + if arch isa GPU + u2 = on_architecture(CPU(), sim.model.velocities.u) + v2 = on_architecture(CPU(), sim.model.velocities.v) + w2 = on_architecture(CPU(), sim.model.velocities.w) + b2 = on_architecture(CPU(), sim.model.tracers.b) + else + u2 = deepcopy(sim.model.velocities.u) + v2 = deepcopy(sim.model.velocities.v) + w2 = deepcopy(sim.model.velocities.w) + b2 = deepcopy(sim.model.tracers.b) + end + end end - end - simulation.output_writers[:checkpointer] = Checkpointer(model, - schedule = TimeInterval(checkpointer_interval), - prefix = filename_checkpointer, - overwrite_existing = true) + add_callback!(simulation, save_at_2) + run!(simulation) + + @test iteration(simulation) == 10 + @test isapprox(time(simulation), stop_time, rtol=1e-3) + + uN, vN, wN, bN = nothing, nothing, nothing, nothing + if arch isa GPU + uN = on_architecture(CPU(), model.velocities.u) + vN = on_architecture(CPU(), model.velocities.v) + wN = on_architecture(CPU(), model.velocities.w) + bN = on_architecture(CPU(), model.tracers.b) + else + uN = model.velocities.u + vN = model.velocities.v + wN = model.velocities.w + bN = model.tracers.b + end - outputs = fields(model) - simulation.output_writers[:fields] = JLD2Writer(model, outputs; - schedule = TimeInterval(save_fields_interval), - filename = filename_output_writer, - verbose = false, - overwrite_existing = true) + if grid == underlying_grid + @info " Restarting simulation from pickup file on conformal cubed sphere grid [$FT, $(typeof(arch)), $cm]..." + else + @info " Restarting simulation from pickup file on immersed boundary conformal cubed sphere grid [$FT, $(typeof(arch)), $cm]..." + end - run!(simulation) + @test ut.times[end] == time(simulation) + @test ut[end] == uN + @test vt[end] == vN + @test wt[end] == wN + @test bt[end] == bN - @test iteration(simulation) == 10 - @test time(simulation) == 10minutes + if grid == underlying_grid + @info " Restarting simulation from pickup file on conformal cubed sphere grid [$FT, $A]..." + else + @info " Restarting simulation from pickup file on immersed boundary conformal cubed sphere grid [$FT, $A]..." + end - u_timeseries = FieldTimeSeries(filename_output_writer * ".jld2", "u"; architecture = CPU()) + new_stop_time = 2stop_time + simulation = Simulation(model; Δt, stop_time=new_stop_time) + + simulation.output_writers[:checkpointer] = Checkpointer(model, + schedule = IterationInterval(4), + prefix = filename_checkpointer, + overwrite_existing = true) + + simulation.output_writers[:fields] = JLD2Writer(model, outputs; + schedule = IterationInterval(2), + filename = filename_output_writer, + verbose = false, + overwrite_existing = false) + + run!(simulation, pickup=true) + + @test time(simulation) ≈ new_stop_time + + ut = FieldTimeSeries(filename_output_writer * ".jld2", "u") + vt = FieldTimeSeries(filename_output_writer * ".jld2", "v") + wt = FieldTimeSeries(filename_output_writer * ".jld2", "w") + bt = FieldTimeSeries(filename_output_writer * ".jld2", "b") + + if arch isa GPU + uN = on_architecture(CPU(), model.velocities.u) + vN = on_architecture(CPU(), model.velocities.v) + wN = on_architecture(CPU(), model.velocities.w) + bN = on_architecture(CPU(), model.tracers.b) + else + uN = model.velocities.u + vN = model.velocities.v + wN = model.velocities.w + bN = model.tracers.b + end - if grid == underlying_grid - @info " Restarting simulation from pickup file on conformal cubed sphere grid [$FT, $(typeof(arch)), $cm]..." - else - @info " Restarting simulation from pickup file on immersed boundary conformal cubed sphere grid [$FT, $(typeof(arch)), $cm]..." + @test ut[end] == uN + @test vt[end] == vN + @test wt[end] == wN + @test bt[end] == bN end - - simulation = Simulation(model, Δt=1minute, stop_time=20minutes) - - simulation.output_writers[:checkpointer] = Checkpointer(model, - schedule = TimeInterval(checkpointer_interval), - prefix = filename_checkpointer, - overwrite_existing = true) - - simulation.output_writers[:fields] = JLD2Writer(model, outputs; - schedule = TimeInterval(save_fields_interval), - filename = filename_output_writer, - verbose = false, - overwrite_existing = true) - - run!(simulation, pickup = true) - - @test iteration(simulation) == 20 - @test time(simulation) == 20minutes - - u_timeseries = FieldTimeSeries(filename_output_writer * ".jld2", "u"; architecture = CPU()) end end end