Add metric calculation in the face grid of the cubed sphere + general clean up of CubedSphere code

[navidcy]

This is a first start on a cleanup on the cubed sphere code.

This PR develops a way to compute the metrics for an OrthogonalSphericalShellGrid using the haversine formula for the central angle between two points on the sphere and also a formula for computing the area of a spherical triangle.

Closes #1586

[navidcy]

021b3f1 made some progress on computing area elements for the conformal cubed sphere grid.

E.g., the area elements at locations (Center, Center, Any) should sum up to 1/6-th of the total surface area of the sphere...

julia> grid = OrthogonalSphericalShellGrid(size=(16, 16, 1), z=(-1, 0), radius=1); abs(sum(grid.Azᶜᶜᵃ) - 4π / 6) / (4π / 6)
0.0009221387464206299

julia> grid = OrthogonalSphericalShellGrid(size=(160, 160, 1), z=(-1, 0), radius=1); abs(sum(grid.Azᶜᶜᵃ) - 4π / 6) / (4π / 6)
8.970705743640503e-6

[navidcy]

julia> grid = OrthogonalSphericalShellGrid(size=(16, 16, 1), z=(-1, 0), radius=1);

julia> abs.((sum(grid.Azᶜᶜᵃ), sum(grid.Azᶠᶜᵃ), sum(grid.Azᶜᶠᵃ), sum(grid.Azᶠᶠᵃ)) .- 4π/6)
(0.00047706407919712035, 0.0006835619385658909, 0.0006835619385658909, 0.10180181116525588)

julia> grid = OrthogonalSphericalShellGrid(size=(256, 256, 1), z=(-1, 0), radius=1);

julia> abs.((sum(grid.Azᶜᶜᵃ), sum(grid.Azᶠᶜᵃ), sum(grid.Azᶜᶠᵃ), sum(grid.Azᶠᶠᵃ)) .- 4π/6)
(2.7476636894618878e-6, 2.8613381126696424e-6, 2.8613381126696424e-6, 0.006253695587488561)

[navidcy]

@glwagner, @simone-silvestri: Question: How is Δx^faa, defined for the F at the edge of the domain?

Is it the distance between C[1] - F[1]? (don't think so...)
Or distance between C[1] - C[0]?

For a curvilinear grid with variable distances (as the OrthogonalSphericalShellGrid we need to compute where the halo faces and centers would be if we actually extended the grid, right? We can't just replicate the C[2] - C[1] distances Hx times...

[simone-silvestri]

Curvilinear and Rectilinear coordinates are done in the same way here

Oceananigans.jl/src/Grids/grid_generation.jl

Lines 100 to 104 in c568df3

	F = range(FT(F₋), FT(F₊), length = TF)
	C = range(FT(C₋), FT(C₊), length = TC)
	F = OffsetArray(F, -H)
	C = OffsetArray(C, -H)

I think for the CubedSphere we need to exchange the metrics as halos between the faces (we have to define a field of metrics and fill_halo_regions! to have to correct metrics in the halos)

[simone-silvestri]

We can do that just as a construction step, then Δxᶜᶜᵃ = Δxᶜᶜᵃ.data[:, :, 1] if we want to have the metrics as a 2D offsetarray

[navidcy]

I think for the CubedSphere we need to exchange the metrics as halos between the faces (we have to define a field of metrics and fill_halo_regions! to have to correct metrics in the halos)

This make sense!

[simone-silvestri]

So this solution means that we have to think about boundary conditions and halo exchange for a Field on a CubedSphere before completing the grid

[navidcy]

So this solution means that we have to think about boundary conditions and halo exchange for a Field on a CubedSphere before completing the grid

Yeap. I'll leave the halo metrics zero for now and think about halo filling across faces ;)

[navidcy]

The MIT grid seems to just copy the metrics from the first cell...

julia> grid = ConformalCubedSphereGrid(cs32_filepath, CPU(), Nz=1, z=(-1, 0));
┌ Warning: ConformalCubedSphereGrid is experimental: use with caution!
└ @ Oceananigans.CubedSpheres ~/Research/OC.jl/src/CubedSpheres/conformal_cubed_sphere_grid.jl:163

julia> fgrid = grid.faces[1]
32×32×1 OrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded} on CPU with 1×1×1 halo and with precomputed metrics
├── longitude: FullyConnected λ ∈ [-45.0, 44.0312) variably spaced with min(Δλ)=80152.3, max(Δλ)=3.23804e5
├── latitude:  FullyConnected φ ∈ [-45.0, 42.2888) variably spaced with min(Δφ)=80152.3, max(Δφ)=3.23804e5
└── z:         Bounded  z ∈ [-1.0, 0.0]            regularly spaced with Δz=1.0

julia> fgrid.Δxᶠᶠᵃ[0:3, 0:3]
4×4 Matrix{Float64}:
 80152.3        80152.3        1.36265e5  170108.0
 80152.3        80152.3        1.36265e5  170108.0
     1.33313e5      1.33313e5  1.51133e5       1.76175e5
     1.69174e5      1.69174e5  1.75781e5       1.90232e5

[navidcy]

doi:10.1002/qj.4230 seems relevant for halo filling. They talk about the gnomonic cubed sphere -- not sure if same properties that they utilize apply to the conformal.

[navidcy]

@simone-silvestri, from symmetry arguments, the metrics on the halo on one face should be the same as the metrics on some other face. E.g., Δxᶠᶠᵃ[0] on one face should be either Δxᶠᶠᵃ[grid.Nx+1] or Δyᶠᶠᵃ[grid.Ny+1] on the other face (or something like that, depending on whether a rotation applies). Right?

So we could hardcode these in from a single face instead of constructing 6 faces and filling halos to fill the halo metrics. It will assume, thought, that a specific face configuration is implied.

[simone-silvestri]

Well, we could also hardcode them, as all six faces should be equal. We just have to mind the rotations and the inter-face divisions

[navidcy]

OK, this is not perfect but there has been a lot of improvements. I think to move forward we need to properly fill halos! I'd like to merge this and then continue over at https://github.com/CliMA/Oceananigans.jl/tree/jmc-ss/cubed_sphere

cc @simone-silvestri @jm-c

[navidcy]

A bunch of utility functions in the orthogonal_spherical_shell_grid.jl are general and perhaps worth moving them to Utils module? E.g., the central_angle or the spherical_triangle_area.

[simone-silvestri]

probably we are going to have to change this when we subdivide the face in multiple regions

[navidcy]

okie!

[simone-silvestri]

A bunch of utility functions in the orthogonal_spherical_shell_grid.jl are general and perhaps worth moving them to Utils module? E.g., the central_angle or the spherical_triangle_area.

Sounds good!

[navidcy]

A bunch of utility functions in the orthogonal_spherical_shell_grid.jl are general and perhaps worth moving them to Utils module? E.g., the central_angle or the spherical_triangle_area.

Sounds good!

Will leave this for a future PR.

[navidcy]

I'm merging. The document failing is irrelevant. It has to do with the size of gh-pages on the OceananigansDocumentation repo which has probably gone overboard once again...