Reformulate hydrostatic model timestepping #4811

simone-silvestri · 2025-09-23T16:09:22Z

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
Clean up the time-stepping code
Introduce transport weights for transport averaging in the SplitExplicitFreeSurface
Implement a conservative, low-storage Runge Kutta method (of arbitrary order) using the Explicit free surface
Implement a conservative, low-storage Runge Kutta method (of arbitrary order) coupled with the SplitExplicitFreeSurface
Implement a conservative, low-storage Runge Kutta method (of arbitrary order) using the (PCG) ImplicitFreeSurface
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

src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl

glwagner · 2025-09-23T16:15:41Z

src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl

+    velocities :: U            # Container for velocity fields `u`, `v`, and `w`
+    transport_velocities :: W  # Container for velocity fields used to transport tracers


the difference between "normal velocities" and "transport velocities" is unclear

Hmm, how should we call them? They are the velocities that transport tracers around. transport is not good because it implies a multiplication times the area...

src/Advection/reconstruction_coefficients.jl

simone-silvestri · 2025-10-02T15:55:10Z

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 · 2025-10-03T12:06:54Z

ext/OceananigansNCDatasetsExt.jl

+    # Synchronize model state if needed
+    for field in fields(model)
+        synchronize_communication!(field)
+    end


how related to this PR?

I am trying to figure out how to distribute operations in time. I don't think this is necessary anymore.

glwagner · 2025-10-03T12:07:24Z

src/Advection/Advection.jl

-# `Centered(order=12`) and `UpwindBiased(order=11)`. The list can be extended in order to
+# `Centered(order=10`) and `UpwindBiased(order=9)`. The list can be extended in order to
 # compile schemes with higher orders; for example `advection_buffers = [1, 2, 3, 4, 5, 6, 8]`
 # will compile schemes for `advection_buffer=8` and thus `Centered(order=16)` and `UpwindBiased(order=15)`.
 # Note that it is not possible to compile schemes for `advection_buffer = 41` or higher.
-const advection_buffers = [1, 2, 3, 4, 5, 6]
+const advection_buffers = [1, 2, 3, 4, 5]


we should compile past the recommended order. If order=9 is recommended we need order=11 to verify this.

Okay, I guess I'll take this suggestion for #4434, which, apart from this change, should be ready.

glwagner · 2025-10-03T12:08:41Z

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

src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl

glwagner · 2025-10-03T12:10:43Z

src/Models/HydrostaticFreeSurfaceModels/z_star_vertical_spacing.jl

-end
-
-@inline function update_grid_scaling!(σᶜᶜⁿ, σᶠᶜⁿ, σᶜᶠⁿ, σᶠᶠⁿ, σᶜᶜ⁻, i, j, grid, ηⁿ)
+@inline function update_grid_scaling!(z, i, j, grid)


Suggested change

@inline function update_grid_scaling!(z, i, j, grid)

@inline function update_grid_scaling!(z_coordinate, i, j, grid)

glwagner · 2025-10-03T12:12:15Z

src/TimeSteppers/runge_kutta_3.jl

    # First stage
    #

-    compute_flux_bc_tendencies!(model)


how are they computed?

I tried putting this back inside the compute_tendencies! given that these cannot be computed in update_state! anymore. I will verify everything is consistent.

glwagner · 2025-10-03T12:13:40Z

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

src/TimeSteppers/split_runge_kutta.jl

glwagner · 2025-10-03T12:18:48Z

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 · 2025-10-03T13:49:13Z

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 · 2025-10-03T13:59:29Z

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!

…nto ss/optimize-weno

…/Oceananigans.jl into ss/reformulate-hydrostatic-model-2

src/TimeSteppers/runge_kutta_3.jl

…tatic-model-2

…/Oceananigans.jl into ss/reformulate-hydrostatic-model-2

glwagner reviewed Sep 23, 2025

View reviewed changes

src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl Outdated Show resolved Hide resolved

glwagner reviewed Sep 23, 2025

View reviewed changes

src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl Outdated Show resolved Hide resolved

glwagner reviewed Sep 23, 2025

View reviewed changes

navidcy added the numerics 🧮 So things don't blow up and boil the lobsters alive label Sep 25, 2025

This was referenced Sep 30, 2025

Create 1/8 degree prototype script CliMA/ClimaOcean.jl#629

