Add ConformalCubedSphere grid via MultiRegion module

[navidcy]

This PR removes the CubedSpheres module and re-implements CubedSphere grid using the MultiRegion module. Furthermore, it improves the OrthogonalSphericalShellGrid.

Implemented in this PR

• uses Distances package which allowed us to clear up various utilities in the grid_utils.jl file regarding great circle distances on the sphere
• fix some bugs that didn't allow an OrthogonalSphericalShellGrid to be constructed with any topology
• better and more accurate show methods that also work on the GPU.
• construct ConformalCubedSphereGrid + tests
• introduces the Connectivity type for MultiRegionGrids + adds a default CubedSphereGrid connectivity
• adds halo filling functionality for tracers + velocities ( + tests)
• introduced CubedSphereField abstraction
• introduces CubedSphereConnectivity + tests

Outstanding issues

(Several issues will be opened as soon as this PR is merged for the following)

• Coordinate and metric halo fillings for the ConformalCubedSphereGrid are hardcoded. We need to re-implement this based on the connectivity of the grid so that it works (i) for any connectivity and (ii) for any CubedSpherePartition.
• Properties ξₗ, ξᵣ, ηₗ, ηᵣ should be taken out from OrthogonalSphericalShellGrid or, even better, be grouped together into a property conformal_cubed_sphere or something. This way, the OrthogonalSphericalShellGrid will be general and not necessarily tied to the conformal cubed sphere.
• Differentiate the OrthogonalSphericalShellGrid constructors from the cubed_sphere_panel constructors
• Shortcut for velocity halo filling
• Alleviate the need for multiple hallo filling passes for velocities (important for performance/scaling)
• Add testing for MultiRegionGrids with YPartition
• Allow uneven x-y partition for ConformalCubedSphere

[navidcy]

So we need to have some "default connectivity" for the 6 panels, e.g., something like

Oceananigans.jl/src/CubedSpheres/conformal_cubed_sphere_grid.jl

Lines 37 to 115 in 0dc0a8a

	function default_panel_connectivity()
	# Adopted from figure 8.4 of https://mitgcm.readthedocs.io/en/latest/phys_pkgs/exch2.html?highlight=cube%20sphere#fig-6tile
	#
	# panel P5 panel P6
	# +----------+----------+
	# | ↑↑ | ↑↑ |
	# | 1W | 1S |
	# |←3N P5 6W→|←5E P6 2S→|
	# | 4N | 4E |
	# panel P3 | ↓↓ | ↓↓ |
	# +----------+----------+----------+
	# | ↑↑ | ↑↑ |
	# | 5W | 5S |
	# |←1N P3 4W→|←3E P4 6S→|
	# | 2N | 2E |
	# | ↓↓ | ↓↓ |
	# +----------+----------+----------+
	# | ↑↑ | ↑↑ | panel P4
	# | 3W | 3S |
	# |←5N P1 2W→|←1E P2 4S→|
	# | 6N | 6E |
	# | ↓↓ | ↓↓ |
	# +----------+----------+
	# panel P1 panel P2
	panel1_connectivity = CubedSpherePanelConnectivity(
	west = CubedSpherePanelConnectivityDetails(5, :north),
	east = CubedSpherePanelConnectivityDetails(2, :west),
	south = CubedSpherePanelConnectivityDetails(6, :north),
	north = CubedSpherePanelConnectivityDetails(3, :west),
	)
	panel2_connectivity = CubedSpherePanelConnectivity(
	west = CubedSpherePanelConnectivityDetails(1, :east),
	east = CubedSpherePanelConnectivityDetails(4, :south),
	south = CubedSpherePanelConnectivityDetails(6, :east),
	north = CubedSpherePanelConnectivityDetails(3, :south),
	)
	panel3_connectivity = CubedSpherePanelConnectivity(
	west = CubedSpherePanelConnectivityDetails(1, :north),
	east = CubedSpherePanelConnectivityDetails(4, :west),
	south = CubedSpherePanelConnectivityDetails(2, :north),
	north = CubedSpherePanelConnectivityDetails(5, :west),
	)
	panel4_connectivity = CubedSpherePanelConnectivity(
	west = CubedSpherePanelConnectivityDetails(3, :east),
	east = CubedSpherePanelConnectivityDetails(6, :south),
	south = CubedSpherePanelConnectivityDetails(2, :east),
	north = CubedSpherePanelConnectivityDetails(5, :south),
	)
	panel5_connectivity = CubedSpherePanelConnectivity(
	west = CubedSpherePanelConnectivityDetails(3, :north),
	east = CubedSpherePanelConnectivityDetails(6, :west),
	south = CubedSpherePanelConnectivityDetails(4, :north),
	north = CubedSpherePanelConnectivityDetails(1, :west),
	)
	panel6_connectivity = CubedSpherePanelConnectivity(
	west = CubedSpherePanelConnectivityDetails(5, :east),
	east = CubedSpherePanelConnectivityDetails(2, :south),
	south = CubedSpherePanelConnectivityDetails(4, :east),
	north = CubedSpherePanelConnectivityDetails(1, :south),
	)
	panel_connectivity = (
	panel1_connectivity,
	panel2_connectivity,
	panel3_connectivity,
	panel4_connectivity,
	panel5_connectivity,
	panel6_connectivity
	)
	return panel_connectivity
	end

