Computing CFL and using TimeStepWizards with ShallowWaterSimulation

[pjrusello]

I have been playing around with the Bickley Jet example and tried adding a TimeStepWizard and an output message with a CFL calculation.

I saw in issue #2615 there are potentially some issues around model.velocities with ShallowWaterModel, so I was just wondering if calculating a CFL like this

advective_cfl = AdvectiveCFL( 1.0 )
timestep_cfl = advective_cfl( simulation.model )

Should work? I get Inf for CFL with this and the TimeStepWizard tries to compensate by setting the time step smaller and smaller, so something isn't working right for me currently.

[pjrusello]

Reading through some of the other issues it seems like there may still be some issues with how model.velocities are calculated for the conservative formulation that probably play into how a CFL calculator will function. It also seems like adding a closure is also still not working for the conservative formulation, so I'll just play around with the grid size and time step and see if I can stabilize my simulation manually. For reference I'm trying to modify the Bickley jet example to simulate some of the other cases in the Poilin, et al. reference.

[pjrusello]

Modified Bickley jet example here

Status message callback is defined at line 125, lifted and modified from the Langmuir turbulence example.
There's an optional TimeStepWizard defined at line 154, but it's commented out right now.

The CFL is always Inf, but the two different methods to access velocity data agree so that's good.

[glwagner]

Can we start with a shorter and more focused piece of code?

[glwagner]

I'm testing this right now:

using Oceananigans
using Oceananigans.Models: ShallowWaterModel
using Printf

grid = RectilinearGrid(size=(4, 4), x=(0, 2π), y=(0, 2π), topology=(Periodic, Bounded, Flat))
model = ShallowWaterModel(; grid, gravitational_acceleration=1)
simulation = Simulation(model, Δt=1.0, stop_iteration=3)
wizard = TimeStepWizard(cfl=0.1)
simulation.callbacks[:wizard] = Callback(wizard)
run!(simulation)

[glwagner]

^ editing the above as I find bugs

[glwagner]

The above runs.

Then we get:

julia> advective_cfl = AdvectiveCFL(1.0)
CFL{Float64, typeof(Oceananigans.Utils.cell_advection_timescale)}(1.0, Oceananigans.Utils.cell_advection_timescale)

julia> advective_cfl(model)
NaN

which is wrong; it should be 0.

[glwagner]

Ok, one issue with the above example is that h=0 everywhere.

Here's a new snippet:

using Oceananigans
using Oceananigans.Models: ShallowWaterModel
using Printf

grid = RectilinearGrid(size=(4, 4), x=(0, 2π), y=(0, 2π), topology=(Periodic, Bounded, Flat))
model = ShallowWaterModel(; grid, gravitational_acceleration=1)
set!(model, h=1)
advective_cfl = AdvectiveCFL(1.0)
@show advective_cfl(model)

This also returns NaN.

[pjrusello]

Yup, I did the same thing at first. model.velocities looks right though, so it seems like something in the CFL calculation isn't referencing the correct field?

[pjrusello]

Could this be because the halo regions contain zeros? I'm assuming a TimeStepWizard would calculate the CFL over the halos as well?

This gives cfl = 0, where the topology is now (Periodic, Periodic, Flat) and the halo regions contains all 1s for h so there's no division by zero.

using Oceananigans
using Oceananigans.Models: ShallowWaterModel
using Printf

grid = RectilinearGrid(size=(4, 4), x=(0, 2π), y=(0, 2π), topology=(Periodic, Periodic, Flat))
model = ShallowWaterModel(; grid, gravitational_acceleration=1)
set!(model, h=1)
advective_cfl = AdvectiveCFL(1.0)
@show advective_cfl(model)

[glwagner]

Well, the CFL calcuation should ignore the halos. The fact that it doesn't is a bug!

Both TimeStepWizard and AdvectiveCFL are using the same function, cell_advection_timescale. The bug lies there; specifically in the version that's used for ShallowWaterModel defined here:

Oceananigans.jl/src/Models/ShallowWaterModels/shallow_water_cell_advection_timescale.jl

Lines 29 to 35 in af21542

	cell_advection_timescale(model::ShallowWaterModel) =
	shallow_water_cell_advection_timescale(
	model.solution[1].data.parent,
	model.solution[2].data.parent,
	model.solution.h.data.parent,
	model.grid, model.formulation
	)

[glwagner]

#2721 might fix it.