Some issue(s) with SplitExplicitFreeSurface

[navidcy]

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

Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl

Lines 305 to 307 in 29a99a0

	if (!isnothing(substeps) && !isnothing(cfl)) || (isnothing(substeps) && isnothing(cfl))
	throw(ArgumentError("either specify a cfl or a number of substeps"))
	end

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

Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface_kernels.jl

Line 324 in 29a99a0

@inline calculate_substeps(substepping::FTS, Δt) = ceil(Int, 2 * Δt / settings.Δtᴮ)

since settings is not provided to calculate_substeps.

cc @djlikesdjs

[glwagner]

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]

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

[navidcy]

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

[glwagner]

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]

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]

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