Implementing AB2 for variable time-step size

[simone-silvestri]

The Adams-Bashforth time stepper relies on an approximate time-integral to calculate the tendencies.
This time integral, when using a constant time step, is simply approximated as

$$G^{n+1} = \left(\frac{3}{2} + \chi \right) G^n - \left( \frac{1}{2} +\chi \right) G^{n-1}$$

The $\chi$ term should be a small deviation from the time integral, added to reduce the noise generated by non-linear terms.
I wanted to open a bit of a discussion about the details of

1. the $\chi$ value
2. the time integral when we use variable time stepping that is a feature we use quite often

By default, we use $\chi = 0.1$, which is quite large (20% of the smaller coefficient). This will introduce quite a bit of implicit dissipation in the model. Generally, there is a tradeoff between the stability obtained by a larger $\chi$ and the dissipation introduced by deviating from the AB2 behavior. I wondered if 0.1 is too high, especially when using diffusive methods.
Did anyone experiment with lower $\chi$, and if yes, with what results?

Would it make sense to exclude $\chi$ from diffusive terms to limit the implicit dissipation of the time-stepping scheme?

Regarding the time integral for variable time-stepping, the correct form of the AB2 with variable step would have to include the time step at $n$ and the time step at $n-1$ to be correct. This might not make a huge impact, but if we want to save time averaged tendencies and the time step changes size every ten iterations or so, the error will compound and the time average will probably be off.
Would it make sense to implement time-step dependent coefficients for the AB2 scheme?

The problem is that if the time step changes size, the tendency terms at $G^{n-1}$ do not cancel, which is what happens with constant time stepping.

$$c^{n+1} = c^{n} + \Delta t (1.5 G^{n} - 0.5 G^{n-1})$$ $$c^{n+2} = c^{n+1} + \Delta t (1.5 G^{n+1} - 0.5 G^{n})$$

In the above, between $c^{n+1}$ and $c^{n+2}$ we have added only one $G^{n}$ because the extra tendency we have added in $c^{n+1}$ is removed in the successive timestep.
With variable time steps, this is not the case! Suppose the time-step changes between $n+1$ and $n+2$

$$c^{n+1} += \Delta t^n (1.5 G^{n} - 0.5 G^{n-1})$$ $$c^{n+2} += \Delta t^{n+1} (1.5 G^{n+1} - 0.5 G^{n})$$

we remain with a spurious $G^n (1.5\Delta t^{n} - 0.5 \Delta t^{n-1}) - G^n$.

Given this reference, the correct formualtion might be

$$G^{n+1} = \frac{1}{2} \left( \left( 2 + \frac{\Delta t^n}{\Delta t^{n - 1}} \right) G^n - \frac{\Delta t^n}{\Delta t^{n - 1}} G^{n-1} \right)$$

However, also in this form, successive tendencies do not cancel out.

[glwagner]

Should this be a discussion or an issue? As an issue we would like to have some source code changes that we can make in order to resolve the issue.

[simone-silvestri]

It would be nice to convert it into a discussion, but I think we want to correct the AB2 for variable time stepping quite soon, because it might be a large source of error.

[glwagner]

Should we have a separate issue for that then, because we want to continue much of this even if that is fixed right?

When the time-step changes we take an euler step:

Oceananigans.jl/src/TimeSteppers/quasi_adams_bashforth_2.jl

Line 96 in b3be950

euler = euler || (Δt != model.clock.last_Δt)

[simone-silvestri]

Oh wow, I hadn't even noticed that! That might cause a lot of stability problems if it happens once every 10 or 20 steps!

[glwagner]

I changed the title if you want to focus on variable time-step but I think you should open another discussion or issue for the other stuff about the size of chi

[glwagner]

the time integral when we use variable time stepping that is a feature we use quite often

What do you mean by this?

[glwagner]

Regarding the time integral for variable time-stepping, the correct form of the AB2 with variable step would have to include the time step at n and the time step at n-1 to be correct.

It seems like a good idea to use this.

Although, I think we basically always recommend using RK3 for NonhydrostaticModel.

Would another path forward be to also implement RK3 for HydrostaticFreeSurfaceModel, and basically deprecate any recommendation to use AB2 at all? It's just so much more convenient to use RK3 and also seems to have better numerical properties, not least in light of this issue...

[simone-silvestri]

the time integral when we use variable time stepping that is a feature we use quite often

What do you mean by this?

Actually, I should probably say that I use the wizard quite often to change the time step. In my opinion, AB2 is a good compromise between accuracy, stability, and performance. RK3 is better only when you can achieve a CFL 3 times larger. We should fix the variable AB2 time stepper or discourage the use of frequent updates of the time step when using AB2. The first option is probably better in my opinion. I ll look into fixing AB2 and what it entails

