diff --git a/docs/conf.py b/docs/conf.py index feb8b5f4c..a6c0339f8 100755 --- a/docs/conf.py +++ b/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 = [ diff --git a/docs/user_guide/examples/tutorial_interaction.ipynb b/docs/user_guide/examples/tutorial_interaction.ipynb index 0888e3c51..7d42fa9d9 100644 --- a/docs/user_guide/examples/tutorial_interaction.ipynb +++ b/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" diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 696aa6544..5276b7b46 100644 --- a/src/parcels/_core/fieldset.py +++ b/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 diff --git a/src/parcels/_core/mesh.py b/src/parcels/_core/mesh.py index 7ed5e4998..d9cd468c5 100644 --- a/src/parcels/_core/mesh.py +++ b/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): diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index cdb2c0800..293425734 100644 --- a/src/parcels/_core/model.py +++ b/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]: ... diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index b1b6b16f6..1090d2e44 100644 --- a/tests/test_fieldset.py +++ b/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"}) diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 9a0a4b709..e43bf3057 100644 --- a/tests/test_mesh.py +++ b/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