Skip to content

Integrate cubed sphere functionality with halo filling and core infrastructure updates #4538

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

Merged
merged 65 commits into from
Jun 10, 2025

Conversation

siddharthabishnu
Copy link
Contributor

@siddharthabishnu siddharthabishnu commented May 22, 2025

This PR integrates the cubed sphere implementation into the main branch of Oceananigans. While grid generation was addressed in a previous PR (#4295, now merged), this PR focuses on accurately filling the halos of grid metrics and relevant fields---the most computationally challenging aspect of the implementation, and introduces numerous enhancements across several core components, including:

  • MultiRegion grid and field infrastructure
  • CubedSphere grid and field infrastructure
  • Boundary condition handling
  • Implicit and split-explicit barotropic solvers for the free surface
  • Difference operators
  • Conjugate gradient solver

and more. The goal is to enable seamless execution of cubed sphere simulations upon merging this PR. Note that I/O functionality (e.g., output writers, checkpointing, field time series) will be addressed in a subsequent PR.

@siddharthabishnu siddharthabishnu changed the title Sb/cubed sphere integration Integrate cubed sphere functionality with halo filling and core infrastructure updates May 22, 2025
@navidcy
Copy link
Member

navidcy commented May 23, 2025

could you resolve the conflicts?

Comment on lines +79 to +82
#=
field[region][1:Nc, Nc+1:Nc+Hc, k] .= reverse(field[region_N][1:Hc, 1:Nc, k], dims=2)'
field[region][1:Nc, 1-Hc:0, k] .= field[region_S][1:Nc, Nc+1-Hc:Nc, k]
=#
Copy link
Member

@navidcy navidcy May 23, 2025

Choose a reason for hiding this comment

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

why we have commented code here (in similar places further up and down?
do we need it?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

You're right---we don't need these commented lines. I’ve kept them because they show the equivalent non-GPU vectorized implementation, which can be useful for visually verifying halo filling against schematics or physical cubed sphere models. But if you'd prefer a cleaner file, I’m happy to remove them.

Copy link
Member

@navidcy navidcy May 23, 2025

Choose a reason for hiding this comment

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

Either delete or add a comment explaining why they are there or what they demonstrate (in the spirit of what you just explained)?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Comments added in commit c3554fd.

@siddharthabishnu
Copy link
Contributor Author

siddharthabishnu commented May 23, 2025

could you resolve the conflicts?

If you're referring to the conflict that arose after fetching and merging the main branch, the only file affected was src/Models/HydrostaticFreeSurfaceModels/update_hydrostatic_free_surface_model_state.jl. I’ve resolved it in commit 184349c.

launch!(arch, grid, :xy, _implicit_free_surface_linear_operation!,
L_ηⁿ⁺¹, grid, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)

return nothing
end

ImplicitFreeSurfaceOperation = typeof(implicit_free_surface_linear_operation!)

auxiliary_actions!(::ImplicitFreeSurfaceOperation, L_ηⁿ⁺¹, ηⁿ⁺¹, args...) = fill_halo_regions!(ηⁿ⁺¹)
Copy link
Collaborator

Choose a reason for hiding this comment

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

This is used in the solver module right? I am not sure this is very clear. Is there a way to maybe write it in a more streamlined way? Like, for example, by removing the @apply_regionally from line 202 of conjugate_gradient_solver.jl and write a function which in the body:

@inline function mylinearoperation(...)
    @apply_regionally my_actual_linear_operation(...)
    fill_halo_regions!()
end

Copy link
Contributor Author

@siddharthabishnu siddharthabishnu May 27, 2025

Choose a reason for hiding this comment

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

Please see the modifications here: f324fa7...4c4f53e.

Comment on lines 200 to 203
auxiliary_actions!(solver.linear_operation!, q, p, args...)

@apply_regionally solver.linear_operation!(q, p, args...)

Copy link
Collaborator

Choose a reason for hiding this comment

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

see my comment in the Module section

Comment on lines 217 to 218
auxiliary_actions!(args...) = nothing

Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
auxiliary_actions!(args...) = nothing

U = Field{Face, Center, Nothing}(free_surface_grid, boundary_conditions=U_bcs)
V = Field{Center, Face, Nothing}(free_surface_grid, boundary_conditions=V_bcs)
U = Field{Face, Center, Nothing}(free_surface_grid, indices = (:, :, 1:1), boundary_conditions=U_bcs)
V = Field{Center, Face, Nothing}(free_surface_grid, indices = (:, :, 1:1), boundary_conditions=V_bcs)
Copy link
Member

Choose a reason for hiding this comment

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

this has no effect, since the third location is already Nothing

Copy link
Member

Choose a reason for hiding this comment

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

and might possibly be regarded as incorrect

Copy link
Contributor Author

@siddharthabishnu siddharthabishnu Jun 6, 2025

Choose a reason for hiding this comment

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

Without specifying indices = (:, :, 1:1), I was running into an out-of-bounds error. After investigating, I identified the root cause: I wasn’t using reduced_dimensions in the fill_halo_regions! methods tailored for the cubed sphere. Without incorporating reduced_dimensions and without the manual indices workaround, fields defined with Nothing in the third dimension on grids with multiple vertical levels would incorrectly attempt to fill halos across all levels, despite the field having none, resulting in the out-of-bounds error. This issue has now been resolved.

# Note: because of `fill_halo_regions!` below, we cannot use `PCGImplicitFreeSurface` on a
# multi-region grid; `fill_halo_regions!` cannot be used within `@apply_regionally`
fill_halo_regions!(ηⁿ⁺¹)

launch!(arch, grid, :xy, _implicit_free_surface_linear_operation!,
Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
launch!(arch, grid, :xy, _implicit_free_surface_linear_operation!,
@apply_regionally launch!(arch, grid, :xy, _implicit_free_surface_linear_operation!,

@siddharthabishnu siddharthabishnu merged commit 6c2944e into main Jun 10, 2025
58 checks passed
@siddharthabishnu siddharthabishnu deleted the sb/cubed-sphere-integration branch June 10, 2025 18:01
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants