Merge branch 'main' into updating_v3_tutorials

Author: erikvansebille

docs/development/unstructured_grid_search.md

@@ -42,12 +42,12 @@ The algorithm has two phases: **initialisation** (once, when the grid is first u
 
 The hash grid is three-dimensional regardless of the source grid:
 
-| Source grid | Mesh type | Hash-grid space                       |
-| ----------- | --------- | ------------------------------------- |
-| `XGrid`     | spherical | Cartesian unit cube (lon/lat → x,y,z) |
-| `XGrid`     | flat      | 2-D lon/lat bounding box (z set to 0) |
-| `UxGrid`    | spherical | Cartesian unit cube                   |
-| `UxGrid`    | flat      | 2-D lon/lat bounding box              |
+| Source grid | Mesh type | Hash-grid space                          |
+| ----------- | --------- | ---------------------------------------- |
+| `XGrid`     | spherical | Cartesian bounding box (lon/lat → x,y,z) |
+| `XGrid`     | flat      | 2-D lon/lat bounding box (z set to 0)    |
+| `UxGrid`    | spherical | Cartesian bounding box (lon/lat → x,y,z) |
+| `UxGrid`    | flat      | 2-D lon/lat bounding box                 |
 
 For spherical grids, node coordinates are converted from degrees to radians and then to Cartesian (x, y, z) on the unit sphere using `parcels._core.index_search._latlon_rad_to_xyz`:
 
@@ -57,7 +57,7 @@ y = sin(lon) * cos(lat)
 z = sin(lat)
 ```
 
-The hash grid then spans the unit cube `[-1, 1]³`. Working in Cartesian space avoids the longitude wrap-around discontinuity that would otherwise cause the bounding boxes of cells crossing the antimeridian to erroneously span the entire domain.
+The hash grid then spans the axis-aligned Cartesian bounding box of the transformed grid nodes. For a global grid this is (nearly) the unit cube `[-1, 1]³`; for a regional grid it is much tighter, so the full quantization resolution (1024 bins per axis) is spent on the region actually covered by the grid rather than the whole sphere. Working in Cartesian space avoids the longitude wrap-around discontinuity that would otherwise cause the bounding boxes of cells crossing the antimeridian to erroneously span the entire domain.
 
 For flat meshes the hash grid simply spans `[lon_min, lon_max] × [lat_min, lat_max]`, with z fixed at 0.
 

docs/user_guide/examples/explanation_kernelloop.md

@@ -59,14 +59,14 @@ import parcels.tutorial
 ds_fields = parcels.tutorial.open_dataset("CopernicusMarine_data_for_Argo_tutorial/data")
 
 # Create an idealised wind field and add it to the dataset
-tdim, ydim, xdim = (len(ds_fields.time),len(ds_fields.latitude), len(ds_fields.longitude))
+ydim, xdim = len(ds_fields.latitude), len(ds_fields.longitude)
 ds_fields["UWind"] = xr.DataArray(
-    data=0.5 * np.ones((tdim, ydim, xdim)) * np.sin(ds_fields.latitude.values - ds_fields.latitude.values.mean())[None, :, None],
-    coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude])
+    data=0.5 * np.ones((ydim, xdim)) * np.sin(ds_fields.latitude.values - ds_fields.latitude.values.mean())[:, None],
+    coords=[ds_fields.latitude, ds_fields.longitude])
 
 ds_fields["VWind"] = xr.DataArray(
-    data=np.zeros((tdim, ydim, xdim)),
-    coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude])
+    data=np.zeros((ydim, xdim)),
+    coords=[ds_fields.latitude, ds_fields.longitude])
 
 fields = {
     "U": ds_fields["uo"],

src/parcels/_core/field.py

@@ -109,6 +109,9 @@ def grid(self):  # TODO PR: Remove in favour of referencing model grid directly
 
     @property
     def time_interval(self):  # TODO PR: Remove in favour of referencing model time_interval directly
+        if "time" not in self.model.data[self.name].dims:
+            # This field does not have a time-dimension
+            return None
         return self.model.time_interval
 
     def __repr__(self):

src/parcels/_core/fieldset.py

@@ -2,6 +2,7 @@
 
 import functools
 import sys
+import warnings
 from collections.abc import Iterable
 from typing import IO, TYPE_CHECKING
 
@@ -21,6 +22,7 @@
 from parcels._core.utils.string import _assert_str_and_python_varname
 from parcels._core.utils.time import get_datetime_type_calendar
 from parcels._core.utils.time import is_compatible as datetime_is_compatible
+from parcels._core.warnings import FieldSetWarning
 from parcels._python import NOTSET, NotSetType
 from parcels._reprs import fieldset_describe
 from parcels.interpolators import (
@@ -73,6 +75,7 @@ def __init__(self, models: list[ModelData]):
         self._fields: dict[str, Field | VectorField] | None = None
         self.reconstruct_fields()
         self.context: dict[str, float] = {}
+        _warn_if_fields_use_different_meshes(self.fields.values())
 
     def __setattr__(self, name, value):
         """Set field attribute by name. If context exists and name in context, raise error to prevent overwriting context variable."""
@@ -151,6 +154,7 @@ def add_field(self, field: Field, name: str | None = None):
             raise ValueError(f"FieldSet already has a Field with name '{name}'")
 
         self.fields[name] = field
+        _warn_if_fields_use_different_meshes(self.fields.values())
 
     def to_windowed_arrays(self, *, max_levels: int | None = None):
         """Wrap dask-backed field data in rolling time-window caches.
@@ -215,6 +219,7 @@ def add_constant_field(self, name: str, value, mesh: ptyping.Mesh = "spherical")
         self.reconstruct_fields()
         field = getattr(self, name)
         field.interp_method = XConstantField()
+        _warn_if_fields_use_different_meshes(self.fields.values())
 
     def add_context(self, name, value):
         """Add context variable to the FieldSet.
@@ -362,6 +367,28 @@ def assert_compatible_fieldsets(left: FieldSet, right: FieldSet) -> None:
         )
 
 
+def _warn_if_fields_use_different_meshes(fields: Iterable[Field | VectorField]):
+    """Warn if multiple fields use different meshes on the underlying grids.
+
+    Parameters
+    ----------
+    fields : Iterable[Field | VectorField]
+        The fields to check for conflicting meshes.
+
+    Warns
+    -----
+    FieldSetWarning
+        If the fields have different meshes on the underlying grids.
+    """
+    meshes = {field.grid._mesh for field in fields}
+    if len(meshes) > 1:
+        warnings.warn(
+            f"FieldSet has multiple different meshes: {meshes}. This may lead to unexpected behavior during execution.",
+            category=FieldSetWarning,
+            stacklevel=3,
+        )
+
+
 class CalendarError(Exception):  # TODO: Move to a parcels errors module
     """Exception raised when the calendar of a field is not compatible with the rest of the Fields. The user should ensure that they only add fields to a FieldSet that have compatible CFtime calendars."""
 

src/parcels/_core/index_search.py

@@ -77,8 +77,8 @@ def _search_time_index(field: Field, time: np.ndarray):
     if field.time_interval is None:
         return {
             "T": {
-                "index": np.zeros(shape=time.shape, dtype=np.int32),
-                "bcoord": np.zeros(shape=time.shape, dtype=np.float32),
+                "index": np.zeros(shape=np.atleast_1d(time).shape, dtype=np.int32),
+                "bcoord": np.zeros(shape=np.atleast_1d(time).shape, dtype=np.float32),
             }
         }