Skip to content

Drop the seam hand-off corner in section-finding, not in transports #11

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
hdrake wants to merge 6 commits into
base: topology-driven-neighbors
Choose a base branch
from
Draft

Drop the seam hand-off corner in section-finding, not in transports #11

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b0286e4
Drop the seam hand-off corner in section-finding, not in transports
hdrake Aug 10, 2026
39d9df6
Check the section-finding invariant instead of re-establishing it
hdrake Aug 10, 2026
3eaf234
Match the invariant check's tolerance comparison to the walker's
hdrake Aug 10, 2026
1523884
Collapse a repeated corner only onto an index that names the same faces
hdrake Aug 11, 2026
5c14a29
Name the fold declaration in the pole-column error
hdrake Aug 11, 2026
c4d3115
Merge topology-driven-neighbors (now carrying #10) into this branch
hdrake Aug 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions CLAUDE.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -146,9 +146,9 @@ The package is organized around a pipeline: define sections → map to grid →

### Core Modules

- **`section.py`** — Section definition and grid path algorithms. `Section` holds named waypoint coordinates; `GriddedSection` extends it with grid index information. `grid_section()` is the main entry point that maps geographic waypoints to grid vorticity-point indices `(i_c, j_c)` by walking the grid between waypoints along the `curve` requested. There are three: `"great circle"` (the default — every segment follows the geodesic), `"latitude circle"` (every segment follows a parallel, and a segment whose endpoints do not share a latitude — meridional ones included — raises), and `"latitude and great circle"` (decided per segment: constant-latitude segments follow the parallel, all others the geodesic). Under all three, every segment follows the **shortest** path between its two waypoints — raw longitudes are never read as a request to go the long way round, so a `"latitude circle"` segment written `0 -> 270` runs 90° *west* — which is why encircling the globe takes intermediate waypoints (e.g. `0 -> 120 -> 240 -> 360`). `_check_segment_span` resolves each segment's curve (via `_is_constant_latitude`, the one classification the metric choice and the legality check share) and raises for ill-posed ones: endpoints written exactly half a circle apart (neither way round is shorter — one error shared by all curves, quoting degrees of longitude along the parallel or degrees of arc, as applicable), or, under `"latitude circle"`, endpoints that do not share a latitude. The walk is deterministic and direction-independent: it admits neighbors strictly closer to the endpoint plus any seam twin of the current cell (so periodic/fold-seam crossings do not depend on floating-point rounding), and breaks ties by index. Grid topology is inferred entirely from `xgcm.Grid` metadata — each axis' `boundary` (periodic wrap, fill/extend wall, or a single-tile bipolar north fold `{"Y": {"fold": ...}}`) and `face_connections` for multi-tile grids — so there is no `topology` keyword. The pathfinder consumes topology-aware neighbor maps built by `gridutils.build_neighbor_maps`.
- **`section.py`** — Section definition and grid path algorithms. `Section` holds named waypoint coordinates; `GriddedSection` extends it with grid index information. `grid_section()` is the main entry point that maps geographic waypoints to grid vorticity-point indices `(i_c, j_c)` by walking the grid between waypoints along the `curve` requested. There are three: `"great circle"` (the default — every segment follows the geodesic), `"latitude circle"` (every segment follows a parallel, and a segment whose endpoints do not share a latitude — meridional ones included — raises), and `"latitude and great circle"` (decided per segment: constant-latitude segments follow the parallel, all others the geodesic). Under all three, every segment follows the **shortest** path between its two waypoints — raw longitudes are never read as a request to go the long way round, so a `"latitude circle"` segment written `0 -> 270` runs 90° *west* — which is why encircling the globe takes intermediate waypoints (e.g. `0 -> 120 -> 240 -> 360`). `_check_segment_span` resolves each segment's curve (via `_is_constant_latitude`, the one classification the metric choice and the legality check share) and raises for ill-posed ones: endpoints written exactly half a circle apart (neither way round is shorter — one error shared by all curves, quoting degrees of longitude along the parallel or degrees of arc, as applicable), or, under `"latitude circle"`, endpoints that do not share a latitude. The walk is deterministic and direction-independent: it admits neighbors strictly closer to the endpoint plus any seam twin of the current cell (so periodic/fold-seam crossings do not depend on floating-point rounding), and breaks ties by index. That seam-twin step visits the same physical corner twice, so it spans no grid cell and is not a velocity face; **`grid_section()` removes the redundant index right after the walk, via `drop_repeated_corners()`** (`create_section_composite()`, the lower-level entry point, does the same when passed `grid=`). Every section `grid_section()` returns therefore satisfies the invariant that **no two consecutive corners are the same physical point** — every consecutive pair is a real velocity face — and, on a multi-tile grid, that **every corner is the canonical native index of its corner node** (`_OuterTopology.node_native`). That is the only place the normalization happens: `transports.uvindices_from_qindices` *checks* the invariant and raises rather than re-establishing it, so hand-built index arrays are the caller's to fix (by calling `drop_repeated_corners` themselves). Which index of a repeated corner survives is not free: it must name the same two velocity faces, with the same transport signs, as the run it replaces — `_collapse_single_tile_runs` checks that with the very arithmetic `uvindices_from_qindices` will apply, preferring the last index (the frame the path continues in, which is what a seam hand-off needs). Where a whole grid column degenerates to one point — a bipolar cap's pole, stored as every corner of a column — a section can arrive along one index and leave along another; then **no single index names both flanking faces and it raises**, naming the corner, rather than emit a face the section does not have. (On a correctly declared fold grid the walk crosses the seam instead and this does not arise; it shows up on a tripolar grid whose fold is left undeclared.) `_check_supported_topology()` additionally rejects, at the front door, a multi-tile grid that registers no cell-center dimension on both horizontal axes: the corner-node topology is rebuilt by fingerprinting the tracer *cells* around each corner, and velocities are staggered onto the center dimensions (U at (X-corner, Y-center), V at (X-center, Y-corner)), so such a grid can carry no velocities and each of its seam corners cannot be resolved to one index. Only the center **dimensions** are needed — `_OuterTopology` numbers cells with `np.arange`, so no tracer longitude/latitude values are required. Single-tile grids are unaffected: their seam twins are detected geometrically from the corner coordinates alone, so a corner-only single-tile grid still traces *and* still enumerates its velocity faces (it just carries no velocity data to put on them). Grid topology is inferred entirely from `xgcm.Grid` metadata — each axis' `boundary` (periodic wrap, fill/extend wall, or a single-tile bipolar north fold `{"Y": {"fold": ...}}`) and `face_connections` for multi-tile grids — so there is no `topology` keyword. The pathfinder consumes topology-aware neighbor maps built by `gridutils.build_neighbor_maps`.

- **`transports.py`** — Transport computation along sections. `uvindices_from_qindices()` converts vorticity-point indices to U/V velocity-point indices using a per-position corner offset (`gridutils.corner_offset`) covering all three C-grid staggerings: 'outer', 'right', and 'left' (see "Corner staggering" below). `convergent_transport()` is the main function: it lazily computes signed normal transports with configurable orientation (positive inward to the polygon defined by the section).
- **`transports.py`** — Transport computation along sections. `uvindices_from_qindices()` converts vorticity-point indices to U/V velocity-point indices using a per-position corner offset (`gridutils.corner_offset`) covering all three C-grid staggerings: 'outer', 'right', and 'left' (see "Corner staggering" below). It emits exactly one face per consecutive pair of corners: it relies on the section-finding invariant above rather than filtering degenerate pairs itself, so an edge whose normal velocity is stored on neither the source nor the destination face raises instead of silently returning a zero-flux placeholder. Single-tile faces come from `_single_tile_face()`, the one implementation of that index arithmetic (shared with `section._collapse_single_tile_runs`); it reads a step of more than one index along an axis as a seam crossing **only on an axis that is actually periodic**, since on a walled axis the same step is a jump between two indices of one physical corner. `convergent_transport()` is the main function: it lazily computes signed normal transports with configurable orientation (positive inward to the polygon defined by the section).

- **`tracers.py`** — `extract_tracer()` interpolates tracer data to U/V points along a section path for cross-section plotting.

Expand Down
135 changes: 36 additions & 99 deletions sectionate/gridutils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -142,9 +142,32 @@ def coord_dict(grid):
dict
Dictionary containing names of "X" and "Y" dimension variables, at both cell 'center'
position and the corner position ('outer', 'right', or 'left'; see `corner_position`).

Raises
------
ValueError
If either horizontal axis registers no 'center' position. Velocities are staggered
onto the center dimensions -- U lives at (X-corner, Y-center) and V at (X-center,
Y-corner) -- so a grid without them cannot locate any velocity, and anything that
derives faces, transports or tracers from a section needs this mapping.
"""
corner_pos = corner_position(grid)

missing = [ax for ax in ("X", "Y") if "center" not in grid.axes[ax].coords]
if missing:
one = len(missing) == 1
raise ValueError(
f"The {' and '.join(missing)} {'axis' if one else 'axes'} of this grid "
f"{'registers' if one else 'register'} no 'center' position, so sectionate "
"cannot locate its velocities: U is stored at (X-corner, Y-center) and V at "
"(X-center, Y-corner). Declare the cell-center dimensions, e.g. "
"coords={'X': {'center': 'xh', 'outer': 'xq'}, "
"'Y': {'center': 'yh', 'outer': 'yq'}}. (A single-tile grid can still be "
"traced by `grid_section` from its corner coordinates alone, and its velocity "
"faces enumerated; a multi-tile grid cannot, because its corner topology is "
"rebuilt from the tracer cells around each corner.)"
)

return {
"X": {
"center": grid.axes["X"].coords["center"],
Expand Down Expand Up @@ -231,7 +254,10 @@ def build_neighbor_maps(grid, geocorners):
(shared-corner) corner topology, resolved to native indices -- see
`_OuterTopology`. Each physical corner appears under a single canonical
native index (a shared 'outer' seam corner is not stepped through twice),
and the returned maps have the same format and native index frame.
and the returned maps have the same format and native index frame. That
construction reads the tracer cells around each corner, so a multi-tile grid
must register a cell-center dimension on both horizontal axes; `coord_dict`
raises by name if it does not.

Parameters
----------
Expand All @@ -253,20 +279,14 @@ def build_neighbor_maps(grid, geocorners):
multitile = facedim is not None

if multitile:
has_centers = all("center" in grid.axes[ax].coords for ax in ("X", "Y"))
if has_centers:
return outer_topology(grid).maps
# Without tracer-center coordinates the cell-identity construction is
# unavailable (and neither are velocities, so only walking is needed):
# fall back to reading xgcm's corner halos directly. Only shared-corner
# ('outer') tilings are safe here -- staggered corner arrays can pad one
# corner off across rotated/reversed seams.
if corner_position(grid) != "outer":
raise ValueError(
"Multi-tile grids with staggered ('left'/'right') corners require "
"tracer-center coordinates to derive their corner topology."
)
return _multitile_padded_maps(grid, geocorners)
# The cell-fingerprint construction identifies each corner by the tracer cells
# around it, so it needs the cell-center *dimensions* (their names and lengths;
# no tracer longitude/latitude values are read). There is no substitute: xgcm's
# halo padding of corner arrays produces one copy of a shared seam corner per
# face -- and, on staggered ('left'/'right') lattices, can land a cross-seam
# neighbor one corner off -- so maps read off it would not step through a single
# canonical index per physical corner.
return outer_topology(grid).maps
da = geocorners["X"]
Ydim, Xdim = da.dims[-2], da.dims[-1]
ny, nx = da.sizes[Ydim], da.sizes[Xdim]
Expand Down Expand Up @@ -310,89 +330,6 @@ def pad(a):
return maps


def _multitile_padded_maps(grid, geocorners):
"""
Fallback multi-tile neighbor maps for shared-corner ('outer') tilings that
carry no tracer-center coordinates: pad index-valued corner arrays with
xgcm's `face_connections` halos and read the neighbors off directly. A
shared seam corner is stored on every face that touches it, so these maps
step through each face's own copy (seam-twin semantics); `_validate_reciprocity`
refuses topologies xgcm's padding cannot represent consistently.
"""
facedim = grid._facedim
da = geocorners["X"]
Ydim, Xdim = da.dims[-2], da.dims[-1]
ny, nx = da.sizes[Ydim], da.sizes[Xdim]
nf = da.sizes[facedim]
dims = (facedim, Ydim, Xdim)
shape = (nf, ny, nx)
iarr = xr.DataArray(np.broadcast_to(np.arange(nx), shape).astype(float), dims=dims)
jarr = xr.DataArray(np.broadcast_to(np.arange(ny)[:, None], shape).astype(float), dims=dims)
farr = xr.DataArray(np.broadcast_to(np.arange(nf)[:, None, None], shape).astype(float), dims=dims)
own_f = np.broadcast_to(np.arange(nf)[:, None, None], shape)
own_j = np.broadcast_to(np.arange(ny)[:, None], shape)
own_i = np.broadcast_to(np.arange(nx), shape)

axes = _pad_axes(grid, dims)
padding = {ax: grid.axes[ax].padding for ax in axes}
padding_width = {ax: (1, 1) for ax in axes}

def pad(a):
return _module_pad(a, grid, padding_width, padding=padding, fill_value=np.nan)

pj, pi, pf = pad(jarr), pad(iarr), pad(farr)

interior = slice(1, -1)
slices = {
"right": (interior, slice(2, None)),
"left": (interior, slice(0, -2)),
"up": (slice(2, None), interior),
"down": (slice(0, -2), interior),
}

maps = {}
for d, (ysl, xsl) in slices.items():
sel = {Ydim: ysl, Xdim: xsl}
jmap = pj.isel(sel).values
imap = pi.isel(sel).values
fmap = pf.isel(sel).values
wall = np.isnan(fmap) | np.isnan(jmap) | np.isnan(imap)
fmap = np.where(wall, own_f, fmap).astype(np.int64)
jmap = np.where(wall, own_j, jmap).astype(np.int64)
imap = np.where(wall, own_i, imap).astype(np.int64)
maps[d] = (fmap, jmap, imap)

_validate_reciprocity(maps, own_f, own_j, own_i)
return maps


def _validate_reciprocity(maps, own_f, own_j, own_i):
"""
Verify that the neighbor maps describe a consistent topology: if B is a
(non-wall) neighbor of A, then A must be one of B's four neighbors. An
inconsistent map would yield silently-wrong sections, so we detect it and
refuse rather than return garbage neighbors.
"""
for d, (fmap, jmap, imap) in maps.items():
same = (fmap == own_f) & (jmap == own_j) & (imap == own_i)
not_wall = ~same
# Gather each neighbor's own four neighbors and look for the point back.
reciprocated = np.zeros(jmap.shape, dtype=bool)
for (f2, j2, i2) in maps.values():
back = (
(f2[fmap, jmap, imap] == own_f)
& (j2[fmap, jmap, imap] == own_j)
& (i2[fmap, jmap, imap] == own_i)
)
reciprocated |= back
if np.any(not_wall & ~reciprocated):
raise NotImplementedError(
"Could not derive a consistent neighbor topology from this grid's "
"`face_connections` metadata (the multi-tile neighbor maps are not "
"reciprocal). Sections on grids with simpler topology are supported."
)


def outer_topology(grid):
"""
The complete corner-point topology of a multi-tile grid, built on its
Expand Down Expand Up @@ -1108,7 +1045,7 @@ def coincident(a, b):
node_reps[n].append((f, J, I))
self.node_reps = node_reps

# No `_validate_reciprocity` here: the maps project an *undirected* node
# The maps need no reciprocity check: they project an *undirected* node
# graph, so node-level reciprocity holds by construction. Index-level
# reciprocity deliberately does not: where one physical point is stored
# natively more than once (e.g. the LLC south-boundary fold, whose rows
Expand Down
Loading