Resolve corner identity across the bipolar fold seam

[hdrake]

Targets topology-driven-neighbors (the branch behind MOM6-community#47), not master.

What & why

@geoffstanley asked on #47 how sectionate handles the ORCA-style duplicated, sign-flipped velocity rows at a bipolar fold, and how it avoids double counting if its list of F-point indices contains no duplicates. Chasing that down turned up a real bug, which this PR fixes.

The two halves of the question have opposite answers.

The sign flip needs no fold-specific handling. Sectionate never compares an index-frame velocity against a geographic direction: on a single-tile grid it compares the stored velocity against the traversal direction in the same index frame, and on multi-tile grids it evaluates both the velocity's positive direction and the section's direction of travel geographically. Both inputs flip together, so a frame reversal cancels. Sections running along the duplicated seam row — including straight through the fold pivot — were already exact and still are.

Corner identity does need handling, and single-tile grids were missing it. A seam gives one physical corner two index representations. Crossing a fold, the walk steps from one of them straight to a corner adjacent to the other: a real physical edge, but not an index-adjacent one. The single-tile arithmetic in uvindices_from_qindices assumes index adjacency and reads the velocity column off the source corner, which is then the wrong representative. It named the mirror image of the face the section actually crossed, with the opposite sign, and on a closed crossing emitted one face twice — exactly the double counting Geoff predicted, for the reason he gave.

On the real MOM6 tripolar grid, a closed 80°N latitude circle whose true transport is zero returned 1.26e-02; a section over the pole was off by 1.8%. In both cases exactly one edge was wrong — the one leaving the seam. Entering it is fine, because there the source corner is the off-seam one and its column is the right representative.

This is not a regression from MOM6-community#47. master's topology="MOM-tripolar" branch takes the same non-adjacent step across the seam (up = (j-1, (nx-1) - i%nx)) and feeds it to the same Xinc/Yinc arithmetic. What MOM6-community#47 changed is where the connectivity comes from, not how the crossing edge is attributed.

Changes

transports._insert_seam_twins splices the seam twin back into the corner chain wherever a step skips it, so every remaining step is index-adjacent. The crossing becomes one ordinary edge plus one zero-length twin edge, and zero-length edges already emit no velocity face — so nothing is counted twice. It runs inside the deterministic corner→face derivation, which means sections reloaded from saved (i_c, j_c) get it too, and it is a no-op unless a step actually skips a twin.

This is the single-tile counterpart of the corner canonicalization _OuterTopology already performs for multi-tile grids, which is why LLC90 and the cubed-sphere were already exact. It is also the same construction Geoff proposed, and structurally what NEMO does when it carries the duplicated points and masks them out of transports with umaskutil/vmaskutil.

A declared fold whose coordinates do not carry the fold symmetry now raises. With no twins to splice, the crossing edge is genuinely undefined; the error says so instead of silently attributing it elsewhere.

grid_section output is unchanged — the walk is untouched, and docs/make_algorithm_animation.py's replay assertion still matches.

Tests

All use a streamfunction fixture: ψ is a single-valued function of physical position, so seam twins get identical ψ, and transports are built from index-space differences of it. That makes the flow exactly non-divergent, makes the expected answer analytic (ψ(end) − ψ(start)), and reproduces the duplicated, sign-flipped seam row of an ORCA-style fold (NEMO 4.2 manual, Appendix E) without hard-coding the flip — the fixture asserts that structure rather than assuming it.

• _pinched_fold_grid — a new synthetic grid that carries the corner-pivot fold identity exactly (a cylinder pinched shut at the top, so the seam row is mirror-symmetric under i ↔ nxh−i), needing no downloaded data. Unlike the existing _fold_grid, whose plain lat/lon coordinates do not carry the fold its metadata declares.
• Crossing the seam is exact, and uses each velocity face exactly once.
• Running along the seam is exact (the sign-flip case; guards against regressing what already worked).
• A fold-crossing on coordinates lacking the fold symmetry raises.
• Real MOM6 tripolar grid: traced sections crossing the seam are exact and count no face twice (skips without the example data, like the existing real-grid fold test).

With the fix reverted, three of these fail — the real-grid one by exactly the 1.26e-02 quoted above — so they genuinely guard the behaviour.

Verification

• Full suite: 81 passed, 0 skipped, stable across PYTHONHASHSEED.
• All five example notebooks re-executed cleanly. No figure and no numerical output changed — none of them computes transports across a single-tile fold seam — only execution metadata and one line number in a captured warning traceback.
• make html SPHINXOPTS="-W" builds clean.
• docs/make_algorithm_animation.py replay still matches grid_section (9 corners, 8 steps); committed artifacts unchanged.

Docs

CLAUDE.md gains a Seam twins (corner identity) concept entry covering both the multi-tile and single-tile paths and why the sign flip needs no special handling. docs/source/algorithm.md gains a paragraph on why the chain is resolved to index-adjacent steps before faces are named.

---

AI-assisted PR, disclosed per the project AI Usage Policy: commits carry Co-Authored-By: trailers. 🤖 Drafted with Claude Code.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[geoffstanley]

Hey @hdrake, nice work fixing this so quick! Glad that the sign flip was handled correctly, but good to catch the double counting. (The double counting seems the more important thing, because once that is fixed the sign flip becomes irrelevant -- at least for the NEMO style grid in which only duplicated points carry a sign flip. I'm not sure if other grids behave differently.)

If I understand _insert_seam_twins correctly, it looks at each pair of indices (ja, ia) and (jb, ib) that are adjacent in the input lists, and then, if these points are not _index_adjacent (ie they differ by 1 in either i or j, accounting for simple periodic wrapping), a new point is inserted between them. The inserted point is

• the same geographic location as (ja, ia) but has (j,i) indices that are _index_adjacent to (jb,ib), or
• the same geographic location as (jb, ib) but has (j,i) indices that are _index_adjacent to (ja,ia).

This will work IF the transect only takes one step to cross the seam, but that is not the only case. Consider for example a transect that zigs along j up into the seam, then zags along i for a few grid points, then zigs up again and gets out of the seam.

In this case where there are multiple twin points, _insert_seam_twins will interleave the twins. For example, if the (j,i) points are
[..., a, b, c, d, E, F, ...]
and there are twins c =C, and d = D,
and the index_adjacent pairs are a ↔ b, b↔c, c↔d, C↔D, D↔E, E↔F
(think lower case points are on the left half of the bifold, and upper case on the right side),
then _insert_seam_twins will give
[..., a, b, c, C, d, D, E, F, ...]
but what we want is
[..., a, b, c, d, C, D, E, F, ...].
The former enters the seam region, then jumps back and forth within the seam, then exits the seam. The transport accumulation algorithm (starting from uvindices_from_qindices, I think) won't work properly.
The latter enters the seam region and goes as far as it can continuously in index space, and then jumps once to the other side of the seam, also jumping backwards in geographic space, to repeat the section in geographic space but with the other set of indices. The transport accumulation algorithm, appropriately masking out duplicate U and V points, ought to work well with that.

Have I got that right?

If so, as I said I did it by running the standard GriddedSection algorithm twice, once forwards and once with the section reversed. I then merged the lists of indices using a function like

def interleave_first_difference(A, B):
    """
    Given two lists A and B, return a list that is A until the first index that A and B
    equal each other after the first index that they differ from each other, followed by
    B from the first index A and B differ from each other.

    Parameters
    ----------
    A, B : lists

    Returns
    -------
    C : list

    Example
    -------
    # A and B the same except for a chunk in the middle
    A = [1, 2, 3, 4, 7, 8]
    B = [1, 2, 5, 6, 7, 8]
    interleave_first_difference(A, B)
    # produces [1, 2, 3, 4, 5, 6, 7, 8]

    # A and B completely different
    A = [1, 2]
    B = [3, 4]
    interleave_first_difference(A, B)
    # produces [1, 2, 3, 4]
    
    # A and B identical
    A = B = [1, 2]
    interleave_first_difference(A, B)
    # produces [1, 2]
    """
    n = len(A)
    if n != len(B):
        # This function is intended for lists of the same length;
        # it could work for lists of different lengths but it becomes more
        # complicated to describe. 
        raise ValueError(f"A has length {n} but B has length {len(B)}")
    i1 = 0
    i2 = n
    for i in range(n):
        if A[i] != B[i]:
            i1 = i
            break
    for i in range(i1, n):
        if A[i] == B[i]:
            i2 = i
            break

    return A[:i2] + B[i1:]

It would be nice not to have to run the basic algorithm twice though.

Instead, we could have a function _swap_twins(j_c, i_c) that iterates through lists of i and j coordinates and replaces any point with its geographic twin, when such a twin exists. Then use interleave_first_difference on our initial list of points and the _swap_twins output.
Specifically, given (j_c, i_c) as a list of points, something like

A = list(zip(j_c, i_c))
B = list(zip(_swap_twins(j_c, i_c)))
C = interleave_first_difference(A, B)
j_c, i_c = zip(*C)

[hdrake]

@geoffstanley , that's a good point about the zigzagging sections and a corner case I hadn't considered.

I'm going to take a step back and think about what the best path forward is for this mega-PR MOM6-community#47

These are the constraints that I am working with, and which are often in conflict with each other:

1. Sectionate should be able to both create sections between adjacent corners and diagnose transports across real grid faces (as opposed to steps between grid indices that correspond to the same physical point), using only xgcm.Grid's grid metadata (i.e. without hard-coding any information about the specific grid, like we currently do with both periodic boundaries and the tripolar seam)
2. These sections should play nicely with my regionate package, which has two core functionalities:
   i. When given a closed section (e.g. as provided by sectionate), it finds the tracer cell masks that this section bounds
   ii. When give a tracer cell mask, it finds the closed section that bounds it
3. All of this should work seamlessly for three main grid topologies:
   i. curvilinear (cartesian/latlon)
   ii. curvilinear but tripolar with northern fold
   iii. multi-tile

Everything seemed pretty straight-forward in my head until I ran up against the reality that the grid topology gets very messy with output from models like ECCOv4r4, which provide outputs on center and right C-grid staggers. The right stagger means that each tile of the ECCOv4r4 grid is missing some of its corners. It gets even worse if two of these tiles are connected, because then there can be corners that don't appear anywhere in the dataset and thus are just implied by the topology. This led me down the rabbit hole of figuring out how to build outer representations of every grid, which gets kind of complicated.

The code currently makes some different choices for whether or how these outer representations are needed, and how sectionates walks through them (e.g. whether or at what stage in the algorithm to drop the repeated points), which is where the problems described in this PR come in.

Any thoughts welcome!

[geoffstanley]

I can see the complexity! I have just been working with NEMO's T-pivot and F-pivot tripolar grid, and not in conjunction with xgcm, so it's been easier for me.

On 1., does xgcm contain carry the umaskutil and vmaskutil masks? I think it does not currently. So then Sectionate needs to carry some information about the grid beyond what xgcm has, in order to correctly compute transports that cross the north fold. Or xgcm needs to evolve to carry those masks. Or another approach is needed, in which there is precisely one (geographically) duplicated F point in the list (i_c,j_c), and the flipped sign of certain U and V grid cells is accounted for; call this Approach 2. But from my experience with NEMO, that still requires umaskutil and vmaskutil.

By the way, I didn't quite understand this:

Sectionate never compares an index-frame velocity against a geographic direction: on a single-tile grid it compares the stored velocity against the traversal direction in the same index frame, and on multi-tile grids it evaluates both the velocity's positive direction and the section's direction of travel geographically

If Sectionate does indeed handle the sign flip properly, then perhaps Approach 2 can work, without *maskutil. Can you point me to the code where this happens?

On 2., would it help to bring regionate into sectionate? This won't fundamentally address anything of course, but it might make it easier to update the Regionate functionality as Sectionate evolves. This is a thought I had when I first saw your two packages, actually.

Your description of the ECCOv4r4 grid sounds like xgcm is not capable of dealing with it completely. If that's so, then taking a step back, either Sectionate relies fully on xgcm for neighbours and it can't handle certain rare cases of the ECCOv4r4 grid (and perhaps others?), or Sectionate needs to do neighbours itself (or perhaps somehow in conjunction with xgcm, just patching as needed). It does seem like the right approach to upgrade xgcm to handle these cases, which you've been working on, but I'm not familiar enough with xgcm to know how feasible that is.