and from that the connectivity amongst all ranks should be inferred, right?

[navidcy]

@simone-silvestri can you resolve the conflicts? I think they came up after merging #2888.

[navidcy]

hm... @simone-silvestri (personal communication) suggested that the culprit was the commented out velocity halo filling in 852381a. But after bringing them back, still tests fail and code snippet in #2867 (comment) gives the same output...

[navidcy]

I admit that I'm not quite satisfied with these:

Oceananigans.jl/test/test_multi_region_cubed_sphere.jl

Lines 9 to 65 in b9c9026

	function get_range_of_indices(operation, index, Nx, Ny)
	if operation == :endpoint && index == :first
	range_x = 1
	range_y = 1
	elseif operation == :endpoint && index == :last
	range_x = Nx
	range_y = Ny
	elseif operation == :subset && index == :first # here index is the index to skip
	range_x = 2:Nx
	range_y = 2:Ny
	elseif operation == :subset && index == :last # here index is the index to skip
	range_x = 1:Nx-1
	range_y = 1:Ny-1
	else
	range_x = 1:Nx
	range_y = 1:Ny
	end
	return range_x, range_y
	end
	function get_halo_data(field, side, k_index=1; operation=nothing, index=:all)
	Nx, Ny, _ = size(field)
	Hx, Hy, _ = halo_size(field.grid)
	range_x, range_y = get_range_of_indices(operation, index, Nx, Ny)
	if side == :west
	return field.data[-Hx+1:0, range_y, k_index]
	elseif side == :east
	return field.data[Nx+1:Nx+Hx, range_y, k_index]
	elseif side == :south
	return field.data[range_x, -Hy+1:0, k_index]
	elseif side == :north
	return field.data[range_x, Ny+1:Ny+Hy, k_index]
	end
	end
	function get_boundary_indices(Nx, Ny, Hx, Hy, side; operation=nothing, index=:all)
	range_x, range_y = get_range_of_indices(operation, index, Nx, Ny)
	if side == :west
	return 1:Hx, range_y
	elseif side == :south
	return range_x, 1:Hy
	elseif side == :east
	return Nx-Hx+1:Nx, range_y
	elseif side == :north
	return range_x, Ny-Hy+1:Ny
	end
	end

tests, mostly because looking at the code I can't understand what's happening -- they are not human-readable. I'll try to fix them. Perhaps @siddharthabishnu could you add docstring with some explanation?

[glwagner]

I admit that I'm not quite satisfied with these:

Oceananigans.jl/test/test_multi_region_cubed_sphere.jl

