-
Notifications
You must be signed in to change notification settings - Fork 191
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
Split-Explicit free-surface solver interface revisit #3311
Conversation
At the moment, julia> using Oceananigans
julia> grid = RectilinearGrid(size = (10, 10), x = (-100, 100), z = (-100, 0), topology = (Periodic, Flat, Bounded))
10×1×10 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── Periodic x ∈ [-100.0, 100.0) regularly spaced with Δx=20.0
├── Flat y
└── Bounded z ∈ [-100.0, 0.0] regularly spaced with Δz=10.0
julia> free_surface = SplitExplicitFreeSurface(; grid, cfl = 0.7)
SplitExplicitFreeSurface with Oceananigans.Models.HydrostaticFreeSurfaceModels.FixedTimeStepSize{Float64, typeof(Oceananigans.Models.HydrostaticFreeSurfaceModels.averaging_shape_function)}(0.44706193949348216, Oceananigans.Models.HydrostaticFreeSurfaceModels.averaging_shape_function) doesn't error! (A show method is in order though...) |
Why here we have Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl Line 276 in afb4661
Can't this result to type instabilities? |
Δx⁻² = topology(grid)[1] == Flat ? 0 : 1 / minimum_xspacing(grid)^2 | ||
Δy⁻² = topology(grid)[2] == Flat ? 0 : 1 / minimum_yspacing(grid)^2 | ||
Δs = sqrt(1 / (Δx⁻² + Δy⁻²)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@simone-silvestri is this the proper way to compute the spacing for the CFL?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we should make this a small function on its own since users might wanna use this in their scripts?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
the function should handle Flat
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
minimum_xspacing
should handle Flat
. You shouldn't need to check manually.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
True!!
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hm. I think we stumbled on this before. What’s the miminum xspacing for Flat x? 0 I presume.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Inf
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
this use case illustrates it perfectly. The neutral value for a minimum calculation is Inf
(everything is less than Inf
).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
At the moment, for Flat dimensions the spacings have the value f 1 in the grids:
Oceananigans.jl/src/Grids/grid_generation.jl
Lines 124 to 126 in 1c2a6f8
# Flat domains | |
generate_coordinate(FT, ::Flat, N, H, coord::Tuple{<:Number, <:Number}, coordinate_name, arch) = | |
FT(1), range(1, 1, length=N), range(1, 1, length=N), FT(1), FT(1) |
Then, when x/y/zminimum_spacing
is called for those Flat dimension the 1 is returned...
We should really have 0 or Inf as the spacing for Flat dimensions, right?
I think the positional argument should be |
I think so too. @simone-silvestri is there any reason it was done this way? |
Sure that is possible. |
It's following the pattern for other constructors (but in those other cases, We haven't written these conventions out. |
Closing this in favor of #3503 |
Closes #3238