Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fully connected cubed sphere #1522

Merged
merged 87 commits into from
Apr 14, 2021
Merged

Fully connected cubed sphere #1522

merged 87 commits into from
Apr 14, 2021

Conversation

ali-ramadhan
Copy link
Member

@ali-ramadhan ali-ramadhan commented Mar 30, 2021

This PR adds a new submodule Oceananigans.CubedSpheres that implements new grid and field types for running HydrostaticFreeSurfaceModel with a fully connected cubed sphere (should eventually also be flexible enough to just use 1 < n < 6 faces but I'm just testing 6 faces for now).

This PR is still a work-in-progress and is quite messy with some pretty horrible code (please don't judge 😭). I'll be adding to it and refactoring it over the next few days.

Compile times have increased significantly with grid::ConformalCubedSphereGrid which has slowed down development but should be pretty close to being able to time step. EDIT: This might only be on Julia 1.6...

cfl = DiffusiveCFL(FT(Δt))
Δt = FT(1.7)

Δx = model.grid.Δx
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you only use each of these grid spacings once, would it be better to just use the full name where they are used?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

True. I think this change is part of #1497 but I'll revise those tests before merging.

Copy link
Collaborator

@francispoulin francispoulin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did look over all the changes and there certianly is a lot there. I don't feel like I understand all the different components enough to approve but I am curious to see the results and how they compare with predictions

@ali-ramadhan
Copy link
Member Author

Coriolis seems fine now that grid metric halos are filled:

cubed_sphere_surface_gravity_waves_face1_with_coriolis.mp4

@ali-ramadhan ali-ramadhan marked this pull request as ready for review April 13, 2021 13:25
@ali-ramadhan
Copy link
Member Author

ali-ramadhan commented Apr 14, 2021

Tests are on track to pass so I think this PR is finally ready to be reviewed/merged (with a regrettable amount of technical debt).

It's unfortunately too large to be effectively reviewed but I'm happy to make changes if needed (or go over the code on Zoom if anyone would like to do so).

Copy link
Member

@christophernhill christophernhill left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ali-ramadhan looks good for where we currently are. I left a few minor comments.

Is there somewhere to note that this is an intermediate "caveat emptor" that runs but is likely to change quite a bit in the future?

@@ -68,13 +68,7 @@ A = 1e-5 * H # Amplitude of the perturbation

η′(λ, φ, z) = A * exp(- (λ - λ₀)^2 / Δλ^2) * exp(- (φ - φ₀)^2 / Δφ^2)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ali-ramadhan I am getting unresolved method for η′(λ, φ, z) = .... (see below).
η′(λ, φ,) = .... works.

ERROR: LoadError: MethodError: no method matching (::var"#η′#2"{Float64,Int64,Int64,Int64,Int64})(::Float64, ::Float64)
Closest candidates are:
  η′(::Any, ::Any, ::Any) at /Users/chrishill/projects/oceananigans/onan-cs-review/validation/cubed_sphere_surface_gravity_waves/cubed_sphere_surface_gravity_waves.jl:93
Stacktrace:
 [1] set!(::Oceananigans.Fields.ReducedField{Center,Center,Nothing,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},ConformalCubedSphereFaceGrid{Float64,Periodic,Periodic,Bounded,OffsetArrays.OffsetArray{Float64,2,Array{Float64,2}},OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},1,NamedTuple{(:x, :y, :z),Tuple{CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.CubedSpheres.CubedSphereExchange,Oceananigans.CubedSpheres.CubedSphereExchangeInformation{Int64,Symbol}},BoundaryCondition{Oceananigans.CubedSpheres.CubedSphereExchange,Oceananigans.CubedSpheres.CubedSphereExchangeInformation{Int64,Symbol}}},CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.CubedSpheres.CubedSphereExchange,Oceananigans.CubedSpheres.CubedSphereExchangeInformation{Int64,Symbol}},BoundaryCondition{Oceananigans.CubedSpheres.CubedSphereExchange,Oceananigans.CubedSpheres.CubedSphereExchangeInformation{Int64,Symbol}}},CoordinateBoundaryConditions{Nothing,Nothing}}}}, ::var"#η′#2"{Float64,Int64,Int64,Int64,Int64}) at /Users/chrishill/projects/oceananigans/onan-cs-review/src/CubedSpheres/cubed_sphere_set!.jl:54

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah sorry about this! I forgot to change this since model.free_surface.η became a ReducedField. I'll update the validation experiments.


