From fa34f51d10b9ac3066b56461d04b688ca5b3e5ad Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Tue, 28 Jul 2026 13:57:22 -0400 Subject: [PATCH 1/2] Update unstructured grid search docs to reflect variable bitwidth in construction --- docs/development/unstructured_grid_search.md | 51 +++++++++++++++----- 1 file changed, 39 insertions(+), 12 deletions(-) diff --git a/docs/development/unstructured_grid_search.md b/docs/development/unstructured_grid_search.md index bedc1aa5e..95d36e78e 100644 --- a/docs/development/unstructured_grid_search.md +++ b/docs/development/unstructured_grid_search.md @@ -9,11 +9,11 @@ This page documents the algorithm used in Parcels to locate which grid cell a pa ## Motivation -On a rectilinear grid, finding which cell contains a particle is a trivial O(1) operation — compute an index from the coordinate directly. On curvilinear and unstructured grids, no such shortcut exists; each cell has an arbitrary shape and position. +On a rectilinear grid, finding which cell contains a particle is a trivial O(1) operation. We can compute an index from the coordinate directly. On curvilinear and unstructured grids, no such shortcut exists; each cell has an arbitrary shape and position. The naive approach is an O(N) linear scan: compute the distance from the particle to each face centroid and find the minimum. For large meshes this becomes a bottleneck, and distance minimisation does not guarantee the nearest centroid is actually the containing cell. -Parcels uses **spatial hashing** to reduce this to an O(1) average-case lookup. The key idea is to impose a cheap, regular "hash grid" on top of the irregular source grid, assign each source-grid face to the hash cells its bounding box overlaps, and then use the hash grid to narrow candidate faces for any query point down to a small set before running an exact point-in-cell test. The key data structure that relates the cells in the hash grid to the source grid is the **hash table**. In Parcels, the hash table is a lookup table that takes as input a position and returns a short list of faces in a grid that likely contain the query point. Once a short list of candidate faces is obtained, Parcels then uses point-in-cell checks to locate which face a particle resides in. +Parcels uses **spatial hashing** to reduce this to an O(1) average-case lookup. The key idea is to impose a cheap, regular "hash grid" on top of the irregular source grid, assign each source-grid face to the hash cells its bounding box overlaps, and then use the hash grid to narrow candidate faces for any query point down to a small set before running an exact point-in-cell test. The key data structure that relates the cells in the hash grid to the source grid is the **hash table**. In Parcels, the hash table is a compressed sparse row matrix whose row indices are morton encodings of the positions in the grid and the non-zero values in each row the of matrix are the face id's in a grid that likely contain the query point. Once a short list of candidate faces is obtained, Parcels then uses point-in-cell checks to locate which face a particle resides in. The search methods that exist in the literature — KD-trees, BVH trees, quad-trees — give O(log N) query time. Spatial hashing achieves O(1) on average, which is particularly important when the same grid is queried millions of times per time step (one query per particle, per step). @@ -26,8 +26,9 @@ The algorithm has two phases: **initialisation** (once, when the grid is first u ### Initialisation 1. Determine the bounding box of every face in the source grid. -2. Map each face's bounding box to a set of integer hash cell coordinates using **Morton encoding** (see below). -3. Sort all (face, morton-code) pairs by morton code and store them in a compact CSR-like structure: arrays of unique codes `keys`, their starting positions `starts`, hit counts `counts`, and the corresponding face indices `i`, `j`. +2. Choose the hash-grid resolution (the quantization `bitwidth`) so that the total table size stays within a fixed memory budget (see **Hash-grid Resolution and the Entry Budget** below). +3. Map each face's bounding box to a set of integer hash cell coordinates using **Morton encoding** (see below). +4. Sort all (face, morton-code) pairs by morton code and store them in a compact CSR-like structure: arrays of unique codes `keys`, their starting positions `starts`, hit counts `counts`, and a single flat array `faces` of the (raveled) face indices. ### Query @@ -79,9 +80,9 @@ Each coordinate axis is mapped from its floating-point range `[min, max]` to an xq = clip( floor( (x - xmin) / (xmax - xmin) * bitwidth ), 0, bitwidth ) ``` -With `bitwidth = 1023` (the default), each axis gets 1024 distinct integer levels. This is equivalent to overlaying a regular 1024×1024×1024 grid — the "hash grid" — on the domain. Each integer triple `(xq, yq, zq)` identifies one cell of this hash grid. +With `bitwidth = 1023` (the maximum), each axis gets 1024 distinct integer levels. This is equivalent to overlaying a regular 1024×1024×1024 grid — the "hash grid" — on the domain. Each integer triple `(xq, yq, zq)` identifies one cell of this hash grid. The `bitwidth` actually used may be lower: it is capped at build time to keep the table within a memory budget. -**Choosing bitwidth:** The default of 1023 was chosen because it fits in 10 bits per axis, and 3 × 10 = 30 bits, which fits exactly in a `uint32`. Increasing bitwidth beyond 1023 would require a `uint64` Morton code. In practice, 1023 provides more than enough resolution — for a global ocean model with ~1 million faces, the hash grid cells are already smaller than individual grid cells. +**Choosing bitwidth:** The ceiling of 1023 was chosen because it fits in 10 bits per axis, and 3 × 10 = 30 bits, which fits exactly in a `uint32`. Increasing bitwidth beyond 1023 would require a `uint64` Morton code. In practice, 1023 provides more than enough resolution — for a global ocean model with ~1 million faces, the hash grid cells are already smaller than individual grid cells — so the build starts at 1023 and only coarsens if the entry budget requires it. **Effect of domain extent:** If the quantization bounds `[xmin, xmax]` exactly match the source grid's coordinate range, then at `bitwidth = 1023` most hash cells contain at most one or two source faces. Expanding the bounds beyond the data range is equivalent to coarsening the hash grid, increasing the average number of candidates per query. @@ -124,7 +125,27 @@ The resulting bit layout (relative to the least significant bit) is: x0, y0, z0, x1, y1, z1, ..., x9, y9, z9 ``` -The combined encode function `parcels._core.spatialhash._encode_morton3d` chains all three steps (quantize → dilate → interleave) and is the function used during queries. +The combined encode function `parcels._core.spatialhash._encode_morton3d` chains all three steps (quantize, dilate, interleave) and is the function used during queries. + +--- + +## Hash-grid Resolution and the Entry Budget + +The CSR table stores one entry per (face, overlapping-hash-cell) pair. A finer hash grid (higher `bitwidth`) means smaller hash cells and fewer candidate faces per query — but each face's bounding box then overlaps more hash cells, so the table holds more entries and costs more memory to build. To keep this bounded, the resolution is chosen at build time against a fixed entry budget rather than always using the maximum 10-bit resolution. + +The relevant constants live at the top of `spatialhash.py`: + +- `_HASH_MAX_BITWIDTH = 1023` — the maximum (and initial) quantization level, 10 bits per axis. +- `_HASH_ENTRIES_PER_FACE = 16` — target average number of hash-cell overlaps per face. +- `_HASH_ENTRY_BUDGET_MIN = 2**22` — a floor so that small grids are never forced onto an over-coarse hash grid. + +The budget on the total number of (face, hash-cell) entries is: + +``` +budget = max(_HASH_ENTRIES_PER_FACE * nfaces, _HASH_ENTRY_BUDGET_MIN) +``` + +`SpatialHash._total_hash_entries(bitwidth)` computes what the table size *would* be at a given resolution: it quantizes the low and high corner of every face bounding box and sums `(nx · ny · nz)` over all faces. If the table at `bitwidth = 1023` already fits the budget, that resolution is kept. Otherwise the constructor binary-searches for the largest bitwidth whose table fits the budget and uses that, coarsening the hash grid just enough to stay within memory. --- @@ -139,11 +160,13 @@ The hash table is built by iterating over every face in the source grid: 3. **Generate one Morton code per (face, hash cell) pair.** For each face, every hash cell within its bounding box gets a Morton code. This is done vectorised: for each face `f` spanning a `nx × ny × nz` block of hash cells, `nx*ny*nz` Morton codes are generated by iterating over all `(xi, yi, zi)` offsets within the block. The result is a flat array `morton_codes` (one entry per face-hashcell overlap) and a parallel array `face_ids`. -4. **Sort by Morton code and compress.** `morton_codes` and `face_ids` are sorted together by code. `np.unique` then finds the unique codes and the start/count for each group, producing a CSR (Compressed Sparse Row) structure: +4. **Sort by Morton code and compress.** Each `(morton_code, face_id)` pair is fused into a single `uint64` — the Morton code in the high 32 bits, the face id in the low 32 bits — and the fused array is sorted in place. Unsigned comparison then orders entries by Morton code, with ties broken by ascending face id, without allocating a separate `argsort` permutation or gather. The code and face id are split back out afterwards. Runs of equal codes are located by comparing adjacent entries (rather than `np.unique`), producing a CSR (Compressed Sparse Row) structure: - `keys` — sorted unique Morton codes - - `starts` — index into `face_sorted` where each unique code's entries begin + - `starts` — index into `faces` where each unique code's entries begin - `counts` — number of face entries for each unique code - - `i`, `j` — the face indices (column, row) for each entry, in sorted order + - `faces` — the flat (raveled) face index for each entry, in sorted order + + Only the flat face id is stored (4 bytes per entry). The `(j, i)` pair is *not* precomputed and stored; instead `query` unravels just the gathered candidate faces on demand (see below). This avoids holding two `int64` index arrays (16 bytes per entry) for the lifetime of the grid. This CSR layout enables O(1) hash table lookup via binary search on `keys`. @@ -161,7 +184,7 @@ Given arrays of particle latitudes `y` and longitudes `x`: 3. **Binary search.** Use `np.searchsorted(keys, query_codes)` to find the position of each query code in the sorted unique-codes array. Queries whose code does not appear exactly in `keys` return no candidates (the hash cell the particle falls in has no registered faces). -4. **Gather candidates.** For each query with a valid hit, use `starts[pos]` and `counts[pos]` to gather the candidate face indices `(j_all, i_all)` from the sorted arrays. This is done fully vectorised using a CSR traversal with `np.repeat` and cumulative sums. At the end of this stage, there are potentially multiple faces to check with a point-in-cell test. +4. **Gather candidates.** For each query with a valid hit, use `starts[pos]` and `counts[pos]` to gather the candidate face ids from the flat `faces` array. This is done fully vectorised using a CSR traversal with `np.repeat` and cumulative sums. The gathered flat face ids are then unravelled to `(j_all, i_all)` with `np.unravel_index` — only the gathered candidates are unravelled, not the whole table. At the end of this stage, there are potentially multiple faces to check with a point-in-cell test. 5. **Point-in-cell test.** Call `self._point_in_cell` (either `curvilinear_point_in_cell` or `uxgrid_point_in_cell`) on all candidates simultaneously. The provided point-in-cell tests for both unstructured and structured grids guarantee that at most candidate cell is found containing each particle. @@ -227,7 +250,11 @@ This two-stage approach (cheap guess check, then hash lookup) means that once pa ## Design Notes and Limitations -**Hash grid resolution vs. memory.** The CSR structure stores one entry per (face, overlapping-hash-cell) pair. For a mesh with highly variable face sizes, large faces can overlap many hash cells, inflating the table size. The 10-bit quantisation (1024 levels per axis) bounds this: a face whose bounding box spans the whole domain contributes at most 1024³ ≈ 10⁹ entries, though in practice face bounding boxes are far smaller. +**Hash grid resolution vs. memory.** The CSR structure stores one entry per (face, overlapping-hash-cell) pair. For a mesh with highly variable face sizes, large faces overlap many hash cells and inflate the table. Rather than relying on the 10-bit quantisation alone to bound this, the build caps the total entry count against a fixed budget and lowers the `bitwidth` as needed (see **Hash-grid Resolution and the Entry Budget** above). The trade-off is moved from memory to query time: a coarser hash grid returns more candidate faces per query. + +**Degenerate faces (XGrid).** On curvilinear spherical grids, a face whose corner nodes are undefined (e.g. land points with lon/lat masked to `0.0`) can span a huge portion of the domain and, on its own, dominate the table size. During construction the curvilinear spherical path runs `_find_degenerate_xgrid_faces`, which flags any cell whose longest edge — the maximum great-circle chord across its four edges and two diagonals — exceeds `threshold_factor` (default 10) times the 99th-percentile edge length, and emits a `FieldSetWarning` naming the first few offending `(j, i)` locations. This detection currently runs **only for the curvilinear (`XGrid`) spherical path**; there is no equivalent check on the `UxGrid` path. + +**Periodic boundaries.** `SpatialHash` does not support queries on periodic elements — there is no longitudinal wrapping or remapping of query points. A particle that crosses a periodic seam is not folded back into the domain by the search. For spherical grids the antimeridian is handled only *implicitly*, by working in Cartesian space (`XGrid`) or projecting onto the face plane (`UxGrid`), not by treating the domain as periodic. **Spherical geometry degeneracy.** Near the poles, lon/lat cells become highly elongated in lon/lat space. Working in Cartesian space mitigates this: all cells have similar extents in the unit cube, so the hash grid resolution is more uniform. From db2a1ab9066a8588d7787de8418b40fe5942eaae Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 28 Jul 2026 17:58:43 +0000 Subject: [PATCH 2/2] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- docs/development/unstructured_grid_search.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/development/unstructured_grid_search.md b/docs/development/unstructured_grid_search.md index 95d36e78e..002102ff0 100644 --- a/docs/development/unstructured_grid_search.md +++ b/docs/development/unstructured_grid_search.md @@ -145,7 +145,7 @@ The budget on the total number of (face, hash-cell) entries is: budget = max(_HASH_ENTRIES_PER_FACE * nfaces, _HASH_ENTRY_BUDGET_MIN) ``` -`SpatialHash._total_hash_entries(bitwidth)` computes what the table size *would* be at a given resolution: it quantizes the low and high corner of every face bounding box and sums `(nx · ny · nz)` over all faces. If the table at `bitwidth = 1023` already fits the budget, that resolution is kept. Otherwise the constructor binary-searches for the largest bitwidth whose table fits the budget and uses that, coarsening the hash grid just enough to stay within memory. +`SpatialHash._total_hash_entries(bitwidth)` computes what the table size _would_ be at a given resolution: it quantizes the low and high corner of every face bounding box and sums `(nx · ny · nz)` over all faces. If the table at `bitwidth = 1023` already fits the budget, that resolution is kept. Otherwise the constructor binary-searches for the largest bitwidth whose table fits the budget and uses that, coarsening the hash grid just enough to stay within memory. --- @@ -166,7 +166,7 @@ The hash table is built by iterating over every face in the source grid: - `counts` — number of face entries for each unique code - `faces` — the flat (raveled) face index for each entry, in sorted order - Only the flat face id is stored (4 bytes per entry). The `(j, i)` pair is *not* precomputed and stored; instead `query` unravels just the gathered candidate faces on demand (see below). This avoids holding two `int64` index arrays (16 bytes per entry) for the lifetime of the grid. + Only the flat face id is stored (4 bytes per entry). The `(j, i)` pair is _not_ precomputed and stored; instead `query` unravels just the gathered candidate faces on demand (see below). This avoids holding two `int64` index arrays (16 bytes per entry) for the lifetime of the grid. This CSR layout enables O(1) hash table lookup via binary search on `keys`. @@ -254,7 +254,7 @@ This two-stage approach (cheap guess check, then hash lookup) means that once pa **Degenerate faces (XGrid).** On curvilinear spherical grids, a face whose corner nodes are undefined (e.g. land points with lon/lat masked to `0.0`) can span a huge portion of the domain and, on its own, dominate the table size. During construction the curvilinear spherical path runs `_find_degenerate_xgrid_faces`, which flags any cell whose longest edge — the maximum great-circle chord across its four edges and two diagonals — exceeds `threshold_factor` (default 10) times the 99th-percentile edge length, and emits a `FieldSetWarning` naming the first few offending `(j, i)` locations. This detection currently runs **only for the curvilinear (`XGrid`) spherical path**; there is no equivalent check on the `UxGrid` path. -**Periodic boundaries.** `SpatialHash` does not support queries on periodic elements — there is no longitudinal wrapping or remapping of query points. A particle that crosses a periodic seam is not folded back into the domain by the search. For spherical grids the antimeridian is handled only *implicitly*, by working in Cartesian space (`XGrid`) or projecting onto the face plane (`UxGrid`), not by treating the domain as periodic. +**Periodic boundaries.** `SpatialHash` does not support queries on periodic elements — there is no longitudinal wrapping or remapping of query points. A particle that crosses a periodic seam is not folded back into the domain by the search. For spherical grids the antimeridian is handled only _implicitly_, by working in Cartesian space (`XGrid`) or projecting onto the face plane (`UxGrid`), not by treating the domain as periodic. **Spherical geometry degeneracy.** Near the poles, lon/lat cells become highly elongated in lon/lat space. Working in Cartesian space mitigates this: all cells have similar extents in the unit cube, so the hash grid resolution is more uniform.