Should the time step wizard take into account of background velocity in calculation of CFL number?

[zhihua-zheng]

Hi,

I'm running a case similar to this example, which has a background velocity from the thermal wind shear. However, I find this setup usually has some numerical instabilities at the beginning. I was wondering if that's related to the time step wizard's adjustment of Δt being too large. From the calculation of cell_advection_timescale here,

Oceananigans.jl/src/Advection/cell_advection_timescale.jl

Lines 13 to 27 in d0b7ec8

	function cell_advection_timescale(grid, velocities)
	u, v, w = velocities
	τ = KernelFunctionOperation{Center, Center, Center}(cell_advection_timescaleᶜᶜᶜ, grid, u, v, w)
	return minimum(τ)
	end
	@inline function cell_advection_timescaleᶜᶜᶜ(i, j, k, grid, u, v, w)
	Δx = Δxᶠᶜᶜ(i, j, k, grid)
	Δy = Δyᶜᶠᶜ(i, j, k, grid)
	Δz = Δzᶜᶜᶠ(i, j, k, grid)
	inverse_timescale = @inbounds abs(u[i, j, k]) / Δx + abs(v[i, j, k]) / Δy + abs(w[i, j, k]) / Δz
	return 1 / inverse_timescale
	end

it seems like it does not use the background velocity.

For this use case, should the time step wizard include the background velocity to calculate a new Δt?

[simone-silvestri]

You are right! Indeed, the CFL calculation should account for additional velocities. Seems like we have to update the TimeStepWizard to account for additional velocities.

It should be possible by changing the model-specific cell_advection_timescale functions

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

Line 61 in d0b7ec8

cell_advection_timescale(model::HydrostaticFreeSurfaceModel) = cell_advection_timescale(model.grid, model.velocities)

and

Oceananigans.jl/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl

Line 51 in d0b7ec8

cell_advection_timescale(model::NonhydrostaticModel) = cell_advection_timescale(model.grid, model.velocities)

to pass total velocities as done here

Oceananigans.jl/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl

Lines 65 to 67 in d0b7ec8

	total_velocities = (u = SumOfArrays{2}(velocities.u, background_fields.velocities.u),
	v = SumOfArrays{2}(velocities.v, background_fields.velocities.v),
	w = SumOfArrays{2}(velocities.w, background_fields.velocities.w))

[glwagner]

For sure and just to clarify, this is not an issue with TimeStepWizard, but rather with the definition of cell_advection_timescale(model), which is a function of each model individually.

I believe background velocities are only supported by NonhydrostaticModel, right?

[zhihua-zheng]

I believe background velocities are only supported by NonhydrostaticModel, right?

Yes.

[glwagner]

@zhihua-zheng this is a bona fide issue, do you mind if I convert to an issue? Would you like to try to fix it?

[zhihua-zheng]

No problem with converting to an issue. I have tried some modifications according to @simone-silvestri 's suggestion, and it seems to work. No experience with fixing issues before, but I think it would be interesting to start doing it.

[glwagner]

All you gotta do is submit your new code to be merged into Oceananigans using a pull request!

[jagoosw]

I think all that needs to be change is this one line:

Oceananigans.jl/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl

Line 51 in cc51ff0

cell_advection_timescale(model::NonhydrostaticModel) = cell_advection_timescale(model.grid, model.velocities)

which can be changed to:

cell_advection_timescale(model::NonhydrostaticModel) = cell_advection_timescale(model.grid, total_velocities(model)) 

since we defined total_velocities here:

Oceananigans.jl/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl

Lines 224 to 227 in cc51ff0

	@inline total_velocities(m::NonhydrostaticModel) =
	(u = SumOfArrays{2}(m.velocities.u, m.background_fields.velocities.u),
	v = SumOfArrays{2}(m.velocities.v, m.background_fields.velocities.v),
	w = SumOfArrays{2}(m.velocities.w, m.background_fields.velocities.w))

(plus using Oceananigans.Models: total_velocities at the start of the file).

[jagoosw]

I suppose arguably we need to include biogeochemical_drift_velocity but then we would have to compute the advection timescale for every tracer because we couldn't know how the drift velocity would effect the timescale

Oceananigans.jl/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl

Lines 250 to 254 in cc51ff0

	biogeochemical_velocities = biogeochemical_drift_velocity(biogeochemistry, val_tracer_name)
	total_velocities = (u = SumOfArrays{3}(velocities.u, background_fields.velocities.u, biogeochemical_velocities.u),
	v = SumOfArrays{3}(velocities.v, background_fields.velocities.v, biogeochemical_velocities.v),
	w = SumOfArrays{3}(velocities.w, background_fields.velocities.w, biogeochemical_velocities.w))

[glwagner]

(plus using Oceananigans.Models: total_velocities at the start of the file).

You won't need that because cell_advection_timescale is defined in the same module as total_velocities

[glwagner]

Another possibility is

function cell_advection_timescale(model::NonhydrostaticModel)
    grid = model.grid
    velocities = total_velocities(model)
    return cell_advection_timescale(grid, velocities)
end

[zhihua-zheng]

I think these are all good solutions! Just make sure I understand the source code here: Both your and @jagoosw's modification of cell_advection_timescale in NonhydrostaticModels.jl is to redefine the method of cell_advection_timescale in NonhydrostaticModel. This method is being used in

Oceananigans.jl/src/Simulations/time_step_wizard.jl

Line 102 in cc51ff0

advective_Δt = wizard.cfl * wizard.cell_advection_timescale(model)

with model as the input, to compute a new time step. Is this interpretation correct?

In addition to @jagoosw's comment about biogeochemical_drift_velocity: does the velocity in AdvectiveForcing also play a role here?

[zhihua-zheng]

You won't need that because cell_advection_timescale is defined in the same module as total_velocities

Because total_velocities that defined in

Oceananigans.jl/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl

Lines 223 to 227 in cc51ff0

	# return the total advective velocities
	@inline total_velocities(m::NonhydrostaticModel) =
	(u = SumOfArrays{2}(m.velocities.u, m.background_fields.velocities.u),
	v = SumOfArrays{2}(m.velocities.v, m.background_fields.velocities.v),
	w = SumOfArrays{2}(m.velocities.w, m.background_fields.velocities.w))

is included by this line below? Sorry, I'm new to Julia and the code structure seems a bit convoluted.

Oceananigans.jl/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl

Line 42 in cc51ff0

include("nonhydrostatic_model.jl")

[glwagner]

The module I was referring to is NonhydrostaticModels. The top-level file NonhydrostaticModels.jl includes all the code for the module NonhydrostaticModels. The directory containing all the code for the module NonhydrostaticModels is called NonhydrostaticModels:

We are definitely open to concrete suggestions or PRs that improve the code structure and make it less convoluted, which is highly undesirable.