Add steady-state test case for the cubed sphere #3302

siddharthabishnu · 2023-10-02T01:41:18Z

This PR introduces a new steady-state test case tailored for the cubed sphere. It initializes surface elevation and velocities using a steady-state solution, aiming to study the temporal evolution arising from the numerical discretization errors. Smooth spatial and temporal error increments are expected. However, if the errors escalate at an unprecedented rate between panels with non-trivial connectivities, it indicates a potential bug, such as improperly filled halos.

Added the validation/multi_region/cubed_sphere_steady_state.jl script. Future commits will initialize surface elevation and velocities using a steady-state solution, aiming to study the temporal evolution arising from the numerical discretization errors. Smooth spatial and temporal error increments are expected. However, if the errors escalate at an unprecedented rate between panels with non-trivial connectivities, it indicates a potential bug, such as improperly filled halos.

Successfully initialized the surface elevation and velocities, but encountered a bug arising from the incorrect implementation of the _fill_south_halos! function.

siddharthabishnu · 2023-10-02T03:32:11Z

Executing the validation/multi_region/cubed_sphere_steady_state.jl file yields an error partially shown below:

julia> include("validation/multi_region/cubed_sphere_steady_state.jl")
[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 41.133 ms, wall time: 0 seconds
[ Info:     ... simulation initialization complete (820.774 ms)
[ Info: Executing initial time step...
ERROR: LoadError: MethodError: no method matching _fill_south_halo!(::Int64, ::Int64, ::Oceananigans.Grids.ZRegOrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded, OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Float64, NamedTuple{(:ξ, :η, :rotation), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}, Rotations.RotXY{Float64}}}, CPU}, ::OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, ::BoundaryCondition{Oceananigans.BoundaryConditions.MultiRegionCommunication, Oceananigans.MultiRegion.CubedSphereRegionalConnectivity{Oceananigans.MultiRegion.South, Oceananigans.MultiRegion.North, Nothing}}, ::Tuple{Center, Center, Face})

Closest candidates are:
  _fill_south_halo!(::Any, ::Any, ::Any, ::Any, ::BoundaryCondition{<:Oceananigans.BoundaryConditions.Open}, ::Any, ::Any...)
   @ Oceananigans /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/cubed-sphere-steady-state/src/BoundaryConditions/fill_halo_regions_open.jl:36
  _fill_south_halo!(::Any, ::Any, ::Any, ::Any, ::BoundaryCondition{<:Oceananigans.BoundaryConditions.Flux}, ::Any...)
   @ Oceananigans /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/cubed-sphere-steady-state/src/BoundaryConditions/fill_halo_regions_flux.jl:32
  _fill_south_halo!(::Any, ::Any, ::Any, ::Any, ::Union{BoundaryCondition{<:Oceananigans.BoundaryConditions.Value}, BoundaryCondition{<:Oceananigans.BoundaryConditions.Gradient}}, ::Any, ::Any...)
   @ Oceananigans /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/cubed-sphere-steady-state/src/BoundaryConditions/fill_halo_regions_value_gradient.jl:44
  ...

Stacktrace:
  [1] macro expansion
    @ /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/cubed-sphere-steady-state/src/BoundaryConditions/fill_halo_regions.jl:174 [inlined]
  [2] cpu__fill_south_and_north_halo!
    @ ~/.julia/packages/KernelAbstractions/cWlFz/src/macros.jl:276 [inlined]
  [3] cpu__fill_south_and_north_halo!(__ctx__::KernelAbstractions.CompilerMetadata{Oceananigans.Utils.OffsetStaticSize{(1:30, 2:2)}, KernelAbstractions.NDIteration.NoDynamicCheck, CartesianIndex{2}, Nothing, KernelAbstractions.NDIteration.NDRange{2, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, KernelAbstractions.NDIteration.StaticSize{(30, 1)}, Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, c::OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, south_bc::BoundaryCondition{Oceananigans.BoundaryConditions.MultiRegionCommunication, Oceananigans.MultiRegion.CubedSphereRegionalConnectivity{Oceananigans.MultiRegion.South, Oceananigans.MultiRegion.North, Nothing}}, north_bc::BoundaryCondition{Oceananigans.BoundaryConditions.MultiRegionCommunication, Oceananigans.MultiRegion.CubedSphereRegionalConnectivity{Oceananigans.MultiRegion.North, Oceananigans.MultiRegion.West, Oceananigans.MultiRegion.↺}}, loc::Tuple{Center, Center, Face}, grid::Oceananigans.Grids.ZRegOrthogonalSphericalShellGrid{Float64, FullyConnected, FullyConnected, Bounded, OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Float64, NamedTuple{(:ξ, :η, :rotation), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}, Rotations.RotXY{Float64}}}, CPU}, args::Tuple{})
    @ Oceananigans.BoundaryConditions ./none:0

