diff --git a/docs/development/unstructured_grid_search.md b/docs/development/unstructured_grid_search.md index bedc1aa5e..c4f519b54 100644 --- a/docs/development/unstructured_grid_search.md +++ b/docs/development/unstructured_grid_search.md @@ -5,6 +5,8 @@ This page documents the algorithm used in Parcels to locate which grid cell a pa - `src/parcels/_core/spatialhash.py` — `SpatialHash` class and Morton encoding utilities - `src/parcels/_core/index_search.py` — point-in-cell tests and the high-level search dispatch +For debugging information regarding the spatial hash grid that underlies a curvilinear or unstructured grid, call `SpatialHash.describe()`, which prints a summary of the hash table's statistics. This table includes the total number of occupied hash cells, the percentage of hash cells that are occupied, the minimum, mean, and maximum number of mesh faces in a hash cell, and more. + --- ## Motivation diff --git a/src/parcels/_core/spatialhash.py b/src/parcels/_core/spatialhash.py index 208cbfe95..b3e6aa677 100644 --- a/src/parcels/_core/spatialhash.py +++ b/src/parcels/_core/spatialhash.py @@ -1,4 +1,6 @@ +import sys import warnings +from typing import IO import numpy as np @@ -10,6 +12,7 @@ ) from parcels._core.warnings import FieldSetWarning from parcels._python import isinstance_noimport +from parcels._reprs import spatialhash_describe # Budget on the total number of (face, hash cell) pairs in the hash table: # max(_HASH_ENTRIES_PER_FACE * nfaces, _HASH_ENTRY_BUDGET_MIN). @@ -524,6 +527,22 @@ def query(self, y, x): coords_best.reshape((num_queries, coordinates.shape[1])), ) + def describe(self, buf: IO | None = None) -> None: + """ + Summary of the SpatialHash's hash-table statistics (resolution, occupancy, + entry counts). + + Parameters + ---------- + buf : file-like, default: sys.stdout + writable buffer + """ + if buf is None: + buf = sys.stdout + assert buf is not None + + buf.write(spatialhash_describe(self)) + def _dilate_bits(n): """ diff --git a/src/parcels/_reprs.py b/src/parcels/_reprs.py index 5a6ae2944..03323b2a7 100644 --- a/src/parcels/_reprs.py +++ b/src/parcels/_reprs.py @@ -18,6 +18,7 @@ from parcels import Field, FieldSet, ParticleSet from parcels._core.field import VectorField from parcels._core.model import ModelData + from parcels._core.spatialhash import SpatialHash from parcels._core.utils.time import TimeInterval @@ -280,6 +281,36 @@ def fieldset_describe(fieldset: FieldSet) -> str: ) +def spatialhash_describe(spatialhash: SpatialHash) -> str: + grid = spatialhash._source_grid + hash_table = spatialhash._hash_table + counts = hash_table["counts"] + + n_faces = int(np.size(spatialhash._xlow)) + n_entries = int(hash_table["faces"].size) + n_occupied_cells = int(hash_table["keys"].size) + n_total_cells = (spatialhash._bitwidth + 1) ** 3 + + rows = { + "Grid type": type(grid).__name__, + "Mesh": grid._mesh, + "Total mesh faces": f"{n_faces:,d}", + "Bitwidth (current / max)": f"{spatialhash._bitwidth} / 1023 (higher = finer resolution hash grid)", + "Total hash cells": f"{n_total_cells:,d}", + "Occupied hash cells": f"{n_occupied_cells:,d}, {n_occupied_cells / n_total_cells * 100:.4f}%", + "Total (hash cell --> gird face) entries": f"{n_entries:,d}", + "Entries per occupied hash cell (avg)": f"{n_entries / n_occupied_cells:.2f}" if n_occupied_cells else "-", + "Entries per face (avg)": f"{n_entries / n_faces:.2f}" if n_faces else "-", + "Faces per occupied hash cell (min / mean / max)": ( + f"{counts.min():,d} / {counts.mean():.2f} / {counts.max():,d}" if n_occupied_cells else "-" + ), + } + key_width = max(len(k) for k in rows) + table = "\n".join(f"{k.ljust(key_width)} : {v}" for k, v in rows.items()) + + return "Spatial Hash Grid Statistics" + "\n" + table + "\n" + + def _get_parent_model(field: Field | VectorField) -> ModelData: if isinstance_noimport(field, "Field"): return field.model # type:ignore[union-attr] diff --git a/tests/test_spatialhash.py b/tests/test_spatialhash.py index 76cdf88d4..78ae1311b 100644 --- a/tests/test_spatialhash.py +++ b/tests/test_spatialhash.py @@ -1,3 +1,5 @@ +from io import StringIO + import numpy as np from parcels._core.fieldset import FieldSet @@ -20,6 +22,30 @@ def test_spatialhash_init(): assert spatialhash is not None +def test_spatialhash_describe(): + ds = datasets["2d_left_rotated"] + grid = FieldSet.from_sgrid_conventions(ds, mesh="flat").data_g.grid + spatialhash = grid.get_spatial_hash() + + io = StringIO() + expected = """\ +Spatial Hash Grid Statistics +Grid type : XGrid +Mesh : FlatMesh() +Total mesh faces : 1,711 +Bitwidth (current / max) : 1023 / 1023 (higher = finer resolution hash grid) +Total hash cells : 1,073,741,824 +Occupied hash cells : 796,054, 0.0741% +Total (hash cell --> gird face) entries : 1,080,194 +Entries per occupied hash cell (avg) : 1.36 +Entries per face (avg) : 631.32 +Faces per occupied hash cell (min / mean / max) : 1 / 1.36 / 4 +""" + spatialhash.describe(io) + actual = io.getvalue() + assert actual == expected + + def test_invalid_positions(): ds = datasets["2d_left_rotated"] grid = FieldSet.from_sgrid_conventions(ds, mesh="flat").data_g.grid