Fill metrics for halo regions for OrthogonalSphericalShellGrid

[navidcy]

This PR fills the metrics for the halo regions for a OrthogonalSphericalShellGrid.

Also adds a validation example for a splash on an OrthogonalSphericalShellGrid, i.e., a HydrostaticFreeSurface model starting from rest + its free surface lifted up in the form of a Gaussian.

Partially resolves #3198.

[glwagner]

beautiful

[simone-silvestri]

Does the metric not reflect across the boundaries?

metric[1-i, j] = metric[i, j]

[navidcy]

I don't understand exactly what you mean here @simone-silvestri

[simone-silvestri]

what is the rule for the metrics in the halo? I would expect them to reflect across the boundary following this rule:

metric[1-i, j] = metric[i, j]

but here you seem to be implying that all metrics in the boundary are equal to the last metric in the interior

parent(metric)[i, j] = parent(metric)[i+1, j]

what is the correct rule here? I think we need halo metrics to be filled with a fill halo to be sure because neither of the options is always correct.

[navidcy]

I don’t understand what you mean by “reflect”. In general it is ambiguous how one can fill the metrics and thus I chose to fill them according to their nearest interior value. For FullyConnected topologies, of course, the halo regions correspond to the metrics of the interior neighboring region.

[simone-silvestri]

It should be like that also for Periodic which takes from the other side of the face. In that case (and FullyConnected) the metrics are not so ambiguous. The only instance where metrics in the halos can be considered equal to the interior ones (or ambiguous) is for solid walls and value/gradient BC (where we only need 1 halo). If we don't fill them correctly we might have distortions across the boundary which can crash the simulation/kick off numerical instabilities that grow and pollute the solution (see #3091 where metrics where not correctly considered in reconstructions)

[navidcy]

Yeah, I found that also and I tried to replicate that functionality. I think I did it correctly.

[navidcy]

I couldn't use that lower_exterior_Δcoordᶠ method; unfortunately I don't remember now why... it was work done over a long flight which hints that I may not have come to that conclusion at my full brain capacity. :(

I'll revisit and see if I can use that utility here.

[simone-silvestri]

I wouldn't stress too much about this at the moment. It is kind of a low priority because, in every practical application of the OSSG, the metrics should be filled with a fill_halo_regions! (not only forConnected but also for Periodic topologies. This might be akin to "extending" in some cases). I agree with extending for Bounded halo metrics but they matter only for Value and Gradient boundary conditions

[glwagner]

Aren't we planning to replace lat-lon grid with OSSG? What happened to that plan?

[glwagner]

I personally would like to move forward with a fully function OSSG so I can use it for regional simulations (eg in the arctic) that are not restricted to borders that exactly align with geographic coordinates. It will also be nice to rotate domains to follow coastlines, etc

[navidcy]

@simone-silvestri can you have a look now?

[navidcy]

@simone-silvestri I plot here how the metric look for a series of overlapping conformal cubed sphere panels constructed with this PR.

using Oceananigans
using Oceananigans.Grids
using Oceananigans.MultiRegion: getregion
using GLMakie

Nx, Ny = 8, 8
H = 4

range_with_offset(N, H, offset::Int) = (-H + 1 + offset*N):(N + H + offset*N)

function plot_series_ossg(grid, filename)
    j_index = 2

    fig = Figure(resolution=(2200, 600), fontsize=30)
    ax = Axis(fig[1, 1])
    lines!(ax, range_with_offset(Nx, H, 0), parent(grid.Δxᶜᶜᵃ)[:, j_index+H], linewidth=4, color = (:green, 0.9), label="an ossg")
    lines!(ax, range_with_offset(Nx, H, 1), parent(grid.Δxᶜᶜᵃ)[:, j_index+H], linewidth=8, color = (:red, 0.5), label="another ossg")
    lines!(ax, range_with_offset(Nx, H, 2), parent(grid.Δxᶜᶜᵃ)[:, j_index+H], linewidth=4, color = (:blue, 0.5), label="another ossg")
    vlines!(ax, [1, Nx+1, 2Nx+1, 3Nx+1], linewidth=8, color=(:black, 0.3))
    fig[1, 2] = Legend(fig, ax, framevisible = false)
    save(filename, fig)

    return fig
end

grid = conformal_cubed_sphere_panel(size=(Nx, Ny, 1), z=(-1, 0), halo=(H, H, 1))
plot_series_ossg(grid, "metric_test_ossg_bounded.png")

grid = conformal_cubed_sphere_panel(size=(Nx, Ny, 1), z=(-1, 0), halo=(H, H, 1), topology=(Periodic, Periodic, Bounded))
plot_series_ossg(grid, "metric_test_ossg_periodic.png")

For Bounded topology:

For Periodic topology:

Isn't this what you were trying to imply they should look like?

[glwagner]

This looks great and it's consistent with what generate_coordinate does

Oceananigans.jl/src/Grids/grid_generation.jl

Line 29 in b99b467

function generate_coordinate(FT, topo::AT, N, H, coord, arch)

[navidcy]

@glwagner with 8d8678a I changed to using normal indexing (not parent-array indices). The results of the script above remain unchanged!

Is this OK to merge now?