Reformulate hydrostatic model timestepping

[simone-silvestri]

This is a work-in-progress PR that aims to reformulate the timestepping scheme of the HydrostaticFreeSurfaceModel to make sure that three numerical properties are satisfied:

• tracers are globally conserved (the integral of a tracer field remains constant in the absence of forcings)
• tracers are locally conserved (an initially constant tracer remains constant in the absence of forcings)
• The free surface is evolved consistently with the grid thickness.

This is easier to achieve with an explicit discretization of the free surface (ExplicitFreeSurface and SplitExplicitFreeSurface) rather than an implicit free surface discretization and with RK3 rather than AB2. However, this PR reformulates much of the timestepping algorithm, so it will be a lengthy WIP that will need to be thoroughly validated for all combinations of barotropic and baroclinic timestepping schemes.

More changes are performed to correctly abstract the timestepping schemes, which, at the moment, are tailored specifically to the NonhydrostaticModel (for example, the timesteppers include a compute_pressure_correction! or a correct_velocities function, which are meaningful only for a non-hydrostatic Navier Stokes solver)

Edit: another objective I have is to remove any timestepping from the update_state! function. We should be able to call update_state! multiple times in a row without problems and timestepping inside update_state! negates this possibility.

TODO:

• [ ] Remove timestepping from update_state! In another PR
• [x] Clean up the time-stepping code
• [x] Introduce transport weights for transport averaging in the SplitExplicitFreeSurface
• [x] Implement a conservative, low-storage Runge Kutta method (of arbitrary order) using the Explicit free surface
• [x] Implement a conservative, low-storage Runge Kutta method (of arbitrary order) coupled with the SplitExplicitFreeSurface
• [x] Implement a conservative, low-storage Runge Kutta method (of arbitrary order) using the (PCG) ImplicitFreeSurface
• [x] Ensure that the conservation is satisfied on DistributedGrids
• [ ] Implement a semi - conservative method for AB2 In another PR
• [ ] Address precision errors that tend to destroy conservation for constant fields In another PR

[simone-silvestri]

I am at the point here that the RK formulation is conservative both locally and globally. However, there seems to be precision errors accumulating when the initial field is constant, which break conservation. I have gotten some suggestions from @milankl as to how to tackle this issue. I will post a MWE later on when I have addressed some other points (conservation on distributed architectures and on AB2).

[glwagner]

why are there so many changes to WENO advection, etc?

[glwagner]

It's really hard for me to understand how this is an improvement. I think some explanation is needed (maybe docs or a docstring...)

[glwagner]

It wouldn't be a completely crazy idea to start a new docs section on "Time steppers" which outlines how each time stepper works. Esp because there appears to be a kind of interface for using different methods, eg AB2 and RK3.

[simone-silvestri]

why are there so many changes to WENO advection, etc?

I have merged in #4434 because I was expecting to merge it before this one. Maybe that was a mistake

[simone-silvestri]

It wouldn't be a completely crazy idea to start a new docs section on "Time steppers" which outlines how each time stepper works. Esp because there appears to be a kind of interface for using different methods, eg AB2 and RK3.

Yeah that would be nice! After this PR we can also use RK for sea ice, for example, by extending rk_substep!

[glwagner]

why are there so many changes to WENO advection, etc?

I have merged in #4434 because I was expecting to merge it before this one. Maybe that was a mistake

Do you want us to also review the WENO changes then or ignore them?

[simone-silvestri]

Yeah, also some of the code here probably still needs to change as I am fixing the implementation. My master plan is to get everything working properly with the numerics correct for all cases (all free surfaces, distributed, immersed boundaries etc...) so that I have a clear idea of all the changes to be done. Then, if I can keep it simple, go for one PR, while if the mole of the changes is too much for one PR, split the work in more PRs. The advection PR is ready to review, though. I have added back the high-order advection schemes.

[glwagner]

TODO:

• [ ] Remove timestepping from update_state! In another PR
• [x] Clean up the time-stepping code
• [x] Implement a conservative, low-storage Runge Kutta method (of arbitrary order) using the Explicit free surface
• [x] Implement a conservative, low-storage Runge Kutta method (of arbitrary order) coupled with the SplitExplicitFreeSurface
• [x] Implement a conservative, low-storage Runge Kutta method (of arbitrary order) using the (PCG) ImplicitFreeSurface
• [x] Ensure that the conservation is satisfied on DistributedGrids
• [x] Implement a semi - conservative method for AB2
• [ ] Address precision errors that tend to destroy conservation for constant fields In another PR

There are 115 files changed on this PR. Does it make sense to split up some of these objectives? It seems like the first objective ("clean up time-stepping") doesn't need to be part of a PR that implements new methods. Also, why wouldn't it make sense to implement each of the different time-stepping methods in a new PR? This would make it easier and faster to get the PR merged and finished.

[simone-silvestri]

Yeah, I will definitely split the PR; however, most of the changes here are already in #4434. So, until that PR is merged, the scale of this one is not really reflected on the file changes

[francispoulin]

Looks like all the boxes left to be checked are checked. Well done @simone-silvestri ! Too bad some of the checks failed.

[milankl]

I haven't worked again on the low precision conservation issues. I believe this is possible and can try to find a solution earlier if urgent, but could possibly also go into another PR?

[simone-silvestri]

I haven't worked again on the low precision conservation issues. I believe this is possible and can try to find a solution earlier if urgent, but could possibly also go into another PR?

No problem, yeah, I think the important thing to figure out was that the conservation breaking was indeed a precision issue, which at least confirms the correctness of the numerics. We can look at a correction in a new PR!