[glwagner]

I find that RK3 usually allows a time-step about 5x larger, up to CFL of 0.7 vs 0.1-0.2 for AB2. I thought this was the typical experience. Why don't we answer this once and for all, there are not two answers to this question.

[glwagner]

discourage the use of frequent updates of the time step when using AB2

This is just buyer beware, as with anything. There are no examples currently that update the time-step at all while using AB2...

[ali-ramadhan]

Out of curiousity, how much work is it to implement RK3 for HydrostaticFreeSurfaceModel?

Seems like fixing AB2 for variable time steps isn't too difficult if it's just changing how $G^{n+1}$ is computed though.

Why don't we answer this once and for all, there are not two answers to this question.

I thought this was answered quite clearly, albeit for a very idealized case, here: #945 (comment) Is AB2 still technically only first-order accurate?

I know it's quite difficult to extrapolate from the convergence test to global simulations, but my experience has been that RK3 beats out AB2 for non-hydrostatic simulations. That said, it's probably still good to have AB2 especially if it can support variable time steps.

[glwagner]

I thought this was answered quite clearly, albeit for a very idealized case, here: #945 (comment) Is AB2 still technically only first-order accurate?

True, that PR does demonstrate convergence with time step! However we seem to rarely concern ourselves with time-stepping errors (apart from catastrophic instability), which I think is justified since spatial errors dominate in our simulations that are typically marginally-resolved in space by design.

The stickier and harder question regards time to solution, which is what we really care about. It's stickier because it is fundamentally heuristic. For example since RK3 is a 3-stage time-stepper and AB2 is 1-stage, we need to be able to take time-steps that are 3 times longer with RK3. Theoretically, AB2 is stable up to CFL=0.5 while RK3 is stable up to CFL=sqrt(3)=1.7 --- just barely achieving the "3x" needed for RK3 to be superior. But that's just in theory. In practice we find that solutions become unstable at lower CFL, and caution warrants a bit of buffer. It seems like we usually use CFL=0.1-0.2 for AB2 but CFL=0.5-1 for RK3. So at the lower end for AB2, RK3 is indeed providing a decent speed up for time to solution. And I agree with you @ali-ramadhan that I have found in practice that I prefer RK3 for LES.

There is more to the story --- RK3 also has superior stability properties for diffusion. Where there is some interesting room for further research is also to understand how much more stable RK3 might make a turbulence closure scheme like CATKE or k-epsilon. We do find that these can impose stricter time-step restrictions than advection, especially on grids with very fine vertical resolution near the surface.

RK3 also has the crucial advantage that it is self-starting. This means we can change the time-step at will (without resolving this issue I guess) and we don't need tendencies to checkpoint. Also of course RK3 is actually 3rd order whereas QAB2 is first order.

Note that the advantages of RK3 have historically not been useful for ocean simulations because at coarse (non-eddy-resolving) resolutions the momentum time-step is limited by Coriolis, not by CFL / advection. In that case a one-stage time-stepper might be preferred.

Out of curiousity, how much work is it to implement RK3 for HydrostaticFreeSurfaceModel?

It may be a research project to adapt the split-explicit free surface. It's little unclear because naively, we may be able to advance the split-explicit free surface each stage and therefore use the same algorithm we use for AB2. But there has been some work to implement an algorithm that resembles Le and Moin 1991 where the free surface is advanced to the end of the time-step, while baroclinic variables are advanced in stages. That would yield even more speed up.

@simone-silvestri has thought about this a bit and might have more to add...

[ali-ramadhan]

Thanks for all the context!

Good points about the time step being more restricted by Coriolis for coarse global simulations and error being dominated by spatial discretization.

I'd also be very interested how AB2 compares against RK3 for time-to-solution. I guess this could be readily tested with the non-hydrostatic model simulating turbulence. I'd be curious if RK3 is always faster, or if it's case-dependent. Although it is risky to be on the edge of stability CFL-wise.

It may be a research project to adapt the split-explicit free surface.

Ah I didn't realize that RK3 was not really used for global ocean models, especially with a split-explicit free surface :(

However, also in this form, successive tendencies do not cancel out.

Is this neccessary or is this why Quasi AB2 is technically only first-order accurate? I guess right now with Euler steps the tendency terms do cancel out when an AB2 time step is taken.

[glwagner]

Ah I didn't realize that RK3 was not really used for global ocean models, especially with a split-explicit free surface :(

As far as I know it has never been used. However, we are aware of some research efforts to implement RK3 in NEMO. Here is a talk about it:

https://usclivar.org/abstract-ocean-model-2024/gain-efficiency-new-time-scheme-nemo-runge-kutta-3rd-order

I also found this:

https://forge.nemo-ocean.eu/nemo/nemo/-/merge_requests/68