diff --git a/CLAUDE.md b/CLAUDE.md index 01dc28e..4eab778 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -139,7 +139,7 @@ The package is organized around a pipeline: define sections → map to grid → - **`section.py`** — Section definition and grid path algorithms. `Section` holds named waypoint coordinates; `GriddedSection` extends it with grid index information. `grid_section()` is the main entry point that maps geographic waypoints to grid vorticity-point indices `(i_c, j_c)` by walking the grid between waypoints along the `curve` requested. There are three: `"great circle"` (the default — every segment follows the geodesic), `"latitude circle"` (every segment follows a parallel, and a segment whose endpoints do not share a latitude — meridional ones included — raises), and `"latitude and great circle"` (decided per segment: constant-latitude segments follow the parallel, all others the geodesic). Under all three, every segment follows the **shortest** path between its two waypoints — raw longitudes are never read as a request to go the long way round, so a `"latitude circle"` segment written `0 -> 270` runs 90° *west* — which is why encircling the globe takes intermediate waypoints (e.g. `0 -> 120 -> 240 -> 360`). `_check_segment_span` resolves each segment's curve (via `_is_constant_latitude`, the one classification the metric choice and the legality check share) and raises for ill-posed ones: endpoints written exactly half a circle apart (neither way round is shorter — one error shared by all curves, quoting degrees of longitude along the parallel or degrees of arc, as applicable), or, under `"latitude circle"`, endpoints that do not share a latitude. The walk is deterministic and direction-independent: it admits neighbors strictly closer to the endpoint plus any seam twin of the current cell (so periodic/fold-seam crossings do not depend on floating-point rounding), and breaks ties by index. Grid topology is inferred entirely from `xgcm.Grid` metadata — each axis' `boundary` (periodic wrap, fill/extend wall, or a single-tile bipolar north fold `{"Y": {"fold": ...}}`) and `face_connections` for multi-tile grids — so there is no `topology` keyword. The pathfinder consumes topology-aware neighbor maps built by `gridutils.build_neighbor_maps`. -- **`transports.py`** — Transport computation along sections. `uvindices_from_qindices()` converts vorticity-point indices to U/V velocity-point indices using a per-position corner offset (`gridutils.corner_offset`) covering all three C-grid staggerings: 'outer', 'right', and 'left' (see "Corner staggering" below). `convergent_transport()` is the main function: it lazily computes signed normal transports with configurable orientation (positive inward to the polygon defined by the section). +- **`transports.py`** — Transport computation along sections. `uvindices_from_qindices()` converts vorticity-point indices to U/V velocity-point indices using a per-position corner offset (`gridutils.corner_offset`) covering all three C-grid staggerings: 'outer', 'right', and 'left' (see "Corner staggering" below). Before doing so it resolves **corner identity**, since a seam gives one physical corner two indices: multi-tile grids canonicalize each corner through `_OuterTopology`, and single-tile grids run `_insert_seam_twins()`, which splices the seam twin into the corner list wherever a step skips it (see "Seam twins" below). `convergent_transport()` is the main function: it lazily computes signed normal transports with configurable orientation (positive inward to the polygon defined by the section). - **`tracers.py`** — `extract_tracer()` interpolates tracer data to U/V points along a section path for cross-section plotting. @@ -151,6 +151,8 @@ The package is organized around a pipeline: define sections → map to grid → - **Vorticity points (q-points):** Sections are defined as paths through vorticity-point indices `(i_c, j_c)`. Consecutive q-points define velocity faces (either U or V). - **Corner staggering (three positions):** Vorticity sits at one of three xgcm positions: `"outer"` (MOM6 symmetric, M+1×N+1), `"right"` (MOM6 non-symmetric, M×N), or `"left"` (MITgcm/ECCO, incl. the lat-lon-cap, M×N). All three are native; they differ only by a per-position velocity index offset (`gridutils.corner_offset`: outer→0, right→+1, left→0). `"left"` indexes like `"outer"`; it differs only in array length and in that the *high* corner row/column is absent (so a section exactly on the north/east domain wall clips one row inside). +- **Seam twins (corner identity):** A seam gives one physical corner more than one index — the periodic wrap's first/last column, a shared multi-tile boundary corner, and the bipolar fold's mirrored seam-row columns. Transport attribution therefore has to know *which* indices denote the same point, or it will name the mirror image of the face a section crossed (with the sign that goes with it) and can count a face twice. Multi-tile grids get this from `_OuterTopology`, which resolves every corner to one canonical native `(face, j, i)`. Single-tile grids get it from `transports._insert_seam_twins()`: wherever consecutive corners are a physical edge but *not* index-adjacent — which is how the walk crosses a fold seam — it splices in the twin, so the crossing becomes one ordinary edge plus one zero-length twin edge, and zero-length edges emit no velocity face. This is the same construction as the `umaskutil`/`vmaskutil` masking of duplicated points in NEMO. The **sign flip** of an ORCA-style duplicated seam row (NEMO 4.2 manual, Appendix E) needs no handling of its own: the traversal direction and the stored velocity are read in the same index frame (geographically on both sides, for multi-tile), so a frame reversal cancels. A declared fold whose corner coordinates do not carry the fold symmetry has no twins to splice, and a crossing raises rather than silently mis-attributing. +- **Seam representation contract (two axes, and `q`):** A section crossing a seam can be written *seam-explicit* (both indices of the shared corner, the step between them spanning no cell) or *twin-free* (one of them). **Input is liberal:** `uvindices_from_qindices()` accepts either and normalizes internally — physical coincidence is tested *before* index adjacency, because the two indices of a fold-seam corner sit far apart in the index lattice (column `i` against column `nx-i`). This is what lets a section arrive from `grid_section`, from a mask traced on the grid (`regionate.boundaries` deliberately keeps both coincident corners at a seam junction, and closes its loops with a coincident edge), or from saved indices reloaded off disk. **Output is strict:** the returned faces are all real velocity faces. So the two axes have different lengths, on purpose: the **corner** axis holds `i_c`/`j_c`/`f_c`/`lons_c`/`lats_c`, all 1:1 with each other (a repeated corner is *kept* — what is dropped is the face between the pair, never the corner); the **`sect`** axis holds the velocity faces and their own `lon`/`lat`. Relating them is what `uvindices["q"]` is for: `q[k]` is the index into the caller's corner arrays that face `k` starts at, so face `k` spans corners `q[k]` and `q[k]+1`. It is strictly increasing, skips exactly the steps that carried no face, and is carried as a coordinate on `convergent_transport()` and `extract_tracer()` output. - **Sign conventions:** For a *closed* section, `convergent_transport()` determines orientation (clockwise/counterclockwise) using stereographic projection and signed polygon area, then applies sign corrections so positive transport means "inward" (toward the enclosed polygon). For an *open* section there is no enclosing polygon, so `positive_in` is undefined; it instead uses the **left-of-transect** convention — `positive_in=True` makes positive transport point to the left of the section as traversed from the first to the last waypoint — and emits a `UserWarning`. `is_section_counterclockwise()` is only consulted in the closed case. - **xgcm.Grid dependency:** The package relies heavily on `xgcm.Grid` for grid metadata (axis boundaries, coordinate positions, dataset access via `grid._ds`). diff --git a/docs/source/algorithm.md b/docs/source/algorithm.md index e82a811..b8e0feb 100644 --- a/docs/source/algorithm.md +++ b/docs/source/algorithm.md @@ -226,7 +226,8 @@ from sectionate.transports import uvindices_from_qindices uv = uvindices_from_qindices(grid, i_c, j_c) # uv["var"] is "U" or "V" per face; uv["i"], uv["j"] index that velocity; -# uv["Xinc"], uv["Yinc"] record the direction of travel through each face. +# uv["Xinc"], uv["Yinc"] record the direction of travel through each face; +# uv["q"] says which step of (i_c, j_c) each face came from. ``` Whether a face is a `U` or a `V` point falls straight out of which way the step went, though @@ -239,6 +240,29 @@ staggering (`outer`, `right` or `left`) and is read from the grid rather than as "Up to" `N-1`, because a consecutive pair that resolves to the same physical point — the twin corners of a seam again — spans no cell and carries no flux, so it emits no face at all. +Twins matter here for a second reason. Naming the face between two corners assumes they are +neighbours *in index space*, and across a bipolar fold that can fail: the walk may leave one +index representation of a seam corner and land next to the *other*, a real physical edge whose +endpoints are far apart in `i`. Attributing that step from the source corner alone would name +the mirrored column's face, with the sign that goes with it. So before converting the chain, +`uvindices_from_qindices` resolves corner identity — multi-tile grids canonicalise every corner +through the outer-lattice topology, single-tile grids splice the skipped twin back into the +chain — after which every remaining step is an ordinary adjacent one, and the zero-length twin +edge it introduces drops out by the rule above rather than being counted twice. + +Which representation the chain arrives in does not matter. A section may already carry both +indices of a seam corner — that is what the walk returns across a periodic wrap, and what a +boundary traced on a cell mask returns at every seam junction — or only one, which is what it +returns across a fold. Either way the faces come out the same, because coincidence is tested +before adjacency: a pair that is *already* the zero-length twin edge is passed straight to the +drop rather than being mistaken for a step needing repair. + +The consequence is that faces and corner steps are not one-to-one, so `uv["q"]` says which is +which: face `k` spans corners `q[k]` and `q[k]+1`, and the steps `q` skips are exactly the ones +that carried no flux. It rides along on `convergent_transport` and `extract_tracer` output, so +a transport can be put back onto the section's corners — to plot it against them, to slice out +the part belonging to one stretch of the path, or to find a child section inside its parent. + `transports.convergent_transport` then accumulates the signed normal transport through those faces. For a **closed** section it works out the traversal orientation and signs everything so that positive means *into* the enclosed region; for an **open** section there is no inside, so diff --git a/examples/1_creating_an_OSNAP_section.ipynb b/examples/1_creating_an_OSNAP_section.ipynb index 7f58a59..6d09c38 100644 --- a/examples/1_creating_an_OSNAP_section.ipynb +++ b/examples/1_creating_an_OSNAP_section.ipynb @@ -23,10 +23,10 @@ "execution_count": 1, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:42.488689Z", - "iopub.status.busy": "2026-08-04T02:01:42.488449Z", - "iopub.status.idle": "2026-08-04T02:01:45.442158Z", - "shell.execute_reply": "2026-08-04T02:01:45.441763Z" + "iopub.execute_input": "2026-08-11T22:00:01.034760Z", + "iopub.status.busy": "2026-08-11T22:00:01.034675Z", + "iopub.status.idle": "2026-08-11T22:00:05.409204Z", + "shell.execute_reply": "2026-08-11T22:00:05.408646Z" } }, "outputs": [ @@ -34,7 +34,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Sectionate version: 0.4.0rc2.dev3\n" + "Sectionate version: 0.4.0rc3.dev19\n" ] } ], @@ -74,10 +74,10 @@ "execution_count": 2, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:45.456437Z", - "iopub.status.busy": "2026-08-04T02:01:45.456296Z", - "iopub.status.idle": "2026-08-04T02:01:45.489296Z", - "shell.execute_reply": "2026-08-04T02:01:45.488946Z" + "iopub.execute_input": "2026-08-11T22:00:05.430876Z", + "iopub.status.busy": "2026-08-11T22:00:05.430658Z", + "iopub.status.idle": "2026-08-11T22:00:05.475380Z", + "shell.execute_reply": "2026-08-11T22:00:05.474895Z" } }, "outputs": [ @@ -122,10 +122,10 @@ "execution_count": 3, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:45.490519Z", - "iopub.status.busy": "2026-08-04T02:01:45.490459Z", - "iopub.status.idle": "2026-08-04T02:01:45.492106Z", - "shell.execute_reply": "2026-08-04T02:01:45.491756Z" + "iopub.execute_input": "2026-08-11T22:00:05.476668Z", + "iopub.status.busy": "2026-08-11T22:00:05.476575Z", + "iopub.status.idle": "2026-08-11T22:00:05.478541Z", + "shell.execute_reply": "2026-08-11T22:00:05.478091Z" } }, "outputs": [], @@ -147,10 +147,10 @@ "execution_count": 4, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:45.493150Z", - "iopub.status.busy": "2026-08-04T02:01:45.493100Z", - "iopub.status.idle": "2026-08-04T02:01:45.875493Z", - "shell.execute_reply": "2026-08-04T02:01:45.874991Z" + "iopub.execute_input": "2026-08-11T22:00:05.479732Z", + "iopub.status.busy": "2026-08-11T22:00:05.479656Z", + "iopub.status.idle": "2026-08-11T22:00:06.089589Z", + "shell.execute_reply": "2026-08-11T22:00:06.088832Z" } }, "outputs": [], @@ -167,10 +167,10 @@ "execution_count": 5, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:45.876703Z", - "iopub.status.busy": "2026-08-04T02:01:45.876630Z", - "iopub.status.idle": "2026-08-04T02:01:46.650395Z", - "shell.execute_reply": "2026-08-04T02:01:46.649969Z" + "iopub.execute_input": "2026-08-11T22:00:06.091657Z", + "iopub.status.busy": "2026-08-11T22:00:06.091523Z", + "iopub.status.idle": "2026-08-11T22:00:07.364580Z", + "shell.execute_reply": "2026-08-11T22:00:07.364011Z" } }, "outputs": [ @@ -218,10 +218,10 @@ "execution_count": 6, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:46.651743Z", - "iopub.status.busy": "2026-08-04T02:01:46.651652Z", - "iopub.status.idle": "2026-08-04T02:01:46.673811Z", - "shell.execute_reply": "2026-08-04T02:01:46.673484Z" + "iopub.execute_input": "2026-08-11T22:00:07.366904Z", + "iopub.status.busy": "2026-08-11T22:00:07.366796Z", + "iopub.status.idle": "2026-08-11T22:00:07.394288Z", + "shell.execute_reply": "2026-08-11T22:00:07.393714Z" } }, "outputs": [], @@ -249,10 +249,10 @@ "execution_count": 7, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:46.674910Z", - "iopub.status.busy": "2026-08-04T02:01:46.674851Z", - "iopub.status.idle": "2026-08-04T02:01:46.808809Z", - "shell.execute_reply": "2026-08-04T02:01:46.808475Z" + "iopub.execute_input": "2026-08-11T22:00:07.395861Z", + "iopub.status.busy": "2026-08-11T22:00:07.395753Z", + "iopub.status.idle": "2026-08-11T22:00:07.617234Z", + "shell.execute_reply": "2026-08-11T22:00:07.616819Z" } }, "outputs": [ @@ -307,10 +307,10 @@ "execution_count": 8, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:46.810059Z", - "iopub.status.busy": "2026-08-04T02:01:46.809995Z", - "iopub.status.idle": "2026-08-04T02:01:46.912776Z", - "shell.execute_reply": "2026-08-04T02:01:46.912484Z" + "iopub.execute_input": "2026-08-11T22:00:07.618807Z", + "iopub.status.busy": "2026-08-11T22:00:07.618701Z", + "iopub.status.idle": "2026-08-11T22:00:07.781094Z", + "shell.execute_reply": "2026-08-11T22:00:07.780609Z" } }, "outputs": [ @@ -318,7 +318,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/Users/hfdrake/miniforge3/envs/docs_env_sectionate_fix-layer-interface-matching/lib/python3.14/site-packages/xgcm/grid.py:463: UserWarning: The north-fold (tripolar) boundary condition is experimental. Its API and numerical behavior may change in future releases, and it has not yet been validated across the full range of grid configurations. Please review results carefully and report any issues at https://github.com/xgcm/xgcm/issues.\n", + "/opt/anaconda3/envs/docs_env_sectionate_seam-corner-identity/lib/python3.14t/site-packages/xgcm/grid.py:463: UserWarning: The north-fold (tripolar) boundary condition is experimental. Its API and numerical behavior may change in future releases, and it has not yet been validated across the full range of grid configurations. Please review results carefully and report any issues at https://github.com/xgcm/xgcm/issues.\n", " warnings.warn(\n" ] } @@ -343,10 +343,10 @@ "execution_count": 9, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:46.913869Z", - "iopub.status.busy": "2026-08-04T02:01:46.913814Z", - "iopub.status.idle": "2026-08-04T02:01:46.915373Z", - "shell.execute_reply": "2026-08-04T02:01:46.915068Z" + "iopub.execute_input": "2026-08-11T22:00:07.782403Z", + "iopub.status.busy": "2026-08-11T22:00:07.782308Z", + "iopub.status.idle": "2026-08-11T22:00:07.784390Z", + "shell.execute_reply": "2026-08-11T22:00:07.783677Z" } }, "outputs": [], @@ -370,10 +370,10 @@ "execution_count": 10, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:46.916215Z", - "iopub.status.busy": "2026-08-04T02:01:46.916166Z", - "iopub.status.idle": "2026-08-04T02:01:46.932327Z", - "shell.execute_reply": "2026-08-04T02:01:46.931936Z" + "iopub.execute_input": "2026-08-11T22:00:07.785639Z", + "iopub.status.busy": "2026-08-11T22:00:07.785549Z", + "iopub.status.idle": "2026-08-11T22:00:07.807944Z", + "shell.execute_reply": "2026-08-11T22:00:07.807526Z" } }, "outputs": [], @@ -390,10 +390,10 @@ "execution_count": 11, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:46.933351Z", - "iopub.status.busy": "2026-08-04T02:01:46.933294Z", - "iopub.status.idle": "2026-08-04T02:01:47.029043Z", - "shell.execute_reply": "2026-08-04T02:01:47.028616Z" + "iopub.execute_input": "2026-08-11T22:00:07.809303Z", + "iopub.status.busy": "2026-08-11T22:00:07.809220Z", + "iopub.status.idle": "2026-08-11T22:00:07.943970Z", + "shell.execute_reply": "2026-08-11T22:00:07.943581Z" } }, "outputs": [ @@ -433,10 +433,10 @@ "execution_count": 12, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:47.030083Z", - "iopub.status.busy": "2026-08-04T02:01:47.030000Z", - "iopub.status.idle": "2026-08-04T02:01:47.045217Z", - "shell.execute_reply": "2026-08-04T02:01:47.044885Z" + "iopub.execute_input": "2026-08-11T22:00:07.945413Z", + "iopub.status.busy": "2026-08-11T22:00:07.945301Z", + "iopub.status.idle": "2026-08-11T22:00:07.966268Z", + "shell.execute_reply": "2026-08-11T22:00:07.965960Z" } }, "outputs": [], @@ -453,10 +453,10 @@ "execution_count": 13, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:47.046152Z", - "iopub.status.busy": "2026-08-04T02:01:47.046094Z", - "iopub.status.idle": "2026-08-04T02:01:47.111056Z", - "shell.execute_reply": "2026-08-04T02:01:47.110790Z" + "iopub.execute_input": "2026-08-11T22:00:07.967809Z", + "iopub.status.busy": "2026-08-11T22:00:07.967736Z", + "iopub.status.idle": "2026-08-11T22:00:08.066435Z", + "shell.execute_reply": "2026-08-11T22:00:08.065974Z" } }, "outputs": [ @@ -529,10 +529,10 @@ "execution_count": 14, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:47.112255Z", - "iopub.status.busy": "2026-08-04T02:01:47.112186Z", - "iopub.status.idle": "2026-08-04T02:01:47.114192Z", - "shell.execute_reply": "2026-08-04T02:01:47.113925Z" + "iopub.execute_input": "2026-08-11T22:00:08.067976Z", + "iopub.status.busy": "2026-08-11T22:00:08.067864Z", + "iopub.status.idle": "2026-08-11T22:00:08.070432Z", + "shell.execute_reply": "2026-08-11T22:00:08.070004Z" } }, "outputs": [ @@ -593,10 +593,10 @@ "execution_count": 15, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:47.115147Z", - "iopub.status.busy": "2026-08-04T02:01:47.115093Z", - "iopub.status.idle": "2026-08-04T02:01:47.117491Z", - "shell.execute_reply": "2026-08-04T02:01:47.117191Z" + "iopub.execute_input": "2026-08-11T22:00:08.071557Z", + "iopub.status.busy": "2026-08-11T22:00:08.071464Z", + "iopub.status.idle": "2026-08-11T22:00:08.074713Z", + "shell.execute_reply": "2026-08-11T22:00:08.074240Z" } }, "outputs": [ @@ -644,10 +644,10 @@ "execution_count": 16, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:47.118350Z", - "iopub.status.busy": "2026-08-04T02:01:47.118300Z", - "iopub.status.idle": "2026-08-04T02:01:47.120829Z", - "shell.execute_reply": "2026-08-04T02:01:47.120449Z" + "iopub.execute_input": "2026-08-11T22:00:08.075740Z", + "iopub.status.busy": "2026-08-11T22:00:08.075670Z", + "iopub.status.idle": "2026-08-11T22:00:08.078496Z", + "shell.execute_reply": "2026-08-11T22:00:08.078138Z" } }, "outputs": [ @@ -702,10 +702,10 @@ "execution_count": 17, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:47.121682Z", - "iopub.status.busy": "2026-08-04T02:01:47.121634Z", - "iopub.status.idle": "2026-08-04T02:01:47.123385Z", - "shell.execute_reply": "2026-08-04T02:01:47.123069Z" + "iopub.execute_input": "2026-08-11T22:00:08.079615Z", + "iopub.status.busy": "2026-08-11T22:00:08.079541Z", + "iopub.status.idle": "2026-08-11T22:00:08.081692Z", + "shell.execute_reply": "2026-08-11T22:00:08.081314Z" } }, "outputs": [ @@ -745,10 +745,10 @@ "execution_count": 18, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:47.124215Z", - "iopub.status.busy": "2026-08-04T02:01:47.124167Z", - "iopub.status.idle": "2026-08-04T02:01:47.162616Z", - "shell.execute_reply": "2026-08-04T02:01:47.162264Z" + "iopub.execute_input": "2026-08-11T22:00:08.082891Z", + "iopub.status.busy": "2026-08-11T22:00:08.082819Z", + "iopub.status.idle": "2026-08-11T22:00:08.140330Z", + "shell.execute_reply": "2026-08-11T22:00:08.139816Z" } }, "outputs": [], @@ -765,10 +765,10 @@ "execution_count": 19, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:47.163631Z", - "iopub.status.busy": "2026-08-04T02:01:47.163580Z", - "iopub.status.idle": "2026-08-04T02:01:47.255905Z", - "shell.execute_reply": "2026-08-04T02:01:47.255590Z" + "iopub.execute_input": "2026-08-11T22:00:08.141743Z", + "iopub.status.busy": "2026-08-11T22:00:08.141647Z", + "iopub.status.idle": "2026-08-11T22:00:08.281095Z", + "shell.execute_reply": "2026-08-11T22:00:08.280404Z" } }, "outputs": [ @@ -817,10 +817,10 @@ "execution_count": 20, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:47.256958Z", - "iopub.status.busy": "2026-08-04T02:01:47.256900Z", - "iopub.status.idle": "2026-08-04T02:01:47.296309Z", - "shell.execute_reply": "2026-08-04T02:01:47.295955Z" + "iopub.execute_input": "2026-08-11T22:00:08.282524Z", + "iopub.status.busy": "2026-08-11T22:00:08.282405Z", + "iopub.status.idle": "2026-08-11T22:00:08.339928Z", + "shell.execute_reply": "2026-08-11T22:00:08.339479Z" } }, "outputs": [ @@ -868,10 +868,10 @@ "execution_count": 21, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:47.297356Z", - "iopub.status.busy": "2026-08-04T02:01:47.297301Z", - "iopub.status.idle": "2026-08-04T02:01:47.299807Z", - "shell.execute_reply": "2026-08-04T02:01:47.299484Z" + "iopub.execute_input": "2026-08-11T22:00:08.341453Z", + "iopub.status.busy": "2026-08-11T22:00:08.341327Z", + "iopub.status.idle": "2026-08-11T22:00:08.345928Z", + "shell.execute_reply": "2026-08-11T22:00:08.345453Z" } }, "outputs": [], @@ -895,10 +895,10 @@ "execution_count": 22, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:47.300821Z", - "iopub.status.busy": "2026-08-04T02:01:47.300769Z", - "iopub.status.idle": "2026-08-04T02:01:47.393286Z", - "shell.execute_reply": "2026-08-04T02:01:47.392961Z" + "iopub.execute_input": "2026-08-11T22:00:08.347390Z", + "iopub.status.busy": "2026-08-11T22:00:08.347304Z", + "iopub.status.idle": "2026-08-11T22:00:08.483354Z", + "shell.execute_reply": "2026-08-11T22:00:08.482901Z" } }, "outputs": [ diff --git a/examples/2_OSNAP_transports_CM4p25.ipynb b/examples/2_OSNAP_transports_CM4p25.ipynb index a1e02c6..1e99719 100644 --- a/examples/2_OSNAP_transports_CM4p25.ipynb +++ b/examples/2_OSNAP_transports_CM4p25.ipynb @@ -21,10 +21,10 @@ "execution_count": 1, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:48.388420Z", - "iopub.status.busy": "2026-08-04T02:01:48.388248Z", - "iopub.status.idle": "2026-08-04T02:01:49.080762Z", - "shell.execute_reply": "2026-08-04T02:01:49.080434Z" + "iopub.execute_input": "2026-08-11T22:00:09.648479Z", + "iopub.status.busy": "2026-08-11T22:00:09.648379Z", + "iopub.status.idle": "2026-08-11T22:00:10.773702Z", + "shell.execute_reply": "2026-08-11T22:00:10.773142Z" } }, "outputs": [ @@ -32,7 +32,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Sectionate version: 0.4.0rc2.dev3\n" + "Sectionate version: 0.4.0rc3.dev19\n" ] } ], @@ -60,10 +60,10 @@ "execution_count": 2, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.082054Z", - "iopub.status.busy": "2026-08-04T02:01:49.081933Z", - "iopub.status.idle": "2026-08-04T02:01:49.208403Z", - "shell.execute_reply": "2026-08-04T02:01:49.207974Z" + "iopub.execute_input": "2026-08-11T22:00:10.775409Z", + "iopub.status.busy": "2026-08-11T22:00:10.775225Z", + "iopub.status.idle": "2026-08-11T22:00:10.993780Z", + "shell.execute_reply": "2026-08-11T22:00:10.992793Z" } }, "outputs": [ @@ -71,7 +71,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/Users/hfdrake/miniforge3/envs/docs_env_sectionate_fix-layer-interface-matching/lib/python3.14/site-packages/xgcm/grid.py:463: UserWarning: The north-fold (tripolar) boundary condition is experimental. Its API and numerical behavior may change in future releases, and it has not yet been validated across the full range of grid configurations. Please review results carefully and report any issues at https://github.com/xgcm/xgcm/issues.\n", + "/opt/anaconda3/envs/docs_env_sectionate_seam-corner-identity/lib/python3.14t/site-packages/xgcm/grid.py:463: UserWarning: The north-fold (tripolar) boundary condition is experimental. Its API and numerical behavior may change in future releases, and it has not yet been validated across the full range of grid configurations. Please review results carefully and report any issues at https://github.com/xgcm/xgcm/issues.\n", " warnings.warn(\n" ] } @@ -94,10 +94,10 @@ "execution_count": 3, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.209584Z", - "iopub.status.busy": "2026-08-04T02:01:49.209516Z", - "iopub.status.idle": "2026-08-04T02:01:49.211241Z", - "shell.execute_reply": "2026-08-04T02:01:49.210950Z" + "iopub.execute_input": "2026-08-11T22:00:10.995290Z", + "iopub.status.busy": "2026-08-11T22:00:10.995161Z", + "iopub.status.idle": "2026-08-11T22:00:10.997352Z", + "shell.execute_reply": "2026-08-11T22:00:10.996858Z" } }, "outputs": [], @@ -124,10 +124,10 @@ "execution_count": 4, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.212231Z", - "iopub.status.busy": "2026-08-04T02:01:49.212180Z", - "iopub.status.idle": "2026-08-04T02:01:49.226750Z", - "shell.execute_reply": "2026-08-04T02:01:49.226317Z" + "iopub.execute_input": "2026-08-11T22:00:10.998779Z", + "iopub.status.busy": "2026-08-11T22:00:10.998682Z", + "iopub.status.idle": "2026-08-11T22:00:11.024320Z", + "shell.execute_reply": "2026-08-11T22:00:11.023705Z" } }, "outputs": [], @@ -144,10 +144,10 @@ "execution_count": 5, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.227864Z", - "iopub.status.busy": "2026-08-04T02:01:49.227804Z", - "iopub.status.idle": "2026-08-04T02:01:49.343097Z", - "shell.execute_reply": "2026-08-04T02:01:49.342741Z" + "iopub.execute_input": "2026-08-11T22:00:11.026191Z", + "iopub.status.busy": "2026-08-11T22:00:11.026080Z", + "iopub.status.idle": "2026-08-11T22:00:11.202010Z", + "shell.execute_reply": "2026-08-11T22:00:11.201517Z" }, "tags": [] }, @@ -195,10 +195,10 @@ "execution_count": 6, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.344132Z", - "iopub.status.busy": "2026-08-04T02:01:49.344069Z", - "iopub.status.idle": "2026-08-04T02:01:49.364171Z", - "shell.execute_reply": "2026-08-04T02:01:49.363778Z" + "iopub.execute_input": "2026-08-11T22:00:11.203427Z", + "iopub.status.busy": "2026-08-11T22:00:11.203317Z", + "iopub.status.idle": "2026-08-11T22:00:11.257923Z", + "shell.execute_reply": "2026-08-11T22:00:11.257471Z" } }, "outputs": [ @@ -206,7 +206,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/Users/hfdrake/code/sectionate/sectionate/transports.py:565: UserWarning: This section is open (its first and last waypoints are far apart), so it does not enclose a polygon and the inward/outward sense of `positive_in` is undefined. Sectionate applies the left-of-transect convention instead: with `positive_in=True` positive transport points to the LEFT of the section as it is traversed from the first to the last waypoint (set `positive_in=False` to flip to right-of-transect). Verify the sign matches your expectations, e.g. by plotting the section's normal vectors.\n", + "/Users/hfdrake/code/sectionate/sectionate/transports.py:738: UserWarning: This section is open (its first and last waypoints are far apart), so it does not enclose a polygon and the inward/outward sense of `positive_in` is undefined. Sectionate applies the left-of-transect convention instead: with `positive_in=True` positive transport points to the LEFT of the section as it is traversed from the first to the last waypoint (set `positive_in=False` to flip to right-of-transect). Verify the sign matches your expectations, e.g. by plotting the section's normal vectors.\n", " warnings.warn(\n" ] } @@ -220,10 +220,10 @@ "execution_count": 7, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.365281Z", - "iopub.status.busy": "2026-08-04T02:01:49.365213Z", - "iopub.status.idle": "2026-08-04T02:01:49.446445Z", - "shell.execute_reply": "2026-08-04T02:01:49.446021Z" + "iopub.execute_input": "2026-08-11T22:00:11.259483Z", + "iopub.status.busy": "2026-08-11T22:00:11.259387Z", + "iopub.status.idle": "2026-08-11T22:00:11.382951Z", + "shell.execute_reply": "2026-08-11T22:00:11.382488Z" } }, "outputs": [ @@ -255,10 +255,10 @@ "execution_count": 8, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.447595Z", - "iopub.status.busy": "2026-08-04T02:01:49.447529Z", - "iopub.status.idle": "2026-08-04T02:01:49.522301Z", - "shell.execute_reply": "2026-08-04T02:01:49.521934Z" + "iopub.execute_input": "2026-08-11T22:00:11.384532Z", + "iopub.status.busy": "2026-08-11T22:00:11.384421Z", + "iopub.status.idle": "2026-08-11T22:00:11.503329Z", + "shell.execute_reply": "2026-08-11T22:00:11.502388Z" } }, "outputs": [ @@ -317,10 +317,10 @@ "execution_count": 9, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.523408Z", - "iopub.status.busy": "2026-08-04T02:01:49.523349Z", - "iopub.status.idle": "2026-08-04T02:01:49.524914Z", - "shell.execute_reply": "2026-08-04T02:01:49.524639Z" + "iopub.execute_input": "2026-08-11T22:00:11.505047Z", + "iopub.status.busy": "2026-08-11T22:00:11.504930Z", + "iopub.status.idle": "2026-08-11T22:00:11.507143Z", + "shell.execute_reply": "2026-08-11T22:00:11.506647Z" } }, "outputs": [], @@ -340,10 +340,10 @@ "execution_count": 10, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.525806Z", - "iopub.status.busy": "2026-08-04T02:01:49.525746Z", - "iopub.status.idle": "2026-08-04T02:01:49.540458Z", - "shell.execute_reply": "2026-08-04T02:01:49.540125Z" + "iopub.execute_input": "2026-08-11T22:00:11.508555Z", + "iopub.status.busy": "2026-08-11T22:00:11.508457Z", + "iopub.status.idle": "2026-08-11T22:00:11.531318Z", + "shell.execute_reply": "2026-08-11T22:00:11.530823Z" } }, "outputs": [], @@ -360,10 +360,10 @@ "execution_count": 11, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.541502Z", - "iopub.status.busy": "2026-08-04T02:01:49.541448Z", - "iopub.status.idle": "2026-08-04T02:01:49.631238Z", - "shell.execute_reply": "2026-08-04T02:01:49.630940Z" + "iopub.execute_input": "2026-08-11T22:00:11.532780Z", + "iopub.status.busy": "2026-08-11T22:00:11.532695Z", + "iopub.status.idle": "2026-08-11T22:00:11.657138Z", + "shell.execute_reply": "2026-08-11T22:00:11.656666Z" } }, "outputs": [ @@ -410,10 +410,10 @@ "execution_count": 12, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.632524Z", - "iopub.status.busy": "2026-08-04T02:01:49.632463Z", - "iopub.status.idle": "2026-08-04T02:01:49.670014Z", - "shell.execute_reply": "2026-08-04T02:01:49.669657Z" + "iopub.execute_input": "2026-08-11T22:00:11.658538Z", + "iopub.status.busy": "2026-08-11T22:00:11.658425Z", + "iopub.status.idle": "2026-08-11T22:00:11.704483Z", + "shell.execute_reply": "2026-08-11T22:00:11.704123Z" } }, "outputs": [ @@ -421,7 +421,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/Users/hfdrake/code/sectionate/sectionate/transports.py:565: UserWarning: This section is open (its first and last waypoints are far apart), so it does not enclose a polygon and the inward/outward sense of `positive_in` is undefined. Sectionate applies the left-of-transect convention instead: with `positive_in=True` positive transport points to the LEFT of the section as it is traversed from the first to the last waypoint (set `positive_in=False` to flip to right-of-transect). Verify the sign matches your expectations, e.g. by plotting the section's normal vectors.\n", + "/Users/hfdrake/code/sectionate/sectionate/transports.py:738: UserWarning: This section is open (its first and last waypoints are far apart), so it does not enclose a polygon and the inward/outward sense of `positive_in` is undefined. Sectionate applies the left-of-transect convention instead: with `positive_in=True` positive transport points to the LEFT of the section as it is traversed from the first to the last waypoint (set `positive_in=False` to flip to right-of-transect). Verify the sign matches your expectations, e.g. by plotting the section's normal vectors.\n", " warnings.warn(\n" ] } @@ -435,10 +435,10 @@ "execution_count": 13, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.671028Z", - "iopub.status.busy": "2026-08-04T02:01:49.670972Z", - "iopub.status.idle": "2026-08-04T02:01:49.749043Z", - "shell.execute_reply": "2026-08-04T02:01:49.748684Z" + "iopub.execute_input": "2026-08-11T22:00:11.705824Z", + "iopub.status.busy": "2026-08-11T22:00:11.705751Z", + "iopub.status.idle": "2026-08-11T22:00:11.818993Z", + "shell.execute_reply": "2026-08-11T22:00:11.818578Z" } }, "outputs": [ @@ -470,10 +470,10 @@ "execution_count": 14, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.750155Z", - "iopub.status.busy": "2026-08-04T02:01:49.750092Z", - "iopub.status.idle": "2026-08-04T02:01:49.836438Z", - "shell.execute_reply": "2026-08-04T02:01:49.836017Z" + "iopub.execute_input": "2026-08-11T22:00:11.820492Z", + "iopub.status.busy": "2026-08-11T22:00:11.820399Z", + "iopub.status.idle": "2026-08-11T22:00:11.952720Z", + "shell.execute_reply": "2026-08-11T22:00:11.952198Z" } }, "outputs": [ @@ -530,10 +530,10 @@ "execution_count": 15, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.837429Z", - "iopub.status.busy": "2026-08-04T02:01:49.837368Z", - "iopub.status.idle": "2026-08-04T02:01:49.838781Z", - "shell.execute_reply": "2026-08-04T02:01:49.838537Z" + "iopub.execute_input": "2026-08-11T22:00:11.954149Z", + "iopub.status.busy": "2026-08-11T22:00:11.954035Z", + "iopub.status.idle": "2026-08-11T22:00:11.956006Z", + "shell.execute_reply": "2026-08-11T22:00:11.955521Z" } }, "outputs": [], @@ -553,10 +553,10 @@ "execution_count": 16, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.839688Z", - "iopub.status.busy": "2026-08-04T02:01:49.839633Z", - "iopub.status.idle": "2026-08-04T02:01:49.867487Z", - "shell.execute_reply": "2026-08-04T02:01:49.867055Z" + "iopub.execute_input": "2026-08-11T22:00:11.957258Z", + "iopub.status.busy": "2026-08-11T22:00:11.957157Z", + "iopub.status.idle": "2026-08-11T22:00:11.995477Z", + "shell.execute_reply": "2026-08-11T22:00:11.994992Z" } }, "outputs": [], @@ -573,10 +573,10 @@ "execution_count": 17, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.868502Z", - "iopub.status.busy": "2026-08-04T02:01:49.868449Z", - "iopub.status.idle": "2026-08-04T02:01:49.898245Z", - "shell.execute_reply": "2026-08-04T02:01:49.897927Z" + "iopub.execute_input": "2026-08-11T22:00:11.996734Z", + "iopub.status.busy": "2026-08-11T22:00:11.996647Z", + "iopub.status.idle": "2026-08-11T22:00:12.038639Z", + "shell.execute_reply": "2026-08-11T22:00:12.038143Z" } }, "outputs": [ @@ -584,7 +584,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/Users/hfdrake/code/sectionate/sectionate/transports.py:565: UserWarning: This section is open (its first and last waypoints are far apart), so it does not enclose a polygon and the inward/outward sense of `positive_in` is undefined. Sectionate applies the left-of-transect convention instead: with `positive_in=True` positive transport points to the LEFT of the section as it is traversed from the first to the last waypoint (set `positive_in=False` to flip to right-of-transect). Verify the sign matches your expectations, e.g. by plotting the section's normal vectors.\n", + "/Users/hfdrake/code/sectionate/sectionate/transports.py:738: UserWarning: This section is open (its first and last waypoints are far apart), so it does not enclose a polygon and the inward/outward sense of `positive_in` is undefined. Sectionate applies the left-of-transect convention instead: with `positive_in=True` positive transport points to the LEFT of the section as it is traversed from the first to the last waypoint (set `positive_in=False` to flip to right-of-transect). Verify the sign matches your expectations, e.g. by plotting the section's normal vectors.\n", " warnings.warn(\n" ] } @@ -598,10 +598,10 @@ "execution_count": 18, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:49.899254Z", - "iopub.status.busy": "2026-08-04T02:01:49.899198Z", - "iopub.status.idle": "2026-08-04T02:01:50.002350Z", - "shell.execute_reply": "2026-08-04T02:01:50.001979Z" + "iopub.execute_input": "2026-08-11T22:00:12.040004Z", + "iopub.status.busy": "2026-08-11T22:00:12.039923Z", + "iopub.status.idle": "2026-08-11T22:00:12.201133Z", + "shell.execute_reply": "2026-08-11T22:00:12.200689Z" } }, "outputs": [ @@ -658,10 +658,10 @@ "execution_count": 19, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:50.003392Z", - "iopub.status.busy": "2026-08-04T02:01:50.003333Z", - "iopub.status.idle": "2026-08-04T02:01:50.004876Z", - "shell.execute_reply": "2026-08-04T02:01:50.004533Z" + "iopub.execute_input": "2026-08-11T22:00:12.203255Z", + "iopub.status.busy": "2026-08-11T22:00:12.203156Z", + "iopub.status.idle": "2026-08-11T22:00:12.205176Z", + "shell.execute_reply": "2026-08-11T22:00:12.204722Z" } }, "outputs": [], @@ -681,10 +681,10 @@ "execution_count": 20, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:50.005800Z", - "iopub.status.busy": "2026-08-04T02:01:50.005751Z", - "iopub.status.idle": "2026-08-04T02:01:50.008642Z", - "shell.execute_reply": "2026-08-04T02:01:50.008176Z" + "iopub.execute_input": "2026-08-11T22:00:12.206455Z", + "iopub.status.busy": "2026-08-11T22:00:12.206350Z", + "iopub.status.idle": "2026-08-11T22:00:12.210185Z", + "shell.execute_reply": "2026-08-11T22:00:12.209816Z" } }, "outputs": [], @@ -709,10 +709,10 @@ "execution_count": 21, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:50.009563Z", - "iopub.status.busy": "2026-08-04T02:01:50.009509Z", - "iopub.status.idle": "2026-08-04T02:01:50.084867Z", - "shell.execute_reply": "2026-08-04T02:01:50.084458Z" + "iopub.execute_input": "2026-08-11T22:00:12.211415Z", + "iopub.status.busy": "2026-08-11T22:00:12.211336Z", + "iopub.status.idle": "2026-08-11T22:00:12.321927Z", + "shell.execute_reply": "2026-08-11T22:00:12.321499Z" } }, "outputs": [ diff --git a/examples/3_Labrador_convergence_CM4p25.ipynb b/examples/3_Labrador_convergence_CM4p25.ipynb index 2c28820..c075a44 100644 --- a/examples/3_Labrador_convergence_CM4p25.ipynb +++ b/examples/3_Labrador_convergence_CM4p25.ipynb @@ -21,10 +21,10 @@ "execution_count": 1, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:50.972072Z", - "iopub.status.busy": "2026-08-04T02:01:50.971866Z", - "iopub.status.idle": "2026-08-04T02:01:51.710112Z", - "shell.execute_reply": "2026-08-04T02:01:51.709792Z" + "iopub.execute_input": "2026-08-11T22:00:13.415468Z", + "iopub.status.busy": "2026-08-11T22:00:13.415374Z", + "iopub.status.idle": "2026-08-11T22:00:14.460836Z", + "shell.execute_reply": "2026-08-11T22:00:14.460267Z" } }, "outputs": [ @@ -32,7 +32,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Sectionate version: 0.4.0rc2.dev3\n" + "Sectionate version: 0.4.0rc3.dev19\n" ] } ], @@ -60,10 +60,10 @@ "execution_count": 2, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:51.711277Z", - "iopub.status.busy": "2026-08-04T02:01:51.711166Z", - "iopub.status.idle": "2026-08-04T02:01:51.832845Z", - "shell.execute_reply": "2026-08-04T02:01:51.832479Z" + "iopub.execute_input": "2026-08-11T22:00:14.462588Z", + "iopub.status.busy": "2026-08-11T22:00:14.462404Z", + "iopub.status.idle": "2026-08-11T22:00:14.653363Z", + "shell.execute_reply": "2026-08-11T22:00:14.652783Z" } }, "outputs": [ @@ -71,7 +71,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/Users/hfdrake/miniforge3/envs/docs_env_sectionate_fix-layer-interface-matching/lib/python3.14/site-packages/xgcm/grid.py:463: UserWarning: The north-fold (tripolar) boundary condition is experimental. Its API and numerical behavior may change in future releases, and it has not yet been validated across the full range of grid configurations. Please review results carefully and report any issues at https://github.com/xgcm/xgcm/issues.\n", + "/opt/anaconda3/envs/docs_env_sectionate_seam-corner-identity/lib/python3.14t/site-packages/xgcm/grid.py:463: UserWarning: The north-fold (tripolar) boundary condition is experimental. Its API and numerical behavior may change in future releases, and it has not yet been validated across the full range of grid configurations. Please review results carefully and report any issues at https://github.com/xgcm/xgcm/issues.\n", " warnings.warn(\n" ] } @@ -94,10 +94,10 @@ "execution_count": 3, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:51.833990Z", - "iopub.status.busy": "2026-08-04T02:01:51.833928Z", - "iopub.status.idle": "2026-08-04T02:01:51.835664Z", - "shell.execute_reply": "2026-08-04T02:01:51.835344Z" + "iopub.execute_input": "2026-08-11T22:00:14.655023Z", + "iopub.status.busy": "2026-08-11T22:00:14.654910Z", + "iopub.status.idle": "2026-08-11T22:00:14.657048Z", + "shell.execute_reply": "2026-08-11T22:00:14.656602Z" } }, "outputs": [], @@ -120,10 +120,10 @@ "execution_count": 4, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:51.836649Z", - "iopub.status.busy": "2026-08-04T02:01:51.836596Z", - "iopub.status.idle": "2026-08-04T02:01:51.864344Z", - "shell.execute_reply": "2026-08-04T02:01:51.864033Z" + "iopub.execute_input": "2026-08-11T22:00:14.658403Z", + "iopub.status.busy": "2026-08-11T22:00:14.658302Z", + "iopub.status.idle": "2026-08-11T22:00:14.699332Z", + "shell.execute_reply": "2026-08-11T22:00:14.698792Z" } }, "outputs": [], @@ -140,10 +140,10 @@ "execution_count": 5, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:51.865600Z", - "iopub.status.busy": "2026-08-04T02:01:51.865535Z", - "iopub.status.idle": "2026-08-04T02:01:51.980369Z", - "shell.execute_reply": "2026-08-04T02:01:51.979958Z" + "iopub.execute_input": "2026-08-11T22:00:14.700768Z", + "iopub.status.busy": "2026-08-11T22:00:14.700679Z", + "iopub.status.idle": "2026-08-11T22:00:14.863953Z", + "shell.execute_reply": "2026-08-11T22:00:14.863467Z" } }, "outputs": [ @@ -192,10 +192,10 @@ "execution_count": 6, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:51.981486Z", - "iopub.status.busy": "2026-08-04T02:01:51.981410Z", - "iopub.status.idle": "2026-08-04T02:01:52.020015Z", - "shell.execute_reply": "2026-08-04T02:01:52.019586Z" + "iopub.execute_input": "2026-08-11T22:00:14.865477Z", + "iopub.status.busy": "2026-08-11T22:00:14.865380Z", + "iopub.status.idle": "2026-08-11T22:00:14.915602Z", + "shell.execute_reply": "2026-08-11T22:00:14.915138Z" } }, "outputs": [], @@ -208,10 +208,10 @@ "execution_count": 7, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:52.021157Z", - "iopub.status.busy": "2026-08-04T02:01:52.021084Z", - "iopub.status.idle": "2026-08-04T02:01:52.062912Z", - "shell.execute_reply": "2026-08-04T02:01:52.062595Z" + "iopub.execute_input": "2026-08-11T22:00:14.917328Z", + "iopub.status.busy": "2026-08-11T22:00:14.917222Z", + "iopub.status.idle": "2026-08-11T22:00:14.974074Z", + "shell.execute_reply": "2026-08-11T22:00:14.973594Z" } }, "outputs": [], @@ -225,10 +225,10 @@ "execution_count": 8, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:52.064165Z", - "iopub.status.busy": "2026-08-04T02:01:52.064104Z", - "iopub.status.idle": "2026-08-04T02:01:52.141159Z", - "shell.execute_reply": "2026-08-04T02:01:52.140798Z" + "iopub.execute_input": "2026-08-11T22:00:14.975528Z", + "iopub.status.busy": "2026-08-11T22:00:14.975438Z", + "iopub.status.idle": "2026-08-11T22:00:15.085181Z", + "shell.execute_reply": "2026-08-11T22:00:15.084693Z" } }, "outputs": [ @@ -264,10 +264,10 @@ "execution_count": 9, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:52.142231Z", - "iopub.status.busy": "2026-08-04T02:01:52.142172Z", - "iopub.status.idle": "2026-08-04T02:01:52.278672Z", - "shell.execute_reply": "2026-08-04T02:01:52.278340Z" + "iopub.execute_input": "2026-08-11T22:00:15.086576Z", + "iopub.status.busy": "2026-08-11T22:00:15.086460Z", + "iopub.status.idle": "2026-08-11T22:00:15.305204Z", + "shell.execute_reply": "2026-08-11T22:00:15.304744Z" } }, "outputs": [ diff --git a/examples/4_sections_on_global_tripolar_grid.ipynb b/examples/4_sections_on_global_tripolar_grid.ipynb index a5b13f9..56dca7e 100644 --- a/examples/4_sections_on_global_tripolar_grid.ipynb +++ b/examples/4_sections_on_global_tripolar_grid.ipynb @@ -23,10 +23,10 @@ "execution_count": 1, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:53.171667Z", - "iopub.status.busy": "2026-08-04T02:01:53.171402Z", - "iopub.status.idle": "2026-08-04T02:01:55.635990Z", - "shell.execute_reply": "2026-08-04T02:01:55.635261Z" + "iopub.execute_input": "2026-08-11T22:00:16.454308Z", + "iopub.status.busy": "2026-08-11T22:00:16.454225Z", + "iopub.status.idle": "2026-08-11T22:00:21.627208Z", + "shell.execute_reply": "2026-08-11T22:00:21.626761Z" } }, "outputs": [ @@ -34,7 +34,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Sectionate version: 0.4.0rc2.dev3\n" + "Sectionate version: 0.4.0rc3.dev19\n" ] } ], @@ -66,10 +66,10 @@ "execution_count": 2, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:55.638398Z", - "iopub.status.busy": "2026-08-04T02:01:55.638098Z", - "iopub.status.idle": "2026-08-04T02:01:55.646099Z", - "shell.execute_reply": "2026-08-04T02:01:55.645540Z" + "iopub.execute_input": "2026-08-11T22:00:21.629045Z", + "iopub.status.busy": "2026-08-11T22:00:21.628833Z", + "iopub.status.idle": "2026-08-11T22:00:21.634331Z", + "shell.execute_reply": "2026-08-11T22:00:21.633798Z" } }, "outputs": [], @@ -143,10 +143,10 @@ "execution_count": 3, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:55.647677Z", - "iopub.status.busy": "2026-08-04T02:01:55.647564Z", - "iopub.status.idle": "2026-08-04T02:01:55.794379Z", - "shell.execute_reply": "2026-08-04T02:01:55.793983Z" + "iopub.execute_input": "2026-08-11T22:00:21.636087Z", + "iopub.status.busy": "2026-08-11T22:00:21.635977Z", + "iopub.status.idle": "2026-08-11T22:00:21.849657Z", + "shell.execute_reply": "2026-08-11T22:00:21.849113Z" } }, "outputs": [ @@ -154,7 +154,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/Users/hfdrake/miniforge3/envs/docs_env_sectionate_fix-layer-interface-matching/lib/python3.14/site-packages/xgcm/grid.py:463: UserWarning: The north-fold (tripolar) boundary condition is experimental. Its API and numerical behavior may change in future releases, and it has not yet been validated across the full range of grid configurations. Please review results carefully and report any issues at https://github.com/xgcm/xgcm/issues.\n", + "/opt/anaconda3/envs/docs_env_sectionate_seam-corner-identity/lib/python3.14t/site-packages/xgcm/grid.py:463: UserWarning: The north-fold (tripolar) boundary condition is experimental. Its API and numerical behavior may change in future releases, and it has not yet been validated across the full range of grid configurations. Please review results carefully and report any issues at https://github.com/xgcm/xgcm/issues.\n", " warnings.warn(\n" ] } @@ -177,10 +177,10 @@ "execution_count": 4, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:55.795569Z", - "iopub.status.busy": "2026-08-04T02:01:55.795503Z", - "iopub.status.idle": "2026-08-04T02:01:55.798005Z", - "shell.execute_reply": "2026-08-04T02:01:55.797632Z" + "iopub.execute_input": "2026-08-11T22:00:21.851197Z", + "iopub.status.busy": "2026-08-11T22:00:21.851086Z", + "iopub.status.idle": "2026-08-11T22:00:21.854081Z", + "shell.execute_reply": "2026-08-11T22:00:21.853719Z" } }, "outputs": [], @@ -203,10 +203,10 @@ "execution_count": 5, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:55.799062Z", - "iopub.status.busy": "2026-08-04T02:01:55.799000Z", - "iopub.status.idle": "2026-08-04T02:01:56.913038Z", - "shell.execute_reply": "2026-08-04T02:01:56.912614Z" + "iopub.execute_input": "2026-08-11T22:00:21.855391Z", + "iopub.status.busy": "2026-08-11T22:00:21.855292Z", + "iopub.status.idle": "2026-08-11T22:00:23.639865Z", + "shell.execute_reply": "2026-08-11T22:00:23.639437Z" } }, "outputs": [ @@ -251,10 +251,10 @@ "execution_count": 6, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:56.915401Z", - "iopub.status.busy": "2026-08-04T02:01:56.915332Z", - "iopub.status.idle": "2026-08-04T02:01:57.818895Z", - "shell.execute_reply": "2026-08-04T02:01:57.818522Z" + "iopub.execute_input": "2026-08-11T22:00:23.643439Z", + "iopub.status.busy": "2026-08-11T22:00:23.643336Z", + "iopub.status.idle": "2026-08-11T22:00:25.059781Z", + "shell.execute_reply": "2026-08-11T22:00:25.059195Z" } }, "outputs": [ @@ -305,10 +305,10 @@ "execution_count": 7, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:57.821173Z", - "iopub.status.busy": "2026-08-04T02:01:57.821093Z", - "iopub.status.idle": "2026-08-04T02:01:57.822690Z", - "shell.execute_reply": "2026-08-04T02:01:57.822323Z" + "iopub.execute_input": "2026-08-11T22:00:25.063686Z", + "iopub.status.busy": "2026-08-11T22:00:25.063577Z", + "iopub.status.idle": "2026-08-11T22:00:25.065440Z", + "shell.execute_reply": "2026-08-11T22:00:25.064973Z" } }, "outputs": [], @@ -328,10 +328,10 @@ "execution_count": 8, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:57.823644Z", - "iopub.status.busy": "2026-08-04T02:01:57.823593Z", - "iopub.status.idle": "2026-08-04T02:01:58.205365Z", - "shell.execute_reply": "2026-08-04T02:01:58.204922Z" + "iopub.execute_input": "2026-08-11T22:00:25.066724Z", + "iopub.status.busy": "2026-08-11T22:00:25.066629Z", + "iopub.status.idle": "2026-08-11T22:00:25.650199Z", + "shell.execute_reply": "2026-08-11T22:00:25.649650Z" } }, "outputs": [ @@ -370,10 +370,10 @@ "execution_count": 9, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:58.207583Z", - "iopub.status.busy": "2026-08-04T02:01:58.207497Z", - "iopub.status.idle": "2026-08-04T02:01:58.606342Z", - "shell.execute_reply": "2026-08-04T02:01:58.605892Z" + "iopub.execute_input": "2026-08-11T22:00:25.653650Z", + "iopub.status.busy": "2026-08-11T22:00:25.653533Z", + "iopub.status.idle": "2026-08-11T22:00:26.261303Z", + "shell.execute_reply": "2026-08-11T22:00:26.260735Z" } }, "outputs": [ @@ -412,10 +412,10 @@ "execution_count": 10, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:58.608545Z", - "iopub.status.busy": "2026-08-04T02:01:58.608465Z", - "iopub.status.idle": "2026-08-04T02:01:59.012843Z", - "shell.execute_reply": "2026-08-04T02:01:59.012464Z" + "iopub.execute_input": "2026-08-11T22:00:26.264608Z", + "iopub.status.busy": "2026-08-11T22:00:26.264493Z", + "iopub.status.idle": "2026-08-11T22:00:26.896808Z", + "shell.execute_reply": "2026-08-11T22:00:26.896231Z" } }, "outputs": [ @@ -456,10 +456,10 @@ "execution_count": 11, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:59.015039Z", - "iopub.status.busy": "2026-08-04T02:01:59.014959Z", - "iopub.status.idle": "2026-08-04T02:01:59.577602Z", - "shell.execute_reply": "2026-08-04T02:01:59.577160Z" + "iopub.execute_input": "2026-08-11T22:00:26.900306Z", + "iopub.status.busy": "2026-08-11T22:00:26.900201Z", + "iopub.status.idle": "2026-08-11T22:00:27.799804Z", + "shell.execute_reply": "2026-08-11T22:00:27.799205Z" } }, "outputs": [ @@ -499,10 +499,10 @@ "execution_count": 12, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:59.579759Z", - "iopub.status.busy": "2026-08-04T02:01:59.579676Z", - "iopub.status.idle": "2026-08-04T02:01:59.972989Z", - "shell.execute_reply": "2026-08-04T02:01:59.972545Z" + "iopub.execute_input": "2026-08-11T22:00:27.803030Z", + "iopub.status.busy": "2026-08-11T22:00:27.802848Z", + "iopub.status.idle": "2026-08-11T22:00:28.395964Z", + "shell.execute_reply": "2026-08-11T22:00:28.395515Z" } }, "outputs": [ @@ -543,10 +543,10 @@ "execution_count": 13, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:01:59.975278Z", - "iopub.status.busy": "2026-08-04T02:01:59.975201Z", - "iopub.status.idle": "2026-08-04T02:02:00.688345Z", - "shell.execute_reply": "2026-08-04T02:02:00.687998Z" + "iopub.execute_input": "2026-08-11T22:00:28.399468Z", + "iopub.status.busy": "2026-08-11T22:00:28.399364Z", + "iopub.status.idle": "2026-08-11T22:00:29.421742Z", + "shell.execute_reply": "2026-08-11T22:00:29.421063Z" } }, "outputs": [ @@ -588,10 +588,10 @@ "execution_count": 14, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:02:00.690842Z", - "iopub.status.busy": "2026-08-04T02:02:00.690767Z", - "iopub.status.idle": "2026-08-04T02:02:01.046170Z", - "shell.execute_reply": "2026-08-04T02:02:01.045845Z" + "iopub.execute_input": "2026-08-11T22:00:29.425127Z", + "iopub.status.busy": "2026-08-11T22:00:29.425017Z", + "iopub.status.idle": "2026-08-11T22:00:29.931978Z", + "shell.execute_reply": "2026-08-11T22:00:29.931422Z" } }, "outputs": [ @@ -623,10 +623,10 @@ "execution_count": 15, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:02:01.048283Z", - "iopub.status.busy": "2026-08-04T02:02:01.048205Z", - "iopub.status.idle": "2026-08-04T02:02:01.418272Z", - "shell.execute_reply": "2026-08-04T02:02:01.417864Z" + "iopub.execute_input": "2026-08-11T22:00:29.934492Z", + "iopub.status.busy": "2026-08-11T22:00:29.934400Z", + "iopub.status.idle": "2026-08-11T22:00:30.461991Z", + "shell.execute_reply": "2026-08-11T22:00:30.461474Z" } }, "outputs": [ @@ -658,10 +658,10 @@ "execution_count": 16, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:02:01.420269Z", - "iopub.status.busy": "2026-08-04T02:02:01.420194Z", - "iopub.status.idle": "2026-08-04T02:02:01.811776Z", - "shell.execute_reply": "2026-08-04T02:02:01.811376Z" + "iopub.execute_input": "2026-08-11T22:00:30.465247Z", + "iopub.status.busy": "2026-08-11T22:00:30.465149Z", + "iopub.status.idle": "2026-08-11T22:00:31.015005Z", + "shell.execute_reply": "2026-08-11T22:00:31.014464Z" } }, "outputs": [ @@ -701,10 +701,10 @@ "execution_count": 17, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:02:01.813687Z", - "iopub.status.busy": "2026-08-04T02:02:01.813606Z", - "iopub.status.idle": "2026-08-04T02:02:02.168679Z", - "shell.execute_reply": "2026-08-04T02:02:02.168282Z" + "iopub.execute_input": "2026-08-11T22:00:31.018242Z", + "iopub.status.busy": "2026-08-11T22:00:31.018142Z", + "iopub.status.idle": "2026-08-11T22:00:31.522163Z", + "shell.execute_reply": "2026-08-11T22:00:31.521659Z" } }, "outputs": [ diff --git a/examples/5_MOC_transports_ECCOv4r4.ipynb b/examples/5_MOC_transports_ECCOv4r4.ipynb index 37f91c3..b7bc675 100644 --- a/examples/5_MOC_transports_ECCOv4r4.ipynb +++ b/examples/5_MOC_transports_ECCOv4r4.ipynb @@ -33,10 +33,10 @@ "execution_count": 1, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:09:56.567055Z", - "iopub.status.busy": "2026-08-04T02:09:56.566859Z", - "iopub.status.idle": "2026-08-04T02:09:57.455411Z", - "shell.execute_reply": "2026-08-04T02:09:57.455061Z" + "iopub.execute_input": "2026-08-11T22:00:32.833089Z", + "iopub.status.busy": "2026-08-11T22:00:32.832749Z", + "iopub.status.idle": "2026-08-11T22:00:34.026277Z", + "shell.execute_reply": "2026-08-11T22:00:34.025847Z" } }, "outputs": [ @@ -44,7 +44,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Sectionate version: 0.4.0rc2.dev3\n" + "Sectionate version: 0.4.0rc3.dev19\n" ] } ], @@ -77,10 +77,10 @@ "execution_count": 2, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:09:57.469675Z", - "iopub.status.busy": "2026-08-04T02:09:57.469513Z", - "iopub.status.idle": "2026-08-04T02:09:58.262430Z", - "shell.execute_reply": "2026-08-04T02:09:58.262018Z" + "iopub.execute_input": "2026-08-11T22:00:34.027886Z", + "iopub.status.busy": "2026-08-11T22:00:34.027707Z", + "iopub.status.idle": "2026-08-11T22:00:34.796161Z", + "shell.execute_reply": "2026-08-11T22:00:34.795533Z" } }, "outputs": [ @@ -129,10 +129,10 @@ "execution_count": 3, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:09:58.263547Z", - "iopub.status.busy": "2026-08-04T02:09:58.263491Z", - "iopub.status.idle": "2026-08-04T02:09:58.265337Z", - "shell.execute_reply": "2026-08-04T02:09:58.265038Z" + "iopub.execute_input": "2026-08-11T22:00:34.797548Z", + "iopub.status.busy": "2026-08-11T22:00:34.797454Z", + "iopub.status.idle": "2026-08-11T22:00:34.799881Z", + "shell.execute_reply": "2026-08-11T22:00:34.799433Z" } }, "outputs": [], @@ -172,10 +172,10 @@ "execution_count": 4, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:09:58.266314Z", - "iopub.status.busy": "2026-08-04T02:09:58.266260Z", - "iopub.status.idle": "2026-08-04T02:10:02.515464Z", - "shell.execute_reply": "2026-08-04T02:10:02.514970Z" + "iopub.execute_input": "2026-08-11T22:00:34.801213Z", + "iopub.status.busy": "2026-08-11T22:00:34.801121Z", + "iopub.status.idle": "2026-08-11T22:00:40.988516Z", + "shell.execute_reply": "2026-08-11T22:00:40.987554Z" } }, "outputs": [ @@ -245,10 +245,10 @@ "execution_count": 5, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:10:02.518605Z", - "iopub.status.busy": "2026-08-04T02:10:02.518471Z", - "iopub.status.idle": "2026-08-04T02:10:17.559173Z", - "shell.execute_reply": "2026-08-04T02:10:17.558828Z" + "iopub.execute_input": "2026-08-11T22:00:40.995420Z", + "iopub.status.busy": "2026-08-11T22:00:40.994818Z", + "iopub.status.idle": "2026-08-11T22:01:39.771745Z", + "shell.execute_reply": "2026-08-11T22:01:39.770180Z" } }, "outputs": [ @@ -293,10 +293,10 @@ "execution_count": 6, "metadata": { "execution": { - "iopub.execute_input": "2026-08-04T02:10:17.560375Z", - "iopub.status.busy": "2026-08-04T02:10:17.560315Z", - "iopub.status.idle": "2026-08-04T02:10:17.727353Z", - "shell.execute_reply": "2026-08-04T02:10:17.726950Z" + "iopub.execute_input": "2026-08-11T22:01:39.779584Z", + "iopub.status.busy": "2026-08-11T22:01:39.779441Z", + "iopub.status.idle": "2026-08-11T22:01:40.099171Z", + "shell.execute_reply": "2026-08-11T22:01:40.098610Z" } }, "outputs": [ diff --git a/sectionate/tests/test_seam_representations.py b/sectionate/tests/test_seam_representations.py new file mode 100644 index 0000000..50f57ea --- /dev/null +++ b/sectionate/tests/test_seam_representations.py @@ -0,0 +1,208 @@ +""" +Tests for the seam representation contract of `uvindices_from_qindices`. + +A seam gives one physical corner more than one index, so a section crossing it can be +written two ways: *seam-explicit*, carrying both indices (the step between them spans no +grid cell), or *twin-free*, carrying only one. Both describe the same path, and sectionate +accepts both -- which is what lets a section arrive from `grid_section`, from a mask traced +on the grid (`regionate.boundaries`, whose single-tile back-end deliberately keeps both +coincident corners at a seam junction), or from saved indices reloaded off disk. + +The contract these pin down: + +* **input is liberal** -- either representation is accepted and yields the *same* faces; +* **output is strict** -- the returned faces are all real velocity faces, so a step between + two indices of one physical corner emits nothing; +* **"q" relates the two** -- faces and corner steps are therefore not one-to-one, and "q" + says which corner step each face came from. +""" +import warnings + +import numpy as np +import pytest + +from sectionate.section import grid_section, distance_on_unit_sphere, COINCIDENT_TOLERANCE_M +from sectionate.transports import convergent_transport, uvindices_from_qindices +from sectionate.gridutils import get_geo_corners + +from test_section_fold import ( + _pinched_fold_grid, + _add_streamfunction_transports, + _xgcm_supports_fold, +) +from test_section_multitile import _two_face_transport_grid + +pytestmark = pytest.mark.skipif( + not _xgcm_supports_fold(), + reason="installed xgcm lacks the bipolar north-fold boundary", +) + + +def _corner_coords(grid, i_c, j_c, f_c=None): + geo = get_geo_corners(grid) + glon, glat = np.asarray(geo["X"].values), np.asarray(geo["Y"].values) + if f_c is None: + return glon[j_c, i_c], glat[j_c, i_c] + return glon[f_c, j_c, i_c], glat[f_c, j_c, i_c] + + +def _coincident_steps(grid, i_c, j_c, f_c=None): + """Indices k where corners k and k+1 are the same physical point.""" + lon, lat = _corner_coords(grid, i_c, j_c, f_c) + d = distance_on_unit_sphere(lon[:-1], lat[:-1], lon[1:], lat[1:]) + return set(np.flatnonzero(d < COINCIDENT_TOLERANCE_M).tolist()) + + +def _to_twin_free(grid, i_c, j_c): + """Delete the second index of every coincident consecutive pair.""" + drop = {k + 1 for k in _coincident_steps(grid, i_c, j_c)} + keep = [k for k in range(len(i_c)) if k not in drop] + return np.asarray(i_c)[keep], np.asarray(j_c)[keep] + + +def _wrap_section(grid): + """A section crossing the periodic-X seam, which `grid_section` returns seam-explicit.""" + geo = get_geo_corners(grid) + glon, glat = np.asarray(geo["X"].values), np.asarray(geo["Y"].values) + i_c, j_c, _, _ = grid_section( + grid, + [float(glon[1, 1]), float(glon[1, 13])], + [float(glat[1, 1]), float(glat[1, 13])], + ) + return i_c, j_c + + +def test_both_seam_representations_give_the_same_faces(): + """The two ways of writing a seam crossing name identical velocity faces. + + `grid_section` returns the periodic-X crossing seam-explicit, so deleting the duplicate + corner produces the twin-free form of the very same path. Nothing downstream may depend + on which one it was handed. + """ + grid = _pinched_fold_grid() + i_x, j_x = _wrap_section(grid) + assert _coincident_steps(grid, i_x, j_x), "fixture must actually cross the seam" + + i_f, j_f = _to_twin_free(grid, i_x, j_x) + assert len(i_f) < len(i_x) + + explicit = uvindices_from_qindices(grid, i_x, j_x) + free = uvindices_from_qindices(grid, i_f, j_f) + + for key in ("var", "i", "j", "Xinc", "Yinc"): + np.testing.assert_array_equal( + explicit[key], free[key], err_msg=f"{key} differs between representations" + ) + + +def _mirror_pair_boundary(): + """The boundary of two cells that are mirror images across the fold seam. + + This is the shape a mask tracer produces for a region joined only through the fold: the + loop hands off between the seam corner's two indices twice, and each hand-off is a + coincident pair. On `_pinched_fold_grid(16, 6)` the seam row is j=5, mirrored under + i <-> 16-i, so cells (4, 1) and (4, 14) are a mirror pair sharing the seam face. + """ + j_c = np.array([5, 4, 4, 5, 5, 4, 4, 5, 5]) + i_c = np.array([1, 1, 2, 2, 14, 14, 15, 15, 1]) + return i_c, j_c + + +def test_mask_traced_fold_loop_is_accepted_and_counts_each_face_once(): + """A seam-explicit loop straddling the fold resolves to its real faces. + + Before physical coincidence was tested ahead of index adjacency, the hand-off pair + (j=5, i=2) -> (j=5, i=14) -- already a zero-length twin edge, but 12 apart in the index + lattice -- was sent down the twin-splicing path and raised. + """ + grid = _pinched_fold_grid() + i_c, j_c = _mirror_pair_boundary() + assert _coincident_steps(grid, i_c, j_c) == {3, 7} + + uv = uvindices_from_qindices(grid, i_c, j_c) + + # 8 corner steps, 2 of them hand-offs -> 6 faces: three per cell. The cells' shared + # seam face is interior to the pair and is correctly not among them. + assert uv["var"].size == 6 + faces = list(zip(uv["var"].tolist(), uv["i"].tolist(), uv["j"].tolist())) + assert len(set(faces)) == 6, f"a face was named twice: {faces}" + assert all(v in ("U", "V") for v in uv["var"]) + + +def test_mask_traced_fold_loop_carries_zero_net_transport(): + """The loop encloses two cells of an exactly non-divergent flow, so the transport + convergence over it is zero -- the discrete divergence theorem across the fold seam. + Naming the mirror image of a crossed face, or counting one twice, breaks this.""" + grid = _pinched_fold_grid() + _, tgrid = _add_streamfunction_transports(grid) + i_c, j_c = _mirror_pair_boundary() + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + out = convergent_transport(tgrid, i_c, j_c, utr="utest", vtr="vtest") + total = float(out["conv_mass_transport"].sum().values) + + assert abs(total) < 1.0e-12, f"net convergence {total} should vanish" + + +def _assert_q_contract(grid, i_c, j_c, f_c=None): + uv = uvindices_from_qindices(grid, i_c, j_c, f_c=f_c) + q = uv["q"] + n = len(i_c) + + assert q.size == uv["var"].size + assert np.all(np.diff(q) > 0), f"q must be strictly increasing, got {q}" + assert q.min() >= 0 and q.max() + 1 <= n - 1, "face k must span corners q[k], q[k]+1" + + # The steps q skips are exactly the ones that carried no face. + skipped = set(range(n - 1)) - set(q.tolist()) + assert skipped == _coincident_steps(grid, i_c, j_c, f_c), ( + f"q skipped {sorted(skipped)}, but the coincident steps are " + f"{sorted(_coincident_steps(grid, i_c, j_c, f_c))}" + ) + assert q.size == (n - 1) - len(skipped) + return uv + + +def test_q_maps_faces_back_to_corner_steps_across_a_periodic_seam(): + grid = _pinched_fold_grid() + i_c, j_c = _wrap_section(grid) + uv = _assert_q_contract(grid, i_c, j_c) + assert uv["q"].size < len(i_c) - 1, "the seam hand-off must leave a gap in q" + + +def test_q_maps_faces_back_to_corner_steps_on_a_mask_traced_loop(): + grid = _pinched_fold_grid() + i_c, j_c = _mirror_pair_boundary() + uv = _assert_q_contract(grid, i_c, j_c) + assert uv["q"].tolist() == [0, 1, 2, 4, 5, 6] + + +def test_q_is_the_identity_when_no_corner_repeats(): + """A section that touches no seam has one face per corner step, and q says so.""" + grid = _pinched_fold_grid() + geo = get_geo_corners(grid) + glon, glat = np.asarray(geo["X"].values), np.asarray(geo["Y"].values) + i_c, j_c, _, _ = grid_section( + grid, + [float(glon[0, 3]), float(glon[0, 6])], + [float(glat[0, 3]), float(glat[0, 6])], + ) + assert not _coincident_steps(grid, i_c, j_c) + uv = _assert_q_contract(grid, i_c, j_c) + np.testing.assert_array_equal(uv["q"], np.arange(len(i_c) - 1)) + + +def test_q_contract_holds_on_a_multitile_grid(): + """Multi-tile corners resolve through the corner-node graph rather than by coordinate + coincidence, but the reported relation between faces and corner steps is the same.""" + grid = _two_face_transport_grid() + geo = get_geo_corners(grid) + glon, glat = np.asarray(geo["X"].values), np.asarray(geo["Y"].values) + i_c, j_c, f_c, _, _ = grid_section( + grid, + [float(glon[0, 1, 1]), float(glon[1, 1, 2])], + [float(glat[0, 1, 1]), float(glat[1, 1, 2])], + ) + assert len(set(f_c.tolist())) > 1, "fixture must actually cross the face seam" + _assert_q_contract(grid, i_c, j_c, f_c=f_c) diff --git a/sectionate/tests/test_section_fold.py b/sectionate/tests/test_section_fold.py index a6f681d..b46d9a7 100644 --- a/sectionate/tests/test_section_fold.py +++ b/sectionate/tests/test_section_fold.py @@ -5,6 +5,7 @@ seam connectivity from xgcm's padding with no fold-specific code of its own. """ import os +import warnings import numpy as np import pytest @@ -12,7 +13,8 @@ import xgcm from sectionate.gridutils import get_geo_corners, build_neighbor_maps -from sectionate.section import grid_section +from sectionate.section import grid_section, distance_on_unit_sphere +from sectionate.transports import convergent_transport, uvindices_from_qindices def _fold_grid(nx=8, ny=5, pivot="corner"): @@ -111,3 +113,216 @@ def test_fold_section_crosses_arctic_on_real_grid(): assert np.all(np.abs(np.asarray(lats_c) - 80.0) < 1.0) # It reaches the northern fold row. assert j_c.max() == ds.yq.size - 1 + + +# --------------------------------------------------------------------------- +# Transports across the fold seam. +# +# The seam row identifies each corner with its mirrored column, so the same +# physical corner carries two indices. Crossing the seam, the walk steps from one +# of them to a corner adjacent to the *other*, which is a physical edge but not an +# index-adjacent one. The tests below pin down that such a crossing names the +# velocity face the section actually crossed, with the right sign, exactly once. +# +# They use a streamfunction fixture: psi is a single-valued function of *physical* +# position, so seam twins get identical psi, and transports are built from +# index-space differences of psi. That makes the flow exactly non-divergent and the +# answer for any section known analytically -- transport = psi(end) - psi(start) -- +# while automatically reproducing the duplicated, sign-flipped seam row of an +# ORCA-style fold (NEMO 4.2 manual, Appendix E) without hard-coding the flip. +# --------------------------------------------------------------------------- + + +def _pinched_fold_grid(nxh=16, ny=6): + """ + A single-tile C-grid that genuinely carries the corner-pivot fold identity. + + Picture a cylinder (X-periodic) whose cross-section is pinched shut as `j` rises: + at row `j` the corners sit on an ellipse with semi-axes ``(a_j, b_j)``, and + ``b_{ny-1} = 0``, so the top row collapses onto a segment traversed out and back. + Under ``i <-> nxh - i`` that row is therefore *exactly* mirror-symmetric, which is + the corner-pivot fold identity, with its two fold poles at ``i = 0`` and + ``i = nxh/2``. ``a_j`` grows with `j` so the ellipses nest and every corner below + the seam stays distinct. + + Unlike `_fold_grid`, whose plain lat/lon coordinates do not carry the fold its + metadata declares, this one is geometrically consistent with the fold, so + transports across the seam are well defined. + """ + i_q, j_q = np.arange(nxh + 1), np.arange(ny) + i_h, j_h = np.arange(nxh) + 0.5, np.arange(ny - 1) + 0.5 + + def lonlat(jj, ii): + theta = 2.0 * np.pi * np.asarray(ii) / nxh + t = np.asarray(jj) / (ny - 1) + x = (1.0 + 0.5 * t) * np.cos(theta) + y = (1.0 - t) * np.sin(theta) + return 25.0 * x, 45.0 + 18.0 * y + + lon_c, lat_c = lonlat(*np.meshgrid(j_q, i_q, indexing="ij")) + lon_h, lat_h = lonlat(*np.meshgrid(j_h, i_h, indexing="ij")) + + ds = xr.Dataset(coords={ + "xq": i_q, "yq": j_q, "xh": i_h, "yh": j_h, + "geolon_c": (("yq", "xq"), lon_c), "geolat_c": (("yq", "xq"), lat_c), + "geolon": (("yh", "xh"), lon_h), "geolat": (("yh", "xh"), lat_h), + }) + return xgcm.Grid( + ds, + coords={"X": {"center": "xh", "outer": "xq"}, "Y": {"center": "yh", "outer": "yq"}}, + padding={"X": "periodic", "Y": {"fold": "corner"}}, + autoparse_metadata=False, + ) + + +def _streamfunction(lon, lat): + """A degree-1 spherical harmonic: smooth everywhere (poles included) and + single-valued in physical space, so seam twins receive identical values.""" + lo, la = np.deg2rad(lon), np.deg2rad(lat) + return (0.3 * np.cos(la) * np.cos(lo) + - 0.8 * np.cos(la) * np.sin(lo) + + 0.5 * np.sin(la)) + + +def _add_streamfunction_transports(grid, utr="utest", vtr="vtest"): + """Attach exactly non-divergent transports derived from `_streamfunction`. + + Returns the corner-array streamfunction and a grid carrying the transports. + `umo[j,i] = psi[j,i] - psi[j+1,i]` and `vmo[j,i] = psi[j,i+1] - psi[j,i]` make + the transport across any section telescope to `psi(end) - psi(start)`. + """ + geo = get_geo_corners(grid) + psi = _streamfunction(np.asarray(geo["X"].values), np.asarray(geo["Y"].values)) + ydim, xdim = geo["X"].dims[-2], geo["X"].dims[-1] + ycenter = grid.axes["Y"].coords["center"] + xcenter = grid.axes["X"].coords["center"] + + ds = grid._ds.copy() + ds[utr] = xr.DataArray(psi[:-1, :] - psi[1:, :], dims=(ycenter, xdim)) + ds[vtr] = xr.DataArray(psi[:, 1:] - psi[:, :-1], dims=(ydim, xcenter)) + regridded = xgcm.Grid( + ds, + coords={ax: dict(grid.axes[ax].coords) for ax in ("X", "Y")}, + padding={ax: grid.axes[ax].padding for ax in ("X", "Y")}, + autoparse_metadata=False, + ) + return psi, regridded + + +def _section_transport(grid, i_c, j_c, utr="utest", vtr="vtest"): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") # open sections warn about `positive_in` + out = convergent_transport(grid, i_c, j_c, utr=utr, vtr=vtr) + return float(out["conv_mass_transport"].sum().values) + + +def test_pinched_fold_grid_carries_the_fold_symmetry(): + """The fixture really does duplicate its seam row, and the streamfunction + transports really do reproduce the ORCA duplicate-with-sign-flip structure.""" + nxh, ny = 16, 6 + grid = _pinched_fold_grid(nxh, ny) + geo = get_geo_corners(grid) + lon, lat = np.asarray(geo["X"].values), np.asarray(geo["Y"].values) + top = ny - 1 + + # Seam row: corner i and its mirror nxh-i are the same physical point ... + mirrored = distance_on_unit_sphere(lon[top], lat[top], lon[top, ::-1], lat[top, ::-1]) + assert np.nanmax(mirrored) < 1.0e-6 + # ... while the row just below it holds two genuinely distinct halves. + below = distance_on_unit_sphere(lon[top-1], lat[top-1], lon[top-1, ::-1], lat[top-1, ::-1]) + assert below[1:nxh // 2].min() > 1.0e4 + + psi, tgrid = _add_streamfunction_transports(grid) + v = tgrid._ds["vtest"].values + # The seam row's V transports are duplicates carrying a sign flip, and the + # transports are exactly non-divergent -- neither was put in by hand. + assert np.abs(v[top] + v[top, ::-1]).max() < 1.0e-12 * max(np.abs(v[top]).max(), 1.0) + u = tgrid._ds["utest"].values + conv = (u[:, :-1] - u[:, 1:]) + (v[:-1, :] - v[1:, :]) + assert np.abs(conv).max() < 1.0e-12 + + +def test_transport_crossing_the_fold_seam_is_exact(): + """A path that goes up to the seam and back down the mirrored side must transport + exactly `psi(end) - psi(start)`, using each velocity face once. + + Before the seam-twin fix this crossing was attributed to the mirror column's face + with the wrong sign, so the total was wrong by ~2x that face's transport. + """ + nxh, ny = 16, 6 + top = ny - 1 + psi, grid = _add_streamfunction_transports(_pinched_fold_grid(nxh, ny)) + + for i0 in (3, 5, 11): + mirror = nxh - i0 + # up to the seam at column i0, across the fold, back down at the mirror column + j_c = np.array([top - 2, top - 1, top, top - 1, top - 2]) + i_c = np.array([i0, i0, i0, mirror, mirror]) + + got = _section_transport(grid, i_c, j_c) + want = psi[j_c[-1], i_c[-1]] - psi[j_c[0], i_c[0]] + assert got == pytest.approx(want, abs=1.0e-12), f"column {i0}" + + # no velocity face is used twice, and the zero-length twin edge emits none + uv = uvindices_from_qindices(grid, i_c, j_c) + faces = list(zip(uv["var"].tolist(), uv["i"].tolist(), uv["j"].tolist())) + assert len(faces) == len(set(faces)) == 4 + + +def test_transport_along_the_fold_seam_is_exact(): + """Running *along* the duplicated, sign-flipped seam row needs no fold-specific + handling: the traversal direction and the stored velocity share an index frame, + so the flip cancels. This already held before the fix and must keep holding.""" + nxh, ny = 16, 6 + top = ny - 1 + psi, grid = _add_streamfunction_transports(_pinched_fold_grid(nxh, ny)) + + for i0, i1 in ((2, 6), (6, 12), (4, 12)): # last one straddles the fold pivot + i_c = np.concatenate([np.arange(i0, i1 + 1), [i1]]) + j_c = np.concatenate([np.full(i1 - i0 + 1, top), [top - 1]]) + got = _section_transport(grid, i_c, j_c) + want = psi[j_c[-1], i_c[-1]] - psi[j_c[0], i_c[0]] + assert got == pytest.approx(want, abs=1.0e-12), f"columns {i0}..{i1}" + + +def test_fold_crossing_raises_when_coordinates_lack_the_fold_symmetry(): + """`_fold_grid` declares a fold its plain lat/lon coordinates do not carry, so the + seam has no twins and the crossing edge is genuinely undefined. Say so, rather + than silently attributing it to some other face.""" + nxh, ny = 8, 5 + grid = _fold_grid(nxh, ny) + j_c = np.array([ny - 2, ny - 1, ny - 2]) + i_c = np.array([3, 3, nxh - 3]) + with pytest.raises(ValueError, match="fold"): + uvindices_from_qindices(grid, i_c, j_c) + + +@pytest.mark.skipif( + not os.path.exists(_mom6_example_path()), + reason="MOM6 example dataset not present (download via examples/load_example_model_grid.py)", +) +def test_fold_crossing_transport_is_exact_on_real_tripolar_grid(): + """The real-grid counterpart: traced sections that cross the tripolar seam -- + a closed Arctic latitude circle and an open path over the pole -- transport + exactly `psi(end) - psi(start)`, and use no velocity face twice.""" + ds = xr.open_dataset(_mom6_example_path()) + coords = {"X": {"center": "xh", "outer": "xq"}, "Y": {"center": "yh", "outer": "yq"}} + grid = xgcm.Grid(ds, coords=coords, padding={"X": "periodic", "Y": {"fold": "corner"}}, + autoparse_metadata=False) + psi, tgrid = _add_streamfunction_transports(grid) + + sections = { + "latitude circle 80N": (np.arange(0.0, 360.0 + 5.0, 5.0), [80.0] * 73, "latitude circle"), + "meridional over pole": ([100.0, 100.0, 280.0, 280.0], [60.0, 89.0, 89.0, 60.0], "great circle"), + } + for name, (lons, lats, curve) in sections.items(): + i_c, j_c, _, _ = grid_section(grid, lons, lats, curve=curve) + assert j_c.max() == ds.yq.size - 1, f"{name} never reached the seam" + + got = _section_transport(tgrid, i_c, j_c) + want = psi[j_c[-1], i_c[-1]] - psi[j_c[0], i_c[0]] + assert got == pytest.approx(want, abs=1.0e-12), name + + uv = uvindices_from_qindices(grid, i_c, j_c) + faces = list(zip(uv["var"].tolist(), uv["i"].tolist(), uv["j"].tolist())) + assert len(faces) == len(set(faces)), f"{name} counted a velocity face twice" diff --git a/sectionate/tracers.py b/sectionate/tracers.py index 387cfdc..272d5b6 100644 --- a/sectionate/tracers.py +++ b/sectionate/tracers.py @@ -38,7 +38,8 @@ def extract_tracer( RETURNS: -------- - xarray.DataArray with data interpolated to the U and V points along the section. + xarray.DataArray with data interpolated to the U and V points along the section, + carrying a "q" coordinate giving the step of (i_c, j_c) each point came from. """ da=grid._ds[name] @@ -79,4 +80,11 @@ def extract_tracer( tracer.name = da.name tracer.attrs = da.attrs + # Same "q" as `convergent_transport`: the step of the caller's corner arrays each face + # came from, so a tracer section and a transport section stay relatable to the corners + # even though faces and corner steps are not one-to-one. + tracer = tracer.assign_coords({ + "q": xr.DataArray(uvindices["q"], dims=sect_coord) + }) + return tracer diff --git a/sectionate/transports.py b/sectionate/transports.py index 6b3aa7a..6e3dc4d 100644 --- a/sectionate/transports.py +++ b/sectionate/transports.py @@ -138,6 +138,135 @@ def _uv_for_edge(A, B, neighbor_maps, offset, ranges, glon, glat): return "0", 0, 0, int(fB), 0 +def _axis_periodic(grid, axis): + """Whether `axis` wraps, so its first and last index columns are neighbors.""" + ax = grid.axes.get(axis) + return (ax is not None) and (ax.padding == "periodic") + + +def _index_adjacent(ja, ia, jb, ib, ny, nx, y_periodic, x_periodic): + """ + Whether two corners are neighbors in the *raw index* lattice -- the assumption the + single-tile velocity-face arithmetic in `uvindices_from_qindices` is built on. A + wrap across a periodic axis counts, since that arithmetic handles it explicitly. + """ + dj, di = int(jb) - int(ja), int(ib) - int(ia) + if dj == 0 and di == 0: + return True # degenerate; dropped downstream + if dj == 0: + return abs(di) == 1 or (x_periodic and abs(di) == nx - 1) + if di == 0: + return abs(dj) == 1 or (y_periodic and abs(dj) == ny - 1) + return False + + +def _adjacent_candidates(j, i, ny, nx, y_periodic, x_periodic): + """Every corner index adjacent to (j,i) in the raw index lattice, in a fixed order.""" + cand = [(j, i - 1), (j, i + 1), (j - 1, i), (j + 1, i)] + if x_periodic: + cand += [(j, 0), (j, nx - 1)] + if y_periodic: + cand += [(0, i), (ny - 1, i)] + out = [] + for (jj, ii) in cand: + if ((0 <= jj < ny) and (0 <= ii < nx) and ((jj, ii) not in out) + and _index_adjacent(j, i, jj, ii, ny, nx, y_periodic, x_periodic)): + out.append((jj, ii)) + return out + + +def _insert_seam_twins(i_c, j_c, glon, glat, grid): + """ + Splice in seam twins so that every consecutive pair of corners is index-adjacent. + + A seam gives one physical corner two index representations. The bipolar north fold is + the single-tile case: it identifies each seam-row corner with its mirrored column, so + the walk can step from one representation straight to a corner adjacent to the *other*. + That pair is a real physical edge but not an index-adjacent one, and the arithmetic in + `uvindices_from_qindices` reads the velocity column off the source corner -- which is + then the wrong representative. It names the mirror image of the face the section + actually crossed, and (because the mirrored half runs the opposite way in index space) + with the opposite sign. On a closed fold-crossing section that also emits one face + twice, double-counting it. + + Splicing the twin in restores the invariant the arithmetic needs: the crossing becomes + one ordinary edge plus one zero-length twin edge, and zero-length edges are dropped at + the end of `uvindices_from_qindices`, so nothing is counted twice. This is the + single-tile counterpart of the corner canonicalization `_OuterTopology` performs for + multi-tile grids; NEMO reaches the same place from the other direction, carrying the + duplicated points and masking them out of transports with `umaskutil`/`vmaskutil`. + + Runs on every single-tile grid and is a no-op unless a step actually skips the twin, so + grids without seam duplicates (and fold sections that never cross the seam) are untouched. + + A pair that is *already* the two indices of one physical corner needs no splicing: it is + the zero-length twin edge this function otherwise creates, so it is passed straight through + to the drop. That case has to be recognised before index adjacency is tested, because the + two indices of a seam corner are in general far apart in the index lattice -- a bipolar + fold seam identifies column i with column nx-i. Sections traced on a cell mask arrive in + exactly this form (`regionate.boundaries` keeps both coincident corners at a seam + junction), as do sections reloaded from saved indices. + + RETURNS: + -------- + i_c, j_c: np.ndarray + The corner chain with seam twins spliced in. + src: np.ndarray of int + For each returned corner, the index into the *input* arrays that it came from. A + spliced twin is attributed to the corner it was derived from, i.e. to the input step + it belongs to, so that the velocity face between returned corners m and m+1 belongs to + input step `src[m]`. This is what `uvindices_from_qindices` reports as "q". + """ + if len(i_c) < 2: + return np.asarray(i_c), np.asarray(j_c), np.arange(len(i_c)) + ny, nx = glon.shape + y_periodic, x_periodic = _axis_periodic(grid, "Y"), _axis_periodic(grid, "X") + + def coincident(ja, ia, jb, ib): + return distance_on_unit_sphere( + glon[ja, ia], glat[ja, ia], glon[jb, ib], glat[jb, ib] + ) < COINCIDENT_TOLERANCE_M + + out_i, out_j, src = [int(i_c[0])], [int(j_c[0])], [0] + for k in range(len(i_c) - 1): + ja, ia = int(j_c[k]), int(i_c[k]) + jb, ib = int(j_c[k + 1]), int(i_c[k + 1]) + if coincident(ja, ia, jb, ib): + pass # already a zero-length twin edge; emits no velocity face + elif not _index_adjacent(ja, ia, jb, ib, ny, nx, y_periodic, x_periodic): + # Either a twin of A adjacent to B, or a twin of B adjacent to A, turns the + # step into two index-adjacent ones. Which of the two exists depends on + # whether the section is entering or leaving the seam. + twin = next( + (t for t in _adjacent_candidates(jb, ib, ny, nx, y_periodic, x_periodic) + if coincident(*t, ja, ia)), + None, + ) + if twin is None: + twin = next( + (t for t in _adjacent_candidates(ja, ia, ny, nx, y_periodic, x_periodic) + if coincident(*t, jb, ib)), + None, + ) + if twin is None: + raise ValueError( + f"Section steps between corners (j={ja}, i={ia}) and (j={jb}, i={ib}), " + "which are neither adjacent in index space nor separated by a seam twin, " + "so the velocity face between them is undefined. On a grid with a bipolar " + "fold this means the corner coordinates do not carry the fold symmetry the " + "grid declares (the seam row's mirrored columns should hold identical " + "coordinates)." + ) + out_j.append(twin[0]) + out_i.append(twin[1]) + src.append(k) + out_j.append(jb) + out_i.append(ib) + src.append(k + 1) + + return np.asarray(out_i), np.asarray(out_j), np.asarray(src) + + def uvindices_from_qindices(grid, i_c, j_c, f_c=None): """ Find the `grid` indices of the N-1 velocity points defined by the consecutive indices of @@ -165,6 +294,11 @@ def uvindices_from_qindices(grid, i_c, j_c, f_c=None): - "j" : "Y"-dimension index of appropriate "U" or "V" velocity - "Yinc" : True if point was passed through while going in positive "j"-index direction - "Xinc" : True if point was passed through while going in positive "i"-index direction + - "q" : index into the *input* corner arrays of the corner each velocity face starts + at, so face k spans corners `q[k]` and `q[k]+1`. Because a consecutive pair of + corners at the same physical point defines no velocity face (see below), there are + in general fewer faces than corner steps and "q" is what relates the two: it is + strictly increasing, and skips exactly the steps that carried no face. For multi-tile grids (`f_c` given) the dict additionally contains: - "face" : face index of the velocity point - "Lsign" : +1 if the velocity's positive direction points left of the section's @@ -190,10 +324,22 @@ def uvindices_from_qindices(grid, i_c, j_c, f_c=None): nodes = ot.node_id[f_c, j_c + ot.t, i_c + ot.t] if (nodes < 0).any(): raise ValueError("Section contains indices that are not grid corners.") + # Keep the LAST index of each run of repeated corners, not the first. Every index of + # a node rewrites to the same canonical native one, so which survives does not change + # a single emitted corner -- but it does decide which input step each face is + # attributed to in "q", and the real face is the one leaving the run, not entering it. keep = np.ones(nodes.size, dtype=bool) - keep[1:] = nodes[1:] != nodes[:-1] + keep[:-1] = nodes[:-1] != nodes[1:] + src = np.flatnonzero(keep) nat = ot.node_native[nodes[keep]] f_c, j_c, i_c = nat[:, 0].copy(), nat[:, 1].copy(), nat[:, 2].copy() + else: + # Single-tile grids have no face topology to canonicalize corners through, but + # they can still carry seam duplicates: a bipolar north fold identifies every + # seam-row corner with its mirrored column. Splice those twins in so that the + # index arithmetic below -- which assumes consecutive corners are index-adjacent + # -- reads the velocity face the section actually crossed. See `_insert_seam_twins`. + i_c, j_c, src = _insert_seam_twins(i_c, j_c, glon, glat, grid) nsec = i_c.size uvindices = { @@ -201,7 +347,11 @@ def uvindices_from_qindices(grid, i_c, j_c, f_c=None): "i":np.zeros(nsec-1, dtype=np.int64), "j":np.zeros(nsec-1, dtype=np.int64), "Yinc":np.zeros(nsec-1, dtype=bool), - "Xinc":np.zeros(nsec-1, dtype=bool) + "Xinc":np.zeros(nsec-1, dtype=bool), + # Which step of the *caller's* corner array each velocity face came from. Faces are + # dropped below, so this is what ties the returned faces back to the corners that + # produced them once the two are no longer one-to-one. + "q":np.asarray(src[:-1], dtype=np.int64), } if f_c is not None: @@ -247,8 +397,8 @@ def uvindices_from_qindices(grid, i_c, j_c, f_c=None): } uvindex["i"] += (offset if zonal else 0) uvindex["j"] += (offset if not(zonal) else 0) - for (key, v) in uvindices.items(): - v[k] = uvindex[key] + for (key, value) in uvindex.items(): + uvindices[key][k] = value # Drop zero-length faces. A single physical corner can carry two indices -- at a # periodic seam, a shared multi-tile boundary corner, or the bipolar fold seam -- so @@ -641,6 +791,14 @@ def convergent_transport( dsout = dsout.assign_coords({ "sign": orient_fact*sect["Lsign"], + # Which step of the caller's corner arrays each face came from. Faces and corner + # steps are not one-to-one (a step between two indices of the same physical corner + # carries no face), so this is what maps a result back onto the section's corners. + "q": xr.DataArray( + uvindices["q"], + coords=(dsout[sect_coord],), + dims=(sect_coord,) + ), "dir": xr.DataArray( np.array(["U" if u else "V" for u in sect["Umask"]]), coords=(dsout[sect_coord],),