Fix bug for lat-lon grid with Periodic y-dimension

[navidcy]

Closes #3465

[navidcy]

there are other places that this can change -- just wanna see if all tests are OK first

[navidcy]

I believe this PR is ready to review

[navidcy]

Interestingly this breaks the distributed tests!
https://buildkite.com/clima/oceananigans-distributed/builds/1218#018da64e-30be-4ecb-9b10-32cdb3daa89c

[glwagner]

That's annoying. There must be some implicit assumptions that we can't see

[glwagner]

@simone-silvestri can you help?

[simone-silvestri]

this is the assumption

Oceananigans.jl/src/DistributedComputations/distributed_grids.jl

Lines 140 to 160 in c55878e

	# The Latitudinal direction is _special_:
	# precompute_metrics assumes that `length(φᵃᶠᵃ) = length(φᵃᶜᵃ) + 1`, which is always the case in a
	# serial grid because `LatitudeLongitudeGrid` should be always `Bounded`, but it is not true for a
	# partitioned `DistributedGrid` with Ry > 1 (one rank will hold a `RightConnected` topology)
	# An extra point is to precompute the Y-direction metrics in case of only one halo, hence
	# we disregard the topology when constructing the metrics and add a halo point!
	# Furthermore, the `LatitudeLongitudeGrid` requires an extra halo on it's latitudinal coordinate to allow calculating
	# the z-area on halo cells. (see: Az = R^2 * Δλ * (sin(φ[j]) - sin(φ[j-1]))
	Lφ, φᵃᶠᵃ, φᵃᶜᵃ, Δφᵃᶠᵃ, Δφᵃᶜᵃ = generate_coordinate(FT, Bounded(), nφ, Hφ + 1, φl, :latitude, arch.child_architecture)
	preliminary_grid = LatitudeLongitudeGrid{TX, TY, TZ}(arch,
	nλ, nφ, nz,
	Hλ, Hφ, Hz,
	Lλ, Lφ, Lz,
	Δλᶠᵃᵃ, Δλᶜᵃᵃ, λᶠᵃᵃ, λᶜᵃᵃ,
	Δφᵃᶠᵃ, Δφᵃᶜᵃ, φᵃᶠᵃ, φᵃᶜᵃ,
	Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ,
	(nothing for i=1:10)..., convert(FT, radius))
	return !precompute_metrics ? preliminary_grid : with_precomputed_metrics(preliminary_grid)
	end

[glwagner]

I don't understand that text on first read...

I think it's just wrong that the meridional direction is always Bounded. We need Periodic too. How can we implement this?

[simone-silvestri]

Hmmm, I am not sure how you solved this in the serial grid, but I guess in the same way. By removing that assumption and calculating the coordinate accordingly?

[glwagner]

Hmmm, I am not sure how you solved this in the serial grid, but I guess in the same way. By removing that assumption and calculating the coordinate accordingly?

Since we removed the assumption, do we also eliminate the need for this code which avoids the assumption?

[simone-silvestri]

If the metrics are correctly calculated on the boundary for FullyConnected, RightConnected, and LeftConnected topologies yes.

Note that in the latitude direction, the halo metrics for FullyConnected and Periodic are very different!
In the sense that FullyConnected is like a Bounded in the halos (Periodic in latitude is a bit wonky)

[glwagner]

If the metrics are correctly calculated on the boundary for FullyConnected, RightConnected, and LeftConnected topologies yes.

Note that in the latitude direction, the halo metrics for FullyConnected and Periodic are very different! In the sense that FullyConnected is like a Bounded in the halos (Periodic in latitude is a bit wonky)

Periodic in latitude is precisely defined...

[simone-silvestri]

sure it is possible to define it, still it is conceptually wrong, you can easily see it if you visualize the grid:

julia> grid = LatitudeLongitudeGrid(size = (1, 2, 1), topology = (Periodic, Periodic, Bounded), latitude = (50, 52), longitude = (10, 11), z = (0, 1))
1×2×1 LatitudeLongitudeGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── longitude: Periodic λ ∈ [10.0, 11.0) regularly spaced with Δλ=1.0
├── latitude:  Periodic φ ∈ [50.0, 52.0) regularly spaced with Δφ=1.0
└── z:         Bounded  z ∈ [0.0, 1.0]   regularly spaced with Δz=1.0

