TimeStepWizard

[parfenyev]

Hello everyone!

I am trying to do a simulation with an adaptive time step, below is my code:

using Oceananigans, Statistics, Printf

grid = RectilinearGrid(size=(256, 256), extent=(2π, 2π),
                              topology=(Periodic, Periodic, Flat))

# parameters
G = 0.05
k = 9.0
u_forcing(x, y, z, t) = G*cos(k*y)
v_forcing(x, y, z, t) = -G*cos(k*x)
nu = 1e-5

model = NonhydrostaticModel(timestepper = :RungeKutta3,
                              advection = UpwindBiasedFifthOrder(),
                                   grid = grid,
                                forcing = (u=u_forcing, v=v_forcing),
                               buoyancy = nothing,
                                tracers = nothing,
                                closure = ScalarDiffusivity(ν=nu))

set!(model, u=0, v=0)
simulation = Simulation(model, Δt=1e-4, stop_time=40)

wizard = TimeStepWizard(cfl=0.5, max_change=1e-2, max_Δt=1e-1)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(100))

progress_message(sim) = @printf("Iteration: %07d, time: %s, Δt: %s, cfl = %.4e, wall time: %s\n",
                                iteration(sim),
                                prettytime(sim),
                                prettytime(sim.Δt),
                                AdvectiveCFL(sim.Δt)(sim.model),
                                prettytime(sim.run_wall_time))
simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100))

u, v, w = model.velocities
ω = ∂x(v) - ∂y(u)

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; u, v, ω),
                                                      schedule = TimeInterval(0.1),
                                                      filename = "256_example.jld2",
                                                      overwrite_existing = true)

run!(simulation)

The initial cfl is small because the fluid is initially at rest and the flow is driven by an external force. I would expect the time step to increase, but it gradually decreases:

[ Info: Initializing simulation...
Iteration: 0000000, time: 0 seconds, Δt: 50 μs, cfl = 0.0000e+00, wall time: 0 seconds
[ Info:     ... simulation initialization complete (63.188 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (2.852 seconds).
Iteration: 0000100, time: 5.000 ms, Δt: 25 μs, cfl = 2.5463e-07, wall time: 18.488 seconds
Iteration: 0000200, time: 7.500 ms, Δt: 12.500 μs, cfl = 1.9097e-07, wall time: 33.809 seconds
Iteration: 0000300, time: 8.750 ms, Δt: 6.250 μs, cfl = 1.1140e-07, wall time: 49.102 seconds
Iteration: 0000400, time: 9.375 ms, Δt: 3.125 μs, cfl = 5.9678e-08, wall time: 1.073 minutes
Iteration: 0000500, time: 9.687 ms, Δt: 1.562 μs, cfl = 3.0834e-08, wall time: 1.329 minutes
Iteration: 0000600, time: 9.844 ms, Δt: 781.250 ns, cfl = 1.5666e-08, wall time: 1.585 minutes

What am I doing wrong?

[glwagner]

I think max_change has to be great than 1! Basically you want max_change = 1 + 1e-2.

This is interesting though --- I kind of like the idea that "max change" is a small positive number (rather than being a number slightly greater than 1). We could change the API so your code works. It might be more intuitive.

@navidcy @tomchor what do you think about that?

Otherwise we should definitely throw an error if max_change < 1.

[parfenyev]

Thank you!

[navidcy]

tricky one :)

Perhaps max_change kwarg name should be something more descriptive?

At the moment, if we have max_change = 1.1 it implies that a timestep of Δt=10 could change within range 9.09 < Δt 11.0, right?

[glwagner]

Actually I think it means that dt increase to 11.0 at most:

Oceananigans.jl/src/Simulations/time_step_wizard.jl

Line 86 in 64ccdc6

new_Δt = min(wizard.max_change * old_Δt, new_Δt)

There's the corresponding min_change that caps the decrease of the time-step.

I think max_change is the wrong word. I think @parfenyev hits on the right meaning of this term. The way @parfenyev uses it, "maximum change" is a positive percentage (ie, dt can change by 10%).

max_change = 1.1 implying 10% possible change seems weird with this new perspective 💡

Maybe max_change being a percentage that caps how much dt can change in both the positive and negative direction is what we want???

[tomchor]

Maybe max_change being a percentage that caps how much dt can change in both the positive and negative direction is what we want???

I don't know, for me having separate args for the lower bound and upper bound has been helpful in the past. In simulations that for some reason are more unstable (especially in the beginning of some LES with transition to turbulence) I've had to set min_change=0.2, max_change=1.05, for example, and it solved the problem.

But I do agree that the names might not intuitively convey the behavior. I can't, at the moment, think of a better way to do things though.

[francispoulin]

I haven't used this before but will say from what I have just read, more descriptive names woudl be very helpful to the user. max_relative_change is probably not the best but maybe better?

[glwagner]

The more we discuss the more I'm convinced the names don't really make sense: "min change" ?!?

We can still have two arguments to limit how much the time-step increases and decreases separately. Sorry for the confusion --- that was a point that's additional to what we are discussing.

So, maybe max_decrease and max_increase --- both percentages > 0, and additionally max_decrease is less than 1 (because decreasing by 100% decreases the time-step to 0).

[francispoulin]

max_decrease and max_increase are much better than the current names.