Draft

Improve documentation of AbstractModel interface #4797

Open

glwagner reviewed Oct 2, 2025

View reviewed changes

src/Advection/reconstruction_coefficients.jl Outdated Show resolved Hide resolved

glwagner reviewed Oct 3, 2025

View reviewed changes

src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl Outdated Show resolved Hide resolved

glwagner reviewed Oct 3, 2025

View reviewed changes

src/TimeSteppers/split_runge_kutta.jl Show resolved Hide resolved

adding coefficients

c854bf4

simone-silvestri mentioned this pull request Oct 3, 2025

A rework for the advection module that improves WENO performance #4434

Open

simone-silvestri added 9 commits October 3, 2025 16:21

Merge branch 'main' into ss/optimize-weno

91cc41a

adding in the smoothness

579b409

Merge branch 'ss/optimize-weno' of github.com:CliMA/Oceananigans.jl i…

1ce4e4a

…nto ss/optimize-weno

complete 11th order implementation

444f7b6

Merge branch 'main' into ss/optimize-weno

97f8a70

small test

1062702

bugfix

312f620

Merge branch 'ss/optimize-weno' of github.com:CliMA/Oceananigans.jl i…

5b553e3

…nto ss/optimize-weno

another bugfix

31c5d7a

simone-silvestri added 11 commits December 1, 2025 17:48

Merge branch 'main' into ss/reformulate-hydrostatic-model-2

7b08ad6

bugfix

ed48f71

add the xflat grid and yflat grid

536a884

remove the shows

deb91e2

Merge branch 'main' into ss/reformulate-hydrostatic-model-2

f579e62

bugfix

c9d0694

Merge branch 'ss/reformulate-hydrostatic-model-2' of github.com:CliMA…

4357daa

…/Oceananigans.jl into ss/reformulate-hydrostatic-model-2

Merge branch 'main' into ss/reformulate-hydrostatic-model-2

2b20f17

some cleanup

9e331a4

Merge branch 'ss/reformulate-hydrostatic-model-2' of github.com:CliMA…

94ff01f

…/Oceananigans.jl into ss/reformulate-hydrostatic-model-2

start fixing up tests

a2c58f7

giordano reviewed Dec 4, 2025

View reviewed changes

src/TimeSteppers/runge_kutta_3.jl Outdated Show resolved Hide resolved

simone-silvestri added 18 commits December 4, 2025 12:16

only distributed zstar tests remaining

dab4d1f

bugfix the nonhydrostatic model

6430b6d

some updates

a075728

this was not the cause

3ac2181

add the free surface

6669f2d

adding more distributed tests

c981828

bugfix

c8cac57

radius should not be 1 if z goes to -1!

4662431

will this change stuff?

0eb35af

just a check

abc7b4f

not sure why this would pass locally and not on buildkite

f0133cc

Merge branch 'main' into ss/reformulate-hydrostatic-model-2

2082c9e

Merge remote-tracking branch 'origin/main' into ss/reformulate-hydros…

751cfd6

…tatic-model-2

correct imports

49cd565

other fixes

bce48fc

Merge branch 'main' into ss/reformulate-hydrostatic-model-2

6862043

cche current fields

16ffe0e

Merge branch 'ss/reformulate-hydrostatic-model-2' of github.com:CliMA…

1fc8936

…/Oceananigans.jl into ss/reformulate-hydrostatic-model-2

		velocities :: U # Container for velocity fields `u`, `v`, and `w`
		transport_velocities :: W # Container for velocity fields used to transport tracers

	@inline function update_grid_scaling!(z, i, j, grid)
	@inline function update_grid_scaling!(z_coordinate, i, j, grid)

Reformulate hydrostatic model timestepping #4811

Are you sure you want to change the base?

Reformulate hydrostatic model timestepping #4811

Conversation

simone-silvestri commented Sep 23, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

simone-silvestri commented Oct 2, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

glwagner commented Oct 3, 2025

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

glwagner commented Oct 3, 2025

Uh oh!

Uh oh!

glwagner commented Oct 3, 2025

Uh oh!

simone-silvestri commented Oct 3, 2025

Uh oh!

simone-silvestri commented Oct 3, 2025

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

7 participants

simone-silvestri commented Sep 23, 2025 •

edited

Loading