julia> grid.Δxᶜᶠᵃ
8-element OffsetArray(::Vector{Float64}, -2:5) with eltype Float64 with indices -2:5:
     0.0
 74403.92868970236
 72950.43560309989
 71474.721107126
 69977.2347187117
 68458.43258671739
 66918.77735298542
     0.0

Note, the delta x in the halos. They are inconsistent with the periodic assumption. Even if it were, it is impossible to reconcile Δxᶠᵃᵃ at Ny + 1 and Δxᶠᵃᵃ at 1 which should be exactly the same (but they end up being different)! This means that velocities and tracers that you use to "fill" the periodic halo will be exchanged between grid volumes of different size.
Another example, these volumes should be the same

julia> Vᶜᶜᶜ(1, 3, 1, grid)
7.526820531985937e9

julia> Vᶜᶜᶜ(1, 1, 1, grid)
7.864569567312662e9

[glwagner]

sure it is possible to define it, still it is conceptually wrong, you can easily see it if you visualize the grid:

julia> grid = LatitudeLongitudeGrid(size = (1, 2, 1), topology = (Periodic, Periodic, Bounded), latitude = (50, 52), longitude = (10, 11), z = (0, 1))
1×2×1 LatitudeLongitudeGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── longitude: Periodic λ ∈ [10.0, 11.0) regularly spaced with Δλ=1.0
├── latitude:  Periodic φ ∈ [50.0, 52.0) regularly spaced with Δφ=1.0
└── z:         Bounded  z ∈ [0.0, 1.0]   regularly spaced with Δz=1.0

julia> grid.Δxᶜᶠᵃ
8-element OffsetArray(::Vector{Float64}, -2:5) with eltype Float64 with indices -2:5:
     0.0
 74403.92868970236
 72950.43560309989
 71474.721107126
 69977.2347187117
 68458.43258671739
 66918.77735298542
     0.0

Note, the delta x in the halos. They are inconsistent with the periodic assumption. Even if it were, it is impossible to reconcile Δxᶠᵃᵃ at Ny + 1 and Δxᶠᵃᵃ at 1 which should be exactly the same (but they end up being different)! This means that velocities and tracers that you use to "fill" the periodic halo will be exchanged between grid volumes of different size. Another example, these volumes should be the same

julia> Vᶜᶜᶜ(1, 3, 1, grid)
7.526820531985937e9

julia> Vᶜᶜᶜ(1, 1, 1, grid)
7.864569567312662e9

It's valid for a single column and that is exactly why we need this feature. If you want, we can throw an error or warning if we have more than one point in y.

[simone-silvestri]

Shouldn't a single column use a Flat topologies? I think that would probably be more appropriate

[simone-silvestri]

In a periodic single column, again, the problem would be that Δxᶠᵃᵃ[1] is not equal to Δxᶠᵃᵃ[2] which it should.

A consequence is that, since there should be zero derivatives in the y direction, v[1] != v[2] because of the difference in areas.

[glwagner]

Shouldn't a single column use a Flat topologies? I think that would probably be more appropriate

Not when we want to interpolate data onto a particular location.

[glwagner]

In a periodic single column, again, the problem would be that Δxᶠᵃᵃ[1] is not equal to Δxᶠᵃᵃ[2] which it should.

A consequence is that, since there should be zero derivatives in the y direction, v[1] != v[2] because of the difference in areas.

The halo region metrics need to be enforced to be periodic. Just like we do for RectilinearGrid.

[simone-silvestri]

Nono, these are not halo metrics, they are part of the actual grid, periodic assumes that Δxᶠᵃᵃ[1] == Δxᶠᵃᵃ[2] (for a single column) I am not sure how you reconcile the fact that they are not the same. In practice, if the grid is not symmetric about the equator, I think that a single column cell (with associated periodic halos) would look like this:

So I guess the only way to be consistent would be to straighten out that cells

[glwagner]

Where in the code does this crop up? Note that advection is turned off for a single column model...

[glwagner]

This is what we get for Flat:

