Tripolar grid in Oceananigans (only serial)

[simone-silvestri]

Our implementation of the tripolar grid is in https://github.com/CliMA/OrthogonalSphericalShellGrids.jl at the moment.

We have figured out that the implementation of the tripolar grid is very much entangled with Oceananigans because it requires extending a lot of time-stepping specific code.

This leads to delays in development when we have to develop the two packages in tandem.
Therefore, this PR migrates the implementation of the Tripolar grid and its specific methods to Oceananigans.

Edit: This PR implements only the serial version. The distributed TripolarGrid requires a little more work and will come in subsequent PRs

[navidcy]

nice

we should update

Oceananigans.jl/docs/src/grids.md

Lines 86 to 93 in dc4a327

	1. [`RectilinearGrid`](@ref Oceananigans.Grids.RectilinearGrid) can be fashioned into lines, rectangles and boxes.
	2. [`LatitudeLongitudeGrid`](@ref Oceananigans.Grids.LatitudeLongitudeGrid) represents sectors of thin spherical shells, with cells bounded by lines of constant latitude and longitude.
	3. [`OrthogonalSphericalShellGrid`](@ref Oceananigans.Grids.OrthogonalSphericalShellGrid) represents sectors of thin spherical shells divided with mesh lines that intersect at right angles (thus, orthogonal) but are otherwise arbitrary.
	!!! note "OrthogonalSphericalShellGrids.jl"
	See the auxiliary package [`OrthogonalSphericalShellGrids.jl`](https://github.com/CliMA/OrthogonalSphericalShellGrids.jl)
	for recipes that implement some useful `OrthogonalSphericalShellGrid`, including the
	["tripolar" grid](https://www.sciencedirect.com/science/article/abs/pii/S0021999196901369).

[navidcy]

@simone-silvestri

Which device is meant to be used here?

Oceananigans.jl/src/Grids/tripolar_grid.jl

Line 108 in db49f4d

loop! = _compute_tripolar_coordinates!(device(CPU()), (16, 16), (Nλ, Nφ))

WARNING: both Architectures and CUDA export "device"; uses of it in module Grids must be qualified
Unit tests...: Error During Test at /Users/navid/Research/OC11.jl/test/test_tripolar_grid.jl:31
  Got exception outside of a @test
  UndefVarError: `device` not defined
  Stacktrace:
    [1] (TripolarGrid)(arch::CPU, FT::DataType; size::Tuple{Int64, Int64, Int64}, southernmost_latitude::Int64, halo::Tuple{Int64, Int64, Int64}, radius::Float64, z::Tuple{Int64, Int64}, north_poles_latitude::Int64, first_pole_longitude::Int64)
      @ Oceananigans.Grids ~/Research/OC11.jl/src/Grids/tripolar_grid.jl:108

[glwagner]

Why draft?

[navidcy]

Why draft?

I guess because it needs work... Still most test fail. It's not ready for review...

[glwagner]

Huh true. Just thought it would have been a copy/paste job but maybe there's more to it

[navidcy]

It shouldn't be hard to fix... I can give it a go later!

[simone-silvestri]

the main problem is that we use Fields and fill_halo_regions! in the construction of the grid. I ll search a way around it

[navidcy]

This part:

Oceananigans.jl/src/Grids/tripolar_grid.jl

Lines 139 to 148 in 44f62c1

	# Boundary conditions to fill halos of the coordinate and metric terms
	# We need to define them manually because of the convention in the
	# ZipperBoundaryCondition that edge fields need to switch sign (which we definitely do not
	# want for coordinates and metrics)
	default_boundary_conditions = FieldBoundaryConditions(north = ZipperBoundaryCondition(),
	south = nothing, # The south should be `continued`
	west = Oceananigans.PeriodicBoundaryCondition(),
	east = Oceananigans.PeriodicBoundaryCondition(),
	top = nothing,
	bottom = nothing)

requires BoundaryConditions before grids?

[navidcy]

Oh I just saw your comment @simone-silvestri!

[navidcy]

Damn...

Fields and fill_halo_regions! are also used in the cubed sphere grid also, but that's defined later in the MultiRegion module.

[simone-silvestri]

We could move the zipper-filling halo kernels in grids and use them just as kernels.
Then later import them in boundary conditions and use them for the ZipperBoundaryCondition.
It is a bit convoluted though.

[navidcy]

Or we can define the TripolarGrid in a module on its own after Fields/Boundary conditions? It does make sense that it lives in Grids... aaaaargh

[navidcy]

Btw, here we define a haversine

Oceananigans.jl/src/Grids/tripolar_grid.jl

Lines 457 to 467 in 44f62c1

	# Is this the same as in Oceananigans?
	# TODO: check it out
	function haversine(a, b, radius)
	λ₁, φ₁ = a
	λ₂, φ₂ = b
	x₁, y₁, z₁ = lat_lon_to_cartesian(φ₁, λ₁, radius)
	x₂, y₂, z₂ = lat_lon_to_cartesian(φ₂, λ₂, radius)
	return radius * acos(max(-1.0, min((x₁ * x₂ + y₁ * y₂ + z₁ * z₂) / radius^2, 1.0)))
	end

Elsewhere in Oceananigans we use the Distances.haversine; see https://github.com/JuliaStats/Distances.jl/blob/master/src/haversine.jl

[navidcy]

The two are exactly the same; the snippet below demonstrates that:

using Oceananigans.Grids: lat_lon_to_cartesian
using Distances

# Is this the same as in Oceananigans? 
# TODO: check it out
function haversine_simone(a, b, radius)
    λ₁, φ₁ = a
    λ₂, φ₂ = b

    x₁, y₁, z₁ = lat_lon_to_cartesian(φ₁, λ₁, radius)
    x₂, y₂, z₂ = lat_lon_to_cartesian(φ₂, λ₂, radius)

    return radius * acos(max(-1.0, min((x₁ * x₂ + y₁ * y₂ + z₁ * z₂) / radius^2, 1.0)))
end

using Test

for i in 1:100
    (λ1, φ1) = rand()*360, (2rand()-1)*90
    (λ2, φ2) = rand()*360, (2rand()-1)*90

    radius = rand() * 6300e3

    @test haversine_simone((λ1, φ1), (λ2, φ2), radius) ≈ Distances.haversine((λ1, φ1), (λ2, φ2), radius)
end

[simone-silvestri]

The tests seem to pass

[glwagner]

why different notation than usual

[glwagner]

What's a "nonorthogonality_angle"? Is this the deviation from orthogonality? (So it should be small)

[simone-silvestri]

yeah, it should be small. I am not sure what is acceptable, so I am comparing it here with a cubes sphere