CliMA / ClimaCore.jl

CliMA model dycore
https://clima.github.io/ClimaCore.jl/dev
Apache License 2.0
80 stars 8 forks source link

CubeMesh and CubeTopology #141

Closed simonbyrne closed 2 years ago

simonbyrne commented 2 years ago

A CubedMesh should represent a nxn mesh on each facet of a cube, with a point on each face having XYPoint coordinates in the interval [-1,1].

Should match orientation used by Oceananigans.jl

A CubedTopology should match the topology interface, with the addition of the facet function mentioned in #75

Zhengyu-Huang commented 2 years ago

Do you have a branch for this?

simonbyrne commented 2 years ago

No, it's a tracking issue

glwagner commented 2 years ago

Should match orientation used by Oceananigans.jl

By "orientation", do you mean the way each facet is connected to one another?

simonbyrne commented 2 years ago

Yes, is this described anywhere?

glwagner commented 2 years ago

We call this the "face connectivity". The default_face_connectivity is defined here:

https://github.com/CliMA/Oceananigans.jl/blob/1786569778c08eb24ba2948f52205a14b973e060/src/CubedSpheres/conformal_cubed_sphere_grid.jl#L36

Here's the comment at the beginning of that function:

function default_face_connectivity()
    # See figure 8.4 of https://mitgcm.readthedocs.io/en/latest/phys_pkgs/exch2.html?highlight=cube%20sphere#fig-6tile
    #
    #                         face  F5   face  F6
    #                       +----------+----------+
    #                       |    ↑↑    |    ↑↑    |
    #                       |    1W    |    1S    |
    #                       |←3N F5 6W→|←5E F6 2S→|
    #                       |    4N    |    4E    |
    #              face  F3 |    ↓↓    |    ↓↓    |
    #            +----------+----------+----------+
    #            |    ↑↑    |    ↑↑    |
    #            |    5W    |    5S    |
    #            |←1N F3 4W→|←3E F4 6S→|
    #            |    2N    |    2E    |
    #            |    ↓↓    |    ↓↓    |
    # +----------+----------+----------+
    # |    ↑↑    |    ↑↑    | face  F4
    # |    3W    |    3S    |
    # |←5N F1 2W→|←1E F2 4S→|
    # |    6N    |    6E    |
    # |    ↓↓    |    ↓↓    |
    # +----------+----------+
    #   face  F1   face  F2

It references this figure in MITgcm:

https://mitgcm.readthedocs.io/en/latest/phys_pkgs/exch2.html?highlight=cube%20sphere#fig-6tile

I was hoping to generalize the implementation in Oceananigans to permit arbitrary connections between different "faces" or "facets". I also think it might make sense to call these "regions" rather than "face" --- we also use "face" as shorthand for the interface between cells, and things get confusing fast.

We should be able to separate the concerns of having a simulation on a "multi-region grid" or "multi-region mesh" from the specific connectedness associated with a cubed sphere. Also there are different region connectivities that might correspond to the cubed sphere.

simonbyrne commented 2 years ago

yes, i think we agreed in #75 to call them "panels".

1) what's the orientation within each panel: i.e. is x always horizontal, y always vertical in the picture you showed? 2) how do you globally orient the sphere? which is the pole, which is the prime meridian, etc.

glwagner commented 2 years ago

Ah, I didn't see #75 . I'm concerned that "panel" seems more specific to a spherical topology... but maybe its ok...

The main thing to note I think is that this abstraction is not specific to spherical domains. We can have different regions laid out on a plane to do simulations in, for example, a "t-shaped" domain. Or, we can have two grids joined with different resolutions along the parallel axis. It's the metric terms and relative rotation of adjoined domains that produces a spherical topology --- a separate concern from simply passing information between two regions / panels.

EDIT: I expressed my preference for "patch" among the options presented, but I am ok with panel as well since I think it covers the planar case described above.

glwagner commented 2 years ago
  1. what's the orientation within each panel: i.e. is x always horizontal, y always vertical in the picture you showed?

Yes, "x" is horizontal and "y" is vertical. This means that local to a panel/region, "x" can be either west-east or south-north, for example.

how do you globally orient the sphere? which is the pole, which is the prime meridian, etc.

I think @christophernhill might have to answer that. Each face/panel/region defines a local grid and coordinates are contained there. I'm not sure it matters; we may want the ability to shift the panels relative to global coordinates to "hide" corners under continents.

simonbyrne commented 2 years ago

Yes, "x" is horizontal and "y" is vertical. This means that local to a panel/region, "x" can be either west-east or south-north, for example.

I guess the question I'm getting at is: if i was to number the elements from 1 to 6n^2, what is the order you traverse them in?

glwagner commented 2 years ago

Yes, "x" is horizontal and "y" is vertical. This means that local to a panel/region, "x" can be either west-east or south-north, for example.

I guess the question I'm getting at is: if i was to number the elements from 1 to 6n^2, what is the order you traverse them in?

For each panel, the ordering is x, y, z so x is fastest (first index), y is second index, and z is third index. The panels are then tupled together, but we could use a different data structure in the future.

For each panel, the first index (x) varies left-to-right in the ASCII diagram. The second index (y) varies down-to-up. Does that make sense?

I think for the grid we currently have, the pole is located on face 3. @christophernhill can confirm.

simonbyrne commented 2 years ago

Thanks, that makes sense.

glwagner commented 2 years ago

Here's some code to probe the properties of the MITgcm-derived cubed sphere grid:

using Pkg
pkg"add Oceananigans DataDeps"

using Oceananigans, DataDeps

dd32 = DataDep("cubed_sphere_32_grid",
    "Conformal cubed sphere grid with 32×32 grid points on each face",
    "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/cubed_sphere_grids/cubed_sphere_32_grid.jld2",
    "b1dafe4f9142c59a2166458a2def743cd45b20a4ed3a1ae84ad3a530e1eff538" # sha256sum
)

DataDeps.register(dd32)

grid_filepath = datadep"cubed_sphere_32_grid/cubed_sphere_32_grid.jld2"
grid = ConformalCubedSphereGrid(grid_filepath, Nz=1, z=(-4000, 0))

I then find

julia> minimum(grid.faces[3].φᶜᶜᵃ)
0.0

and

julia> minimum(grid.faces[6].φᶜᶜᵃ)
-87.9406638719625
OsKnoth commented 2 years ago

We should follow the work by Sriharsha. The Orientation should than follow A PARALLEL EDGE ORIENTATION ALGORITHM FOR QUADRILATERAL MESHES M. HOMOLYA AND D. A. HAM SIAM J. SCI. COMPUT. 2016, Vol. 38, No. 5, pp. S48–S61

valeriabarra commented 2 years ago

Do you have a branch for this?

I had started a cubed-sphere branch. But at this point it sounds like we are not taking this direction anymore.

szy21 commented 2 years ago

Could we close this?

simonbyrne commented 2 years ago

I think so.