# This field needs BCs otherwise errors happen so I'll assume all faces have
# the same boundary conditions. A very bad assumption...
cubed_sphere_bcs = faces[1].boundary_conditions
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

😸

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah haha, in hindsight I should have went with CubedSphereData and CubedSphereBoundaryConditions structs. But will do it next time we get a chance to refactor.

model = HydrostaticFreeSurfaceModel(
architecture = CPU(),
grid = grid,
momentum_advection = VectorInvariant(),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ali-ramadhan do we want to have this on here - since we know it is broken?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

True. I'll revert back to momentum_advection = nothing for now until VectorInvariant() is correct for the cubed sphere.

model = HydrostaticFreeSurfaceModel(
architecture = CPU(),
grid = grid,
momentum_advection = VectorInvariant(),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ali-ramadhan do we want inertial terms in this setup? We could just do the linear problem for this example?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah I think I was just testing it but will change to momentum_advection = nothing for now, at least until VectorInvariant() is fixed for the cubed sphere.

@ali-ramadhan
Copy link
Member Author

Thanks for taking a look @christophernhill!

Is there somewhere to note that this is an intermediate "caveat emptor" that runs but is likely to change quite a bit in the future?

Yes, we just have a warning that gets printed when you construct a cubed sphere grid saying that's it's an experimental feature:

@warn "ConformalCubedSphereGrid is experimental: use with caution!"

λnode(LX::Face, LY::Center, LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = grid.λᶠᶠᵃ[i, j]
λnode(LX::Center, LY::Face, LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = grid.λᶠᶠᵃ[i, j]
φnode(LX::Face, LY::Center, LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = grid.φᶠᶠᵃ[i, j]
φnode(LX::Center, LY::Face, LZ, i, j, k, grid::ConformalCubedSphereFaceGrid) = grid.φᶠᶠᵃ[i, j]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do these all need @inbounds?

# Not sure how to best represent these so will concatenate along dim 3 for now.

λnodes(LX, LY, LZ, grid::ConformalCubedSphereGrid) = cat(Tuple(λnodes(LX, LY, LZ, grid_face) for grid_face in grid.faces)..., dims=3)
φnodes(LX, LY, LZ, grid::ConformalCubedSphereGrid) = cat(Tuple(φnodes(LX, LY, LZ, grid_face) for grid_face in grid.faces)..., dims=3)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the point is probably to support broadcasting with CubedSphereData...

offset_data = zeros(Nx + 1 + 2Hx, Ny + 1 + 2Hy)
offset_data[1+Hx:Nx+1+Hx, 1+Hy:Ny+1+Hy] .= data
return OffsetArray(offset_data, -Hx, -Hy)
function offset_data_2d(data, loc, topo, N, H)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm... this is really in xy (there are two other two dimensional configurations)

# Launch velocity update kernels

velocities_events = []

for name in (:u, :v)
model.velocities isa PrescribedVelocityFields && break
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function should never be called when using PrescribedVelocityFields. This means there is a bug in the code in prescribed_velocity_fields.jl. This line shouldn't be here...

Copy link
Member

@glwagner glwagner left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I didn't review every file in the PR. Existing code looks largely untouched so I think we are fine. My only comment is that we definitely should not have any code that refers to PrescribedVelocityFields anywhere outside of prescribed_velocity_fields.jl. I'm ok to keep this here in the spirit of moving forward but the underlying issue should be addressed in the future.

@ali-ramadhan
Copy link
Member Author

Thanks for taking a look @glwagner! I'll open a PR to make sure we're working with PrescribedVelocityFields correctly. And I'll also open a second PR to start looking at cubed sphere performance + getting it working on the GPU.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants