search

CliMA / Oceananigans.jl

🌊 Julia software for fast, friendly, flexible, ocean-flavored fluid dynamics on CPUs and GPUs
https://clima.github.io/OceananigansDocumentation/stable
MIT License
991 stars 194 forks source link

Some issue(s) with `SplitExplicitFreeSurface` #3238

Closed navidcy closed 7 months ago

navidcy commented 1 year ago

So initially me and @djlikesdjs were trying to use the SplitExplicitFreeSurface with the adaptive barotropic step based on CFL.

This:

using Oceananigans

grid = RectilinearGrid(size = (10, 10), x = (-100, 100), z = (-100, 0), topology = (Periodic, Flat, Bounded))

free_surface = SplitExplicitFreeSurface(; grid, cfl = 0.7)

errors

ERROR: ArgumentError: either specify a cfl or a number of substeps
Stacktrace:
 [1] Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitSettings(FT::DataType; substeps::Int64, cfl::Float64, grid::RectilinearGrid{Float64, Periodic, Flat, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, fixed_Δt::Nothing, gravitational_acceleration::Float64, averaging_kernel::Function, timestepper::Oceananigans.Models.HydrostaticFreeSurfaceModels.ForwardBackwardScheme)
   @ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/.julia/packages/Oceananigans/pbNSE/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl:306
 [2] SplitExplicitFreeSurface(FT::DataType; gravitational_acceleration::Float64, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:grid, :cfl), Tuple{RectilinearGrid{Float64, Periodic, Flat, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, Float64}}})
   @ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/.julia/packages/Oceananigans/pbNSE/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl:74
 [3] top-level scope
   @ REPL[10]:1

This seems counter-intuitive to me because we had actually prescribed a cfl. We figured that the error was because substeps has a default value of 200 so when we also prescribe a cfl kwarg then

https://github.com/CliMA/Oceananigans.jl/blob/29a99a0c235b2f6bf0cec525f2249125ad254ccc/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl#L305-L307

gets triggered. A solution was to add substeps = nothing. This is a bit counter-intuitive though. Is there a different way we should have done it? If not, perhaps a good idea is to add a note in the SplitExplicitFreeSurface dostring....

Moving on, we got another error down the road. Running this:

using Oceananigans

grid = RectilinearGrid(size = (10, 10), x = (-100, 100), z = (-100, 0), topology = (Periodic, Flat, Bounded))

free_surface = SplitExplicitFreeSurface(; grid, cfl = 0.7, substeps = nothing)

model = HydrostaticFreeSurfaceModel(; grid, free_surface,
                                    buoyancy = BuoyancyTracer(),
                                    tracers = :b,
                                    momentum_advection = WENO(),
                                    tracer_advection = WENO())

simulation = Simulation(model, Δt=10, stop_time=3600)

run!(simulation)

spits out 

ERROR: LoadError: UndefVarError: `settings` not defined
Stacktrace:
  [1] calculate_substeps
    @ ~/.julia/packages/Oceananigans/pbNSE/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface_kernels.jl:324 [inlined]
  [2] iterate_split_explicit!(free_surface::SplitExplicitFreeSurface{Field{Center, Center, Face, Nothing, RectilinearGrid{Float64, Periodic, Flat, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, Tuple{Colon, Colon, UnitRange{Int64}}, OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, Float64, FieldBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{Nothing, Nothing, Nothing, Nothing}}, Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitState{Field{Center, Center, Face, Nothing, RectilinearGrid{Float64, Periodic, Flat, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, Tuple{Colon, Colon, UnitRange{Int64}}, OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, Float64, FieldBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{Nothing, Nothing, Nothing, Nothing}}, Field{Face, Center, Nothing, Nothing, RectilinearGrid{Float64, Periodic, Flat, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, Float64, FieldBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{Nothing, Nothing, Nothing, Nothing}}, Field{Center, Face, Nothing, Nothing, RectilinearGrid{Float64, Periodic, Flat, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, Float64, FieldBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{Nothing, Nothing, Nothing, Nothing}}}, Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitAuxiliaryFields{Field{Center, Face, Nothing, Nothing, RectilinearGrid{Float64, Periodic, Flat, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, Float64, FieldBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{Nothing, Nothing, Nothing, Nothing}}, Field{Face, Center, Nothing, Nothing, RectilinearGrid{Float64, Periodic, Flat, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, Float64, FieldBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{Nothing, Nothing, Nothing, Nothing}}, Field{Center, Center, Nothing, Nothing, RectilinearGrid{Float64, Periodic, Flat, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, Float64, FieldBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, BoundaryCondition{Oceananigans.BoundaryConditions.Flux, Nothing}}, Nothing, Oceananigans.Fields.FieldBoundaryBuffers{Nothing, Nothing, Nothing, Nothing}}, Symbol, Tuple{Int64, Int64}}, Float64, Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitSettings{Oceananigans.Models.HydrostaticFreeSurfaceModels.FixedTimeStepSize{Float64, typeof(Oceananigans.Models.HydrostaticFreeSurfaceModels.averaging_shape_function)}, Oceananigans.Models.HydrostaticFreeSurfaceModels.ForwardBackwardScheme}}, grid::RectilinearGrid{Float64, Periodic, Flat, Bounded, Float64, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, CPU}, Δt::Float64)
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/.julia/packages/Oceananigans/pbNSE/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface_kernels.jl:339

...

 [15] top-level scope
    @ ~/Research/mean_flow_test/mwe.jl:15
 [16] include(fname::String)
    @ Base.MainInclude ./client.jl:478
 [17] top-level scope
    @ REPL[1]:1
in expression starting at /Users/navid/Research/mean_flow_test/mwe.jl:15

Most probably this is due to a bug at

https://github.com/CliMA/Oceananigans.jl/blob/29a99a0c235b2f6bf0cec525f2249125ad254ccc/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface_kernels.jl#L324

since settings is not provided to calculate_substeps.

cc @djlikesdjs

glwagner commented 1 year ago

On the user interface side, @simone-silvestri discussed a refactor that would implement something like:

# To fix the number of substeps at 200
free_surface = SplitExplicitFreeSurface(time_step=FixedNumber(200))

# To fix the number of substeps according to CFL criteria, given the time-step passed to simulation
free_surface = SplitExplicitFreeSurface(time_step=FixedNumber(simulation_Δt=2minutes, cfl=0.7))

# Fixed time-step (variable number of substeps)
free_surface = SplitExplicitFreeSurface(time_step=FixedSize(cfl=0.7))

Do you have any feedback on that? Partly, this was motivated by my own confusion with the API. (Note, we also discussed some internal changes like getting rid of settings.) One detail is that we shouldn't have to pass grid because the free surface already must be "materialized" on the model grid in order to allocate memory for fields. My opinion is that we either should pass grid and cut out the materialization, or we should not pass grid and materialize under the hood. But not both because that's the worst of both words (less convenient for users, complicated under the hood).

We're waiting for #3125 to do this since there are some changes introduced there

navidcy commented 1 year ago

I'd only change the kwarg arg to be barotropic_time_step = FixedNumber(200), etc.

navidcy commented 1 year ago

(@simone-silvestri, until we implement the refactor and/or changes, what's the best way forward using v0.87.1?)

glwagner commented 1 year ago

I don't like barotropic because not everyone knows what that is, not to mention that it can even be confusing to people who know what "barotropic" means, because oceanographers sort of abuse the term. Really, this is the time-step for the free surface solver (further interpretation requires fairly detailed understanding of the physics). What might be a less technical name?

glwagner commented 1 year ago

I thought just saying "time step" was good because the context is already there (it's being specified inside the free surface type, thus we have to delve into the docs for that type to find further information)

simone-silvestri commented 1 year ago

Yeah, there a couple of bugs with the adaptive time step that I meant to fix with #3125 that reworks the internals of the split explicit solver to enable overlayed communication, but that is taking a bit too much, and we want to refactor the constructor (and internal structure) anyway so it's better to do it in a separate (maybe quicker) PR

I can adapt #3125 later