diff --git a/src/parcels/__init__.py b/src/parcels/__init__.py index 7ae1f6928..1378340fd 100644 --- a/src/parcels/__init__.py +++ b/src/parcels/__init__.py @@ -61,6 +61,7 @@ "OutsideTimeInterval", "StatusCode", # Warnings + "FieldEvalWarning", "FieldSetWarning", "FileWarning", "KernelWarning", diff --git a/src/parcels/_core/field.py b/src/parcels/_core/field.py index 5a1b6ed87..057256594 100644 --- a/src/parcels/_core/field.py +++ b/src/parcels/_core/field.py @@ -17,6 +17,7 @@ from parcels._core.utils.string import _assert_str_and_python_varname from parcels._core.utils.time import TimeInterval from parcels._core.uxgrid import UxGrid +from parcels._core.warnings import FieldEvalWarning from parcels._core.xgrid import XGrid, _transpose_xfield_data_to_tzyx, assert_all_field_dims_have_axis from parcels._python import assert_same_function_signature from parcels._reprs import field_repr, vectorfield_repr @@ -196,6 +197,7 @@ def eval(self, time: datetime, z, y, x, particles=None): value = self._interp_method(particle_positions, grid_positions, self) _update_particle_states_interp_value(particles, value) + _mask_outofbounds_values(grid_positions, value) return value @@ -299,6 +301,7 @@ def eval(self, time: datetime, z, y, x, particles=None): for vel in (u, v, w): _update_particle_states_interp_value(particles, vel) + _mask_outofbounds_values(grid_positions, vel) if "3D" in self.vector_type: return (u, v, w) @@ -367,6 +370,20 @@ def _update_particle_states_position(particles, grid_positions: dict): ) +def _mask_outofbounds_values(grid_positions: dict, value): + mask = np.zeros(value.shape, dtype=bool) + for dim in ["X", "Y", "Z", "FACE"]: + if dim in grid_positions: + mask[grid_positions[dim]["index"] < 0] = True + if np.any(mask): + warnings.warn( + "Some interpolated values are out-of-bounds. These values are set to 0. Treat carefully.", + FieldEvalWarning, + stacklevel=2, + ) + value[mask] = 0.0 + + def _update_particle_states_interp_value(particles, value): """Update the particle states based on the interpolated value, but only if state is not an Error already.""" if particles: diff --git a/src/parcels/_core/kernel.py b/src/parcels/_core/kernel.py index 5d0fbe35a..ab1986353 100644 --- a/src/parcels/_core/kernel.py +++ b/src/parcels/_core/kernel.py @@ -16,7 +16,7 @@ _raise_grid_searching_error, _raise_outside_time_interval_error, ) -from parcels._core.warnings import KernelWarning +from parcels._core.warnings import FieldEvalWarning, KernelWarning from parcels._python import assert_same_function_signature from parcels.kernels import ( AdvectionAnalytical, @@ -206,13 +206,16 @@ def execute(self, pset, endtime, dt): # run kernels for all particles that need to be evaluated for f in self._kernels: - f(pset[evaluate_particles], self._fieldset) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", FieldEvalWarning) - # check for particles that have to be repeated - repeat_particles = pset.state == StatusCode.Repeat - while np.any(repeat_particles): - f(pset[repeat_particles], self._fieldset) + f(pset[evaluate_particles], self._fieldset) + + # check for particles that have to be repeated repeat_particles = pset.state == StatusCode.Repeat + while np.any(repeat_particles): + f(pset[repeat_particles], self._fieldset) + repeat_particles = pset.state == StatusCode.Repeat # apply position/time update only to particles still in a normal state # (particles that signalled Stop*/Delete/errors should not have time/position advanced) diff --git a/src/parcels/_core/warnings.py b/src/parcels/_core/warnings.py index 40d953608..d6f34fb93 100644 --- a/src/parcels/_core/warnings.py +++ b/src/parcels/_core/warnings.py @@ -28,6 +28,16 @@ class FileWarning(UserWarning): pass +class FieldEvalWarning(UserWarning): + """Warning that is raised when there are issues during the evaluation of a Field. + + These warnings can be related to out-of-bounds indices during interpolation, + or other issues that arise during the evaluation of a Field at particle positions. + """ + + pass + + class KernelWarning(RuntimeWarning): """Warning that is raised when there are issues with the Kernel. diff --git a/src/parcels/interpolators/_uxinterpolators.py b/src/parcels/interpolators/_uxinterpolators.py index f434ac96e..3fb4e7b6c 100644 --- a/src/parcels/interpolators/_uxinterpolators.py +++ b/src/parcels/interpolators/_uxinterpolators.py @@ -158,5 +158,5 @@ def Ux_Velocity( if "3D" in vectorfield.vector_type: w = vectorfield.W._interp_method(particle_positions, grid_positions, vectorfield.W) else: - w = 0.0 + w = np.zeros_like(u) return u, v, w diff --git a/src/parcels/interpolators/_xinterpolators.py b/src/parcels/interpolators/_xinterpolators.py index 725422d02..573cf3427 100644 --- a/src/parcels/interpolators/_xinterpolators.py +++ b/src/parcels/interpolators/_xinterpolators.py @@ -22,7 +22,7 @@ def ZeroInterpolator( field: Field, ) -> np.float32 | np.float64: """Template function used for the signature check of the lateral interpolation methods.""" - return 0.0 + return np.zeros_like(particle_positions["lon"]) def ZeroInterpolator_Vector( @@ -31,7 +31,7 @@ def ZeroInterpolator_Vector( vectorfield: VectorField, ) -> np.float32 | np.float64: """Template function used for the signature check of the interpolation methods for velocity fields.""" - return 0.0 + return np.zeros_like(particle_positions["lon"]) def _get_corner_data_Agrid( @@ -140,8 +140,8 @@ def XConstantField( grid_positions: dict[XgridAxis, dict[str, int | float | np.ndarray]], field: Field, ): - """Returning the single value of a Constant Field (with a size=(1,1,1,1) array)""" - return field.data[0, 0, 0, 0].values + """Returning the value of a Constant Field (with a size=(1,1,1,1) array)""" + return field.data[0, 0, 0, 0].values * np.ones_like(particle_positions["lon"]) def XLinear_Velocity( @@ -159,7 +159,7 @@ def XLinear_Velocity( if vectorfield.W: w = XLinear(particle_positions, grid_positions, vectorfield.W) else: - w = 0.0 + w = np.zeros_like(u) return u, v, w @@ -489,7 +489,7 @@ def is_land(ti: int, zi: int, yi: int, xi: int): w *= f_w else: - w = None + w = np.zeros_like(u) return u, v, w diff --git a/src/parcels/kernels/_sigmagrids.py b/src/parcels/kernels/_sigmagrids.py index 5da1b81b2..bfa408eed 100644 --- a/src/parcels/kernels/_sigmagrids.py +++ b/src/parcels/kernels/_sigmagrids.py @@ -42,44 +42,44 @@ def AdvectionRK4_3D_CROCO(particles, fieldset): # pragma: no cover It also uses linear interpolation of the W field, which gives much better results than the default C-grid interpolation. """ dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.time, particles.dt) - sigma = particles.z / fieldset.h[particles.time, 0, particles.lat, particles.lon] + sigma = particles.z / fieldset.h[particles.time, np.zeros_like(particles.z), particles.lat, particles.lon] sig = convert_z_to_sigma_croco(fieldset, particles.time, particles.z, particles.lat, particles.lon, particles) (u1, v1) = fieldset.UV[particles.time, sig, particles.lat, particles.lon, particles] w1 = fieldset.W[particles.time, sig, particles.lat, particles.lon, particles] - w1 *= sigma / fieldset.h[particles.time, 0, particles.lat, particles.lon] + w1 *= sigma / fieldset.h[particles.time, np.zeros_like(particles.z), particles.lat, particles.lon] lon1 = particles.lon + u1 * 0.5 * dt lat1 = particles.lat + v1 * 0.5 * dt sig_dep1 = sigma + w1 * 0.5 * dt - dep1 = sig_dep1 * fieldset.h[particles.time, 0, lat1, lon1] + dep1 = sig_dep1 * fieldset.h[particles.time, np.zeros_like(particles.z), lat1, lon1] sig1 = convert_z_to_sigma_croco(fieldset, particles.time + 0.5 * dt, dep1, lat1, lon1, particles) (u2, v2) = fieldset.UV[particles.time + 0.5 * dt, sig1, lat1, lon1, particles] w2 = fieldset.W[particles.time + 0.5 * dt, sig1, lat1, lon1, particles] - w2 *= sig_dep1 / fieldset.h[particles.time, 0, lat1, lon1] + w2 *= sig_dep1 / fieldset.h[particles.time, np.zeros_like(particles.z), lat1, lon1] lon2 = particles.lon + u2 * 0.5 * dt lat2 = particles.lat + v2 * 0.5 * dt sig_dep2 = sigma + w2 * 0.5 * dt - dep2 = sig_dep2 * fieldset.h[particles.time, 0, lat2, lon2] + dep2 = sig_dep2 * fieldset.h[particles.time, np.zeros_like(particles.z), lat2, lon2] sig2 = convert_z_to_sigma_croco(fieldset, particles.time + 0.5 * dt, dep2, lat2, lon2, particles) (u3, v3) = fieldset.UV[particles.time + 0.5 * dt, sig2, lat2, lon2, particles] w3 = fieldset.W[particles.time + 0.5 * dt, sig2, lat2, lon2, particles] - w3 *= sig_dep2 / fieldset.h[particles.time, 0, lat2, lon2] + w3 *= sig_dep2 / fieldset.h[particles.time, np.zeros_like(particles.z), lat2, lon2] lon3 = particles.lon + u3 * dt lat3 = particles.lat + v3 * dt sig_dep3 = sigma + w3 * dt - dep3 = sig_dep3 * fieldset.h[particles.time, 0, lat3, lon3] + dep3 = sig_dep3 * fieldset.h[particles.time, np.zeros_like(particles.z), lat3, lon3] sig3 = convert_z_to_sigma_croco(fieldset, particles.time + dt, dep3, lat3, lon3, particles) (u4, v4) = fieldset.UV[particles.time + dt, sig3, lat3, lon3, particles] w4 = fieldset.W[particles.time + dt, sig3, lat3, lon3, particles] - w4 *= sig_dep3 / fieldset.h[particles.time, 0, lat3, lon3] + w4 *= sig_dep3 / fieldset.h[particles.time, np.zeros_like(particles.z), lat3, lon3] lon4 = particles.lon + u4 * dt lat4 = particles.lat + v4 * dt sig_dep4 = sigma + w4 * dt - dep4 = sig_dep4 * fieldset.h[particles.time, 0, lat4, lon4] + dep4 = sig_dep4 * fieldset.h[particles.time, np.zeros_like(particles.z), lat4, lon4] particles.dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * dt particles.dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * dt particles.dz += ( diff --git a/tests/test_advection.py b/tests/test_advection.py index 64bb429ba..3141c62ac 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np import pandas as pd import pytest @@ -15,6 +17,7 @@ convert, ) from parcels._core.utils.time import timedelta_to_float +from parcels._core.warnings import FieldEvalWarning from parcels._datasets.structured.generated import ( decaying_moving_eddy_dataset, moving_eddy_dataset, @@ -176,7 +179,9 @@ def SubmergeParticle(particles, fieldset): # pragma: no cover kernels.append(DeleteParticle) pset = ParticleSet(fieldset=fieldset, lon=0.5, lat=0.5, z=0.9) - pset.execute(kernels, runtime=np.timedelta64(10, "s"), dt=np.timedelta64(1, "s")) + with warnings.catch_warnings(): + warnings.simplefilter("error", FieldEvalWarning) + pset.execute(kernels, runtime=np.timedelta64(10, "s"), dt=np.timedelta64(1, "s")) if direction == "up" and resubmerge_particle: np.testing.assert_allclose(pset.lon[0], 0.6, atol=1e-5) diff --git a/tests/test_field.py b/tests/test_field.py index 5e57c43cb..da4e0e1e0 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -5,8 +5,11 @@ from parcels import Field, UxGrid, VectorField, XGrid from parcels._core.fieldset import FieldSet +from parcels._core.warnings import FieldEvalWarning +from parcels._datasets.structured.generated import simple_UV_dataset from parcels._datasets.structured.generic import T as T_structured from parcels._datasets.structured.generic import datasets as datasets_structured +from parcels._datasets.unstructured.generic import _ux_constant_flow_face_centered_2D from parcels._datasets.unstructured.generic import datasets as datasets_unstructured from parcels.interpolators import ( UxConstantFaceConstantZC, @@ -270,6 +273,57 @@ def test_field_constant_in_time(): assert np.isclose(P1, P2) +@pytest.mark.parametrize( + "field_name, location, expected", + [ + ("U", (0.0, 0.0, 0.0, 5e6), 0.0), + ("UV", (0.0, 0.0, 0.0, 5e6), [[0.0], [0.0]]), + ("U", (0.0, 0.0, 5e6, 0.0), 0.0), + ("UV", (0.0, 0.0, 5e6, 0.0), [[0.0], [0.0]]), + ("U", (0.0, 5e6, 0.0, 0.0), 0.0), + ("UV", (0.0, 5e6, 0.0, 0.0), [[0.0], [0.0]]), + ], +) +def test_field_eval_out_of_bounds_structured(field_name, location, expected): + """Test that Field.eval returns IndexError when queried outside the grid boundaries.""" + # eval outside of bounds should return 0.0 + ds = simple_UV_dataset(mesh="flat") + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") + fieldset.U.data[:] = 1.0 + fieldset.V.data[:] = 2.0 + field = getattr(fieldset, field_name) + with pytest.warns( + FieldEvalWarning, + match="Some interpolated values are out-of-bounds. These values are set to 0. Treat carefully.", + ): + np.testing.assert_allclose(field.eval(*location), expected) + + +@pytest.mark.parametrize( + "field_name, location, expected", + [ + ("U", (0.0, 0.0, 0.0, 5e6), 0.0), + ("UV", (0.0, 0.0, 0.0, 5e6), [[0.0], [0.0]]), + ("U", (0.0, 0.0, 5e6, 0.0), 0.0), + ("UV", (0.0, 0.0, 5e6, 0.0), [[0.0], [0.0]]), + ("U", (0.0, 5e6, 0.0, 0.0), 0.0), + ("UV", (0.0, 5e6, 0.0, 0.0), [[0.0], [0.0]]), + ], +) +def test_field_eval_out_of_bounds_unstructured(field_name, location, expected): + """Test that Field.eval returns IndexError when queried outside the grid boundaries.""" + ds = _ux_constant_flow_face_centered_2D() + fieldset = FieldSet.from_ugrid_conventions(ds, mesh="flat") + fieldset.U.data[:] = 1.0 + fieldset.V.data[:] = 2.0 + field = getattr(fieldset, field_name) + with pytest.warns( + FieldEvalWarning, + match="Some interpolated values are out-of-bounds. These values are set to 0. Treat carefully.", + ): + np.testing.assert_allclose(field.eval(*location), expected) + + def test_field_unstructured_grid_creation(): ...