From dfbe0b296f4b80411b2c4fa9030a51578cbf81c1 Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Mon, 10 Aug 2026 15:23:12 -0700 Subject: [PATCH 1/5] Address PR #47 review: docs wording, ECCO record, topology check Documentation and comment fixes from review 4899521952, plus two small behaviour changes that came out of it. - examples/load_example_ECCO_grid.py: point at the current Zenodo release (v3.0.0, record 21479854). The module docstring now names only the version-independent concept DOI and defers to ZENODO_RECORD_ID for the exact version, instead of restating a version DOI in prose. v3.0.0 concatenates each collection into one file per collection, so the twelve monthly volume-flux entries collapse to a single 2010 file; its published MD5 is recorded and the geometry file's checksum is unchanged. - gridutils.py: drop the "clip-to-edge" jargon in favour of plain language about closed boundaries (both occurrences), mark where the single-tile branch of build_neighbor_maps begins, define "grid cut" in the _OuterTopology docstring, and reframe that class's LLC90-specific remarks as illustrations of general statements. Measured tolerances are kept, now explicitly labelled as LLC90 measurements. Comments only. - transports.py: explain in _in_velocity_range what the check is for -- it is how _uv_for_edge decides whether a seam edge's velocity is stored on the source face or the destination face -- and handle "U" explicitly, so an unexpected component raises instead of silently taking the U branch. - section.py: _check_supported_topology no longer describes the north fold in terms of an xgcm face_connections representation that was never released. It now checks generically that the multi-tile topology is one sectionate can trace: every neighbour face exists, every connection is mutual, and no face is glued to itself. Verified that a self-glued face is not caught anywhere downstream (it silently traces a section as if the seam were absent), and that unknown/one-sided connections previously failed with a KeyError/TypeError from inside xgcm's padding. - CLAUDE.md: apply the requested edits to the project overview and the development-setup section. Tests: pin the self-glued, unknown-face and one-sided-connection failures, check that supported topologies (rotated seam, LLC90's 13-tile connections) still pass, and cover _in_velocity_range including the new raise. Co-Authored-By: Claude Opus 5 (1M context) --- CLAUDE.md | 27 ++-- examples/load_example_ECCO_grid.py | 46 +++---- sectionate/gridutils.py | 141 +++++++++++++-------- sectionate/section.py | 60 +++++++-- sectionate/tests/test_section_multitile.py | 68 +++++++++- sectionate/transports.py | 22 +++- 6 files changed, 248 insertions(+), 116 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 3d90752..5b226f3 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -4,14 +4,12 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co ## Project Overview -Sectionate is a Python package for sampling grid-consistent hydrographic sections from structured ocean model outputs (e.g. MOM6). It traces paths along C-grid velocity faces between geographic waypoints and computes transports/tracer values along those sections. Only structured grids are supported. +Sectionate is a Python package for sampling grid-consistent hydrographic sections from structured ocean model outputs. It traces paths along C-grid velocity faces between geographic waypoints and computes transports/tracer values along those sections. It supports any structured model whose grid can be described by an `xgcm.Grid` object. ## Development Setup **One conda environment per branch/worktree, named `docs_env_sectionate_`.** -Branches diverge in their dependencies — the xgcm floor in particular — so a single -shared environment silently tests the wrong versions. Create it from -`docs/environment.yml`, overriding the baked-in `name:` with `-n`: +Create it from `docs/environment.yml`, overriding the baked-in `name:` with `-n`: ```bash ENV="docs_env_sectionate_$(git rev-parse --abbrev-ref HEAD)" @@ -28,21 +26,14 @@ with `conda env remove -n "$ENV"` once the branch or worktree is gone. (`ci/environment.yml` is the minimal CI environment — pytest only, no plotting or notebook stack. It is not sufficient for the notebooks.) -The environment must carry **xgcm >= 0.10.1**, which -introduces the bipolar north fold (`padding={"Y": {"fold": ...}}`, formerly #711), -the face-connection padding fix (#712), and the bare-`DataArray` `other_component` -vector-pad fix (#749). It is a plain PyPI release, so `pip install -e .` pulls it in -automatically (sectionate requires `xgcm>=0.10.1`). The ECCO example (notebook 5) -pulls its data subset from Zenodo (concept DOI `10.5281/zenodo.21051424`) via stdlib -`urllib` with MD5 checks — no `earthaccess`/Earthdata login needed. - A correctly set up environment runs the full suite with **0 skips**. Skips mean one of -two things: fold and multi-tile tests `skip` on `xgcm < 0.10.1`, and the ECCO LLC90 and -MOM6-fold tests `skip` when their example data is absent from `data/`. In a fresh -worktree the latter is the usual cause — `data/` holds ~1.4 GB of downloaded input that -each checkout otherwise re-fetches. Symlink it from a checkout that already has it -rather than downloading again (`data/*.nc` is gitignored, but `MOM5_global_example_grid.nc` -is tracked, so link the individual files, not the directory): +two things: fold and multi-tile tests `skip` on an xgcm older than the floor in +`pyproject.toml`, and the ECCO LLC90 and MOM6-fold tests `skip` when their example data +is absent from `data/`. In a fresh worktree the latter is the usual cause — `data/` +holds ~1.4 GB of downloaded input that each checkout otherwise re-fetches. Symlink it +from a checkout that already has it rather than downloading again (`data/*.nc` is +gitignored, but `MOM5_global_example_grid.nc` is tracked, so link the individual files, +not the directory): ```bash cd data && for f in /path/to/other/checkout/data/*.nc; do diff --git a/examples/load_example_ECCO_grid.py b/examples/load_example_ECCO_grid.py index 8bf75fb..2393296 100644 --- a/examples/load_example_ECCO_grid.py +++ b/examples/load_example_ECCO_grid.py @@ -2,13 +2,15 @@ fields needed for the meridional overturning streamfunction -- as an ``xgcm.Grid`` that sectionate can consume directly. -The required files are a small subset of NASA's ECCO V4r4 state estimate -(geometry plus the twelve monthly volume-flux files for 2010), redistributed on -Zenodo (concept DOI `10.5281/zenodo.21051424`, this version `10.5281/zenodo.21051920`) -purely to make this example reproducible without a NASA Earthdata Login. Any file -not already present under ``../data/`` is downloaded from Zenodo and checked against -its published MD5; the original PO.DAAC datasets and their DOIs are listed in the -Zenodo record's README. Please cite the original NASA sources (see that README). +The required files are a small subset of NASA's ECCO V4r4 state estimate (the grid +geometry plus the 2010 monthly volume fluxes), redistributed on Zenodo purely to make +this example reproducible without a NASA Earthdata Login. Cite it by its concept DOI, +`10.5281/zenodo.21051424`, which is version-independent and always resolves to the +latest release; the particular version this module downloads is pinned by +``ZENODO_RECORD_ID`` below, which is the single place that version is recorded. Any +file not already present under ``../data/`` is downloaded from that record and checked +against its published MD5; the original PO.DAAC datasets and their DOIs are listed in +the Zenodo record's README. Please cite the original NASA sources (see that README). The native grid is MITgcm-staggered: vorticity points sit on the SW ('left') corner and the coordinates are named ``XC/YC`` (centers) and ``XG/YG`` (corners). @@ -23,30 +25,22 @@ import xarray as xr import xgcm -# Zenodo record holding the redistributed ECCO V4r4 subset. ZENODO_RECORD_ID is the -# version-specific record (v2.0.0); ZENODO_CONCEPT_DOI is the version-independent -# (citeable) DOI that always resolves to the latest version. -ZENODO_RECORD_ID = "21051920" +# Zenodo record holding the redistributed ECCO V4r4 subset. ZENODO_RECORD_ID pins the +# exact version that is downloaded, and is the only place that version is recorded; +# ZENODO_CONCEPT_DOI is the version-independent (citeable) DOI that always resolves to +# the latest version. Bumping ZENODO_RECORD_ID means re-checking ECCO_FILE_MD5 below, +# since each version publishes its own file list and checksums. +ZENODO_RECORD_ID = "21479854" ZENODO_CONCEPT_DOI = "10.5281/zenodo.21051424" ECCO_GEOMETRY_FILE = "GRID_GEOMETRY_ECCO_V4r4_native_llc0090.nc" # Published MD5 checksums (from the Zenodo record) for the files this example needs: -# the grid geometry and the twelve 2010 monthly volume-flux files. +# the grid geometry and the 2010 monthly volume fluxes (one file per collection, +# concatenated along `time`). ECCO_FILE_MD5 = { "GRID_GEOMETRY_ECCO_V4r4_native_llc0090.nc": "2663a7e86d7a0e6f7ddf84124c8376a6", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-01_ECCO_V4r4_native_llc0090.nc": "9d4371c969b2887a6ec61bd32d8f94e9", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-02_ECCO_V4r4_native_llc0090.nc": "1cad72f3f1030995e6e2b07b41c7d01a", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-03_ECCO_V4r4_native_llc0090.nc": "2bd83ca7a7b4e96b91e6179e683d88e8", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-04_ECCO_V4r4_native_llc0090.nc": "dffc775c76c3418f7fa0aa919a3854a5", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-05_ECCO_V4r4_native_llc0090.nc": "b46303fe1ae85a31b57cd5213a0928b2", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-06_ECCO_V4r4_native_llc0090.nc": "1540a87364c33614af29f97b8bca0eec", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-07_ECCO_V4r4_native_llc0090.nc": "226f7bdc68bff510189f2a8242d7193e", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-08_ECCO_V4r4_native_llc0090.nc": "fbf24a53e650a9c07052aa5369472b04", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-09_ECCO_V4r4_native_llc0090.nc": "e1a4ac8847c8948317d018445d49556e", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-10_ECCO_V4r4_native_llc0090.nc": "17845a2cfeb1f53d4ea7f97150bb9b5d", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-11_ECCO_V4r4_native_llc0090.nc": "55086ea0467dd7bd3b21e17c5afb097d", - "OCEAN_3D_VOLUME_FLUX_mon_mean_2010-12_ECCO_V4r4_native_llc0090.nc": "1487f31096153b5dfe188755aceb6553", + "OCEAN_3D_VOLUME_FLUX_mon_mean_2010_ECCO_V4r4_native_llc0090.nc": "a957bd1d133156b0ad6dd611bd39af7a", } ECCO_VOLUME_FLUX_FILES = sorted( @@ -120,8 +114,8 @@ def download_ECCO_geometry(data_dir="../data"): def download_ECCO_volume_flux(data_dir="../data"): - """Return a sorted list of the twelve 2010 monthly volume-flux files, fetching - any that are missing from Zenodo. + """Return a sorted list of the 2010 monthly volume-flux files, fetching any that + are missing from Zenodo. These files hold the native 'left'-staggered mass-weighted velocities ``UVELMASS`` (on the cell 'u'/west face, dims ``(k, tile, j, i_g)``) and diff --git a/sectionate/gridutils.py b/sectionate/gridutils.py index 9accc51..77ae492 100644 --- a/sectionate/gridutils.py +++ b/sectionate/gridutils.py @@ -181,9 +181,11 @@ def check_outer(grid): # hard-coding how to step across periodic boundaries, multi-tile face seams, or # the bipolar north fold, we precompute -- for every corner point -- the # ([face,] j, i) index of each of its four neighbors ("right", "left", "up", -# "down"). Walls (no neighbor, e.g. a "fill"/"extend" edge) are represented as -# the point itself, so the pathfinder simply never finds them closer to the -# target -- reproducing clip-to-edge. +# "down"). Where there is no neighbor -- a closed domain edge, i.e. a "fill" or +# "extend" boundary -- the point is recorded as its own neighbor. A step in that +# direction therefore stays put, so it never gets the walk closer to the target +# and is never taken: the section stops at the boundary instead of running off +# the array. # # For single-tile grids and multi-tile grids with shared ('outer') corners, the # topology logic lives upstream in `xgcm`: we pad index-valued arrays with the @@ -267,6 +269,13 @@ def build_neighbor_maps(grid, geocorners): "tracer-center coordinates to derive their corner topology." ) return _multitile_padded_maps(grid, geocorners) + + # --- single-tile grids (no face dimension) --- + # One corner array of dims (Y, X), whose only topology is each axis' own + # boundary condition: a periodic wrap, a closed ("fill"/"extend") edge, or a + # bipolar north fold. All three are what xgcm's halo padding already encodes, + # so the neighbor maps are read straight off the padded index arrays below and + # the returned maps carry `fmap = None`. da = geocorners["X"] Ydim, Xdim = da.dims[-2], da.dims[-1] ny, nx = da.sizes[Ydim], da.sizes[Xdim] @@ -300,8 +309,9 @@ def pad(a): sel = {Ydim: ysl, Xdim: xsl} jmap = pj.isel(sel).values imap = pi.isel(sel).values - # A NaN halo means "no neighbor" (an unconnected/fill edge) -- represent - # the wall as the point itself, matching the clip-to-edge behavior. + # A NaN halo means "no neighbor" (a closed, unconnected/fill edge): record + # the point as its own neighbor, so a step in that direction stays put and + # the walk stops at the boundary rather than running off the array. wall = np.isnan(jmap) | np.isnan(imap) jmap = np.where(wall, own_j, jmap).astype(np.int64) imap = np.where(wall, own_i, imap).astype(np.int64) @@ -433,14 +443,26 @@ class _OuterTopology: cubed-sphere -- a single face's corner array is incomplete along two of its edges, and a physical vorticity point may be stored **once** (the common case), **not at all** (a junction that lives on no face -- a cube vertex, or - an LLC tile junction), or **more than once** (an 'outer' seam twin; a - bit-reversed copy across a self-folded boundary; the coincident lips of a - grid cut). Any topology derived by padding the per-face corner lattice - therefore breaks down exactly where regions are hardest to trace: rotated - seams, cube-vertex / tile junctions, and the polar grid cut. This class - instead reconstructs, per face, the full `(Nyc+1, Nxc+1)` outer lattice of - corner *slots* and resolves each slot to the native corner that stores its - physical point: + a point where three tiles meet), or **more than once** (an 'outer' seam twin; + a bit-reversed copy across a self-folded boundary; the coincident lips of a + *grid cut*). + + **Grid cut** (used throughout this class): a seam along which the mesh has + been slit open, so that one line of physical corner points appears twice in + the lattice, as the two coincident "lips" of the cut. A cut differs from a + face seam in that it is *not* declared in `face_connections`: neither lip + enters the other's halo, so nothing combinatorial reveals that the two lips + are the same points -- only their coordinates do. Either lip may be stored + natively or not stored at all, and the two need not agree along the length of + the cut. For example, the ECCO LLC90 grid is slit under Antarctica along the + 65E/115W great circle, leaving coincident, reversed-index lips (tiles 0/3 on + one side, tiles 9/12 on the other). + + Any topology derived by padding the per-face corner lattice therefore breaks + down exactly where regions are hardest to trace: rotated seams, cube-vertex / + tile junctions, and grid cuts. This class instead reconstructs, per face, the + full `(Nyc+1, Nxc+1)` outer lattice of corner *slots* and resolves each slot + to the native corner that stores its physical point: 1. **Cell-identity fill (topological).** Tracer *cells* pad reliably across any seam (rotation or reversal included), so each slot is keyed by the @@ -458,20 +480,22 @@ class _OuterTopology: 2. **Coordinate fallback (only for the under-determined residue).** A few slots are left: points stored on NO face (open walls; the two coincident lips of a grid cut), and slots with only two usable cells that DO have a - native storage (the LLC90 polar-cut corners, where a face's edge wraps onto - the pole and loses a cell). Two cells cannot fix a corner topologically, so - these last are matched to their stored native corner by extrapolating the - slot's position and snapping to the nearest native corner within a fraction - of the local spacing. This is the ONLY identity inferred from coordinates, - confined to the genuinely under-determined residue. Wall / cut slots get an - extrapolated position but no identity; those that snap to nothing keep it. + native storage (on LLC90, for instance, the corners along its polar cut, + where a face's edge wraps onto the pole and loses a cell). Two cells cannot + fix a corner topologically, so these last are matched to their stored native + corner by extrapolating the slot's position and snapping to the nearest + native corner within a fraction of the local spacing. This is the ONLY + identity inferred from coordinates, confined to the genuinely + under-determined residue. Wall / cut slots get an extrapolated position but + no identity; those that snap to nothing keep it. 3. **Nodes.** Slots are merged into physical corner *nodes* by shared native identity and by exact physical coincidence (a unit-sphere position key, which also collapses the pole's degenerate longitudes). Node adjacency is read off the outer lattices: two nodes are neighbours iff some face holds index-adjacent slots for them. Every physical - velocity face is such an edge of exactly this graph, so cube-vertex - junctions are ordinary 4-valent (Arctic cap: 3-valent) nodes. + velocity face is such an edge of exactly this graph, so a junction where + several faces meet is an ordinary node of the usual valence (4 at a + cube vertex; 3 at the corners of LLC90's Arctic cap, for example). 4. **Native maps.** The graph is projected back onto native corner indices as the standard `NEIGHBOR_DIRECTIONS` maps (`build_neighbor_maps` format): each native corner's neighbour in a direction is the adjacent @@ -506,14 +530,14 @@ class _OuterTopology: # Coordinate-fallback acceptance threshold. Native identity is resolved # topologically wherever the tracer-cell fingerprint determines it; only a # small under-determined residue (slots with two usable cells that still have - # a native storage -- the LLC90 polar-cut corners) is matched by coordinates, - # by snapping the extrapolated slot to a native corner within this fraction of - # the local corner spacing. It is also the tolerance of the `loose` merge that - # unites the coincident lips of a grid cut. Extrapolation/round-off error is - # O(curvature * spacing**2), far below this; distinct corners are a full - # spacing apart, far above it -- though this margin narrows near coordinate - # singularities (the pole/cut), which is why identity is snapped only where no - # topological fingerprint exists. + # a native storage -- on LLC90, the corners along its polar cut) is matched by + # coordinates, by snapping the extrapolated slot to a native corner within this + # fraction of the local corner spacing. It is also the tolerance of the `loose` + # merge that unites the coincident lips of a grid cut. Extrapolation/round-off + # error is O(curvature * spacing**2), far below this; distinct corners are a + # full spacing apart, far above it -- though this margin narrows near + # coordinate singularities (a pole, or a cut), which is why identity is snapped + # only where no topological fingerprint exists. SNAP_FRACTION = 0.35 def __init__(self, grid): @@ -694,9 +718,9 @@ def _fill_diagonal(f_, Jp, Ip, jcell, icell, ga, gb): # fingerprint identifies the slot exactly -- keyed against every native # corner's usable-cell 3-subsets (a subset is shared by at most one native # corner, so any match is unambiguous). This resolves such junctions with - # no coordinates. A junction stored on NO face (a cube vertex; the LLC - # un-stored polar corners) matches nothing here and is left for the - # coordinate-free `by_junction` merge below. + # no coordinates. A junction stored on NO face (for example a cube's + # un-stored vertex, or LLC90's un-stored polar-cut corners) matches nothing + # here and is left for the coordinate-free `by_junction` merge below. cell3_to_native = {} for f in range(nf): for j in range(nyq): @@ -740,14 +764,15 @@ def _fill_diagonal(f_, Jp, Ip, jcell, icell, ga, gb): # and 3-cell junction match) resolve every slot whose native identity is # combinatorially determined. What can be left are: # (a) slots with only TWO usable cells that DO have a native storage -- - # chiefly the LLC90 polar-cut corners where a face's edge wraps onto - # the 65E/115W great circle right at the pole and loses a cell. Two - # cells do not fix a corner, so there is no topological fingerprint; - # these are matched to their stored native corner by snapping the - # extrapolated position to the nearest native corner within a small - # fraction of the local spacing (`SNAP_FRACTION`). This is the only - # place identity is inferred from coordinates, and only for this - # genuinely under-determined residue (30 corners on LLC90). + # corners where a face's edge runs into a grid cut or a pole and so + # loses a cell. (On LLC90 these are its polar-cut corners, where a + # face's edge wraps onto the 65E/115W great circle right at the pole: + # 30 corners in all.) Two cells do not fix a corner, so there is no + # topological fingerprint; these are matched to their stored native + # corner by snapping the extrapolated position to the nearest native + # corner within a small fraction of the local spacing + # (`SNAP_FRACTION`). This is the only place identity is inferred from + # coordinates, and only for this genuinely under-determined residue. # (b) points stored on NO face -- open-wall corners and the two coincident # lips of a grid cut -- which get an extrapolated position (no # identity) so the cut's coincident lips can later be merged by @@ -756,11 +781,12 @@ def _fill_diagonal(f_, Jp, Ip, jcell, icell, ga, gb): nat_xyz = _lonlat_to_xyz(nat_lon, nat_lat).reshape(-1, 3) nat_ok = np.isfinite(nat_xyz).all(axis=1) - # A genuine, physically distinct corner sits at least ~half a spacing away - # (measured minimum on LLC90: 0.5); a true coincidence lands within - # SNAP_FRACTION (measured maximum: 0.06 -- a ~6x margin). The band between - # is a "no man's land" that a well-posed grid never puts a two-cell fallback - # corner in: `_snap` returns the nearest native corner's index and, for a + # A genuine, physically distinct corner sits at least ~half a spacing away, + # while a true coincidence lands well within SNAP_FRACTION. (Measured on + # LLC90: the closest distinct corner is 0.5 spacings away and the farthest + # true coincidence 0.06 -- a ~6x margin.) The band between is a "no man's + # land" that a well-posed grid never puts a two-cell fallback corner in: + # `_snap` returns the nearest native corner's index and, for a # two-usable-cell slot, its distance ratio so the caller can raise if that # ratio is unexpectedly in the ambiguous band (see the raise below). DISTINCT_FRACTION = 0.5 @@ -823,8 +849,8 @@ def _snap(p, h): # instead lands in the ambiguous band -- too far to match a # native corner, too close to be a distinct one -- the # geometry is pathological and silently demoting it to a wall - # could drop real transport. Raise instead (this never fires - # on LLC90; the fallback corners snap with a ~6x margin). + # could drop real transport. Raise instead (on LLC90 this + # never fires: its fallback corners snap with a ~6x margin). if (usable4[f, J, I].sum() == 2 and self.SNAP_FRACTION <= best_fail_ratio < DISTINCT_FRACTION): raise ValueError( @@ -948,15 +974,18 @@ def coincident(a, b): by_cells3[ck] = a # Identity-less slots (points stored on no face) at the coincident lips of - # a grid cut are merged by *approximate* coincidence. The LLC90 southern - # domain is slit along the 65E/115W polar great circle: the mesh is cut - # open along that meridian pair under Antarctica, leaving coincident, - # reversed-index lips. Because the two lips are unconnected in - # `face_connections`, neither tile's cells enter the other's halo, so the - # coincidence carries NO shared tracer-cell signal -- it can only be seen - # from coordinates. + # a grid cut are merged by *approximate* coincidence. A cut is not declared + # in `face_connections` (see the class docstring), so neither lip's cells + # enter the other's halo: the coincidence carries NO shared tracer-cell + # signal and can only be seen from coordinates. A cut's two lips need not + # be uniform along its length -- part of a lip may have a native storage + # (and so real transport) while the rest is stored on no face at all -- and + # only the identity-less part is merged here; the stored part was already + # resolved above, by exact identity or by the two-cell coordinate fallback. # - # Along the 115W (tiles 9/12) side the lip is NOT uniform: + # LLC90 illustrates both. Its southern domain is slit under Antarctica + # along the 65E/115W polar great circle, leaving coincident, + # reversed-index lips. On the 115W (tiles 9/12) side: # - Equatorward (lat ~ -88 to -80), tile-9 and tile-12 corners are each # stored on no face. This step merges those coincident pairs, so the # cut's two sides read as one wall for a region tracer. These are the diff --git a/sectionate/section.py b/sectionate/section.py index b730346..2aba28b 100644 --- a/sectionate/section.py +++ b/sectionate/section.py @@ -300,27 +300,63 @@ def grid_section(grid, lons, lats, curve="great circle"): def _check_supported_topology(grid): """ - Raise if the multi-tile `grid` requires topology features sectionate does not yet support. - - A tripolar/bipolar north fold expressed as a `face_connections` self-connection (a face that - connects to itself along the "Y" axis) is not supported. Use the single-tile bipolar-fold - padding instead -- ``padding={"X": "periodic", "Y": {"fold": "corner"}}`` (xgcm >= the - bipolar-fold release) -- which sectionate handles natively via xgcm's fold padding. + Raise if the multi-tile `grid` describes a topology sectionate cannot trace a section on. + + A section is traced by walking a graph of corner points, and a step across a tile seam is + recognised by its two ends having different face indices (see `transports._uv_for_edge`). + Two things are therefore required of `face_connections`, neither of which xgcm enforces + when the grid is built: + + * Every face named as a neighbour must exist, and every connection must be mutual: the + named neighbour must name this face back, along the axis it was reached by. Otherwise + the walk either steps into a face that is not there or into one with no way back -- + today an obscure failure from deep inside xgcm's padding. + * No face may be glued to itself. Both sides of such a seam carry the same face index, so + a crossing is indistinguishable from an ordinary step within the face, and the seam's + velocity is read from the wrong side -- silently, since nothing about the indices looks + out of place. Topologies of that shape belong on the axis rather than in + `face_connections`: a zonally periodic axis as ``padding="periodic"``, a bipolar/tripolar + north fold as ``padding={"Y": {"fold": ...}}`` on a single-tile grid. sectionate handles + both of those natively. + + A grid can still prove unsupported once its corner graph is actually built -- a corner that + ends up with more than four neighbours, or a staggered grid whose seam corners are stored on + no face, for instance. Those are rejected there, each with its own reason; see + `gridutils.build_neighbor_maps` and `gridutils._OuterTopology`. """ facedim = grid._facedim + faces = set(range(grid._ds.sizes[facedim])) connections = (getattr(grid, "_face_connections", None) or {}).get(facedim, {}) for face, axis_sides in connections.items(): - for sides in axis_sides.values(): + for axis, sides in axis_sides.items(): for side in sides: if side is None: continue - neighbor_face = side[0] + neighbor_face, neighbor_axis = side[0], side[1] if neighbor_face == face: raise NotImplementedError( - "Grids that represent the tripolar/bipolar north fold as a face that " - "connects to itself (a `face_connections` self-connection) are not " - "supported. Build the grid as a single tile with the fold expressed as " - "padding instead: padding={'X':'periodic','Y':{'fold':'corner'}}." + f"Face {face} of this grid's `face_connections` is glued to itself " + f"(along its own '{axis}' axis). sectionate tells a seam crossing from " + "an ordinary step by the face index changing, so it cannot trace a face " + "joined to itself. Express that topology on the axis instead, on a " + "single-tile grid: padding='periodic' for a periodic axis, or " + "padding={'Y': {'fold': 'corner'}} for a bipolar/tripolar north fold." + ) + if neighbor_face not in faces: + raise ValueError( + f"Face {face} of this grid's `face_connections` names a neighbour face " + f"{neighbor_face} along '{axis}', but the grid's '{facedim}' dimension " + f"only has {len(faces)} faces." + ) + back = connections.get(neighbor_face, {}).get(neighbor_axis, ()) + if not any( + s is not None and s[0] == face and s[1] == axis for s in back + ): + raise ValueError( + f"This grid's `face_connections` are not mutual: face {face} connects to " + f"face {neighbor_face} along '{axis}', but face {neighbor_face} does not " + f"connect back to face {face} along its '{neighbor_axis}' axis. Every " + "seam must be declared from both sides." ) def create_section_composite( diff --git a/sectionate/tests/test_section_multitile.py b/sectionate/tests/test_section_multitile.py index e42f766..4bb926b 100644 --- a/sectionate/tests/test_section_multitile.py +++ b/sectionate/tests/test_section_multitile.py @@ -511,15 +511,79 @@ def test_rotated_seam_transport_streamfunction(): assert np.isclose(abs(conv), abs(dpsi), rtol=1e-9) -def test_north_fold_self_connection_raises(): +def test_face_glued_to_itself_raises(): + """A face connected to itself is refused up front. Both sides of such a seam carry + the same face index, so a crossing looks exactly like an ordinary step within the + face -- nothing further down the pipeline notices, and a section is silently traced + as if the seam were not there. (A periodic axis or a north fold belongs in the axis + `padding` of a single-tile grid, which sectionate does support.)""" ng = 5 lon = np.zeros((1, ng, ng)); lat = np.zeros((1, ng, ng)) fc = {"face": {0: {"Y": ((0, "Y", True), (0, "Y", True))}}} grid = _make_grid(lon, lat, fc) - with pytest.raises(NotImplementedError): + with pytest.raises(NotImplementedError, match="glued to itself"): grid_section(grid, [0., 1.], [0., 1.]) +def test_unknown_neighbour_face_raises(): + """A `face_connections` entry naming a face the grid does not have is refused with a + clear message, rather than failing obscurely inside xgcm's padding later on.""" + grid = two_face_x_to_x() + grid._face_connections["face"][0]["X"] = (None, (7, "X", False)) + with pytest.raises(ValueError, match="only has 2 faces"): + grid_section(grid, [10., 100.], [0., 0.]) + + +def test_one_sided_face_connection_raises(): + """A seam declared from only one side is refused: the walk could cross it but never + step back, so the topology is not one sectionate can trace.""" + grid = two_face_x_to_x() + grid._face_connections["face"][1]["X"] = (None, None) + with pytest.raises(ValueError, match="not mutual"): + grid_section(grid, [10., 100.], [0., 0.]) + + +def test_supported_multitile_topologies_pass_the_check(): + """The front-door check must not reject topologies sectionate really does support -- + including rotated seams (face 0's X axis glued to face 1's Y axis) and the real + 13-tile lat-lon-cap connections, whose seams are mutual across different axes.""" + import os, sys + from sectionate.section import _check_supported_topology + + _check_supported_topology(two_face_x_to_x()) + _check_supported_topology(rotated_two_face_streamfunction()) + + examples = os.path.abspath( + os.path.join(os.path.dirname(__file__), "..", "..", "examples") + ) + if examples not in sys.path: + sys.path.insert(0, examples) + from load_example_ECCO_grid import LLC90_FACE_CONNECTIONS # metadata only, no data + + ng = 4 + zeros = np.zeros((13, ng, ng)) + llc90 = _make_grid(zeros, zeros.copy(), {"face": LLC90_FACE_CONNECTIONS["tile"]}) + _check_supported_topology(llc90) + + +def test_in_velocity_range_rejects_unknown_component(): + """`_in_velocity_range` decides whether a velocity index exists in a face's own + arrays: in range for an edge stored on that face, out of range for one stored on the + face across the seam. It only ever sees the "U"/"V" components `_anchor_velocity` + produces, so anything else is a caller bug and is reported as one rather than + silently falling through to the "U" branch.""" + from sectionate.transports import _in_velocity_range + ranges = {"Xc": 4, "Yc": 4, "Xq": 5, "Yq": 5} + + assert _in_velocity_range("V", 3, 4, ranges) # (X-center, Y-corner): both in range + assert not _in_velocity_range("V", 4, 4, ranges) # past the end of the X-center axis + assert _in_velocity_range("U", 4, 3, ranges) # (X-corner, Y-center): both in range + assert not _in_velocity_range("U", 4, 4, ranges) # past the end of the Y-center axis + assert not _in_velocity_range("U", -1, 0, ranges) # off the low end + with pytest.raises(ValueError, match="'U' or 'V'"): + _in_velocity_range("0", 0, 0, ranges) + + def test_save_load_roundtrip_preserves_face_indices(tmp_path): """A multi-tile GriddedSection round-trips through save/load with its face indices intact. (f_c was previously not persisted, so a reloaded multi-tile section silently diff --git a/sectionate/transports.py b/sectionate/transports.py index 5e24a91..5ec1502 100644 --- a/sectionate/transports.py +++ b/sectionate/transports.py @@ -33,10 +33,28 @@ def _edge_direction(A, neighbor_maps): def _in_velocity_range(var, vi, vj, ranges): - """Whether a velocity index is a real point on its face (vs. an off-face artifact).""" + """ + Whether the velocity index (vi, vj) is one that a face's own arrays actually hold. + + `_anchor_velocity` builds a section edge's velocity index by offsetting the corner the + edge leaves from. Read in the frame of the face that corner belongs to, an edge that + stays inside that face always lands on a stored velocity. An edge that crosses a seam, + however, is stored on only ONE of the two faces it touches: read in the other face's + frame the same edge produces an index that falls outside that face's array -- past its + high end, or negative. `_uv_for_edge` uses this check to tell the two apart, taking the + source face's index when it is in range and otherwise the destination face's, so a + velocity is always read in the frame of the face that stores it and never has to be + rotated across the seam. + + `ranges` holds a face's array lengths along each axis: "Xc"/"Yc" for the tracer-center + axes and "Xq"/"Yq" for the corner axes. A "V" (Y-direction) velocity sits at + (X-center, Y-corner) and a "U" (X-direction) velocity at (X-corner, Y-center). + """ if var == "V": # vmo at (X-center, Y-corner) return (0 <= vi < ranges["Xc"]) and (0 <= vj < ranges["Yq"]) - return (0 <= vi < ranges["Xq"]) and (0 <= vj < ranges["Yc"]) # umo at (X-corner, Y-center) + if var == "U": # umo at (X-corner, Y-center) + return (0 <= vi < ranges["Xq"]) and (0 <= vj < ranges["Yc"]) + raise ValueError(f"velocity component must be 'U' or 'V', got {var!r}") def _anchor_velocity(d, f, j, i, offset): From 568ac7b12edc5466c7e5ccbfd7a89d83badad431 Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Mon, 10 Aug 2026 15:32:41 -0700 Subject: [PATCH 2/5] Correct two overclaiming comments surfaced by the review questions Answering the reviewer's two questions about `gridutils.py` turned up documentation that did not survive measurement. `_OuterTopology`'s corner/outer-lattice guard now says what it actually catches. It is unreachable on a well-formed grid -- corner dims are Nc+1 for 'outer' and Nc for 'left'/'right', so the test is an equality for 'outer' and 'right' and leaves one slot of slack for 'left'. It fires only when the declared corner position contradicts the corner arrays' length, or when centers and corners come from different grids; xgcm does not validate dim lengths, so such a grid builds without complaint. Its message was also wrong for the commonest trigger -- a symmetric array declared 'right' is not "larger than the outer lattice", it just cannot sit at the +1 offset 'right' implies -- so it now prints both shapes and names the likely cause. `padded_transports` claimed that no halo pad can supply a value for edges stored on no face. Measured against the pinned xgcm that is not what happens: on ECCO LLC90 it and a dict-form vector pad agree bit for bit, both summing cell convergence to exactly 0.0 globally, because the slots it zeroes are exactly the ones the pad fills with a `fill_value` that is already 0 there. The docstring now claims only what holds: independence from the axis `fill_value` (on the cubed-sphere fixture, `fill_value=nan`, a vector pad leaves 23 of 128 cells' convergence NaN where a wall's transport is 0 by definition), independence from halo correctness across rotated seams, and that it takes plain arrays. CLAUDE.md's matching claim that the dict form of `xgcm.pad` "was always exact" is softened for the same reason: xgcm#712 was a general `_pad_face_connections` bug that could corrupt any face-connection pad, dict vector form included, while #749 fixed only the bare-DataArray path. Co-Authored-By: Claude Opus 5 (1M context) --- CLAUDE.md | 2 +- sectionate/gridutils.py | 37 ++++++++++++++++++++++++++++++------- 2 files changed, 31 insertions(+), 8 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 5b226f3..f8a6208 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -119,7 +119,7 @@ The package is organized around a pipeline: define sections → map to grid → - **`tracers.py`** — `extract_tracer()` interpolates tracer data to U/V points along a section path for cross-section plotting. -- **`gridutils.py`** — Grid introspection utilities: `corner_position()` returns the shared vorticity corner position (`"outer"`, `"right"`, or `"left"`) and `corner_offset()` the corresponding velocity index shift used throughout transports/tracers; `check_outer()` is a thin wrapper (True iff `"outer"`). (Package source names positions only by these xgcm labels; their correspondence to MOM6/MITgcm/ECCO conventions is documented under "Corner staggering" below and in the example notebooks.) `coord_dict()` and `get_geo_corners()` extract coordinate/dimension names from `xgcm.Grid` metadata; `build_neighbor_maps()` builds the topology-aware neighbor maps the pathfinder consumes. For single-tile grids these come from xgcm halo padding; for every multi-tile grid they are projected from **`outer_topology(grid)`** (cached `_OuterTopology`), which reconstructs each face's full 'outer' (shared-corner) lattice, resolves every extended corner slot to the native corner storing that physical point (topological matching of the surrounding tracer cells — including a diagonal-cell recovery for normal-seam face corners and a 3-cell fingerprint for 3-tile junctions with native storage — with a small coordinate fallback only for the polar-cut corners no cell fingerprint can determine), and merges slots into physical corner *nodes* — so rotated/reversed seams, cube-vertex/tile junctions, boundary folds (e.g. LLC's south rows), and the polar grid cut are ordinary graph nodes, each resolved to its canonical native `(face, j, i)` (a physical corner may be stored once, on no face, or more than once). `_OuterTopology.padded_transports(u, v)` extends native transports to the outer lattice by reading each missing edge slot's *stored* twin through the node graph — a topology-exact, xgcm-independent alternative to `xgcm.pad(..., other_component=...)` (cell convergence built from `padded_transports` sums to exactly zero globally on ECCO LLC90). Even a correct vector pad (xgcm#749 fixes the bare-DataArray path; the dict form was always exact) cannot supply a halo for edges stored on *no* face, so `padded_transports` remains the exact path on grids with cuts/caps. Points stored on **no** face (a cubed-sphere stores `6N²` corners for `6N²+2` points; LLC90's 4th Arctic-cap vertex; the unstored lip of its 65°E/115°W polar grid cut under Antarctica) become edge-less nodes: sections and traced boundaries cannot pass through them and raise an error rather than return invalid indices. +- **`gridutils.py`** — Grid introspection utilities: `corner_position()` returns the shared vorticity corner position (`"outer"`, `"right"`, or `"left"`) and `corner_offset()` the corresponding velocity index shift used throughout transports/tracers; `check_outer()` is a thin wrapper (True iff `"outer"`). (Package source names positions only by these xgcm labels; their correspondence to MOM6/MITgcm/ECCO conventions is documented under "Corner staggering" below and in the example notebooks.) `coord_dict()` and `get_geo_corners()` extract coordinate/dimension names from `xgcm.Grid` metadata; `build_neighbor_maps()` builds the topology-aware neighbor maps the pathfinder consumes. For single-tile grids these come from xgcm halo padding; for every multi-tile grid they are projected from **`outer_topology(grid)`** (cached `_OuterTopology`), which reconstructs each face's full 'outer' (shared-corner) lattice, resolves every extended corner slot to the native corner storing that physical point (topological matching of the surrounding tracer cells — including a diagonal-cell recovery for normal-seam face corners and a 3-cell fingerprint for 3-tile junctions with native storage — with a small coordinate fallback only for the polar-cut corners no cell fingerprint can determine), and merges slots into physical corner *nodes* — so rotated/reversed seams, cube-vertex/tile junctions, boundary folds (e.g. LLC's south rows), and the polar grid cut are ordinary graph nodes, each resolved to its canonical native `(face, j, i)` (a physical corner may be stored once, on no face, or more than once). `_OuterTopology.padded_transports(u, v)` extends native transports to the outer lattice by reading each missing edge slot's *stored* twin through the node graph — a topology-exact, xgcm-independent alternative to `xgcm.pad(..., other_component=...)` (cell convergence built from `padded_transports` sums to exactly zero globally on ECCO LLC90, where it agrees bit-for-bit with the dict-form vector pad against the pinned xgcm). What it does not depend on is the face-connection halo being right (xgcm#712 was a general `_pad_face_connections` bug that could corrupt any face-connection pad, dict vector form included; the separate #749 fixed only the bare-`DataArray` dispatch), and it does not depend on the axis `fill_value`: edges stored on *no* face resolve to zero, a wall's true transport, rather than to whatever `fill_value` the axis declares — which is what makes it exact on grids whose `fill_value` is not 0 (on the cubed-sphere fixture, `fill_value=np.nan`, a vector pad leaves 23 of 128 cells' convergence NaN). Points stored on **no** face (a cubed-sphere stores `6N²` corners for `6N²+2` points; LLC90's 4th Arctic-cap vertex; the unstored lip of its 65°E/115°W polar grid cut under Antarctica) become edge-less nodes: sections and traced boundaries cannot pass through them and raise an error rather than return invalid indices. - **`utils.py`** — Section catalog I/O. Sections can be loaded by name from JSON catalog files in `sectionate/catalog/`. Also provides `save_gridded_section()`/`load_gridded_section()` for persisting gridded sections. diff --git a/sectionate/gridutils.py b/sectionate/gridutils.py index 77ae492..c7476f9 100644 --- a/sectionate/gridutils.py +++ b/sectionate/gridutils.py @@ -574,10 +574,25 @@ def __init__(self, grid): # native corner (j, i) lives at outer slot (j+t, i+t) t = 1 if pos == "right" else 0 nqy, nqx = Nyc + 1, Nxc + 1 + # Unreachable on a well-formed grid: corner dims are Nc+1 for 'outer' and Nc + # for 'left'/'right', so this is an equality for 'outer' and 'right' (no slack + # at all) and leaves exactly one slot of slack for 'left'. What it really + # catches is a grid whose declared corner *position* contradicts the length of + # its corner arrays -- e.g. a symmetric (N+1)-length corner array declared + # 'right', which xgcm builds without complaint (`Axis.__init__` validates + # position names and duplicate dims, never dim lengths) -- or centers and + # corners taken from different grids. Worth keeping: without it the failure + # surfaces ~100 lines below as an opaque NumPy broadcast error out of the + # `ident[f, t:t+nyq, t:t+nxq, 1] = jn` fill. It is one-sided, though: corner + # arrays that are too *small* for the declared position (a 'left'-sized array + # declared 'outer') fit inside the lattice and are not caught here. if (t + nyq > nqy) or (t + nxq > nqx): raise ValueError( - "Corner arrays are larger than the grid's outer lattice; the " - "corner/center coordinates are inconsistent." + f"Corner dims ({nyq}, {nxq}) do not fit the outer lattice " + f"({nqy}, {nqx}) at the offset implied by '{pos}' staggering. Check " + "that the corner coordinates are declared at the right xgcm position " + "(a symmetric, (N+1)-length corner array is 'outer', not " + "'right'/'left') and that centers and corners come from the same grid." ) self.facedim, self.nf, self.t = facedim, nf, t self.nqy, self.nqx = nqy, nqx @@ -1188,13 +1203,21 @@ def padded_transports(self, u, v): each missing edge slot's value from the native storage of that physical edge (with the sign of the receiving face's own axis direction). - This is the topology-exact, xgcm-independent analogue of + This is a topology-exact, xgcm-independent analogue of ``xgcm.pad(..., other_component=...)``: every added slot is resolved through the corner-node graph to the *stored* velocity of the same - physical face. Unlike a halo pad it needs no vector rotation and, more - importantly, resolves edges that are stored on *no* face (open walls, - grid cuts, a cap's un-stored vertex) to zero -- a wall carries no - transport -- which no halo pad can supply a value for. + physical face. That buys three things over a halo pad -- it does not + depend on the halo being right across a rotated or reversed seam (cf. + xgcm#712, #749), it needs no vector rotation, and it takes plain arrays. + It is also independent of the axis `fill_value`: an edge stored on *no* + face (an open wall, a grid cut, a cap's un-stored vertex) is resolved to + zero, because a wall carries no transport, where a halo pad can only + insert whatever `fill_value` the axis declares. Where that `fill_value` + is already 0 the two agree exactly -- on ECCO LLC90 this and a dict-form + vector pad match bit for bit, both summing cell convergence to exactly + 0.0 globally -- but where it is not they diverge: on the cubed-sphere + test fixture, whose `fill_value` is NaN, a vector pad leaves 23 of 128 + cells' convergence NaN where this returns the correct 0. Parameters ---------- From c566d7da5080f6a1a55e920b6f64720211544757 Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Tue, 11 Aug 2026 10:31:48 -0700 Subject: [PATCH 3/5] Correct false claims in the docs added by the review follow-up An audit of this branch found several statements that are not true of the code, plus a topology check that rejects a legal grid. Fix all of them. `_check_supported_topology` claimed that face existence and connection mutuality are things "neither of which xgcm enforces when the grid is built". xgcm enforces both, at `Grid.__init__`: naming a face the grid does not have raises KeyError, and a one-sided connection raises TypeError, each while it builds its face-connection table. Only the self-gluing rejection is load bearing -- xgcm accepts a face glued to itself and nothing downstream notices, so a section is traced as though the seam were absent. The docstring now says that, and labels the other two branches as the defence in depth they are. The two tests that fired them by corrupting `grid._face_connections` after construction are replaced by one that asserts xgcm rejects those grids at construction, which is where the real enforcement lives. The same function built its set of known faces as `range(n_faces)`, so a grid labelled `face = [1, 2]` -- legal, and accepted by xgcm -- was rejected with "only has 2 faces" even though face 2 exists. Read the labels instead, and report them in the message. Covered by a new test. `padded_transports` claimed to be independent of "the halo being right across a rotated or reversed seam". It is not: `_OuterTopology` builds its whole node graph by padding tracer-cell ids with the grid's own `face_connections` via `xgcm.padding.pad`, and `padded_transports` reads that graph. What it really buys over a *vector* pad is that it needs no vector rotation, takes plain arrays, and is independent of the axis `fill_value`. The list also promised three things and gave four. Separately, a zeroed cut lip was justified as "a wall carries no transport" -- true for a wall, but for a cut lip or an un-stored cap vertex the edge exists and zero is the conservative choice, not a truth. Say so. The corner-array guard rendered as "Corner dims (9, 9) do not fit the outer lattice (9, 9)" for the very case it was written for -- two identical shapes and "do not fit" -- and leaked "outer lattice", a term that appears in no user-facing doc and is defined only inside the private `_OuterTopology`. Rewritten in terms of the corner and centre array shapes and the offset the declared staggering implies. Its twelve-line comment above a two-line `if` is trimmed. `_in_velocity_range` said a seam edge "is stored on only ONE of the two faces it touches"; `_uv_for_edge`'s third branch exists precisely for edges stored on neither. Corrected to "at most one", with that case named. The LLC90 generalisation begun on this branch stopped at the class body: the `gridutils.py` module header and the CLAUDE.md bullet still stated polar-cut specifics as general truth. Generalised, but not over-generalised -- a corner loses a cell at a cut that reaches a pole, not at "a grid cut or a pole", since a bipolar fold alone loses none. Per project policy, no source file references a merged PR or a closed issue outside a TODO; the four such references introduced here are replaced by descriptions of the behaviour or bug. (`test_section_fold.py` is left alone; its references are handled elsewhere.) Nits: define "fingerprint" and "seam twin" at first use, "bit-reversed" -> "index-reversed", trim the ten-line grid-cut definition, and fix the leftovers in `examples/load_example_ECCO_grid.py` from the twelve-files-to-one change -- including `load_ECCO_MOC_grid` claiming its transports are "time-averaged over the twelve months of 2010" when `utr`/`vtr` keep a length-12 `time` dimension. Full suite: 81 passed, 0 skipped. Co-Authored-By: Claude Opus 5 (1M context) --- CLAUDE.md | 2 +- examples/load_example_ECCO_grid.py | 15 +- sectionate/gridutils.py | 178 ++++++++++++--------- sectionate/section.py | 38 +++-- sectionate/tests/test_cube_left_grid.py | 16 +- sectionate/tests/test_ecco_llc90.py | 9 +- sectionate/tests/test_section_multitile.py | 63 ++++++-- sectionate/transports.py | 7 +- 8 files changed, 199 insertions(+), 129 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index f8a6208..4f59a70 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -119,7 +119,7 @@ The package is organized around a pipeline: define sections → map to grid → - **`tracers.py`** — `extract_tracer()` interpolates tracer data to U/V points along a section path for cross-section plotting. -- **`gridutils.py`** — Grid introspection utilities: `corner_position()` returns the shared vorticity corner position (`"outer"`, `"right"`, or `"left"`) and `corner_offset()` the corresponding velocity index shift used throughout transports/tracers; `check_outer()` is a thin wrapper (True iff `"outer"`). (Package source names positions only by these xgcm labels; their correspondence to MOM6/MITgcm/ECCO conventions is documented under "Corner staggering" below and in the example notebooks.) `coord_dict()` and `get_geo_corners()` extract coordinate/dimension names from `xgcm.Grid` metadata; `build_neighbor_maps()` builds the topology-aware neighbor maps the pathfinder consumes. For single-tile grids these come from xgcm halo padding; for every multi-tile grid they are projected from **`outer_topology(grid)`** (cached `_OuterTopology`), which reconstructs each face's full 'outer' (shared-corner) lattice, resolves every extended corner slot to the native corner storing that physical point (topological matching of the surrounding tracer cells — including a diagonal-cell recovery for normal-seam face corners and a 3-cell fingerprint for 3-tile junctions with native storage — with a small coordinate fallback only for the polar-cut corners no cell fingerprint can determine), and merges slots into physical corner *nodes* — so rotated/reversed seams, cube-vertex/tile junctions, boundary folds (e.g. LLC's south rows), and the polar grid cut are ordinary graph nodes, each resolved to its canonical native `(face, j, i)` (a physical corner may be stored once, on no face, or more than once). `_OuterTopology.padded_transports(u, v)` extends native transports to the outer lattice by reading each missing edge slot's *stored* twin through the node graph — a topology-exact, xgcm-independent alternative to `xgcm.pad(..., other_component=...)` (cell convergence built from `padded_transports` sums to exactly zero globally on ECCO LLC90, where it agrees bit-for-bit with the dict-form vector pad against the pinned xgcm). What it does not depend on is the face-connection halo being right (xgcm#712 was a general `_pad_face_connections` bug that could corrupt any face-connection pad, dict vector form included; the separate #749 fixed only the bare-`DataArray` dispatch), and it does not depend on the axis `fill_value`: edges stored on *no* face resolve to zero, a wall's true transport, rather than to whatever `fill_value` the axis declares — which is what makes it exact on grids whose `fill_value` is not 0 (on the cubed-sphere fixture, `fill_value=np.nan`, a vector pad leaves 23 of 128 cells' convergence NaN). Points stored on **no** face (a cubed-sphere stores `6N²` corners for `6N²+2` points; LLC90's 4th Arctic-cap vertex; the unstored lip of its 65°E/115°W polar grid cut under Antarctica) become edge-less nodes: sections and traced boundaries cannot pass through them and raise an error rather than return invalid indices. +- **`gridutils.py`** — Grid introspection utilities: `corner_position()` returns the shared vorticity corner position (`"outer"`, `"right"`, or `"left"`) and `corner_offset()` the corresponding velocity index shift used throughout transports/tracers; `check_outer()` is a thin wrapper (True iff `"outer"`). (Package source names positions only by these xgcm labels; their correspondence to MOM6/MITgcm/ECCO conventions is documented under "Corner staggering" below and in the example notebooks.) `coord_dict()` and `get_geo_corners()` extract coordinate/dimension names from `xgcm.Grid` metadata; `build_neighbor_maps()` builds the topology-aware neighbor maps the pathfinder consumes. For single-tile grids these come from xgcm halo padding; for every multi-tile grid they are projected from **`outer_topology(grid)`** (cached `_OuterTopology`), which reconstructs each face's full 'outer' (shared-corner) lattice, resolves every extended corner slot to the native corner storing that physical point (topological matching of the surrounding tracer cells — a corner's *fingerprint* — including a diagonal-cell recovery for normal-seam face corners and a 3-cell fingerprint for 3-tile junctions with native storage, with a small coordinate fallback only where too few of those cells survive to identify a corner at all, as where a face's edge runs into a grid cut that reaches a pole), and merges slots into physical corner *nodes* — so rotated/reversed seams, cube-vertex/tile junctions, boundary folds (e.g. LLC's south rows), and grid cuts are ordinary graph nodes, each resolved to its canonical native `(face, j, i)` (a physical corner may be stored once, on no face, or more than once). `_OuterTopology.padded_transports(u, v)` extends native transports to the outer lattice by reading each missing edge slot's *stored* twin through the node graph — a topology-exact alternative to `xgcm.pad(..., other_component=...)` (cell convergence built from `padded_transports` sums to exactly zero globally on ECCO LLC90, where it agrees bit-for-bit with the dict-form vector pad against the pinned xgcm). What it buys over a **vector** pad is three things: it needs no vector rotation (every slot is filled with a value already stored in the frame it is read in, so nothing can be wrong about how a rotated or reversed seam maps `u` onto `v`); it takes plain arrays instead of an `{axis: component}` mapping; and it is independent of the axis `fill_value` — an edge stored on *no* face resolves to 0.0, where a halo pad can only insert the declared `fill_value` (on the cubed-sphere fixture, `fill_value=np.nan`, a vector pad leaves 23 of 128 cells' convergence NaN). For an open wall 0.0 is the truth; for a cut lip or an un-stored cap vertex the edge exists but is stored nowhere, so 0.0 is the conservative choice rather than a truth. It is *not* independent of xgcm's **scalar** padding: `_OuterTopology` builds the node graph in the first place by padding tracer-cell ids with the grid's own `face_connections` via `xgcm.padding.pad`. Points stored on **no** face (a cubed-sphere stores `6N²` corners for `6N²+2` points; LLC90's 4th Arctic-cap vertex; the unstored lip of its 65°E/115°W polar grid cut under Antarctica) become edge-less nodes: sections and traced boundaries cannot pass through them and raise an error rather than return invalid indices. - **`utils.py`** — Section catalog I/O. Sections can be loaded by name from JSON catalog files in `sectionate/catalog/`. Also provides `save_gridded_section()`/`load_gridded_section()` for persisting gridded sections. diff --git a/examples/load_example_ECCO_grid.py b/examples/load_example_ECCO_grid.py index 2393296..0ee4fdc 100644 --- a/examples/load_example_ECCO_grid.py +++ b/examples/load_example_ECCO_grid.py @@ -36,8 +36,8 @@ ECCO_GEOMETRY_FILE = "GRID_GEOMETRY_ECCO_V4r4_native_llc0090.nc" # Published MD5 checksums (from the Zenodo record) for the files this example needs: -# the grid geometry and the 2010 monthly volume fluxes (one file per collection, -# concatenated along `time`). +# the grid geometry, and a single file holding all twelve 2010 monthly volume fluxes +# along a `time` dimension. ECCO_FILE_MD5 = { "GRID_GEOMETRY_ECCO_V4r4_native_llc0090.nc": "2663a7e86d7a0e6f7ddf84124c8376a6", "OCEAN_3D_VOLUME_FLUX_mon_mean_2010_ECCO_V4r4_native_llc0090.nc": "a957bd1d133156b0ad6dd611bd39af7a", @@ -114,8 +114,8 @@ def download_ECCO_geometry(data_dir="../data"): def download_ECCO_volume_flux(data_dir="../data"): - """Return a sorted list of the 2010 monthly volume-flux files, fetching any that - are missing from Zenodo. + """Return the local paths of the 2010 monthly volume-flux data (currently a single + file covering all twelve months), fetching from Zenodo any that are missing. These files hold the native 'left'-staggered mass-weighted velocities ``UVELMASS`` (on the cell 'u'/west face, dims ``(k, tile, j, i_g)``) and @@ -181,8 +181,9 @@ def load_ECCO_MOC_grid(data_dir="../data"): Same native ('left') multi-tile grid as ``load_ECCO_LLC90_grid``, but the dataset additionally carries ``utr``/``vtr`` -- the volume transports (m^3/s) across the U (west) and V (south) cell faces, with vertical dimension ``k`` - (depth coordinate ``Z``), time-averaged over the twelve months of 2010 -- ready - to pass to ``sectionate.transports.convergent_transport`` as - ``utr="utr", vtr="vtr"``. + (depth coordinate ``Z``) and a length-12 ``time`` dimension holding the twelve + 2010 monthly means -- ready to pass to + ``sectionate.transports.convergent_transport`` as ``utr="utr", vtr="vtr"`` + (average over ``time`` first for an annual-mean overturning). """ return _ecco_grid(_ecco_dataset(data_dir=data_dir, with_transports=True)) diff --git a/sectionate/gridutils.py b/sectionate/gridutils.py index c7476f9..4dcd0d5 100644 --- a/sectionate/gridutils.py +++ b/sectionate/gridutils.py @@ -201,10 +201,10 @@ def check_outer(grid): # lattice: each face is extended by its missing corner row+column, and every # extended slot is resolved to the native corner that stores that physical point # by matching the surrounding tracer cells, which pad reliably (with a small -# coordinate fallback only for the polar-cut corners no cell fingerprint can -# determine). The resulting global corner graph -- on which face seams, -# cube-vertex / tile junctions, and even the polar grid cut are ordinary nodes -- -# is projected back onto native ([face,] j, i) indices. See `_OuterTopology`. +# coordinate fallback only where too few of those cells survive to identify a +# corner at all). The resulting global corner graph -- on which face seams, tile +# junctions and grid cuts are all ordinary nodes -- is projected back onto native +# ([face,] j, i) indices. See `_OuterTopology`. # --------------------------------------------------------------------------- NEIGHBOR_DIRECTIONS = ("right", "left", "up", "down") @@ -229,11 +229,13 @@ def build_neighbor_maps(grid, geocorners): Multi-tile grids do not go through xgcm's halo padding at all: on staggered ('left'/'right') corner lattices padding can land one corner off across rotated/reversed seams, and on shared-corner ('outer') tilings it produces - seam-twin duplicates. Their maps are instead derived from the grid's outer - (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. + *seam twins* -- one physical corner stored once on each of the two faces that + share a seam, so it carries two different native indices. Their maps are + instead derived from the grid's outer (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. Parameters ---------- @@ -443,20 +445,18 @@ class _OuterTopology: cubed-sphere -- a single face's corner array is incomplete along two of its edges, and a physical vorticity point may be stored **once** (the common case), **not at all** (a junction that lives on no face -- a cube vertex, or - a point where three tiles meet), or **more than once** (an 'outer' seam twin; - a bit-reversed copy across a self-folded boundary; the coincident lips of a - *grid cut*). - - **Grid cut** (used throughout this class): a seam along which the mesh has - been slit open, so that one line of physical corner points appears twice in - the lattice, as the two coincident "lips" of the cut. A cut differs from a - face seam in that it is *not* declared in `face_connections`: neither lip - enters the other's halo, so nothing combinatorial reveals that the two lips - are the same points -- only their coordinates do. Either lip may be stored - natively or not stored at all, and the two need not agree along the length of - the cut. For example, the ECCO LLC90 grid is slit under Antarctica along the - 65E/115W great circle, leaving coincident, reversed-index lips (tiles 0/3 on - one side, tiles 9/12 on the other). + a point where three tiles meet), or **more than once** (a *seam twin* -- the + same physical corner stored once on each of two faces sharing an 'outer' + seam; an index-reversed copy across a self-folded boundary; the coincident + lips of a *grid cut*). + + **Grid cut**: a seam along which the mesh has been slit open, so one line of + physical corners appears twice, as the cut's two coincident "lips". Unlike a + face seam a cut is not declared in `face_connections`, so neither lip enters + the other's halo and only their coordinates reveal that they are the same + points. Either lip may be stored natively or not at all, and the two need not + agree along the cut's length. (ECCO LLC90 is slit under Antarctica along the + 65E/115W great circle, leaving reversed-index lips: tiles 0/3 against 9/12.) Any topology derived by padding the per-face corner lattice therefore breaks down exactly where regions are hardest to trace: rotated seams, cube-vertex / @@ -466,28 +466,32 @@ class _OuterTopology: 1. **Cell-identity fill (topological).** Tracer *cells* pad reliably across any seam (rotation or reversal included), so each slot is keyed by the - up-to-four global cell ids around it; a slot on a face seam has the same - cells as its native twin on the neighbouring face, giving an exact, purely - combinatorial identification (no coordinate tolerance, no assumptions - about xgcm's corner-halo alignment). The four *diagonally* padded corner - cells are pads of pads (unreliable across two seams), so they are first - **recovered topologically** as the unique common edge-neighbour of the two - reliably-padded cells flanking each face corner -- giving normal-seam face - corners a full four-cell fingerprint. A genuine 3-face junction has no such - diagonal but still meets three cells at a unique point, so if it has a - native storage it is matched by that three-cell fingerprint (a 3-cell key - is shared by at most one native corner). All of this is coordinate-free. + up-to-four global cell ids around it -- its *fingerprint*, the set of + tracer cells that meet at that corner. A fingerprint is a property of the + physical point, not of any face's indexing, so a slot on a face seam has + the same one as its native twin on the neighbouring face, giving an exact, + purely combinatorial identification (no coordinate tolerance, no + assumptions about xgcm's corner-halo alignment). The four *diagonally* + padded corner cells are pads of pads (unreliable across two seams), so + they are first **recovered topologically** as the unique common + edge-neighbour of the two reliably-padded cells flanking each face corner + -- giving normal-seam face corners a full four-cell fingerprint. A genuine + 3-face junction has no such diagonal but still meets three cells at a + unique point, so if it has a native storage it is matched by that + three-cell fingerprint (a 3-cell key is shared by at most one native + corner). All of this is coordinate-free. 2. **Coordinate fallback (only for the under-determined residue).** A few slots are left: points stored on NO face (open walls; the two coincident lips of a grid cut), and slots with only two usable cells that DO have a - native storage (on LLC90, for instance, the corners along its polar cut, - where a face's edge wraps onto the pole and loses a cell). Two cells cannot - fix a corner topologically, so these last are matched to their stored native - corner by extrapolating the slot's position and snapping to the nearest - native corner within a fraction of the local spacing. This is the ONLY - identity inferred from coordinates, confined to the genuinely - under-determined residue. Wall / cut slots get an extrapolated position but - no identity; those that snap to nothing keep it. + native storage -- corners where the surrounding cells are themselves + incomplete, as where a face's edge runs into a cut that reaches a pole (on + LLC90, the corners along its polar cut). Two cells cannot fix a corner + topologically, so these last are matched to their stored native corner by + extrapolating the slot's position and snapping to the nearest native + corner within a fraction of the local spacing. This is the ONLY identity + inferred from coordinates, confined to the genuinely under-determined + residue. Wall / cut slots get an extrapolated position but no identity; + those that snap to nothing keep it. 3. **Nodes.** Slots are merged into physical corner *nodes* by shared native identity and by exact physical coincidence (a unit-sphere position key, which also collapses the pole's degenerate longitudes). @@ -574,25 +578,24 @@ def __init__(self, grid): # native corner (j, i) lives at outer slot (j+t, i+t) t = 1 if pos == "right" else 0 nqy, nqx = Nyc + 1, Nxc + 1 - # Unreachable on a well-formed grid: corner dims are Nc+1 for 'outer' and Nc - # for 'left'/'right', so this is an equality for 'outer' and 'right' (no slack - # at all) and leaves exactly one slot of slack for 'left'. What it really - # catches is a grid whose declared corner *position* contradicts the length of - # its corner arrays -- e.g. a symmetric (N+1)-length corner array declared - # 'right', which xgcm builds without complaint (`Axis.__init__` validates - # position names and duplicate dims, never dim lengths) -- or centers and - # corners taken from different grids. Worth keeping: without it the failure - # surfaces ~100 lines below as an opaque NumPy broadcast error out of the - # `ident[f, t:t+nyq, t:t+nxq, 1] = jn` fill. It is one-sided, though: corner - # arrays that are too *small* for the declared position (a 'left'-sized array - # declared 'outer') fit inside the lattice and are not caught here. + # Catches a grid whose declared corner *position* contradicts the length of its + # corner arrays -- e.g. a symmetric (N+1)-length corner array declared 'right', + # which xgcm builds without complaint (`Axis.__init__` validates position names + # and duplicate dims, never dim lengths) -- or centers and corners taken from + # different grids. Without it the failure surfaces ~100 lines below as an opaque + # NumPy broadcast error. It is one-sided: corner arrays too *small* for the + # declared position (a 'left'-sized array declared 'outer') still fit and pass. if (t + nyq > nqy) or (t + nxq > nqx): raise ValueError( - f"Corner dims ({nyq}, {nxq}) do not fit the outer lattice " - f"({nqy}, {nqx}) at the offset implied by '{pos}' staggering. Check " - "that the corner coordinates are declared at the right xgcm position " - "(a symmetric, (N+1)-length corner array is 'outer', not " - "'right'/'left') and that centers and corners come from the same grid." + f"This grid's corner coordinates are {nyq} x {nxq}, too large to sit at " + f"the index offset that '{pos}' corner staggering implies. Its {Nyc} x " + f"{Nxc} tracer cells have {nqy} x {nqx} corners in all, and at '{pos}' " + f"the stored corner arrays sit at index offset {t} within that grid of " + f"corners, so they can be at most {nqy - t} x {nqx - t}. " + "Check that the corner coordinates are declared at the xgcm position " + "they are really stored at -- a symmetric grid, whose corner arrays are " + "one longer than its centers along each axis, is 'outer', not 'right' " + "or 'left' -- and that the centers and corners come from the same grid." ) self.facedim, self.nf, self.t = facedim, nf, t self.nqy, self.nqx = nqy, nqx @@ -779,10 +782,11 @@ def _fill_diagonal(f_, Jp, Ip, jcell, icell, ga, gb): # and 3-cell junction match) resolve every slot whose native identity is # combinatorially determined. What can be left are: # (a) slots with only TWO usable cells that DO have a native storage -- - # corners where a face's edge runs into a grid cut or a pole and so - # loses a cell. (On LLC90 these are its polar-cut corners, where a - # face's edge wraps onto the 65E/115W great circle right at the pole: - # 30 corners in all.) Two cells do not fix a corner, so there is no + # corners where a face's edge runs into a cut that reaches a pole, so + # that the cells which would complete the fingerprint are themselves + # missing. (On LLC90 these are its polar-cut corners, where a face's + # edge wraps onto the 65E/115W great circle right at the pole: 30 + # corners in all.) Two cells do not fix a corner, so there is no # topological fingerprint; these are matched to their stored native # corner by snapping the extrapolated position to the nearest native # corner within a small fraction of the local spacing @@ -1203,21 +1207,39 @@ def padded_transports(self, u, v): each missing edge slot's value from the native storage of that physical edge (with the sign of the receiving face's own axis direction). - This is a topology-exact, xgcm-independent analogue of - ``xgcm.pad(..., other_component=...)``: every added slot is resolved - through the corner-node graph to the *stored* velocity of the same - physical face. That buys three things over a halo pad -- it does not - depend on the halo being right across a rotated or reversed seam (cf. - xgcm#712, #749), it needs no vector rotation, and it takes plain arrays. - It is also independent of the axis `fill_value`: an edge stored on *no* - face (an open wall, a grid cut, a cap's un-stored vertex) is resolved to - zero, because a wall carries no transport, where a halo pad can only - insert whatever `fill_value` the axis declares. Where that `fill_value` - is already 0 the two agree exactly -- on ECCO LLC90 this and a dict-form - vector pad match bit for bit, both summing cell convergence to exactly - 0.0 globally -- but where it is not they diverge: on the cubed-sphere - test fixture, whose `fill_value` is NaN, a vector pad leaves 23 of 128 - cells' convergence NaN where this returns the correct 0. + This is a topology-exact analogue of ``xgcm.pad(..., other_component=...)``: + every added slot is resolved through the corner-node graph to the *stored* + velocity of the same physical face. That buys three things over a vector + halo pad: + + * It needs no vector rotation. Each slot is filled with a value already + stored in the frame it is read in, so nothing here can be wrong about + how a rotated or reversed seam maps `u` onto `v`. + * It takes plain arrays, rather than the ``{axis: component}`` mapping a + vector pad needs to know which component is which. + * It is independent of the axis `fill_value`. An edge stored on *no* + face -- an open wall, a lip of a grid cut, a cap's un-stored vertex -- + is resolved to 0.0, where a halo pad can only insert whatever + `fill_value` the axis declares. Where that `fill_value` is already 0 + the two agree exactly (on ECCO LLC90 this and a dict-form vector pad + match bit for bit, both summing cell convergence to exactly 0.0 + globally); where it is not they diverge -- on the cubed-sphere test + fixture, whose `fill_value` is NaN, a vector pad leaves 23 of 128 + cells' convergence NaN where this returns 0. + + For an open wall 0.0 is the true transport. For a cut lip or an un-stored + cap vertex the physical edge exists but nothing in the dataset stores its + value, so 0.0 is the conservative choice rather than a truth: it keeps + cell convergence exactly conservative and confines the unknown flux to + those edges, instead of letting a `fill_value` spread through the sum. + + This is not independent of xgcm's *scalar* padding: `_OuterTopology` + builds the node graph in the first place by padding tracer-cell ids with + the grid's own `face_connections` (`xgcm.padding.pad`), and this method + reads that graph. What it removes is the dependence on the vector pad -- + the part that has to rotate and re-sign components across a seam -- and + on the pad supplying a halo at all for edges no face stores, which it + cannot do. Parameters ---------- diff --git a/sectionate/section.py b/sectionate/section.py index 2aba28b..9692931 100644 --- a/sectionate/section.py +++ b/sectionate/section.py @@ -304,20 +304,21 @@ def _check_supported_topology(grid): A section is traced by walking a graph of corner points, and a step across a tile seam is recognised by its two ends having different face indices (see `transports._uv_for_edge`). - Two things are therefore required of `face_connections`, neither of which xgcm enforces - when the grid is built: - - * Every face named as a neighbour must exist, and every connection must be mutual: the - named neighbour must name this face back, along the axis it was reached by. Otherwise - the walk either steps into a face that is not there or into one with no way back -- - today an obscure failure from deep inside xgcm's padding. - * No face may be glued to itself. Both sides of such a seam carry the same face index, so - a crossing is indistinguishable from an ordinary step within the face, and the seam's - velocity is read from the wrong side -- silently, since nothing about the indices looks - out of place. Topologies of that shape belong on the axis rather than in - `face_connections`: a zonally periodic axis as ``padding="periodic"``, a bipolar/tripolar - north fold as ``padding={"Y": {"fold": ...}}`` on a single-tile grid. sectionate handles - both of those natively. + A face glued to *itself* therefore cannot be traced: both sides of such a seam carry the + same face index, so a crossing is indistinguishable from an ordinary step within the face, + and the seam's velocity is read from the wrong side -- silently, since nothing about the + indices looks out of place. xgcm builds such a grid without complaint and nothing further + down the pipeline notices, so this is the only place it is caught. Topologies of that shape + belong on the axis rather than in `face_connections`: a zonally periodic axis as + ``padding="periodic"``, a bipolar/tripolar north fold as ``padding={"Y": {"fold": ...}}`` + on a single-tile grid. sectionate handles both of those natively. + + The other two checks below -- that every face named as a neighbour exists, and that every + connection is mutual -- are defence in depth, not new requirements: `xgcm.Grid.__init__` + already rejects both while it builds its face-connection table (an unknown neighbour raises + `KeyError`, a one-sided connection `TypeError`), so any grid that reaches here has passed + them. They are kept because sectionate reads the `_face_connections` mapping directly, and + they say in one sentence what xgcm reports from deep inside its own construction. A grid can still prove unsupported once its corner graph is actually built -- a corner that ends up with more than four neighbours, or a staggered grid whose seam corners are stored on @@ -325,7 +326,9 @@ def _check_supported_topology(grid): `gridutils.build_neighbor_maps` and `gridutils._OuterTopology`. """ facedim = grid._facedim - faces = set(range(grid._ds.sizes[facedim])) + # Face *labels*, not a 0..N-1 range: a grid may label its faces however it likes + # (e.g. `face = [1, 2]`), and `face_connections` is keyed by those labels. + faces = set(grid._ds[facedim].values.tolist()) connections = (getattr(grid, "_face_connections", None) or {}).get(facedim, {}) for face, axis_sides in connections.items(): for axis, sides in axis_sides.items(): @@ -342,11 +345,14 @@ def _check_supported_topology(grid): "single-tile grid: padding='periodic' for a periodic axis, or " "padding={'Y': {'fold': 'corner'}} for a bipolar/tripolar north fold." ) + # Defence in depth (see the docstring): xgcm already rejects both of the + # cases below at `Grid.__init__`, so neither fires on a grid built the + # ordinary way. if neighbor_face not in faces: raise ValueError( f"Face {face} of this grid's `face_connections` names a neighbour face " f"{neighbor_face} along '{axis}', but the grid's '{facedim}' dimension " - f"only has {len(faces)} faces." + f"has no such face; its faces are {sorted(faces)}." ) back = connections.get(neighbor_face, {}).get(neighbor_axis, ()) if not any( diff --git a/sectionate/tests/test_cube_left_grid.py b/sectionate/tests/test_cube_left_grid.py index f1025a7..d60ef2a 100644 --- a/sectionate/tests/test_cube_left_grid.py +++ b/sectionate/tests/test_cube_left_grid.py @@ -263,14 +263,14 @@ def _coordinate_jittered_cube_left_grid(eps=2.0): def test_disagreeing_extrapolation_junction_merges_and_conserves(): - """Regression for MOM6-community/sectionate#49 (the LLC90 tiles 2/6/10 - rotated three-tile corner junction). When the three faces meeting at a corner - stored on no native face extrapolate *disagreeing* corner coordinates, the - coordinate merges cannot see that they are one physical point, so without the - topological (surrounding-cell) merge the junction splits into separate - degree-0 nodes and the stored radial velocity faces around it read as zero - from every neighbour's frame -- introducing a spurious convergence into the - junction cells. + """Regression for a three-tile corner junction that split into separate nodes + (as the LLC90 tiles 2/6/10 junction once did). When the three faces meeting + at a corner stored on no native face extrapolate *disagreeing* corner + coordinates, the coordinate merges cannot see that they are one physical + point, so without the topological (surrounding-cell) merge the junction + splits into separate degree-0 nodes and the stored radial velocity faces + around it read as zero from every neighbour's frame -- introducing a spurious + convergence into the junction cells. Here the exact cube's two un-stored vertices are turned into that failure mode by jittering each face's corner coordinates. The topological merge must diff --git a/sectionate/tests/test_ecco_llc90.py b/sectionate/tests/test_ecco_llc90.py index 7a08a25..a868a40 100644 --- a/sectionate/tests/test_ecco_llc90.py +++ b/sectionate/tests/test_ecco_llc90.py @@ -52,8 +52,13 @@ def test_llc90_outer_topology_builds_and_is_well_formed(): def test_llc90_three_tile_junction_resolves_to_one_node(): """The tiles-2/6/10 three-tile junction (a corner where three rotated tiles meet, stored on no face) must resolve to a SINGLE node whose representations - span all three tiles -- the topological `by_junction` merge, not three split - degree-0 nodes that would leak convergence (regression for #49).""" + span all three tiles -- the topological `by_junction` merge. + + Regression: when the three tiles' extrapolated coordinates for that corner + disagreed, the coordinate merges could not see it was one physical point and + the junction split into three separate degree-0 nodes. The velocity faces + around it then read as zero from every neighbour's frame, leaking a spurious + convergence into the junction's cells.""" from sectionate.gridutils import outer_topology ot = outer_topology(_load_grid()) diff --git a/sectionate/tests/test_section_multitile.py b/sectionate/tests/test_section_multitile.py index 4bb926b..e9c47da 100644 --- a/sectionate/tests/test_section_multitile.py +++ b/sectionate/tests/test_section_multitile.py @@ -525,22 +525,25 @@ def test_face_glued_to_itself_raises(): grid_section(grid, [0., 1.], [0., 1.]) -def test_unknown_neighbour_face_raises(): - """A `face_connections` entry naming a face the grid does not have is refused with a - clear message, rather than failing obscurely inside xgcm's padding later on.""" - grid = two_face_x_to_x() - grid._face_connections["face"][0]["X"] = (None, (7, "X", False)) - with pytest.raises(ValueError, match="only has 2 faces"): - grid_section(grid, [10., 100.], [0., 0.]) - - -def test_one_sided_face_connection_raises(): - """A seam declared from only one side is refused: the walk could cross it but never - step back, so the topology is not one sectionate can trace.""" - grid = two_face_x_to_x() - grid._face_connections["face"][1]["X"] = (None, None) - with pytest.raises(ValueError, match="not mutual"): - grid_section(grid, [10., 100.], [0., 0.]) +@pytest.mark.parametrize("fc, exc", [ + # a neighbour face the grid does not have + ({"face": {0: {"X": (None, (7, "X", False))}, + 1: {"X": ((0, "X", False), None)}}}, KeyError), + # a seam declared from only one side + ({"face": {0: {"X": (None, (1, "X", False))}, + 1: {"X": (None, None)}}}, TypeError), +]) +def test_xgcm_rejects_broken_face_connections_at_construction(fc, exc): + """xgcm -- not sectionate -- is what enforces that every face named as a neighbour + exists and that every connection is mutual: `xgcm.Grid.__init__` aborts on both while + it builds its face-connection table, so no such grid can ever reach `grid_section`. + `_check_supported_topology` restates both as defence in depth; this test pins where + the real enforcement lives, so that if xgcm ever stops doing it the backstop's status + is revisited rather than silently relied on.""" + ng = 6 + lon = np.zeros((2, ng, ng)); lat = np.zeros((2, ng, ng)) + with pytest.raises(exc): + _make_grid(lon, lat, fc) def test_supported_multitile_topologies_pass_the_check(): @@ -566,6 +569,34 @@ def test_supported_multitile_topologies_pass_the_check(): _check_supported_topology(llc90) +def test_non_zero_based_face_labels_pass_the_check(): + """Face labels need not be 0..N-1. A grid labelled `face = [1, 2]` is perfectly legal + -- xgcm accepts it and keys `face_connections` by those labels -- so the check must + read the labels themselves, not assume a range. (It previously built the set of known + faces as `range(n_faces)`, and rejected this grid for naming face 2 as a neighbour.)""" + from sectionate.section import _check_supported_topology + + ng = 6 + lon = np.zeros((2, ng, ng)); lat = np.zeros((2, ng, ng)) + lon[0] = np.broadcast_to(np.linspace(0, 90, ng), (ng, ng)) + lon[1] = np.broadcast_to(np.linspace(90, 180, ng), (ng, ng)) + ds = xr.Dataset({}, coords={ + "xg": (("xg",), np.arange(ng)), + "yg": (("yg",), np.arange(ng)), + "face": (("face",), np.array([1, 2])), + "geolon_c": (("face", "yg", "xg"), lon), + "geolat_c": (("face", "yg", "xg"), lat), + }) + grid = xgcm.Grid( + ds, coords={"X": {"outer": "xg"}, "Y": {"outer": "yg"}}, + padding="fill", fill_value=np.nan, + face_connections={"face": {1: {"X": (None, (2, "X", False))}, + 2: {"X": ((1, "X", False), None)}}}, + autoparse_metadata=False, + ) + _check_supported_topology(grid) + + def test_in_velocity_range_rejects_unknown_component(): """`_in_velocity_range` decides whether a velocity index exists in a face's own arrays: in range for an edge stored on that face, out of range for one stored on the diff --git a/sectionate/transports.py b/sectionate/transports.py index 5ec1502..6b3aa7a 100644 --- a/sectionate/transports.py +++ b/sectionate/transports.py @@ -39,13 +39,18 @@ def _in_velocity_range(var, vi, vj, ranges): `_anchor_velocity` builds a section edge's velocity index by offsetting the corner the edge leaves from. Read in the frame of the face that corner belongs to, an edge that stays inside that face always lands on a stored velocity. An edge that crosses a seam, - however, is stored on only ONE of the two faces it touches: read in the other face's + however, is stored on at most ONE of the two faces it touches: read in the other face's frame the same edge produces an index that falls outside that face's array -- past its high end, or negative. `_uv_for_edge` uses this check to tell the two apart, taking the source face's index when it is in range and otherwise the destination face's, so a velocity is always read in the frame of the face that stores it and never has to be rotated across the seam. + "At most one" because a staggered tiling can also leave a seam edge stored on NEITHER + face -- a crossing through a corner the faces share but neither one stores. Such an edge + is out of range in both frames, which is `_uv_for_edge`'s third branch: it is degenerate + and carries no flux. + `ranges` holds a face's array lengths along each axis: "Xc"/"Yc" for the tracer-center axes and "Xq"/"Yq" for the corner axes. A "V" (Y-direction) velocity sits at (X-center, Y-corner) and a "U" (X-direction) velocity at (X-corner, Y-center). From 2269bb0ee787aebec70746098e86d8782b02d29f Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Tue, 11 Aug 2026 13:14:32 -0700 Subject: [PATCH 4/5] Drop the face-existence and mutuality checks; trust xgcm `_check_supported_topology` re-checked two things `xgcm.Grid.__init__` already enforces: that every face named as a neighbour exists (an unknown one raises KeyError) and that every connection is mutual (a one-sided one raises TypeError). Both abort while xgcm builds its face-connection table, so no grid reaching this function could ever fail them. Keeping them as a backstop only guards against a `_face_connections` dict corrupted after construction, which is the caller's problem, not something to carry code for. Remove both branches, along with the `faces` label set that existed only to serve the existence check and the `neighbor_axis` unpacking that existed only to serve the mutuality check. What remains is the self-gluing rejection, which is load bearing: xgcm accepts a face glued to itself, nothing further downstream notices, and a section is then traced as though the seam were absent. The docstring now describes a function that enforces exactly one rule and records why that one is sectionate's to catch while the others are not. The test that pinned the face-label indexing bug would have passed vacuously once the existence check was gone, so it is reframed rather than dropped: it now asserts both that a legal `face = [1, 2]` grid passes and that the same grid with face 2 glued to itself is still rejected and named as face 2. That keeps a real regression guard on the surviving check comparing labels rather than positions along the face dimension. The test asserting xgcm rejects the malformed grids at construction is kept and its docstring updated -- that upstream contract is now the entire reason not to check here. Full suite: 81 passed, 0 skipped. Co-Authored-By: Claude Opus 5 (1M context) --- sectionate/section.py | 48 +++++-------------- sectionate/tests/test_section_multitile.py | 55 +++++++++++++++------- 2 files changed, 49 insertions(+), 54 deletions(-) diff --git a/sectionate/section.py b/sectionate/section.py index 9692931..cf51eca 100644 --- a/sectionate/section.py +++ b/sectionate/section.py @@ -307,18 +307,17 @@ def _check_supported_topology(grid): A face glued to *itself* therefore cannot be traced: both sides of such a seam carry the same face index, so a crossing is indistinguishable from an ordinary step within the face, and the seam's velocity is read from the wrong side -- silently, since nothing about the - indices looks out of place. xgcm builds such a grid without complaint and nothing further - down the pipeline notices, so this is the only place it is caught. Topologies of that shape - belong on the axis rather than in `face_connections`: a zonally periodic axis as - ``padding="periodic"``, a bipolar/tripolar north fold as ``padding={"Y": {"fold": ...}}`` - on a single-tile grid. sectionate handles both of those natively. - - The other two checks below -- that every face named as a neighbour exists, and that every - connection is mutual -- are defence in depth, not new requirements: `xgcm.Grid.__init__` - already rejects both while it builds its face-connection table (an unknown neighbour raises - `KeyError`, a one-sided connection `TypeError`), so any grid that reaches here has passed - them. They are kept because sectionate reads the `_face_connections` mapping directly, and - they say in one sentence what xgcm reports from deep inside its own construction. + indices looks out of place. Topologies of that shape belong on the axis rather than in + `face_connections`: a zonally periodic axis as ``padding="periodic"``, a bipolar/tripolar + north fold as ``padding={"Y": {"fold": ...}}`` on a single-tile grid. sectionate handles + both of those natively. + + That is the only thing checked here, because it is the only one left to check. xgcm builds + a self-glued grid without complaint, and nothing further downstream notices, so sectionate + has to catch it. The other ways `face_connections` can be malformed -- naming a face the + grid does not have, or declaring a seam from one side only -- already abort in + `xgcm.Grid.__init__` while it builds its face-connection table, so no such grid ever + reaches this function and re-checking them here would be dead code. A grid can still prove unsupported once its corner graph is actually built -- a corner that ends up with more than four neighbours, or a staggered grid whose seam corners are stored on @@ -326,17 +325,13 @@ def _check_supported_topology(grid): `gridutils.build_neighbor_maps` and `gridutils._OuterTopology`. """ facedim = grid._facedim - # Face *labels*, not a 0..N-1 range: a grid may label its faces however it likes - # (e.g. `face = [1, 2]`), and `face_connections` is keyed by those labels. - faces = set(grid._ds[facedim].values.tolist()) connections = (getattr(grid, "_face_connections", None) or {}).get(facedim, {}) for face, axis_sides in connections.items(): for axis, sides in axis_sides.items(): for side in sides: if side is None: continue - neighbor_face, neighbor_axis = side[0], side[1] - if neighbor_face == face: + if side[0] == face: raise NotImplementedError( f"Face {face} of this grid's `face_connections` is glued to itself " f"(along its own '{axis}' axis). sectionate tells a seam crossing from " @@ -345,25 +340,6 @@ def _check_supported_topology(grid): "single-tile grid: padding='periodic' for a periodic axis, or " "padding={'Y': {'fold': 'corner'}} for a bipolar/tripolar north fold." ) - # Defence in depth (see the docstring): xgcm already rejects both of the - # cases below at `Grid.__init__`, so neither fires on a grid built the - # ordinary way. - if neighbor_face not in faces: - raise ValueError( - f"Face {face} of this grid's `face_connections` names a neighbour face " - f"{neighbor_face} along '{axis}', but the grid's '{facedim}' dimension " - f"has no such face; its faces are {sorted(faces)}." - ) - back = connections.get(neighbor_face, {}).get(neighbor_axis, ()) - if not any( - s is not None and s[0] == face and s[1] == axis for s in back - ): - raise ValueError( - f"This grid's `face_connections` are not mutual: face {face} connects to " - f"face {neighbor_face} along '{axis}', but face {neighbor_face} does not " - f"connect back to face {face} along its '{neighbor_axis}' axis. Every " - "seam must be declared from both sides." - ) def create_section_composite( gridlon, diff --git a/sectionate/tests/test_section_multitile.py b/sectionate/tests/test_section_multitile.py index e9c47da..430c188 100644 --- a/sectionate/tests/test_section_multitile.py +++ b/sectionate/tests/test_section_multitile.py @@ -534,12 +534,13 @@ def test_face_glued_to_itself_raises(): 1: {"X": (None, None)}}}, TypeError), ]) def test_xgcm_rejects_broken_face_connections_at_construction(fc, exc): - """xgcm -- not sectionate -- is what enforces that every face named as a neighbour - exists and that every connection is mutual: `xgcm.Grid.__init__` aborts on both while - it builds its face-connection table, so no such grid can ever reach `grid_section`. - `_check_supported_topology` restates both as defence in depth; this test pins where - the real enforcement lives, so that if xgcm ever stops doing it the backstop's status - is revisited rather than silently relied on.""" + """xgcm -- not sectionate -- enforces that every face named as a neighbour exists and + that every connection is mutual: `xgcm.Grid.__init__` aborts on both while it builds + its face-connection table, so no such grid can reach `grid_section` in the first place. + + This is why `_check_supported_topology` does not re-check either one. The test pins + that upstream contract, so that if xgcm ever stops enforcing it we find out here rather + than by tracing a section through a topology nothing validated.""" ng = 6 lon = np.zeros((2, ng, ng)); lat = np.zeros((2, ng, ng)) with pytest.raises(exc): @@ -569,13 +570,8 @@ def test_supported_multitile_topologies_pass_the_check(): _check_supported_topology(llc90) -def test_non_zero_based_face_labels_pass_the_check(): - """Face labels need not be 0..N-1. A grid labelled `face = [1, 2]` is perfectly legal - -- xgcm accepts it and keys `face_connections` by those labels -- so the check must - read the labels themselves, not assume a range. (It previously built the set of known - faces as `range(n_faces)`, and rejected this grid for naming face 2 as a neighbour.)""" - from sectionate.section import _check_supported_topology - +def _labelled_grid(face_connections, labels=(1, 2)): + """A two-face 'outer' grid whose face coordinate carries arbitrary labels.""" ng = 6 lon = np.zeros((2, ng, ng)); lat = np.zeros((2, ng, ng)) lon[0] = np.broadcast_to(np.linspace(0, 90, ng), (ng, ng)) @@ -583,18 +579,41 @@ def test_non_zero_based_face_labels_pass_the_check(): ds = xr.Dataset({}, coords={ "xg": (("xg",), np.arange(ng)), "yg": (("yg",), np.arange(ng)), - "face": (("face",), np.array([1, 2])), + "face": (("face",), np.array(labels)), "geolon_c": (("face", "yg", "xg"), lon), "geolat_c": (("face", "yg", "xg"), lat), }) - grid = xgcm.Grid( + return xgcm.Grid( ds, coords={"X": {"outer": "xg"}, "Y": {"outer": "yg"}}, padding="fill", fill_value=np.nan, - face_connections={"face": {1: {"X": (None, (2, "X", False))}, - 2: {"X": ((1, "X", False), None)}}}, + face_connections={"face": face_connections}, autoparse_metadata=False, ) - _check_supported_topology(grid) + + +def test_check_reads_face_labels_not_positions(): + """Face labels need not be 0..N-1: a grid labelled `face = [1, 2]` is legal, and xgcm + keys `face_connections` by those labels. The check must therefore compare labels, never + positions along the face dimension. + + Both directions are pinned. A legal `[1, 2]` grid must pass -- an earlier version built + its set of known faces as `range(n_faces)` and rejected this grid for naming a face 2 + it thought did not exist. And a `[1, 2]` grid whose face 2 is glued to itself must still + be caught, and named as face 2, which is what would break if the surviving comparison + ever drifted back to positions.""" + from sectionate.section import _check_supported_topology + + _check_supported_topology(_labelled_grid( + {1: {"X": (None, (2, "X", False))}, + 2: {"X": ((1, "X", False), None)}} + )) + + glued = _labelled_grid( + {1: {"X": (None, (2, "X", False))}, + 2: {"X": ((1, "X", False), None), "Y": ((2, "Y", True), (2, "Y", True))}} + ) + with pytest.raises(NotImplementedError, match="Face 2 .* glued to itself"): + _check_supported_topology(glued) def test_in_velocity_range_rejects_unknown_component(): From 9fc40cd696c6057bf94672b567bf6de9c61c028f Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Tue, 11 Aug 2026 13:24:05 -0700 Subject: [PATCH 5/5] Refactor comments in section tracing logic Removed outdated comments regarding self-glued grids and clarified the purpose of the section tracing function. --- sectionate/section.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/sectionate/section.py b/sectionate/section.py index cf51eca..a856004 100644 --- a/sectionate/section.py +++ b/sectionate/section.py @@ -7,7 +7,7 @@ build_neighbor_maps, ) -# Two corner indices map to the same physical point on seams that fold or wrap (e.g. +# Two corner indices map to the same physical point on seams that fold Tor wrap (e.g. # the bipolar north fold). Geodesic distances below this many metres are treated as # the same point: far below any real grid spacing, far above float round-trip error. COINCIDENT_TOLERANCE_M = 1.e-3 @@ -304,6 +304,7 @@ def _check_supported_topology(grid): A section is traced by walking a graph of corner points, and a step across a tile seam is recognised by its two ends having different face indices (see `transports._uv_for_edge`). + A face glued to *itself* therefore cannot be traced: both sides of such a seam carry the same face index, so a crossing is indistinguishable from an ordinary step within the face, and the seam's velocity is read from the wrong side -- silently, since nothing about the @@ -312,13 +313,6 @@ def _check_supported_topology(grid): north fold as ``padding={"Y": {"fold": ...}}`` on a single-tile grid. sectionate handles both of those natively. - That is the only thing checked here, because it is the only one left to check. xgcm builds - a self-glued grid without complaint, and nothing further downstream notices, so sectionate - has to catch it. The other ways `face_connections` can be malformed -- naming a face the - grid does not have, or declaring a seam from one side only -- already abort in - `xgcm.Grid.__init__` while it builds its face-connection table, so no such grid ever - reaches this function and re-checking them here would be dead code. - A grid can still prove unsupported once its corner graph is actually built -- a corner that ends up with more than four neighbours, or a staggered grid whose seam corners are stored on no face, for instance. Those are rejected there, each with its own reason; see