Skip to content

Open boundary conditions using SplitExplicitFreeSurface#5351

Open
simone-silvestri wants to merge 70 commits intomainfrom
ss/open-boundary-conditions
Open

Open boundary conditions using SplitExplicitFreeSurface#5351
simone-silvestri wants to merge 70 commits intomainfrom
ss/open-boundary-conditions

Conversation

@simone-silvestri
Copy link
Copy Markdown
Collaborator

@simone-silvestri simone-silvestri commented Mar 2, 2026

This PR implements open boundary conditions for the hydrostatic free surface model compatible with the split explicit free surface. The boundary conditions follow what discussed in #5229, meaning Flather boundary conditions for the barotropic component and "Orlanski" (Radiation) boundary conditions for the three-dimensional evolution.

I need to clean up a bit the PR (for example Flather now stores the external values inside its own struct, but it would be probably better to store it in the condition of the boundary_condition struct.

Since the flow of the algorithm has to change, because the barotropic correction needs to be done after filling up the boundary conditions for the 3d velocities, such that the barotropic components of the boundaries is consistent with the 2d evolution, some refactor is necessary and I would like to open a PR specifically for this (this is actually a positive because it means we can compute the entire w velocity in the update_state in distributed computations).

There is an example that showcases the use of these boundary conditions which we can remove later on.
A part from all the changes that need to be done to clean up the code, this feature should be already fully functional and testable, so if you want to test it @vtamsitt you can give it a try

cc @tomchor

  • Include tests
  • Clean up Flather BC

@vtamsitt
Copy link
Copy Markdown

vtamsitt commented Mar 2, 2026

This PR implements open boundary conditions for the hydrostatic free surface model compatible with the split explicit free surface. The boundary conditions follow what discussed in #5229, meaning Flather boundary conditions for the barotropic component and "Orlanski" (Radiation) boundary conditions for the three-dimensional evolution.

I need to clean up a bit the PR (for example Flather now stores the external values inside its own struct, but it would be probably better to store it in the condition of the boundary_condition struct.

Since the flow of the algorithm has to change, because the barotropic correction needs to be done after filling up the boundary conditions for the 3d velocities, such that the barotropic components of the boundaries is consistent with the 2d evolution, some refactor is necessary and I would like to open a PR specifically for this (this is actually a positive because it means we can compute the entire w velocity in the update_state in distributed computations).

There is an example that showcases the use of these boundary conditions which we can remove later on. A part from all the changes that need to be done to clean up the code, this feature should be already fully functional and testable, so if you want to test it @vtamsitt you can give it a try

cc @tomchor

  • Include tests
  • Clean up Flather BC

@simone-silvestri that was fast! I will see if I can get to this at the end of the week, or @rafmudaf may be able to test first with our current Med config.

BoundaryCondition,
FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition, OpenBoundaryCondition,
PerturbationAdvection,
PerturbationAdvection, Flather, Radiation,
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's the difference between these three? Isn't Radiation a subset of PerturbationAdvection?

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 think it's the opposite. PerturbationAdvection requires specifying a phase velocity at the boundary while radiation computes it (from what I understood of perturbation advection)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I think you could make the phase velocity computation the condition and then use perturbation advection?

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.

Heh, I mean in either case --- whether Radiation is a subset of PerturbationAdvection, or whether PerturbationAdvection is a subset of Radiation --- in either of those cases, we should not introduce a new type but instead reuse the same type to express both schemes.

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.

we should not introduce a new type but instead reuse the same type to express both schemes.

I'm not sure. Radiation schemes can differ widely in their internals. Having different scheme types for each one may be useful just for the sake of simplicity. That said, the name Radiation is pretty generic and I think most people would call PertubationAdvection, Orlanski, the oblique scheme, and even Flather, radiation schemes.

@simone-silvestri from talking to you at the conference it seems like what you're implementing here for the baroclinic flow is an Orkanski scheme, no?

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.

I'd say just replace the existing periodic one with this. We don't need both.

…ith `SplitExplicitFreeSurface` (#5378)

* Validate any bc on nothing-dimension fields

* Add U and V to field location map

* Fix missing u, v reference

* Add missing b_bg helper function for sponge layer

* Maintain barotropic velocities for free surface

Co-authored-by: Simone Silvestri <[email protected]>

* Add sea-surface height to field location map

Co-authored-by: Simone Silvestri <[email protected]>

* Remove extra validation dispatch

---------

Co-authored-by: Simone Silvestri <[email protected]>
glwagner and others added 30 commits March 21, 2026 06:53
XLA folds Δt + zero(Δt) → Δt, so the aliasing between last_Δt and
last_stage_Δt persists. Use Reactant.Ops.optimization_barrier which
XLA cannot fold, ensuring distinct buffers for the while loop carry.

Co-Authored-By: Claude Opus 4.6 (1M context) <[email protected]>
Add `materialize_clock!(clock, timestepper)` called in HFSM and NH model
constructors. For QuasiAdamsBashforth2TimeStepper, this sets
`clock.last_Δt = clock.last_stage_Δt` to ensure they are distinct objects.
This is needed for Reactant, where aliased ConcreteRNumber fields cause
XLA buffer donation errors in compiled loops.

With materialize_clock! breaking the alias at construction time, the
optimization_barrier workaround in the Reactant tick! override is no
longer needed and is removed.

Co-Authored-By: Claude Opus 4.6 (1M context) <[email protected]>
materialize_clock! in the constructor doesn't prevent MLIR-level aliasing
that occurs at runtime inside tick! when both last_Δt and last_stage_Δt
are assigned the same Δt value. The optimization_barrier is still needed
until Reactant provides a proper buffer copy op.

See #5389 (comment)

Co-Authored-By: Claude Opus 4.6 (1M context) <[email protected]>
materialize_clock! in the model constructor ensures last_Δt and
last_stage_Δt are distinct ConcreteRNumber objects. Assigning the same
value to both in tick! via .mlir_data does not re-alias the objects.

Co-Authored-By: Claude Opus 4.6 (1M context) <[email protected]>
For QAB2, last_Δt and last_stage_Δt always hold the same value after
tick!. In the Reactant extension, materialize_clock! aliases them via
setfield! (bypassing the ReactantClock setproperty! override) so that
Reactant's tracer sees one buffer, avoiding XLA buffer donation errors.

The src/ QAB2 method is removed since aliasing is only needed for
Reactant's ConcreteRNumber fields (normal Float64 fields are immutable
values with no aliasing concern).

Co-Authored-By: Claude Opus 4.6 (1M context) <[email protected]>
…ananigans.jl into ss/prepare-for-open-boundaries
…liMA/Oceananigans.jl into glw/remove-initialization-update-state
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

6 participants