Fix spurious sea-surface height from p-star partial-cell snapping - #674
Merged
Conversation
xylar
force-pushed
the
fix-p-star-partial-cells
branch
from
July 28, 2026 11:55
f3ce936 to
1228295
Compare
6 tasks
Collaborator
Author
|
@katsmith133, this bug fix relates to the realistic, global initialization. I know it's not really in your wheelhouse yet but could I ask you to review nonetheless? There's not really anyone else who is more familiar than you at the moment anyway that I can turn to. If you just give the code a look over and let me know if things are clearly wrong or confusing, that would be a help. |
xylar
force-pushed
the
fix-p-star-partial-cells
branch
from
July 29, 2026 10:16
1228295 to
03627a1
Compare
The p-star init iteration matched the geometric water-column thickness to the target bathymetry and pinned bottomDepth to the raw target, letting ssh float as a diagnostic. For cells whose seafloor lands in the partial-cell snap band, the snap freezes the pseudo bottom depth so the geometric column cannot equal the target, and the entire quantization residual surfaced as a nonzero ssh (up to ~24 m on ~2% of cells). Flip the invariant to match z-star partial cells: ssh is prescribed (sea_surface_height) and bottomDepth is diagnosed from the column. The geometric column is shifted so its top interface lands on the prescribed ssh; the snap residual then adjusts bottomDepth (the representable bathymetry) instead of ssh. For representable cells the shift is zero and bottomDepth still equals the target to machine precision; ssh is now exactly its prescribed value everywhere. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Being at rest (zero velocity) does not imply ssh = 0: the sea surface is depressed by surface loading (atmospheric pressure, sea ice, ice shelves) even at rest. The zero-surface-pressure coordinate is the reference coordinate, not the resting coordinate. Remove the resting-vs-zero-ssh conflation from comments and docstrings in the p-star vertical coordinate and design doc. No behavior change. Cherry-picked from add-realistic-global-ocean-general-forward (9f60956), dropping that commit's polaris/tasks/ocean/realistic_global/ hunks, which do not exist on main. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The message said "full-cell snap" when partial-cell snapping triggers it too, and framed an unreachable bathymetry as a failure. Since the column is now anchored at the prescribed sea surface, snapping to the nearest representable depth is the intended behaviour: ssh stays exact and bottomDepth absorbs the residual. Reword accordingly, drop it from warning to info, and report how far the sea floor had to move. Note that this branch does not fire often: the check requires every column to be frozen at once ((...).all()), so on a multi-column mesh a handful of snapped columns among many converging ones never reaches it. Making the report per-column would be a separate change. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The stagnation branch requires *every* column to be frozen at once, so on a mesh where some columns are still converging it never fires and the loop simply runs out of iterations, silently. A post-loop `info` message now reports how many columns cell snapping could not place on the requested bathymetry and how far `bottomDepth` moved, however the loop ended.
Polaris and Omega build the geometric column from opposite ends: PStarInitStep anchors at the prescribed sea surface and diagnoses bottomDepth, while VertCoord::computeGeomZHeight anchors at BottomGeomDepth and diagnoses ssh. They agree only because run_pstar_init writes back the shifted seafloor rather than the raw target bathymetry, which until now was guaranteed by a comment. Add a test that reimplements Omega's upward accumulation in NumPy -- deliberately not calling geom_height_from_pseudo_height, since re-anchoring the written column with the helper that built it is an identity -- and checks it reproduces GeomZInterface and ssh from the written bottomDepth, SpecVol and PseudoThickness. The check only has content where the surface-anchoring shift is nonzero, so the snapped cases use partial cells with a target bathymetry between layer interfaces (455 m snaps up to 450 m, 460 m down to 465 m) and assert the shift is nonzero before asserting anything else. A converged, unsnapped control marks the boundary of what those cases add. Writing the raw pre-shift bathymetry instead of the shifted one breaks the reconstruction by 5 m (the snap residual) on all three snapped cases and leaves the control at 6e-14, confirming the test fails against the defect it exists to catch. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
xylar
force-pushed
the
fix-p-star-partial-cells
branch
from
July 31, 2026 15:45
03627a1 to
fa10bff
Compare
Collaborator
Author
|
@katsmith133, another ping that I could use your review on this. Could you let me know when you might have time? |
Collaborator
Author
|
I'll report my testing shortly... |
Contributor
|
Apologies, I don't remember getting pinged the first time. I just looked through and it all makes sense to me. I think if the testing looks good, I am good with approving it. |
Collaborator
Author
TestingPolaris
|
xylar
commented
Jul 31, 2026
Comment on lines
-158
to
+160
| the resting surface-pressure depression for a reference-density fluid | ||
| the sea-surface depression that balances the surface pressure for a | ||
| reference-density fluid (being at rest does not imply $\eta_0 = 0$; surface | ||
| loading depresses the sea surface even at rest) |
Collaborator
Author
There was a problem hiding this comment.
We don't usually retroactively update design docs but this seems fine for clarity.
Collaborator
Author
|
Thanks @katsmith133. Take a quick look when you can, please. |
katsmith133
approved these changes
Aug 1, 2026
katsmith133
left a comment
katsmith133
left a comment
Contributor
There was a problem hiding this comment.
Approved based upon visual inspection and testing by @xylar
Collaborator
Author
|
Thanks @katsmith133! |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
With
coord_type = p-starandpartial_cell_type = partial, a resting, unforced column with zero surface pressure could be initialized with a non-zero sea-surface height, which is a spurious barotropic pressure gradient.Cause: partial-cell snapping quantizes the achievable pseudo bottom depth (
min_pc_fractionforbids thin bottom cells by design), so some target bathymetries are unreachable and therun_pstar_inititeration stagnates on them.geom_height_from_pseudo_heightanchored the column at the requested seafloor, sobottomDepthcame out exactly on target and the whole quantization residual surfaced asssh.Fix: flip the invariant to match z-star partial cells —
sshis prescribed,bottomDepthis diagnosed. The converged column is shifted so its top interface lands onsea_surface_height, andbottomDepthfollows from the shifted bottom interface. Where the iteration converges the shift is zero andbottomDepthstill equals the target to machine precision; where it cannot, the residual moves the seafloor (the representable bathymetry, which is whatpartial_cell_typeis documented to alter) instead of the sea surface.The first two commits here are cherry-picked from #643, where this was independently hit and fixed on realistic global setups; 8427a63 drops that branch's
polaris/tasks/ocean/realistic_global/hunks.Checklist
Testingcomment in the PR documents testing used to verify the changes