RK3 for the hydrostatic free surface model

[simone-silvestri]

paired with a splitting of the barotropic- baroclinic mode described in
https://www.sciencedirect.com/science/article/abs/pii/S0021999122001127?fr=RR-2&ref=pdf_download&rr=8d81dd3bab0dba9a which seems to ensure consistency in the forcing of the barotropic and baroclinic mode.

The algorithm follows:

substep 1

• $u^{*, n+1/2} = u^{n} + \Delta t G(u^{n})$
• $U^{n+1/2} = \text{barotropic-step}\left( u^n, \eta^n, G(u^{n}) \right)$
• $u^{n+1/2} = \text{barotropic-correction}\left(u^{*, n+1/2}, U^{n+1/2} \right)$

substep 2

• $u^{*, n+1/3} = \frac{3}{4}u^{n} + \frac{1}{4}\left(u^{n+1/2} + \Delta t G(u^{n+1/2}) \right)$
• $U^{*, n+1/3} = \text{barotropic-step}\left(u^{n+1/2}, \eta^{n+1/2}, G(u^{n+1/2}) \right)$
• $U^{n+1/3} = \frac{3}{4}U^{n} + \frac{1}{4}U^{*, n+1/3}$
• $u^{n+1/3} = \text{barotropic-correction}\left(u^{*, n+1/3}, U^{n+1/3} \right)$

substep 3

• $u^{*, n+1} = \frac{1}{3}u^{n} + \frac{2}{3}\left(u^{n+1/3} + \Delta t G(u^{n+1/3})\right)$
• $\overline{G}^{n+1} = \frac{1}{6}G(u^{n}) + \frac{1}{6}G(u^{n+1/2}) + \frac{2}{3}G(u^{n+1/3})$
• $U^{n+1} = \text{barotropic-step}\left(u^{n}, \eta^{n}, \overline{G}^{n+1} \right)$
• $u^{n+1} = \text{barotropic-correction}\left(u^{*, n+1}, U^{n+1} \right)$

The proof of concept in this PR is working and allowing a much larger CFL than the AB2 code we use (with the added benefit of reducing the implicit diffusion of the AB2 scheme and being self-starting from time-step n).

This is just a PR thrown together to check this numerical scheme works. I would like to do a couple of improvements on the SplitExplicit code before bringing this PR to a mergeable state. The suggested improvements are described in issue #3873

[simone-silvestri]

For the moment this works only with the SplitExplicitFreeSurface but I think it can be extended to work also with an implicit free surface following the above mentioned algorithm where $\text{barotropic-step}$ and $\text{barotropic-correction}$ are substituted with a solve of the linear system (starting from interpolated velocities) and a pressure correction step with interpolated eta, i.e.:

substep 1:

• $u^{*, n+1/2} = u^{n} + \Delta t G(u^{n})$
• $\eta^{n+1/2} = \text{barotropic-solve}\left(\eta^n, \int_z u^{*, n+1/2} \right)$
• $u^{n+1/2} = \text{pressure-correction}\left(u^{*, n+1/2}, \eta^{n+1/2} \right)$

substep 2

• $u^{*, n+1/3} = \frac{3}{4}u^{n} + \frac{1}{4}\left(u^{n+1/2} + \Delta t G(u^{n+1/2}) \right)$
• $\eta^* = \text{barotropic-solve}\left(\eta^{n+1/2}, \int_z u^{*, n+1/3} \right)$
• $\eta^{n+1/3} = \frac{3}{4}\eta^{n} + \frac{1}{4}\eta^*$
• $u^{n+1/3} = \text{pressure-correction}\left(u^{*, n+1/3}, \eta^{n+1/3} \right)$

substep 3

• $u^{*, n+1} = \frac{1}{3}u^{n} + \frac{2}{3}\left(u^{n+1/3} + \Delta t G(u^{n+1/3})\right)$
• $\int_z u ^{n+1} = \frac{1}{6}\int_z u^{n} + \frac{1}{6}\int_z u^{n+1/2} + \frac{2}{3}\int_z u^{*, n+1}$
• $\eta^{n+1} = \text{barotropic-solve}\left(\eta^{n}, \int_z u ^{n+1} \right)$
• $u^{n+1} = \text{pressure-correction}\left(u^{*, n+1}, \eta^{n+1}\right)$

[simone-silvestri]

It is indeed equivalent in terms of results: when testing tracer advection with a fixed velocity it looks like the profiles are identical for this timestepper in the HydrostaticModel and a NonhydrostaticModel (so probably it is good to have just one rk3 timestepper type)

There is some difference in the nomenclature which I think allows a better understanding of how to split the barotropic and baroclinic modes, because the previous old field are stored and used to restart the substeps instead of having two different tendencies that are averaged at each substep.

[glwagner]

There are some differences between how this algorithm is implemented and how we do the stepping for the NonhydrostaticModel --- do you have any comment about that? We could try implementing this kind of algorithm for the NonhydrostaticModel as well. The implicit free surface solve can be interpreted as a 2D pressure solve, which in the nonhydrostatic context is replaced by a 3D pressure solve.

In particular, it seems that the intermediate predictor velocities $u^\star$ are computed differently and different coefficients are used (1/4, 3/4, 1/3, 2/3) than what we currently use for RK3 in the NonhydrostaticModel:

Oceananigans.jl/src/TimeSteppers/runge_kutta_3.jl

Lines 57 to 62 in d66ed9b

	γ¹ = 8 // 15
	γ² = 5 // 12
	γ³ = 3 // 4
	ζ² = -17 // 60
	ζ³ = -5 // 12

PS ^\star is superior to *!

[glwagner]

It is indeed equivalent in terms of results: when testing tracer advection with a fixed velocity it looks like the profiles are identical for this timestepper in the HydrostaticModel and a NonhydrostaticModel (so probably it is good to have just one rk3 timestepper type)

There is some difference in the nomenclature which I think allows a better understanding of how to split the barotropic and baroclinic modes, because the previous old field are stored and used to restart the substeps instead of having two different tendencies that are averaged at each substep.

Ah sorry after I wrote that (a now deleted comment) I realized there were some subtle differences which I noted above.

[glwagner]

Perhaps worth stating here that I think RK3 is vastly superior to AB2 not only because of the faster time-to-solution but also because it is "self-starting" (does not depend on tendencies from previous time-steps), checkpointing is cheaper and simplified, and the chance of bugs (which cost us nearly 6 months during parameterization development in the past) wherein "restarted" simulations can fail when the prior tendencies have NaN even if we attempt to start with an Euler time-step. So RK3 is both faster and simpler.

[simone-silvestri]

This PR is a little stale, I ll open a version two that is cleaner and close this one.

[simone-silvestri]

closed in favor of #3930