Skip to content

(0.104.0) Reformulate hydrostatic model timestepping#4811

Merged
simone-silvestri merged 756 commits intomainfrom
ss/reformulate-hydrostatic-model-2
Jan 17, 2026
Merged

(0.104.0) Reformulate hydrostatic model timestepping#4811
simone-silvestri merged 756 commits intomainfrom
ss/reformulate-hydrostatic-model-2

Conversation

@simone-silvestri
Copy link
Copy Markdown
Collaborator

@simone-silvestri simone-silvestri commented Sep 23, 2025

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

Comment thread src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl Outdated
Comment thread src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl Outdated
Comment on lines +49 to +50
velocities :: U # Container for velocity fields `u`, `v`, and `w`
transport_velocities :: W # Container for velocity fields used to transport tracers
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what are the units of the transport velocities?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it is also m/s so the same as the velocities

@navidcy navidcy added the numerics 🧮 So things don't blow up and boil the lobsters alive label Sep 25, 2025
Comment thread src/Advection/reconstruction_coefficients.jl Outdated
@simone-silvestri
Copy link
Copy Markdown
Collaborator Author

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

Comment thread ext/OceananigansNCDatasetsExt.jl Outdated
Comment thread src/Advection/Advection.jl Outdated
@glwagner
Copy link
Copy Markdown
Member

glwagner commented Oct 3, 2025

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

Comment thread src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl Outdated
Comment thread src/Models/HydrostaticFreeSurfaceModels/z_star_vertical_spacing.jl Outdated
Comment thread src/TimeSteppers/runge_kutta_3.jl
@glwagner
Copy link
Copy Markdown
Member

glwagner commented Oct 3, 2025

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

Comment thread src/TimeSteppers/split_runge_kutta.jl
@glwagner
Copy link
Copy Markdown
Member

glwagner commented Oct 3, 2025

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
Copy link
Copy Markdown
Collaborator Author

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
Copy link
Copy Markdown
Collaborator Author

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
Copy link
Copy Markdown
Member

glwagner commented Oct 3, 2025

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
Copy link
Copy Markdown
Collaborator Author

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.

@simone-silvestri simone-silvestri added the benchmark performance runs preconfigured benchamarks and spits out timing label Oct 14, 2025
@glwagner
Copy link
Copy Markdown
Member

TODO:

  • Remove timestepping from update_state! In another PR
  • Clean up the time-stepping code
  • 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
  • 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
Copy link
Copy Markdown
Collaborator Author

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

Comment thread src/TimeSteppers/quasi_adams_bashforth_2.jl Outdated
for stages in 2:5
@eval TimeStepper(::Val{Symbol(:SplitRungeKutta, $stages)}, args...; kwargs...) =
SplitRungeKuttaTimeStepper(args...; coefficients=tuple(collect($stages:-1:1)...), kwargs...)
end
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

huh?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

to be able to do

timestepper = :SplitRungeKutta2 

of

timestepper = :SplitRungeKutta4

we discussed above with @ali-ramadhan that we do not want to support schemes with higher number of stages since they are not tested and not be beneficial

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

a github discussion is good, but that fact is not documented! also you could consider simply requiring users to build the scheme manually rather than generating constructors for the symbols (up to you)

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can add something in the docs to say that we can build a timestepper appending the number of stages (:SplitRungeKutta3 for example).

@navidcy navidcy self-requested a review January 17, 2026 19:01
@simone-silvestri simone-silvestri merged commit 5fc355c into main Jan 17, 2026
78 checks passed
@simone-silvestri simone-silvestri deleted the ss/reformulate-hydrostatic-model-2 branch January 17, 2026 23:19
Comment on lines +64 to +65
function assemble_advective_dissipation!(P, grid, ts::RungeKuttaScheme, substep, Fⁿ, Fⁿ⁻¹, Uⁿ, Uⁿ⁻¹, cⁿ⁺¹, cⁿ)
if substep == ts.Nstages
Copy link
Copy Markdown
Collaborator

@giordano giordano Jan 19, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From JETLS:

FieldError: type Oceananigans.TimeSteppers.RungeKutta3TimeStepper has no field Nstages, available fields: γ¹, γ², γ³, ζ², ζ³, Gⁿ, G⁻, implicit_solver

struct RungeKutta3TimeStepper{FT, TG, TI} <: AbstractTimeStepper
γ¹ :: FT
γ² :: FT
γ³ :: FT
ζ² :: FT
ζ³ :: FT
Gⁿ :: TG
G⁻ :: TG
implicit_solver :: TI
end

RungeKuttaScheme is

const RungeKuttaScheme = Union{RungeKutta3TimeStepper, SplitRungeKuttaTimeStepper}
so this method effectively only works with SplitRungeKuttaTimeStepper:
struct SplitRungeKuttaTimeStepper{B, TG, PF, TI} <: AbstractTimeStepper
Nstages :: Int
β :: B
Gⁿ :: TG
Ψ⁻ :: PF # prognostic state at the previous timestep
implicit_solver :: TI
end

Same for assemble_diffusive_dissipation! below.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right we can probable remove that alias, for RK3 this method is not supported. I have started in #4551 but then did not get to complete it

@simone-silvestri simone-silvestri restored the ss/reformulate-hydrostatic-model-2 branch March 9, 2026 09:40
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

breaking change 💔 Concerning a change which breaks the API numerics 🧮 So things don't blow up and boil the lobsters alive

Projects

None yet

Development

Successfully merging this pull request may close these issues.

9 participants