Drop the seam hand-off corner in section-finding, not in transports - #11
Drop the seam hand-off corner in section-finding, not in transports#11hdrake wants to merge 6 commits into
Conversation
A physical vorticity corner can be stored under more than one index: the
duplicated seam column of a symmetric ('outer') periodic grid, a corner shared
by two tiles of a multi-tile grid, or the two lips of a grid cut. The walk in
`infer_grid_path_from_geo` deliberately steps through both identities -- it
admits the "seam twin" of the corner it is standing on -- so that seam crossings
are decided by the grid's topology rather than by floating-point rounding. That
hand-off step spans no grid cell, so the pair of corners is not a velocity face.
Until now the redundant identity survived into every public section object and
was filtered three separate times inside `transports.uvindices_from_qindices`:
a node-canonicalisation plus dedup at the top, a `var="0"` zero-flux placeholder
in `_uv_for_edge`, and a trailing geographic zero-length-face mask at the bottom.
Move the removal to where the path is built instead, so the public API never
surfaces the artificial pair in the first place.
`section.drop_repeated_corners` normalises a traced path so that
* no two consecutive corners are the same physical point -- every consecutive
pair is a real velocity face; and
* on a multi-tile grid, every corner is the canonical native index of its
corner node (`_OuterTopology.node_native`), which is what lets the face
attribution read a corner in either the source or destination face's frame.
It runs from `grid_section`, from `GriddedSection.__init__` when indices are
supplied directly, and from `utils.load_gridded_section`, so the invariant holds
at every entry point -- including for sections saved before the drop existed,
which still contain both identities and are normalised on load.
Of the two identities the later one is kept -- the frame the path continues in --
so the faces on either side of the hand-off keep exactly the indices they had
before: the face past the seam is untouched, and the face before it becomes a
wrapped step. The single-tile index arithmetic gains the Y counterpart of its
existing X wrap rule, and wraps a seam crossing's cell-center index back into
range; both are inert for paths that never hand off across a seam.
With the invariant established up front, all three filters in `transports.py`
are unreachable and are removed. `_uv_for_edge` now raises instead of returning
a zero-flux placeholder, so an edge whose normal velocity is stored on neither
the source nor the destination face is reported rather than silently dropping
real transport.
Numerics are unchanged: the ECCOv4r4 LLC90 overturning streamfunction is
bit-identical at every latitude tested, over the same number of velocity faces.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
hdrake
left a comment
hdrake
left a comment
There was a problem hiding this comment.
These concerns need to be addressed.
| if f_c is not None and not _has_corner_topology(grid): | ||
| # A multi-tile grid without tracer-center coordinates has no corner-node table | ||
| # to canonicalise against (`build_neighbor_maps` falls back to reading xgcm's | ||
| # corner halos for these). Its two identities of a shared boundary corner are | ||
| # separate indices on separate faces with no common index to collapse to, and | ||
| # dropping one would leave a step that is not a neighbor step. Such a grid | ||
| # carries no velocities either -- `uvindices_from_qindices` needs the centers, | ||
| # so no face is ever derived from it -- so the path is left as traced. | ||
| return (i_c, j_c, f_c, | ||
| np.asarray(glon[f_c, j_c, i_c]), np.asarray(glat[f_c, j_c, i_c])) |
There was a problem hiding this comment.
I don't understand this. Does this also apply to multi-tile grids with that are built with outer axes? I don't really understand why those would need tracer coordinates. (By tracer coordinates, do you mean center dimensions? I understand that velocity faces have one dimension/coordinate corresponding to the tracer grid, but where do those coordinates actually get used in sectionate? Do they need to be, or was that an unnecessary restriction?
In any case, any grid for which we cannot find velocities should raise an error as soon as it is passed to sectionate.
There was a problem hiding this comment.
Taking your questions in order.
Does it apply to multi-tile outer grids? It did, but only to those that also declare no center dimensions — the branch in build_neighbor_maps keys off centers, not off corner position. A multi-tile outer grid with center dimensions goes through outer_topology and is fully normalised. Only a center-less one fell back to _multitile_padded_maps, and those maps step onto each face's own copy of a shared seam corner, so a section traced on them keeps the zero-length hand-off step and cannot be normalised away.
Yes, I meant center dimensions, and the old wording overstated the requirement — I've fixed it here and in build_neighbor_maps. Only the dimensions are needed, not tracer coordinate values: _OuterTopology numbers cells with np.arange(nf*Nyc*Nxc) and never reads a tracer longitude or latitude. A center entry naming a bare dataset dimension is enough.
Where they actually get used, and whether the restriction is real. Three places. The corner-node topology is rebuilt by fingerprinting the tracer cells around each corner — cells pad reliably across any seam, corner arrays don't — so it needs the center dimension names and lengths. _in_velocity_range needs the center axis lengths to decide whether a V point at (X-center, Y-corner) is on-face. And convergent_transport/extract_tracer index the velocities on those dimensions. That last one is really your own framing back at you: you noted velocity faces have one dimension on the tracer grid, and that is exactly the point — no center dimensions means no velocities can exist. So the requirement is real, not artificial. What was unnecessary was implying that tracer coordinate values were needed.
Acting on your last sentence. grid_section now refuses a multi-tile grid with no center dimension in _check_supported_topology, naming the offending axes and showing the coords={...} to add, and the degraded fallback in drop_repeated_corners is gone, so the invariant holds unconditionally. Where velocities are actually looked up, coord_dict now raises a named error instead of a bare KeyError: 'center'.
One judgement call I'd rather you made than me. I stopped short of the fully literal reading. Single-tile corner-only grids also have no velocities, but they are not degraded — their seam twins are found geometrically from the corner coordinates, so those sections normalise correctly; they simply can't produce transports. Refusing them at grid_section would break test_section_cornercases.py and test_section_class.py, whose corner-only grid predates this PR, so it's an API change I didn't want to make unilaterally here. Happy to extend the check to all grids if that's what you want.
Related: with the front door closed, _multitile_padded_maps is now reachable only by calling build_neighbor_maps directly. Its three unit tests still pass and I've left it alone, but it's a candidate for removal if you'd rather not carry it.
Written by Claude Code.
There was a problem hiding this comment.
Remove build_neighbor_maps
There was a problem hiding this comment.
Removed _multitile_padded_maps — I read this as the fallback rather than build_neighbor_maps
itself, since grid_section and uvindices_from_qindices both call the latter and CLAUDE.md
documents it as a core utility. Say the word if you meant something broader.
With the front door closed to centre-less multi-tile grids it was unreachable anyway, so it and its
three unit tests are gone, along with _validate_reciprocity, which was only ever called from it.
build_neighbor_maps now always returns maps that step through a single canonical index per
physical corner, which is what everything downstream assumes; its docstring says so directly rather
than describing a branch that no longer exists.
Written by Claude Code.
Review follow-up on #11. Drop the two backwards-compatibility normalisation calls. `master` is at v0.3.4, where `save_gridded_section` does not persist `f_c` at all, there is no multi-tile support, and a symmetric periodic grid had its duplicated seam column stripped before the walk (`gridlon[:,:-1]`), so no released version could have written a saved section containing a repeated corner. `GriddedSection.__init__` and `load_gridded_section` are restored byte-for-byte to their previous form. Removing them would have left one path that silently derives a wrong face: the single-tile index arithmetic turns a repeated-corner pair into a spurious duplicate of a neighbouring face, double-counting its flux. So `uvindices_from_qindices` now *checks* the invariant it documents and raises, naming `drop_repeated_corners` as the remedy for hand-built indices. Checking is a better contract than silently fixing up: normalisation stays in exactly one place, section-finding. Also act on "any grid for which we cannot find velocities should raise an error as soon as it is passed to sectionate". Velocities are staggered onto the cell-center dimensions -- U at (X-corner, Y-center), V at (X-center, Y-corner) -- and the multi-tile corner-node topology is rebuilt by fingerprinting the tracer cells around each corner, so a multi-tile grid with no center dimensions can neither carry velocities nor resolve its seam corners to one identity. `_check_supported_topology` now rejects it at the front door instead of letting `drop_repeated_corners` fall back to returning an unnormalised path; that fallback is deleted. Where velocities are actually looked up, `coord_dict` replaces a bare `KeyError: 'center'` with a message naming what is missing and how to declare it. The requirement is narrower than the old wording implied, and the wording is fixed throughout: only the center *dimensions* are needed, not tracer longitudes and latitudes -- `_OuterTopology` numbers cells with `np.arange` and reads no tracer coordinate values. Single-tile grids are unaffected: their seam twins are found geometrically from the corner coordinates alone, so a corner-only single-tile grid still traces, exactly as it has since before this PR. Numerics are unchanged: the ECCOv4r4 LLC90 overturning streamfunction is still bit-identical to the unmodified base at every latitude tested, over the same number of velocity faces. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The new guard in `uvindices_from_qindices` used `<= COINCIDENT_TOLERANCE_M` where `drop_repeated_corners`, the walker's stop test and its seam-twin test all use `<`. At exactly the tolerance the two disagreed: normalisation would keep the pair and the guard would reject it, so a section straight from `grid_section` could fail the very invariant `grid_section` establishes. Unreachable in practice at 1e-3 m, but the two definitions should be the same one. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Dropping a repeated corner at trace time turned two verified sections on the
MOM6 tripolar example grid into wrong answers, because that grid has three
degenerate corner columns -- the bipolar cap's two poles at 60E and 120W, and
60E's periodic twin -- where 41 corner rows are one physical point:
* A meridional section at 120W (50N -> 85N) collapsed 41 rows into a step of
41 along "Y". The wrap rule added here read that as a periodic seam
crossing, flipped the direction, and put the U face on `j = ny`, one past
the last cell center -- `convergent_transport` then raised IndexError.
Neither wrap rule was gated on its own axis' periodicity; a bipolar-fold
grid's "Y" axis does not wrap. Both are now gated (X too, which is a
pre-existing behaviour change and has its own test), so doubly-periodic
idealized configurations keep them and everything else does not.
* An oblique section entering the 60E pole along corner row 180 and leaving
it along row 140 collapsed to the last index of the run, which is adjacent
to the corner after it but not to the one before it -- turning a V face
into a U face. Which index survives is no longer a guess: each candidate is
checked against the faces the raw path implied, using the very arithmetic
`uvindices_from_qindices` applies, and only one reproducing both flanking
faces and their transport signs is kept (last first, which is what a seam
hand-off needs). Where the path enters a degenerate column on one side and
leaves on another, no index can name both faces; that now raises and names
the corner instead of emitting a face the section does not have. It only
arises on a tripolar grid whose fold is left undeclared -- with the fold
declared the walk crosses the seam and never enters the column.
Also from review:
* `create_section_composite` is public and notebook 4 uses it as the
lower-level entry point, but it returned un-normalised corners that
`uvindices_from_qindices` now rejects. It takes an optional `grid=` and
normalises when given one; `grid_section` passes it, so the normalisation
lives in one place.
* A corner-only single-tile grid traces and enumerates its faces again. The
single-tile face derivation looked up the cell-center dimensions up front,
where it needs them only to wrap a periodic seam crossing back into the
centers' range.
* `_multitile_padded_maps` is gone, with its three unit tests and the
reciprocity check that only it used. With `grid_section` refusing a
centre-less multi-tile grid at the front door it was unreachable except by
calling `build_neighbor_maps` directly, and it produced seam-twin maps that
cannot satisfy the invariant.
* The single-tile face arithmetic is factored into `_single_tile_face`, so
section-finding and transports cannot drift apart.
* "identities", for the several indices of one physical corner, is replaced
throughout by "repeated corner" / "one corner stored under two indices".
`_uv_for_edge`'s error names its corners as (face, j, i) and suggests
rerouting; `coord_dict`'s error agrees in number and no longer claims that
tracing always works without cell centers; and the fold tests describe the
xgcm behaviour they need instead of citing a merged PR.
Transports are unchanged: every face of a 77-section sweep over the MOM6 fold
grid matches the previous behaviour exactly, and the only difference in the
independent transport check is the intended `n_corners` reduction.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The error already says a section cannot pass through the degenerate column at
a bipolar cap's pole and names the corner it entered and left by, but it only
offered "route around the pole". Measured over 200 random sections on the MOM6
tripolar grid, this case never arises when the fold is declared
(padding={"Y": {"fold": "corner"}}): 200/200 trace, with face lists identical
to before the drop. With the fold left as "extend" it arises in 7. Declaring
the fold is therefore the actual fix in almost every case, so say so first and
keep rerouting as the fallback.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Two regressions found and fixed, plus the review findingsAn independent reviewer audited this PR and found two regressions the 83-test suite did not 1.
Both wrap rules — the new Y one and the pre-existing X one — are now gated on whether that axis is 2. A V face silently became a U face. Where a collapsed run changes both indices, The collapse rule is now verified rather than positional: each candidate index is scored by Where no candidate qualifies, it now raises rather than emitting a wrong face. That is a
So a correctly-declared tripolar grid sees no change at all. The losses are confined to grids 3. Corner-only single-tile grids regressed — Also in this round: Numerics unchanged: reference values captured from the unmodified base still differ in exactly One thing this does not fix, reported rather than papered over: Written by Claude Code. |
Two conflicts, both additive documentation: - `create_section_composite`'s parameter list: #10 documented `curve`, this branch documented `grid`. Both kept. - `CLAUDE.md`'s `section.py` bullet: the base now states the three-curve contract from #10; this branch states the repeated-corner invariant. The base's wording is authoritative on curves and carries hdrake's own edits, so the invariant sentences are spliced into it. Two tests on this branch encoded #10-era path lengths that the corner dropping changes: `test_latitude_circle_takes_shortest_path_west` (the seam vertex is one corner under two indices, so the path starts at index 6 rather than carrying both) and `test_parallel_is_held_where_the_geodesic_bows` (the geodesic crosses the periodic seam diagonally, 81 points rather than 83). Expectations updated and the reasons recorded in each docstring. Suite: 96 passed, 0 skips. Transports re-verified against the unmodified base: still exactly one recorded difference, `mom6_fold.zonal_periodic.n_corners: 308 -> 307`. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
1 participant
This is the large change singled out in the review of MOM6-community#47 (MOM6-community#47): "changing at which step in the
sectionateworkflow the artificial cross-tile faces with zero physical length get dropped". Addressesr3753380399,r3753419093,r3753431811,r3753439048. Targetstopology-driven-neighbors.The problem
A physical vorticity corner can be stored under more than one index: the duplicated seam column of a symmetric (
'outer') periodic grid, a corner shared by two tiles, or the two lips of a grid cut. The walk deliberately steps through both — that is what makes a seam crossing depend on the grid's topology rather than on floating-point rounding, and it is unchanged here.But that hand-off step spans no grid cell, so the pair is not a velocity face. It used to survive into every public section object and be filtered three separate times inside
transports.uvindices_from_qindices.The change
Normalisation happens in exactly one place — section-finding. Everywhere else verifies.
section.drop_repeated_cornersruns immediately after the walk, so every sectiongrid_sectionreturns satisfies: no two consecutive corners are the same physical point, and on a multi-tile grid every corner is the canonical native index of its corner node.uvindices_from_qindiceschecks that invariant and raises rather than re-establishing it — without the check, a hand-built path repeating a corner quietly became a spurious duplicate of the neighbouring face and double-counted its flux.All three filters are gone, plus the dead NaN branch in
uvcoords_from_uvindices._uv_for_edgeraises instead of returning a zero-flux placeholder, because silently contributing zero would drop real transport rather than report a problem.create_section_compositetakes an optionalgrid, so the public lower-level entry point can produce a section the rest of the API accepts.Which index a run collapses onto
Choosing "the last" is not safe in general: where a collapsed run changes both indices,
zonalflips and a V face becomes a U face. The rule is now verified rather than positional — each candidate index is scored by running the same face arithmeticuvindices_from_qindicesapplies, on both flanking steps, and kept only if it reproduces the raw path's(var, i, j)and the direction flag carrying that face's sign. Collapsing therefore provably changes no face and no sign.Where no candidate qualifies it raises, naming the corner and the indices the section entered and left by. This happens at a bipolar cap's pole, where a whole column of corner rows degenerates to a single point. Over 200 random sections on the MOM6 tripolar grid:
Y: {"fold": "corner"}Y: "extend"A correctly-declared tripolar grid sees no change at all. The losses are confined to grids whose metadata does not declare the fold, where the walk was passing straight through 41 stored copies of the pole point — so the error points at declaring the fold first, and rerouting second.
Regressions found by an independent review and fixed
Full detail in a comment on this PR.
IndexErroron a plain meridional section at 120°W, because a newly-added Y wrap rule read the 41-row pole-column collapse as a periodic seam crossing and took the velocity index from the wrong end (max_j(U)180 against an axis of size 180). Both wrap rules are now gated on whether that axis is actually periodic — kept, because doubly-periodic idealized configurations are real, but a jump of more than one index can only mean a seam when the axis wraps.coord_dictwas called unconditionally at the top of the single-tile loop, so a grid with no centre dimensions raised where base returned faces. Now looked up only where used.Also:
_multitile_padded_mapsremoved (with_validate_reciprocity, reachable only from it) perr3754203594; "identities" replaced with "repeated corner" throughout, matching the publicdrop_repeated_corners;xgcm#711references removed per the no-links-in-source rule;_uv_for_edgeandcoord_dictmessages made actionable.Numerics
Reference values captured from the unmodified base and compared under two
PYTHONPATHs — corner and face counts,psi_min/psi_max, the first five overturning levels, U/V counts, index heads and theLsignsum — across six ECCO LLC90 latitude circles and three MOM6 tripolar-fold sections. Exactly one recorded value changes:which is the point of the PR. Everything derived from it is identical.
One thing not fixed, reported rather than papered over:
uvindices["Xinc"]can differ from base at one position on a U face. Transports are provably unaffected —Usignuses~Yinc,VsignusesXinc, combined asUsign*Umask + Vsign*Vmask, so the stale flag is masked to zero — butXincis a documented output and is inaccurate in that narrow case.Tests
87 passed, 0 skips, indocs_env_sectionate_pr47-drop-zero-length-faces.Notebook 4's cell 10 calls
create_section_compositewithoutgrid=; adding it makes that cell agree withgrid_section(34 corners rather than 35). The cell only plots, so it executes either way.🤖 Generated with Claude Code