From 23f96058f1ba25083485b67dca5d2a7f42ffcc28 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Thu, 27 Feb 2025 11:49:59 +0100 Subject: [PATCH 1/3] Removing particle.ti Variable --- parcels/particledata.py | 4 ++-- parcels/particleset.py | 3 --- tests/test_fieldset.py | 2 +- 3 files changed, 3 insertions(+), 6 deletions(-) diff --git a/parcels/particledata.py b/parcels/particledata.py index 6ef36fb59a..39809530d7 100644 --- a/parcels/particledata.py +++ b/parcels/particledata.py @@ -122,8 +122,8 @@ def __init__(self, pclass, lon, lat, depth, time, lonlatdepth_dtype, pid_orig, n self._ncount = len(lon) for v in self.ptype.variables: - if v.name in ["ei", "ti"]: - self._data[v.name] = np.empty((len(lon), ngrid), dtype=v.dtype) + if v.name == "ei": + self._data[v.name] = np.empty((len(lon), ngrid), dtype=v.dtype) # TODO len(lon) can be self._ncount? else: self._data[v.name] = np.empty(self._ncount, dtype=v.dtype) diff --git a/parcels/particleset.py b/parcels/particleset.py index 0fa834b907..e6946e2d8c 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -128,13 +128,11 @@ def ArrayClass_init(self, *args, **kwargs): self.ngrids = type(self).ngrids.initial if self.ngrids >= 0: self.ei = np.zeros(self.ngrids, dtype=np.int32) - self.ti = -1 * np.ones(self.ngrids, dtype=np.int32) super(type(self), self).__init__(*args, **kwargs) array_class_vdict = { "ngrids": Variable("ngrids", dtype=np.int32, to_write=False, initial=-1), "ei": Variable("ei", dtype=np.int32, to_write=False), - "ti": Variable("ti", dtype=np.int32, to_write=False, initial=-1), "__init__": ArrayClass_init, } array_class = type(class_name, (pclass,), array_class_vdict) @@ -719,7 +717,6 @@ def from_particlefile( v.name not in [ "ei", - "ti", "dt", "depth", "id", diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index bab9acaeb3..2e31060a9b 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -571,7 +571,7 @@ def test_fieldset_write(tmp_zarrfile): def UpdateU(particle, fieldset, time): # pragma: no cover tmp1, tmp2 = fieldset.UV[particle] _, yi, xi = fieldset.U.unravel_index(particle.ei) - fieldset.U.data[particle.ti, yi, xi] += 1 + fieldset.U.data[0, yi, xi] += 1 fieldset.U.grid.time[0] = time pset = ParticleSet(fieldset, pclass=Particle, lon=5, lat=5) From 17a8d3c7a16f2b09cec91f29050b47742155005c Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Thu, 27 Feb 2025 13:12:54 +0100 Subject: [PATCH 2/3] Cleanign up of search_indices method Not setting to default ti=-1 anymore (as not necessary, and confusing) Also reordering arguments list to align with (time, depth, lat, lon) order throughout Parcels --- parcels/_index_search.py | 4 ++-- parcels/application_kernels/advection.py | 2 +- parcels/field.py | 22 +++++++++++----------- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/parcels/_index_search.py b/parcels/_index_search.py index 1b7e77884e..a52baf56ff 100644 --- a/parcels/_index_search.py +++ b/parcels/_index_search.py @@ -137,7 +137,7 @@ def search_indices_vertical_s( def _search_indices_rectilinear( - field: Field, time: float, z: float, y: float, x: float, ti=-1, particle=None, search2D=False + field: Field, time: float, z: float, y: float, x: float, ti: int, particle=None, search2D=False ): grid = field.grid @@ -223,7 +223,7 @@ def _search_indices_rectilinear( return (zeta, eta, xsi, zi, yi, xi) -def _search_indices_curvilinear(field: Field, time, z, y, x, ti=-1, particle=None, search2D=False): +def _search_indices_curvilinear(field: Field, time, z, y, x, ti, particle=None, search2D=False): if particle: zi, yi, xi = field.unravel_index(particle.ei) else: diff --git a/parcels/application_kernels/advection.py b/parcels/application_kernels/advection.py index e7893e45fa..86ae5c3ba7 100644 --- a/parcels/application_kernels/advection.py +++ b/parcels/application_kernels/advection.py @@ -185,7 +185,7 @@ def AdvectionAnalytical(particle, fieldset, time): # pragma: no cover ds_t = min(ds_t, time_i[np.where(time - fieldset.U.grid.time[ti] < time_i)[0][0]]) zeta, eta, xsi, zi, yi, xi = fieldset.U._search_indices( - -1, particle.depth, particle.lat, particle.lon, particle=particle + time, particle.depth, particle.lat, particle.lon, ti, particle=particle ) if withW: if abs(xsi - 1) < tol: diff --git a/parcels/field.py b/parcels/field.py index deac1a8181..e3e0a5d032 100644 --- a/parcels/field.py +++ b/parcels/field.py @@ -877,15 +877,15 @@ def cell_areas(self): """ return _calc_cell_areas(self.grid) - def _search_indices(self, time, z, y, x, ti=-1, particle=None, search2D=False): + def _search_indices(self, time, z, y, x, ti, particle=None, search2D=False): if self.grid._gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]: return _search_indices_rectilinear(self, time, z, y, x, ti, particle=particle, search2D=search2D) else: return _search_indices_curvilinear(self, time, z, y, x, ti, particle=particle, search2D=search2D) - def _interpolator2D(self, ti, z, y, x, particle=None): + def _interpolator2D(self, time, z, y, x, ti, particle=None): """Impelement 2D interpolation with coordinate transformations as seen in Delandmeter and Van Sebille (2019), 10.5194/gmd-12-3571-2019..""" - (_, eta, xsi, _, yi, xi) = self._search_indices(-1, z, y, x, particle=particle) + (_, eta, xsi, _, yi, xi) = self._search_indices(time, z, y, x, ti, particle=particle) ctx = InterpolationContext2D(self.data, eta, xsi, ti, yi, xi) try: @@ -899,7 +899,7 @@ def _interpolator2D(self, ti, z, y, x, particle=None): raise RuntimeError(self.interp_method + " is not implemented for 2D grids") return f(ctx) - def _interpolator3D(self, ti, z, y, x, time, particle=None): + def _interpolator3D(self, time, z, y, x, ti, particle=None): """Impelement 3D interpolation with coordinate transformations as seen in Delandmeter and Van Sebille (2019), 10.5194/gmd-12-3571-2019..""" (zeta, eta, xsi, zi, yi, xi) = self._search_indices(time, z, y, x, ti, particle=particle) ctx = InterpolationContext3D(self.data, zeta, eta, xsi, ti, zi, yi, xi, self.gridindexingtype) @@ -931,13 +931,13 @@ def temporal_interpolate_fullfield(self, ti, time): f1 = self.data[ti + 1, :] return f0 + (f1 - f0) * ((time - t0) / (t1 - t0)) - def _spatial_interpolation(self, ti, z, y, x, time, particle=None): - """Interpolate horizontal field values using a SciPy interpolator.""" + def _spatial_interpolation(self, time, z, y, x, ti, particle=None): + """Interpolate spatial field values.""" try: if self.grid.zdim == 1: - val = self._interpolator2D(ti, z, y, x, particle=particle) + val = self._interpolator2D(time, z, y, x, ti, particle=particle) else: - val = self._interpolator3D(ti, z, y, x, time, particle=particle) + val = self._interpolator3D(time, z, y, x, ti, particle=particle) if np.isnan(val): # Detect Out-of-bounds sampling and raise exception @@ -1002,8 +1002,8 @@ def eval(self, time, z, y, x, particle=None, applyConversion=True): if self.gridindexingtype == "croco" and self not in [self.fieldset.H, self.fieldset.Zeta]: z = _croco_from_z_to_sigma_scipy(self.fieldset, time, z, y, x, particle=particle) if ti < self.grid.tdim - 1 and time > self.grid.time[ti]: - f0 = self._spatial_interpolation(ti, z, y, x, time, particle=particle) - f1 = self._spatial_interpolation(ti + 1, z, y, x, time, particle=particle) + f0 = self._spatial_interpolation(time, z, y, x, ti, particle=particle) + f1 = self._spatial_interpolation(time, z, y, x, ti + 1, particle=particle) t0 = self.grid.time[ti] t1 = self.grid.time[ti + 1] value = f0 + (f1 - f0) * ((time - t0) / (t1 - t0)) @@ -1011,7 +1011,7 @@ def eval(self, time, z, y, x, particle=None, applyConversion=True): # Skip temporal interpolation if time is outside # of the defined time range or if we have hit an # exact value in the time array. - value = self._spatial_interpolation(ti, z, y, x, self.grid.time[ti], particle=particle) + value = self._spatial_interpolation(self.grid.time[ti], z, y, x, ti, particle=particle) if applyConversion: return self.units.to_target(value, z, y, x) From 3abf522ec4ee10facdfa73dccd25e2990571e275 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Thu, 27 Feb 2025 14:57:48 +0100 Subject: [PATCH 3/3] Calculating ti in test Kernel --- tests/test_fieldset.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 2e31060a9b..4ab303b949 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -571,7 +571,8 @@ def test_fieldset_write(tmp_zarrfile): def UpdateU(particle, fieldset, time): # pragma: no cover tmp1, tmp2 = fieldset.UV[particle] _, yi, xi = fieldset.U.unravel_index(particle.ei) - fieldset.U.data[0, yi, xi] += 1 + ti = fieldset.U._time_index(time) + fieldset.U.data[ti, yi, xi] += 1 fieldset.U.grid.time[0] = time pset = ParticleSet(fieldset, pclass=Particle, lon=5, lat=5)