julia> flat_grid = LatitudeLongitudeGrid(size=grid.Nz, longitude=(215, 215.25), latitude=(50, 50.25), z=(-400, 0), topology=(Flat, Flat, Bounded))
1×1×100 LatitudeLongitudeGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo and without precomputed metrics
├── longitude: Flat λ
├── latitude:  Flat φ
└── z:         Bounded  z ∈ [-400.0, 0.0]    regularly spaced with Δz=4.0

julia> flat_grid.λᶜᵃᵃ
StepRangeLen(1.0, 0.0, 1)

[glwagner]

Changes in both grid construction and also interpolation would be needed to support Flat, Flat, Bounded with non-trivial location:

julia> flat_grid = LatitudeLongitudeGrid(size=grid.Nz, longitude=(215, 215.25), latitude=(50, 50.25), z=(-400, 0), topology=(Flat, Flat, Bounded))
1×1×100 LatitudeLongitudeGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo and without precomputed metrics
├── longitude: Flat λ
├── latitude:  Flat φ
└── z:         Bounded  z ∈ [-400.0, 0.0]    regularly spaced with Δz=4.0

julia> flat_grid.λᶜᵃᵃ
StepRangeLen(1.0, 0.0, 1)

julia> non_flat_grid = LatitudeLongitudeGrid(size=(8, 8, grid.Nz), longitude=(214, 216), latitude=(49, 51), z=(-400, 0), topology=(Bounded, Bounded, Bounded))
8×8×100 LatitudeLongitudeGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── longitude: Bounded  λ ∈ [214.0, 216.0] regularly spaced with Δλ=0.25
├── latitude:  Bounded  φ ∈ [49.0, 51.0]   regularly spaced with Δφ=0.25
└── z:         Bounded  z ∈ [-400.0, 0.0]  regularly spaced with Δz=4.0

julia> c_src = CenterField(non_flat_grid)
8×8×100 Field{Center, Center, Center} on LatitudeLongitudeGrid on CPU
├── grid: 8×8×100 LatitudeLongitudeGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: ZeroFlux, east: ZeroFlux, south: ZeroFlux, north: ZeroFlux, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 14×14×106 OffsetArray(::Array{Float64, 3}, -2:11, -2:11, -2:103) with eltype Float64 with indices -2:11×-2:11×-2:103
    └── max=0.0, min=0.0, mean=0.0

julia> set!(c_src, (λ, φ, z) -> rand())
8×8×100 Field{Center, Center, Center} on LatitudeLongitudeGrid on CPU
├── grid: 8×8×100 LatitudeLongitudeGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: ZeroFlux, east: ZeroFlux, south: ZeroFlux, north: ZeroFlux, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 14×14×106 OffsetArray(::Array{Float64, 3}, -2:11, -2:11, -2:103) with eltype Float64 with indices -2:11×-2:11×-2:103
    └── max=0.999998, min=1.2255e-5, mean=0.497663

julia> c_dst = CenterField(flat_grid)
1×1×100 Field{Center, Center, Center} on LatitudeLongitudeGrid on CPU
├── grid: 1×1×100 LatitudeLongitudeGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo and without precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 1×1×106 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:103) with eltype Float64 with indices 1:1×1:1×-2:103
    └── max=0.0, min=0.0, mean=0.0

julia> Oceananigans.Fields.interpolate!(c_dst, c_src)
ERROR: BoundsError: attempt to access Tuple{Float64} at index [2]
Stacktrace:
  [1] indexed_iterate
    @ ./tuple.jl:92 [inlined]
  [2] _fractional_indices
    @ ~/Projects/Oceananigans.jl/src/Fields/interpolate.jl:163 [inlined]
  [3] fractional_indices
    @ ~/Projects/Oceananigans.jl/src/Fields/interpolate.jl:152 [inlined]
  [4] interpolate
    @ ~/Projects/Oceananigans.jl/src/Fields/interpolate.jl:222 [inlined]
  [5] macro expansion
    @ ~/Projects/Oceananigans.jl/src/Fields/interpolate.jl:302 [inlined]

[glwagner]

Another possibility is to build two grids, then use one for interpolation and then manually copy the data onto a new field generated on the Flat grid. But this is complicated.

Another question is whether this should work for grids that are symmetric about the equator. It seems like it should.

[glwagner]

