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...