Upon inspection, the root of the problem appears to stem from the type specification of BoundaryCondition in the methods _fill_west_halo!, _fill_east_halo!, _fill_south_halo!, and _fill_north_halo!. I attempted to resolve the issue by adding MCBC to the union of type specifications in the pertinent methods. However, this did not rectify the error. It became evident that numerous other locations required modifications.

Moreover, I realized we might be invoking the incorrect methods. Specifically, the methods fill_west_halo!, fill_east_halo!, fill_south_halo!, and fill_north_halo! defined in src/MultiRegion/multi_region_boundary_conditions.jl should be referenced. Instead, we are calling their underscore-prefixed counterparts from src/BoundaryConditions scripts.

@navidcy, @glwagner and @simone-silvestri, what is the most effective strategy to address and resolve this issue?

validation/multi_region/cubed_sphere_steady_state.jl

navidcy · 2023-10-02T10:09:18Z

I'm trying to debug the situation and the code hangs long time at model construction. Is this happening for you @siddharthabishnu?

Even this simpler:

using Oceananigans

grid = ConformalCubedSphereGrid(; panel_size = (10, 10, 2),
                                  z = (-1, 0),
                                  horizontal_direction_halo = 3)

model = HydrostaticFreeSurfaceModel(; grid, momentum_advection = VectorInvariant())

took 10 min to run.

navidcy · 2023-10-02T11:29:40Z

I followed the stack trace and it got triggered from this:

Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl

Lines 135 to 136 in cf4c7e2

    
           # REMEMBER!!! This is going to create problems!!!! 
        
           fill_halo_regions!(ηⁿ⁺¹)

Not sure what the comment above refers to...

Why don't you try split-explicit or some other free surface solver?

siddharthabishnu · 2023-10-02T12:38:58Z

I'm trying to debug the situation and the code hangs long time at model construction. Is this happening for you @siddharthabishnu?

Even this simpler:
using Oceananigans

grid = ConformalCubedSphereGrid(; panel_size = (10, 10, 2),
                                  z = (-1, 0),
                                  horizontal_direction_halo = 3)

model = HydrostaticFreeSurfaceModel(; grid, momentum_advection = VectorInvariant())
took 10 min to run.

Yes, it happened for me too.

navidcy · 2023-10-02T13:37:40Z

Looking at the git blame it turns out that I wrote that comment "REMEMBER!!! This is going to create problems!!!!". But I don't remember why... @simone-silvestri do you?

glwagner · 2023-10-02T14:14:00Z

I followed the stack trace and it got triggered from this:

Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl

Lines 135 to 136 in cf4c7e2

# REMEMBER!!! This is going to create problems!!!!

fill_halo_regions!(ηⁿ⁺¹)

Not sure what the comment above refers to...

Why don't you try split-explicit or some other free surface solver?

HAHA. Remember... trying to remember... rememberr what?

navidcy · 2023-10-02T14:18:22Z

@siddharthabishnu this remark is relevant:

Oceananigans.jl/src/Fields/field.jl

Lines 704 to 720 in 7e411f6

    
           # To correctly fill the halo regions of fields with non-default indices, we'd have to 
        
           # offset indices in the fill halo regions kernels. 
        
           # For now we punt and don't support filling halo regions on windowed fields. 
        
           # Note that `FieldBoundaryConditions` _can_ filter boundary conditions in 
        
           # windowed directions: 
        
           # 
        
           #   filtered_bcs = FieldBoundaryConditions(field.indices, field.boundary_conditions) 
        
           #   
        
           # which will be useful for implementing halo filling for windowed fields in the future. 
        
           fill_halo_regions!(field.data, 
        
                              field.boundary_conditions, 
        
                              field.indices, 
        
                              instantiated_location(field), 
        
                              field.grid, 
        
                              args...; 
        
                              reduced_dimensions = reduced_dims, 
        
                              kwargs...)

since the error you came across comes from filling the halo region of the free surface which is actually a field with non-trivial indices.

glwagner · 2023-10-02T14:21:10Z

True, but free surface only has non-trivial indices in the vertical... is that relevant for cubed sphere?

navidcy · 2023-10-02T14:29:53Z

True, but free surface only has non-trivial indices in the vertical... is that relevant for cubed sphere?

