FEAT: Add mesh properties, fail on mismatched meshes

docs/conf.py

@@ -522,6 +522,9 @@ def linkcode_resolve(domain, info):
 nb_execution_excludepatterns = ["jupyter_execute", ".jupyter_cache"]
 nb_execution_raise_on_error = True
 nb_execution_timeout = 75
+nb_execution_excludepatterns = [
+    "user_guide/examples/tutorial_nestedgrids.ipynb"
+]  # TODO: Remove once https://github.com/Parcels-code/Parcels/issues/2878 is fixed
 suppress_warnings = ["mystnb.unknown_mime_type"]
 nitpicky = True
 nitpick_ignore_regex = [

docs/user_guide/examples/tutorial_interaction.ipynb

@@ -46,6 +46,7 @@
     "from matplotlib.animation import FuncAnimation\n",
     "\n",
     "import parcels\n",
+    "from parcels._datasets.structured.generated import simple_UV_dataset\n",
     "\n",
     "# for interactive display of animations\n",
     "plt.rcParams[\"animation.html\"] = \"jshtml\""
@@ -98,7 +99,8 @@
    "source": [
     "def DiffusionFieldSet():\n",
     "    \"\"\"Define a fieldset with only diffusion\"\"\"\n",
-    "    fieldset = parcels.FieldSet([])\n",
+    "    ds = simple_UV_dataset(dims=(1, 1, 1, 1), mesh=\"flat\")\n",
+    "    fieldset = parcels.FieldSet.from_sgrid_conventions(ds, mesh=\"flat\")\n",
     "    fieldset.add_constant_field(\"Kh_zonal\", 0.0005, mesh=\"flat\")\n",
     "    fieldset.add_constant_field(\"Kh_meridional\", 0.0005, mesh=\"flat\")\n",
     "    return fieldset"

src/parcels/_core/fieldset.py

@@ -1,7 +1,6 @@
 from __future__ import annotations
 
 import sys
-import warnings
 from collections.abc import Iterable
 from typing import IO, TYPE_CHECKING
 
@@ -11,6 +10,7 @@
 
 import parcels._typing as ptyping
 from parcels._core.field import Field, VectorField
+from parcels._core.mesh import FlatMesh, SphericalMesh
 from parcels._core.model import (
     ModelData,
     StructuredModelData,
@@ -20,7 +20,6 @@
 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._repr_utils import fieldset_describe
 from parcels.interpolators import (
@@ -66,17 +65,23 @@ class FieldSet:
     """
 
     def __init__(self, models: list[ModelData]):
+        if models == []:
+            raise ValueError("List of models can't be empty.")
         for model in models:
             if not isinstance(model, ModelData):
                 raise ValueError(f"Expected `model` to be a ModelData object. Got {model}")
         # assert_compatible_calendars(fields)
 
-        self.models = list(models)
+        self.models = models
         self.constant_model: StructuredModelData | None = None
         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())
+        assert_models_have_same_mesh(self.models)
+
+    @property
+    def mesh(self) -> FlatMesh | SphericalMesh:
+        return self.models[0].mesh
 
     def __setattr__(self, name, value):
         """Set field attribute by name. If context exists and name in context, raise error to prevent overwriting context variable."""
@@ -214,7 +219,7 @@ def add_constant_field(self, name: str, value, mesh: ptyping.TMesh = "spherical"
         self.reconstruct_fields()
         field = getattr(self, name)
         field.interp_method = XConstantField()
-        _warn_if_fields_use_different_meshes(self.fields.values())
+        assert_models_have_same_mesh(self.models)
 
     def add_context(self, name, value):
         """Add context variable to the FieldSet.
@@ -367,26 +372,23 @@ 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.
+class IncompatibleMeshesException(Exception): ...
 
-    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,
-        )
+def assert_models_have_same_mesh(models: list[ModelData]):
+    if models == []:
+        return
+
+    first_mesh = None
+    for i, model in enumerate(models):
+        if first_mesh is None:
+            first_mesh = model.mesh
+            continue
+
+        if model.mesh != first_mesh:
+            raise IncompatibleMeshesException(
+                f"All ModelData objects must have the same meshes. ModelData at index 0 has a mesh of {first_mesh!r} while ModelData at index {i} has mesh {model.mesh!r} "
+            )
 
 
 class CalendarError(Exception):  # TODO: Move to a Parcels errors module

src/parcels/_core/mesh.py

@@ -3,7 +3,7 @@
 
 import numpy as np
 
-EARTH_RADIUS = 6366707.019493707
+EARTH_RADIUS = 6_366_707.019493707  # m
 
 
 class BaseMesh(ABC):

src/parcels/_core/model.py

@@ -15,6 +15,7 @@
 from parcels._core._windowed_array import maybe_windowed
 from parcels._core.basegrid import BaseGrid
 from parcels._core.field import Field, VectorField
+from parcels._core.mesh import FlatMesh, SphericalMesh
 from parcels._core.utils.time import TimeInterval
 from parcels._core.uxgrid import UxGrid
 from parcels._core.xgrid import (
@@ -46,6 +47,10 @@ class ModelData(ABC):
     field_to_interpolator: dict[str, ScalarInterpolator | VectorInterpolator]
     vector_field_components: ptyping.VectorFields
 
+    @property
+    def mesh(self) -> FlatMesh | SphericalMesh:
+        return self.grid._mesh
+
     @abstractmethod
     def construct_fields(self) -> list[Field | VectorField]: ...
 

tests/test_fieldset.py

@@ -11,7 +11,8 @@
 import parcels.tutorial
 import tests
 from parcels import ParticleFile, ParticleSet, convert, open_raw_zarr
-from parcels._core.fieldset import FieldSet, _datetime_to_msg
+from parcels._core.fieldset import FieldSet, IncompatibleMeshesException, _datetime_to_msg
+from parcels._core.mesh import SphericalMesh
 from parcels._core.model import _default_vector_field_components
 from parcels._datasets.structured.generic import datasets as datasets_structured
 from parcels._datasets.structured.generic import datasets_sgrid
@@ -259,14 +260,6 @@ def test_multi_model_nonoverlapping_time_interval():
     assert fieldset.time_interval is None
 
 
-def test_fieldset_time_interval_constant_fields():
-    fieldset = FieldSet([])
-    fieldset.add_constant_field("constant_field", 1.0)
-    fieldset.add_constant_field("constant_field2", 2.0)
-
-    assert fieldset.time_interval is None
-
-
 def test_fieldset_add_incompatible_calendars():
     # tests the adding of fieldsets that have incompatible calendars
     ...
@@ -386,6 +379,26 @@ def test_fieldset_add():
     assert set(fields_before) == set(fset.fields.keys())
 
 
+def test_fieldset_add_different_meshes():
+    ds1 = datasets_structured["ds_2d_left"][["U_A_grid", "V_A_grid", "grid"]].rename({"U_A_grid": "U", "V_A_grid": "V"})
+    ds2 = datasets_structured["ds_2d_left"][["U_A_grid", "V_A_grid", "grid"]].rename(
+        {"U_A_grid": "U_wind", "V_A_grid": "V_wind"}
+    )
+
+    fset1 = FieldSet.from_sgrid_conventions(
+        ds1,
+        mesh=SphericalMesh(71_492_000),  # Jupiter
+    )
+    fset2 = FieldSet.from_sgrid_conventions(
+        ds2,
+        mesh="spherical",  # earth
+        vector_fields={"UV_wind": ("U_wind", "V_wind")},
+    )
+
+    with pytest.raises(IncompatibleMeshesException, match="All ModelData objects must have the same meshes."):
+        _ = fset1 + fset2
+
+
 def test_vectorfields_without_time():
     """Test that vector fields without a time dimension can be evaluated."""
     ds1 = datasets_structured["ds_2d_left"][["U_A_grid", "V_A_grid", "grid"]].rename({"U_A_grid": "U", "V_A_grid": "V"})

tests/test_mesh.py

@@ -1,8 +1,8 @@
 import numpy as np
 import pytest
 
-from parcels import FieldSet, ParticleSet, SphericalMesh
-from parcels._core.mesh import EARTH_RADIUS
+from parcels import FieldSet, ParticleSet
+from parcels._core.mesh import EARTH_RADIUS, FlatMesh, SphericalMesh
 from parcels._datasets.structured.generated import simple_UV_dataset
 from parcels.kernels import AdvectionRK4
 
@@ -82,3 +82,22 @@ def test_spherical_mesh_rejects_non_numeric_radius(bad_radius):
 def test_spherical_mesh_rejects_nonpos_radius(bad_radius):
     with pytest.raises(ValueError):
         SphericalMesh(radius=bad_radius)
+
+
+@pytest.mark.parametrize(
+    "lhs, op, rhs",
+    [
+        (FlatMesh(), "==", FlatMesh()),
+        (FlatMesh(), "!=", SphericalMesh()),
+        (SphericalMesh(), "==", SphericalMesh()),
+        (SphericalMesh(20_000), "!=", SphericalMesh(30_000)),
+    ],
+)
+def test_mesh_equality(lhs, op, rhs):
+    match op:
+        case "==":
+            assert lhs == rhs
+        case "!=":
+            assert lhs != rhs
+        case _:
+            raise NotImplementedError