I think we can close this PR and #3465 because the fix was already merged in #3450 (sorry for putting that in silently, I was just trying to get the code working). As @simone-silvestri points out, simulations that are Periodic-in-latitude may be physically incorrect. However, the code will run and does not necessary blow up immediately, for example:

toroidal_world.mp4

I think we could consider adding a warning for this case in a future PR, but we can discuss that elsewhere.

[glwagner]

Code for the above here:

using Oceananigans
using Oceananigans.Units
using Oceananigans.BoundaryConditions: fill_halo_regions!
using GLMakie

Δ = 30
δ = 10

grid = LatitudeLongitudeGrid(size = (128, 128, 1),
                             longitude = (-Δ, Δ),
                             latitude = (-Δ, Δ),
                             z = (-1, 0),
                             topology = (Periodic, Periodic, Bounded))

model = HydrostaticFreeSurfaceModel(; grid,
                                    coriolis = HydrostaticSphericalCoriolis(),
                                    free_surface = SplitExplicitFreeSurface(cfl=0.7; grid),
                                    tracer_advection = nothing,
                                    momentum_advection = WENO())

Δx = minimum_xspacing(grid)
ψᵢ(λ, φ, z) = exp(-(λ^2 + φ^2) / (2δ^2))
ψ = Field{Face, Face, Center}(grid)
set!(ψ, ψᵢ)
fill_halo_regions!(ψ)
set!(model, u=-∂y(ψ), v=∂x(ψ))

#uᵢ(λ, φ, z) = exp(-(λ^2 + φ^2) / (2δ^2))
#set!(model, u=uᵢ)

simulation = Simulation(model, Δt=20minutes, stop_iteration=1000)

ut = []
vt = []
t = []

function save_uvt(sim)
    push!(ut, deepcopy(interior(sim.model.velocities.u, :, :, 1)))
    push!(vt, deepcopy(interior(sim.model.velocities.v, :, :, 1)))
    push!(t, time(sim))
    @info string("Iter: ", iteration(sim), ", time: ", prettytime(sim))
    return nothing
end

add_callback!(simulation, save_uvt, IterationInterval(10))

run!(simulation)

fig = Figure(size=(800, 400))
axu = Axis(fig[1, 1])
axv = Axis(fig[1, 2])
Nt = length(t)
n = Observable(1)
un = @lift ut[$n]
vn = @lift vt[$n]

heatmap!(axu, un)
heatmap!(axv, vn)

record(fig, "toroidal_world.mp4", 1:Nt) do nn
    @info "Drawing frame $nn of $Nt..."
    n[] = nn
end

[glwagner]

An inertial oscillation solution at a midlatitude location also illustrates potential concerns about physical validity:

generated by this code:

using Oceananigans
using Oceananigans.Units
using GLMakie

φ = 50
δ = 0.01

grid = LatitudeLongitudeGrid(size = (1, 1, 1),
                             longitude = (-δ, δ),
                             latitude = (φ-δ, φ+δ),
                             z = (-1, 0),
                             topology = (Periodic, Periodic, Bounded))

model = HydrostaticFreeSurfaceModel(; grid,
                                    tracers = (),
                                    buoyancy = nothing,
                                    coriolis = HydrostaticSphericalCoriolis(),
                                    tracer_advection = nothing,
                                    momentum_advection = nothing)

@show model

set!(model, u=1)

simulation = Simulation(model, Δt=20minutes, stop_iteration=200)

ut = Float64[]
vt = Float64[]
t = Float64[]

function save_uvt(sim)
    push!(ut, first(sim.model.velocities.u))
    push!(vt, first(sim.model.velocities.v))
    push!(t, time(sim))
    return nothing
end

add_callback!(simulation, save_uvt)

run!(simulation)

Ω = model.coriolis.rotation_rate
f₀ = 2Ω * sind(φ)

fig = Figure()
ax = Axis(fig[1, 1])
lines!(ax, t, ut, linestyle=:dash, label="u")
lines!(ax, t, vt, linestyle=:dash, label="v")
lines!(ax, t, cos.(f₀ * t), linewidth=4, label="cos(ft)")
axislegend(ax)

though a bit more thinking about the math would be needed for a complete interpretation. Nevertheless, it's not an f-plane.

[navidcy]

Since #3465 was closed, do we still need this PR open?

[glwagner]

no we don't