It's not relevant for halo filling... I'm not sure why @siddharthabishnu bumps on this tho.
Question: @siddharthabishnu would you get the same with a MultiRegionGrid?

siddharthabishnu · 2023-10-02T15:49:00Z

@navidcy thanks for the commits. I am now testing the script against a free-explicit free surface and a MultiRegionGrid.

siddharthabishnu · 2023-10-02T17:47:20Z

@navidcy thanks for the commits. I am now testing the script against a free-explicit free surface and a MultiRegionGrid.

@navidcy, as you may have anticipated, a MultiRegionGrid does not result in this error.

Code Modification:

#=
grid = ConformalCubedSphereGrid(; panel_size = (Nx, Ny, Nz),
                                z = (-1, 0),
                                radius = R,
                                horizontal_direction_halo = 6,
                                partition = CubedSpherePartition(; R = 1))
=#

latlongrid = LatitudeLongitudeGrid(size=(Nx, Ny, Nz),
                                   longitude = (-90, 90),
                                   latitude = (-45, 45),
                                   z = (-1, 0))
grid = MultiRegionGrid(latlongrid, partition = XPartition(2))

Terminal Output:

julia> include("validation/multi_region/cubed_sphere_steady_state.jl")
┌ Warning: MultiRegion functionalities are experimental: help the development by reporting bugs or non-implemented features!
└ @ Oceananigans.MultiRegion /Users/Sid/Library/CloudStorage/Dropbox/StudyFolder/PostDocMITDesktop/Codes/Oceananigans/cubed-sphere-steady-state/src/MultiRegion/multi_region_grid.jl:102
[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 7.722 days, wall time: 0 seconds
[ Info:     ... simulation initialization complete (111.288 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (994.424 ms).
[ Info: Simulation is stopping after running for 0 seconds.
[ Info: Simulation time 628.319 ms equals or exceeds stop time 628.319 ms.
[ Info: Making an animation from the saved data...

When utilizing a split-explicit free surface, additional errors arise. For instance, for MultiRegionGrid and ConformalCubedSphereGrid (when specifying the number of substeps), we encounter:

ERROR: LoadError: UndefVarError: `settings` not defined

When specifying both grid and cfl for ConformalCubedSphereGrid, the following error occurs:

ERROR: LoadError: type OrthogonalSphericalShellGrid has no field Lz

Introduced an option to toggle between the conformal cubed sphere grid and a multi-region grid. Interestingly, the latter eliminates the existing bug. However, the lines of code defining the multi-region grid are commented out for now.

Added some lines to specify a split-explicit free surface but commented them out.

glwagner · 2023-10-03T13:27:03Z

When utilizing a split-explicit free surface, additional errors arise. For instance, for MultiRegionGrid and ConformalCubedSphereGrid (when specifying the number of substeps), we encounter:
ERROR: LoadError: UndefVarError: `settings` not defined
When specifying both grid and cfl for ConformalCubedSphereGrid, the following error occurs:
ERROR: LoadError: type OrthogonalSphericalShellGrid has no field Lz

Do you need help fixing these?

siddharthabishnu · 2023-10-04T15:40:12Z

When utilizing a split-explicit free surface, additional errors arise. For instance, for MultiRegionGrid and ConformalCubedSphereGrid (when specifying the number of substeps), we encounter:
ERROR: LoadError: UndefVarError: `settings` not defined
When specifying both grid and cfl for ConformalCubedSphereGrid, the following error occurs:
ERROR: LoadError: type OrthogonalSphericalShellGrid has no field Lz
Do you need help fixing these?

@glwagner thanks for your help on Zoom yesterday. @navidcy thanks for creating a new PR #3305 to fix the second issue.

navidcy · 2023-10-04T17:34:38Z

@siddharthabishnu, perhaps it's simpler to try explicit free surface. Have a look, e.g., at

Oceananigans.jl/validation/multi_region/cubed_sphere_momentum_dynamics.jl

Line 63 in c9aa1b3

free_surface = ExplicitFreeSurface(; gravitational_acceleration),

Major modifications include: (a) Switching to fully explicit free surface; (b) Computing the surface elevation amplitude to set the limits for plots; (c) Extending the simulation time to 10 times longer; (d) Saving the output every time step. The numerical solution blows up due to artifacts at the panel corners.

siddharthabishnu · 2023-10-04T21:34:59Z

@siddharthabishnu, perhaps it's simpler to try explicit free surface. Have a look, e.g., at

Oceananigans.jl/validation/multi_region/cubed_sphere_momentum_dynamics.jl

Line 63 in c9aa1b3

free_surface = ExplicitFreeSurface(; gravitational_acceleration),

Thanks @navidcy. Switching to fully explicit free surface from the default implicit one resolved the boundary condition issue I mentioned in the beginning of this conversation.

Merged the remote-tracking branch 'origin/main' into cubed-sphere-steady-state.

siddharthabishnu added 2 commits September 27, 2023 21:32

Initialized successfully but faced fill_halo bug

cf4c7e2

Successfully initialized the surface elevation and velocities, but encountered a bug arising from the incorrect implementation of the _fill_south_halos! function.

siddharthabishnu force-pushed the sb-jmc/cubed-sphere-steady-state branch from 4f4ace4 to cf4c7e2 Compare October 2, 2023 02:06

navidcy reviewed Oct 2, 2023

View reviewed changes

validation/multi_region/cubed_sphere_steady_state.jl Outdated Show resolved Hide resolved

navidcy added the cubed sphere 🧊🌎 label Oct 2, 2023

navidcy added this to the 🧊 Global simulations on the cubed sphere milestone Oct 2, 2023

navidcy assigned siddharthabishnu Oct 2, 2023

navidcy added 3 commits October 2, 2023 16:49

η is ZFaceField

622456b

cleaner indices for η

a78e935

remove duplicate

7e411f6

siddharthabishnu added 2 commits October 2, 2023 13:59

Added option to use MultiRegionGrid

705d4a6

Introduced an option to toggle between the conformal cubed sphere grid and a multi-region grid. Interestingly, the latter eliminates the existing bug. However, the lines of code defining the multi-region grid are commented out for now.

Added lines to specify split-explicit free surface

cbcd0de

Added some lines to specify a split-explicit free surface but commented them out.

This was referenced Oct 4, 2023

Add more clarification in a comment #3304

Merged

Add Lz property for OrthogonalSphericalShellGrid #3305

Merged

navidcy mentioned this pull request Oct 5, 2023

Momentum dynamics validation on the CubedSphere #3306

Open

siddharthabishnu added a commit that referenced this pull request Oct 5, 2023

Addded steady-state test case from PR #3302

aeb39ed

siddharthabishnu and others added 2 commits October 11, 2023 18:34

Merged main into cubed-sphere-steady-state

caaf3a3

Merged the remote-tracking branch 'origin/main' into cubed-sphere-steady-state.

Merge branch 'main' into sb-jmc/cubed-sphere-steady-state

a81fa0b

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add steady-state test case for the cubed sphere #3302

Add steady-state test case for the cubed sphere #3302

siddharthabishnu commented Oct 2, 2023

siddharthabishnu commented Oct 2, 2023 •

edited

Loading

navidcy commented Oct 2, 2023

navidcy commented Oct 2, 2023

siddharthabishnu commented Oct 2, 2023

navidcy commented Oct 2, 2023

glwagner commented Oct 2, 2023

navidcy commented Oct 2, 2023

glwagner commented Oct 2, 2023

navidcy commented Oct 2, 2023

siddharthabishnu commented Oct 2, 2023 •

edited

Loading

siddharthabishnu commented Oct 2, 2023 •

edited

Loading

glwagner commented Oct 3, 2023

siddharthabishnu commented Oct 4, 2023

navidcy commented Oct 4, 2023

siddharthabishnu commented Oct 4, 2023

Add steady-state test case for the cubed sphere #3302

Are you sure you want to change the base?

Add steady-state test case for the cubed sphere #3302

Conversation

siddharthabishnu commented Oct 2, 2023

siddharthabishnu commented Oct 2, 2023 • edited Loading

navidcy commented Oct 2, 2023

navidcy commented Oct 2, 2023

siddharthabishnu commented Oct 2, 2023

navidcy commented Oct 2, 2023

glwagner commented Oct 2, 2023

navidcy commented Oct 2, 2023

glwagner commented Oct 2, 2023

navidcy commented Oct 2, 2023

siddharthabishnu commented Oct 2, 2023 • edited Loading

siddharthabishnu commented Oct 2, 2023 • edited Loading

glwagner commented Oct 3, 2023

siddharthabishnu commented Oct 4, 2023

navidcy commented Oct 4, 2023

siddharthabishnu commented Oct 4, 2023

siddharthabishnu commented Oct 2, 2023 •

edited

Loading

siddharthabishnu commented Oct 2, 2023 •

edited

Loading

siddharthabishnu commented Oct 2, 2023 •

edited

Loading