From 2ebc4584604d82be428b552c28dee89218da8d65 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 22 Jul 2024 15:14:21 +0200 Subject: [PATCH 1/8] Removing redundant kernel attributes --- parcels/kernel.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/parcels/kernel.py b/parcels/kernel.py index 1c9758440..47e0e53aa 100644 --- a/parcels/kernel.py +++ b/parcels/kernel.py @@ -55,8 +55,6 @@ def __init__(self, fieldset, ptype, pyfunc=None, funcname=None, funccode=None, p self._ptype = ptype self._lib = None self.delete_cfiles = delete_cfiles - self._cleanup_files = None - self._cleanup_lib = None self._c_include = c_include # Derive meta information from pyfunc, if not given From 55071fef7d18d22be71aac97c59f9c02a207449e Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 22 Jul 2024 15:25:10 +0200 Subject: [PATCH 2/8] Using particle_ddepth in analyticaladvection --- parcels/application_kernels/advection.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/application_kernels/advection.py b/parcels/application_kernels/advection.py index c661b4610..d90c8c3c0 100644 --- a/parcels/application_kernels/advection.py +++ b/parcels/application_kernels/advection.py @@ -265,7 +265,7 @@ def compute_rs(ds, r, B, delta, s_min): if withW: rs_z = compute_rs(ds_z, zeta, B_z, delta_z, s_min) - particle.depth = (1.-rs_z) * pz[0] + rs_z * pz[1] + particle_ddepth = (1.-rs_z) * pz[0] + rs_z * pz[1] - particle.depth # noqa if particle.dt > 0: particle.dt = max(direction * s_min * (dxdy * dz), 1e-7) From e621047e185e91475233880b6689152867ca0237 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 22 Jul 2024 15:53:42 +0200 Subject: [PATCH 3/8] Fixing way to assign xiyizi in codegenerator --- parcels/compilation/codegenerator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/compilation/codegenerator.py b/parcels/compilation/codegenerator.py index 1081ff98a..9678993f8 100644 --- a/parcels/compilation/codegenerator.py +++ b/parcels/compilation/codegenerator.py @@ -657,7 +657,7 @@ def visit_Subscript(self, node): if isinstance(node.value, FieldNode) or isinstance(node.value, VectorFieldNode): node.ccode = node.value.__getitem__(node.slice.ccode).ccode elif isinstance(node.value, ParticleXiYiZiTiAttributeNode): - node.ccode = f"{node.value.obj}->{node.value.attr}[pnum, {node.slice.ccode}]" + node.ccode = f"{node.value.obj}->{node.value.attr}[{node.slice.ccode}]" elif isinstance(node.value, IntrinsicNode): raise NotImplementedError(f"Subscript not implemented for object type {type(node.value).__name__}") else: From fb689ee59db44e736432b252afe9d967768dab1e Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 23 Jul 2024 16:42:55 +0200 Subject: [PATCH 4/8] Removing redundant argument from helper function in AnalyticalAdvection --- parcels/application_kernels/advection.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/application_kernels/advection.py b/parcels/application_kernels/advection.py index d90c8c3c0..79a64e86b 100644 --- a/parcels/application_kernels/advection.py +++ b/parcels/application_kernels/advection.py @@ -251,7 +251,7 @@ def compute_ds(F0, F1, r, direction, tol): s_min = min(abs(ds_x), abs(ds_y), abs(ds_z), abs(ds_t / (dxdy * dz))) # calculate end position in time s_min - def compute_rs(ds, r, B, delta, s_min): + def compute_rs(r, B, delta, s_min): if abs(B) < tol: return -delta * s_min + r else: From f3f54b7b8f5c30d76016e1cfc6e48b6e7bd94ffd Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Thu, 25 Jul 2024 08:30:11 +0200 Subject: [PATCH 5/8] Fixing xi sampling code, and making unit test more subtle --- parcels/compilation/codegenerator.py | 3 ++- tests/test_particlefile.py | 35 +++++++++++++++------------- 2 files changed, 21 insertions(+), 17 deletions(-) mode change 100644 => 100755 tests/test_particlefile.py diff --git a/parcels/compilation/codegenerator.py b/parcels/compilation/codegenerator.py index 9678993f8..197ec1279 100644 --- a/parcels/compilation/codegenerator.py +++ b/parcels/compilation/codegenerator.py @@ -657,7 +657,8 @@ def visit_Subscript(self, node): if isinstance(node.value, FieldNode) or isinstance(node.value, VectorFieldNode): node.ccode = node.value.__getitem__(node.slice.ccode).ccode elif isinstance(node.value, ParticleXiYiZiTiAttributeNode): - node.ccode = f"{node.value.obj}->{node.value.attr}[{node.slice.ccode}]" + ngrid = str(self.fieldset.gridset.size if self.fieldset is not None else 1) + node.ccode = f"{node.value.obj}->{node.value.attr}[pnum*{ngrid}+{node.slice.ccode}]" elif isinstance(node.value, IntrinsicNode): raise NotImplementedError(f"Subscript not implemented for object type {type(node.value).__name__}") else: diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py old mode 100644 new mode 100755 index 22f2c35e7..22e05f8dc --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -284,7 +284,7 @@ def Update_lon(particle, fieldset, time): def test_write_xiyi(fieldset, mode, tmpdir): outfilepath = tmpdir.join("pfile_xiyi.zarr") fieldset.U.data[:] = 1 # set a non-zero zonal velocity - fieldset.add_field(Field(name='P', data=np.zeros((2, 20)), lon=np.linspace(0, 1, 20), lat=[0, 2])) + fieldset.add_field(Field(name='P', data=np.zeros((3, 20)), lon=np.linspace(0, 1, 20), lat=[-2, 0, 2])) dt = 3600 XiYiParticle = ptype[mode].add_variables([ @@ -306,24 +306,27 @@ def SampleP(particle, fieldset, time): if time > 5*3600: tmp = fieldset.P[particle] # noqa - pset = ParticleSet(fieldset, pclass=XiYiParticle, lon=[0], lat=[0.2], lonlatdepth_dtype=np.float64) + pset = ParticleSet(fieldset, pclass=XiYiParticle, lon=[0, 0.2], lat=[0.2, 1], lonlatdepth_dtype=np.float64) pfile = pset.ParticleFile(name=outfilepath, outputdt=dt) - pset.execute([Get_XiYi, SampleP, AdvectionRK4], endtime=10*dt, dt=dt, output_file=pfile) + pset.execute([SampleP, Get_XiYi, AdvectionRK4], endtime=10*dt, dt=dt, output_file=pfile) ds = xr.open_zarr(outfilepath) - pxi0 = ds['pxi0'][:].values[0].astype(np.int32) - pxi1 = ds['pxi1'][:].values[0].astype(np.int32) - lons = ds['lon'][:].values[0] - pyi = ds['pyi'][:].values[0].astype(np.int32) - lats = ds['lat'][:].values[0] - - assert (pxi0[0] == 0) and (pxi0[-1] == 11) # check that particle has moved - assert np.all(pxi1[:7] == 0) # check that particle has not been sampled on grid 1 until time 6 - assert np.all(pxi1[7:] > 0) # check that particle has not been sampled on grid 1 after time 6 - for xi, lon in zip(pxi0[1:], lons[1:]): - assert fieldset.U.grid.lon[xi] <= lon < fieldset.U.grid.lon[xi+1] - for yi, lat in zip(pyi[1:], lats[1:]): - assert fieldset.U.grid.lat[yi] <= lat < fieldset.U.grid.lat[yi+1] + pxi0 = ds['pxi0'][:].values.astype(np.int32) + pxi1 = ds['pxi1'][:].values.astype(np.int32) + lons = ds['lon'][:].values + pyi = ds['pyi'][:].values.astype(np.int32) + lats = ds['lat'][:].values + + for p in range(pyi.shape[0]): + assert (pxi0[p, 0] == 0) and (pxi0[p, -1] == pset[p].pxi0) # check that particle has moved + assert np.all(pxi1[p, :6] == 0) # check that particle has not been sampled on grid 1 until time 6 + assert np.all(pxi1[p, 6:] > 0) # check that particle has not been sampled on grid 1 after time 6 + for xi, lon in zip(pxi0[p, 1:], lons[p, 1:]): + assert fieldset.U.grid.lon[xi] <= lon < fieldset.U.grid.lon[xi+1] + for xi, lon in zip(pxi1[p, 6:], lons[p, 6:]): + assert fieldset.P.grid.lon[xi] <= lon < fieldset.P.grid.lon[xi+1] + for yi, lat in zip(pyi[p, 1:], lats[p, 1:]): + assert fieldset.U.grid.lat[yi] <= lat < fieldset.U.grid.lat[yi+1] ds.close() From 0257a4023e5097191e11bcfb9a7daf3b8138b376 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 26 Jul 2024 10:25:26 +0200 Subject: [PATCH 6/8] FIxing Scipy AnalyticalAdvection to add to dlon,dlat,ddepth instead of overwrite --- parcels/application_kernels/advection.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/parcels/application_kernels/advection.py b/parcels/application_kernels/advection.py index 79a64e86b..a4c5b7753 100644 --- a/parcels/application_kernels/advection.py +++ b/parcels/application_kernels/advection.py @@ -260,12 +260,12 @@ def compute_rs(r, B, delta, s_min): rs_x = compute_rs(ds_x, xsi, B_x, delta_x, s_min) rs_y = compute_rs(ds_y, eta, B_y, delta_y, s_min) - particle_dlon = (1.-rs_x)*(1.-rs_y) * px[0] + rs_x * (1.-rs_y) * px[1] + rs_x * rs_y * px[2] + (1.-rs_x)*rs_y * px[3] - particle.lon # noqa - particle_dlat = (1.-rs_x)*(1.-rs_y) * py[0] + rs_x * (1.-rs_y) * py[1] + rs_x * rs_y * py[2] + (1.-rs_x)*rs_y * py[3] - particle.lat # noqa + particle_dlon += (1.-rs_x)*(1.-rs_y) * px[0] + rs_x * (1.-rs_y) * px[1] + rs_x * rs_y * px[2] + (1.-rs_x)*rs_y * px[3] - particle.lon # noqa + particle_dlat += (1.-rs_x)*(1.-rs_y) * py[0] + rs_x * (1.-rs_y) * py[1] + rs_x * rs_y * py[2] + (1.-rs_x)*rs_y * py[3] - particle.lat # noqa if withW: rs_z = compute_rs(ds_z, zeta, B_z, delta_z, s_min) - particle_ddepth = (1.-rs_z) * pz[0] + rs_z * pz[1] - particle.depth # noqa + particle_ddepth += (1.-rs_z) * pz[0] + rs_z * pz[1] - particle.depth # noqa if particle.dt > 0: particle.dt = max(direction * s_min * (dxdy * dz), 1e-7) From 104fbc4dd2ee0ae5d92cdf4359fe528044b0d12e Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 26 Jul 2024 11:32:36 +0200 Subject: [PATCH 7/8] Fixing calls to compute_rs now that ds is removed as argument --- parcels/application_kernels/advection.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/parcels/application_kernels/advection.py b/parcels/application_kernels/advection.py index a4c5b7753..239d755fe 100644 --- a/parcels/application_kernels/advection.py +++ b/parcels/application_kernels/advection.py @@ -257,14 +257,14 @@ def compute_rs(r, B, delta, s_min): else: return (r + delta / B) * math.exp(-B * s_min) - delta / B - rs_x = compute_rs(ds_x, xsi, B_x, delta_x, s_min) - rs_y = compute_rs(ds_y, eta, B_y, delta_y, s_min) + rs_x = compute_rs(xsi, B_x, delta_x, s_min) + rs_y = compute_rs(eta, B_y, delta_y, s_min) particle_dlon += (1.-rs_x)*(1.-rs_y) * px[0] + rs_x * (1.-rs_y) * px[1] + rs_x * rs_y * px[2] + (1.-rs_x)*rs_y * px[3] - particle.lon # noqa particle_dlat += (1.-rs_x)*(1.-rs_y) * py[0] + rs_x * (1.-rs_y) * py[1] + rs_x * rs_y * py[2] + (1.-rs_x)*rs_y * py[3] - particle.lat # noqa if withW: - rs_z = compute_rs(ds_z, zeta, B_z, delta_z, s_min) + rs_z = compute_rs(zeta, B_z, delta_z, s_min) particle_ddepth += (1.-rs_z) * pz[0] + rs_z * pz[1] - particle.depth # noqa if particle.dt > 0: From c66c28551e05e12704e9fae02de1e564490cb57e Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 29 Jul 2024 08:18:05 +0200 Subject: [PATCH 8/8] Changing tmp variable to _ As suggested by @VeckoTheGecko --- tests/test_particlefile.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index 22e05f8dc..cf19dd10c 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -304,7 +304,7 @@ def Get_XiYi(particle, fieldset, time): def SampleP(particle, fieldset, time): if time > 5*3600: - tmp = fieldset.P[particle] # noqa + _ = fieldset.P[particle] # To trigger sampling of the P field pset = ParticleSet(fieldset, pclass=XiYiParticle, lon=[0, 0.2], lat=[0.2, 1], lonlatdepth_dtype=np.float64) pfile = pset.ParticleFile(name=outfilepath, outputdt=dt)