Lines 9 to 65 in b9c9026

	function get_range_of_indices(operation, index, Nx, Ny)
	if operation == :endpoint && index == :first
	range_x = 1
	range_y = 1
	elseif operation == :endpoint && index == :last
	range_x = Nx
	range_y = Ny
	elseif operation == :subset && index == :first # here index is the index to skip
	range_x = 2:Nx
	range_y = 2:Ny
	elseif operation == :subset && index == :last # here index is the index to skip
	range_x = 1:Nx-1
	range_y = 1:Ny-1
	else
	range_x = 1:Nx
	range_y = 1:Ny
	end
	return range_x, range_y
	end
	function get_halo_data(field, side, k_index=1; operation=nothing, index=:all)
	Nx, Ny, _ = size(field)
	Hx, Hy, _ = halo_size(field.grid)
	range_x, range_y = get_range_of_indices(operation, index, Nx, Ny)
	if side == :west
	return field.data[-Hx+1:0, range_y, k_index]
	elseif side == :east
	return field.data[Nx+1:Nx+Hx, range_y, k_index]
	elseif side == :south
	return field.data[range_x, -Hy+1:0, k_index]
	elseif side == :north
	return field.data[range_x, Ny+1:Ny+Hy, k_index]
	end
	end
	function get_boundary_indices(Nx, Ny, Hx, Hy, side; operation=nothing, index=:all)
	range_x, range_y = get_range_of_indices(operation, index, Nx, Ny)
	if side == :west
	return 1:Hx, range_y
	elseif side == :south
	return range_x, 1:Hy
	elseif side == :east
	return Nx-Hx+1:Nx, range_y
	elseif side == :north
	return range_x, Ny-Hy+1:Ny
	end
	end

tests, mostly because looking at the code I can't understand what's happening -- they are not human-readable. I'll try to fix them. Perhaps @siddharthabishnu could you add docstring with some explanation?

Shouldn't this be implemented with multiple dispatch? Chains of if-statements are the red flag.

[navidcy]

I admit that I'm not quite satisfied with these:

Oceananigans.jl/test/test_multi_region_cubed_sphere.jl

Lines 9 to 65 in b9c9026

	function get_range_of_indices(operation, index, Nx, Ny)
	if operation == :endpoint && index == :first
	range_x = 1
	range_y = 1
	elseif operation == :endpoint && index == :last
	range_x = Nx
	range_y = Ny
	elseif operation == :subset && index == :first # here index is the index to skip
	range_x = 2:Nx
	range_y = 2:Ny
	elseif operation == :subset && index == :last # here index is the index to skip
	range_x = 1:Nx-1
	range_y = 1:Ny-1
	else
	range_x = 1:Nx
	range_y = 1:Ny
	end
	return range_x, range_y
	end
	function get_halo_data(field, side, k_index=1; operation=nothing, index=:all)
	Nx, Ny, _ = size(field)
	Hx, Hy, _ = halo_size(field.grid)
	range_x, range_y = get_range_of_indices(operation, index, Nx, Ny)
	if side == :west
	return field.data[-Hx+1:0, range_y, k_index]
	elseif side == :east
	return field.data[Nx+1:Nx+Hx, range_y, k_index]
	elseif side == :south
	return field.data[range_x, -Hy+1:0, k_index]
	elseif side == :north
	return field.data[range_x, Ny+1:Ny+Hy, k_index]
	end
	end
	function get_boundary_indices(Nx, Ny, Hx, Hy, side; operation=nothing, index=:all)
	range_x, range_y = get_range_of_indices(operation, index, Nx, Ny)
	if side == :west
	return 1:Hx, range_y
	elseif side == :south
	return range_x, 1:Hy
	elseif side == :east
	return Nx-Hx+1:Nx, range_y
	elseif side == :north
	return range_x, Ny-Hy+1:Ny
	end
	end

tests, mostly because looking at the code I can't understand what's happening -- they are not human-readable. I'll try to fix them. Perhaps @siddharthabishnu could you add docstring with some explanation?

Shouldn't this be implemented with multiple dispatch? Chains of if-statements are the red flag.

Yea, another issue of mine.

[siddharthabishnu]

I will add a docstring with some explanation, and replace some of the if-statements with multiple dispatch.

[navidcy]

I will add a docstring with some explanation, and replace some of the if-statements with multiple dispatch.

What will you dispatch on though?
A docstring is for sure needed. Thanks!

[siddharthabishnu]

I will add a docstring with some explanation, and replace some of the if-statements with multiple dispatch.

What will you dispatch on though? A docstring is for sure needed. Thanks!

I see your point. For dispatch, the argument types need to be different, which is not the case here. I used so many if statements to minimize the number of functions. I think I will just rewrite them in more readable way.

[navidcy]

issues: #3198, #3199, #3200, #3201, #3202, #3203 were opened as a response of this PR.