diff --git a/CLAUDE.md b/CLAUDE.md index 8271b14..01dc28e 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 @@ -152,7 +143,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 — 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 8bf75fb..0ee4fdc 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 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-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 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 @@ -187,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 9accc51..4dcd0d5 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 @@ -199,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") @@ -227,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 ---------- @@ -267,6 +271,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 +311,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,45 +445,61 @@ 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** (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 / + 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 - 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 (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 -- 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). 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 +534,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): @@ -550,10 +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 + # 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( - "Corner arrays are larger than the grid's outer lattice; the " - "corner/center coordinates are inconsistent." + 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 @@ -694,9 +736,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 +782,16 @@ 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 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 + # (`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 +800,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 +868,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 +993,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 @@ -1159,13 +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 the 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. + 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 dcbffab..d63b243 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 @@ -335,27 +335,39 @@ 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. + Raise if the multi-tile `grid` describes a topology sectionate cannot trace a section on. - 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. + 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 + 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 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] - if neighbor_face == face: + if side[0] == 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." ) def create_section_composite( 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 27a764c..691f6fd 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 e42f766..430c188 100644 --- a/sectionate/tests/test_section_multitile.py +++ b/sectionate/tests/test_section_multitile.py @@ -511,15 +511,129 @@ 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.]) +@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 -- 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): + _make_grid(lon, lat, fc) + + +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 _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)) + 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(labels)), + "geolon_c": (("face", "yg", "xg"), lon), + "geolat_c": (("face", "yg", "xg"), lat), + }) + return xgcm.Grid( + ds, coords={"X": {"outer": "xg"}, "Y": {"outer": "yg"}}, + padding="fill", fill_value=np.nan, + face_connections={"face": face_connections}, + autoparse_metadata=False, + ) + + +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(): + """`_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..6b3aa7a 100644 --- a/sectionate/transports.py +++ b/sectionate/transports.py @@ -33,10 +33,33 @@ 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 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). + """ 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):