From 2ebc4584604d82be428b552c28dee89218da8d65 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 22 Jul 2024 15:14:21 +0200 Subject: [PATCH 01/11] Removing redundant kernel attributes --- parcels/kernel.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/parcels/kernel.py b/parcels/kernel.py index 1c9758440a..47e0e53aa1 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 02/11] 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 c661b4610d..d90c8c3c0e 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 03/11] 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 1081ff98ab..9678993f8e 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 04/11] 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 d90c8c3c0e..79a64e86be 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 05/11] 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 9678993f8e..197ec12798 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 22f2c35e71..22e05f8dc9 --- 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 06/11] 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 79a64e86be..a4c5b77532 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 07/11] 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 a4c5b77532..239d755fe0 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 32ac7e87ec0771c8ad1d8a48a2f4ce86af1d3fa9 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 26 Jul 2024 13:24:56 +0200 Subject: [PATCH 08/11] Refactoring GridCode to be called GridType The terms gridtype and gridcode were used interchangeably in the Parcels codebase, leading to confusion. Therefore, this refactoring changes all gridcode to gridtype --- parcels/field.py | 26 +++++----- parcels/grid.py | 12 ++--- parcels/include/index_search.h | 18 +++---- parcels/include/parcels.h | 72 ++++++++++++++-------------- parcels/kernel.py | 4 +- parcels/particleset.py | 4 +- tests/test_data/create_testfields.py | 6 +-- 7 files changed, 71 insertions(+), 71 deletions(-) diff --git a/parcels/field.py b/parcels/field.py index ece0f9c53b..7af24caf10 100644 --- a/parcels/field.py +++ b/parcels/field.py @@ -31,7 +31,7 @@ DeferredNetcdfFileBuffer, NetcdfFileBuffer, ) -from .grid import CGrid, Grid, GridCode +from .grid import CGrid, Grid, GridType __all__ = ['Field', 'VectorField', 'NestedField'] @@ -178,7 +178,7 @@ def __init__(self, name, data, lon=None, lat=None, depth=None, time=None, grid=N self.interp_method = interp_method self.gridindexingtype = gridindexingtype if self.interp_method in ['bgrid_velocity', 'bgrid_w_velocity', 'bgrid_tracer'] and \ - self.grid.gtype in [GridCode.RectilinearSGrid, GridCode.CurvilinearSGrid]: + self.grid.gtype in [GridType.RectilinearSGrid, GridType.CurvilinearSGrid]: logger.warning_once('General s-levels are not supported in B-grid. RectilinearSGrid and CurvilinearSGrid can still be used to deal with shaved cells, but the levels must be horizontal.') self.fieldset = None @@ -687,7 +687,7 @@ def calc_cell_edge_sizes(self): Currently only works for Rectilinear Grids """ if not self.grid.cell_edge_sizes: - if self.grid.gtype in (GridCode.RectilinearZGrid, GridCode.RectilinearSGrid): + if self.grid.gtype in (GridType.RectilinearZGrid, GridType.RectilinearSGrid): self.grid.cell_edge_sizes['x'] = np.zeros((self.grid.ydim, self.grid.xdim), dtype=np.float32) self.grid.cell_edge_sizes['y'] = np.zeros((self.grid.ydim, self.grid.xdim), dtype=np.float32) @@ -877,7 +877,7 @@ def search_indices_rectilinear(self, x, y, z, ti=-1, time=-1, particle=None, sea yi, eta = -1, 0 if grid.zdim > 1 and not search2D: - if grid.gtype == GridCode.RectilinearZGrid: + if grid.gtype == GridType.RectilinearZGrid: # Never passes here, because in this case, we work with scipy try: (zi, zeta) = self.search_indices_vertical_z(z) @@ -885,7 +885,7 @@ def search_indices_rectilinear(self, x, y, z, ti=-1, time=-1, particle=None, sea raise FieldOutOfBoundError(x, y, z, field=self) except FieldOutOfBoundSurfaceError: raise FieldOutOfBoundSurfaceError(x, y, z, field=self) - elif grid.gtype == GridCode.RectilinearSGrid: + elif grid.gtype == GridType.RectilinearSGrid: (zi, zeta) = self.search_indices_vertical_s(x, y, z, xi, yi, xsi, eta, ti, time) else: zi, zeta = -1, 0 @@ -973,12 +973,12 @@ def search_indices_curvilinear(self, x, y, z, ti=-1, time=-1, particle=None, sea eta = min(1., eta) if grid.zdim > 1 and not search2D: - if grid.gtype == GridCode.CurvilinearZGrid: + if grid.gtype == GridType.CurvilinearZGrid: try: (zi, zeta) = self.search_indices_vertical_z(z) except FieldOutOfBoundError: raise FieldOutOfBoundError(x, y, z, field=self) - elif grid.gtype == GridCode.CurvilinearSGrid: + elif grid.gtype == GridType.CurvilinearSGrid: (zi, zeta) = self.search_indices_vertical_s(x, y, z, xi, yi, xsi, eta, ti, time) else: zi = -1 @@ -995,7 +995,7 @@ def search_indices_curvilinear(self, x, y, z, ti=-1, time=-1, particle=None, sea return (xsi, eta, zeta, xi, yi, zi) def search_indices(self, x, y, z, ti=-1, time=-1, particle=None, search2D=False): - if self.grid.gtype in [GridCode.RectilinearSGrid, GridCode.RectilinearZGrid]: + if self.grid.gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]: return self.search_indices_rectilinear(x, y, z, ti, time, particle=particle, search2D=search2D) else: return self.search_indices_curvilinear(x, y, z, ti, time, particle=particle, search2D=search2D) @@ -1399,12 +1399,12 @@ def write(self, filename, varname=None): vname_depth = 'depth%s' % self.name.lower() # Create DataArray objects for file I/O - if self.grid.gtype == GridCode.RectilinearZGrid: + if self.grid.gtype == GridType.RectilinearZGrid: nav_lon = xr.DataArray(self.grid.lon + np.zeros((self.grid.ydim, self.grid.xdim), dtype=np.float32), coords=[('y', self.grid.lat), ('x', self.grid.lon)]) nav_lat = xr.DataArray(self.grid.lat.reshape(self.grid.ydim, 1) + np.zeros(self.grid.xdim, dtype=np.float32), coords=[('y', self.grid.lat), ('x', self.grid.lon)]) - elif self.grid.gtype == GridCode.CurvilinearZGrid: + elif self.grid.gtype == GridType.CurvilinearZGrid: nav_lon = xr.DataArray(self.grid.lon, coords=[('y', range(self.grid.ydim)), ('x', range(self.grid.xdim))]) nav_lat = xr.DataArray(self.grid.lat, coords=[('y', range(self.grid.ydim)), @@ -1553,7 +1553,7 @@ def spatial_c_grid_interpolation2D(self, ti, z, y, x, time, particle=None, apply grid = self.U.grid (xsi, eta, zeta, xi, yi, zi) = self.U.search_indices(x, y, z, ti, time, particle=particle) - if grid.gtype in [GridCode.RectilinearSGrid, GridCode.RectilinearZGrid]: + if grid.gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]: px = np.array([grid.lon[xi], grid.lon[xi+1], grid.lon[xi+1], grid.lon[xi]]) py = np.array([grid.lat[yi], grid.lat[yi], grid.lat[yi+1], grid.lat[yi+1]]) else: @@ -1621,7 +1621,7 @@ def spatial_c_grid_interpolation3D_full(self, ti, z, y, x, time, particle=None): grid = self.U.grid (xsi, eta, zet, xi, yi, zi) = self.U.search_indices(x, y, z, ti, time, particle=particle) - if grid.gtype in [GridCode.RectilinearSGrid, GridCode.RectilinearZGrid]: + if grid.gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]: px = np.array([grid.lon[xi], grid.lon[xi+1], grid.lon[xi+1], grid.lon[xi]]) py = np.array([grid.lat[yi], grid.lat[yi], grid.lat[yi+1], grid.lat[yi+1]]) else: @@ -1737,7 +1737,7 @@ def spatial_c_grid_interpolation3D(self, ti, z, y, x, time, particle=None, apply interpolating linearly V depending on the latitude coordinate. Curvilinear grids are treated properly, since the element is projected to a rectilinear parent element. """ - if self.U.grid.gtype in [GridCode.RectilinearSGrid, GridCode.CurvilinearSGrid]: + if self.U.grid.gtype in [GridType.RectilinearSGrid, GridType.CurvilinearSGrid]: (u, v, w) = self.spatial_c_grid_interpolation3D_full(ti, z, y, x, time, particle=particle) else: (u, v) = self.spatial_c_grid_interpolation2D(ti, z, y, x, time, particle=particle) diff --git a/parcels/grid.py b/parcels/grid.py index cc9bcaf2da..5693a09264 100644 --- a/parcels/grid.py +++ b/parcels/grid.py @@ -7,10 +7,10 @@ from parcels.tools.converters import TimeConverter from parcels.tools.loggers import logger -__all__ = ['GridCode', 'RectilinearZGrid', 'RectilinearSGrid', 'CurvilinearZGrid', 'CurvilinearSGrid', 'CGrid', 'Grid'] +__all__ = ['GridType', 'RectilinearZGrid', 'RectilinearSGrid', 'CurvilinearZGrid', 'CurvilinearSGrid', 'CGrid', 'Grid'] -class GridCode(IntEnum): +class GridType(IntEnum): RectilinearZGrid = 0 RectilinearSGrid = 1 CurvilinearZGrid = 2 @@ -326,7 +326,7 @@ def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh='flat if isinstance(depth, np.ndarray): assert (len(depth.shape) <= 1), 'depth is not a vector' - self.gtype = GridCode.RectilinearZGrid + self.gtype = GridType.RectilinearZGrid self.depth = np.zeros(1, dtype=np.float32) if depth is None else depth if not self.depth.flags['C_CONTIGUOUS']: self.depth = np.array(self.depth, order='C') @@ -371,7 +371,7 @@ def __init__(self, lon, lat, depth, time=None, time_origin=None, mesh='flat'): super().__init__(lon, lat, time, time_origin, mesh) assert (isinstance(depth, np.ndarray) and len(depth.shape) in [3, 4]), 'depth is not a 3D or 4D numpy array' - self.gtype = GridCode.RectilinearSGrid + self.gtype = GridType.RectilinearSGrid self.depth = depth if not self.depth.flags['C_CONTIGUOUS']: self.depth = np.array(self.depth, order='C') @@ -481,7 +481,7 @@ def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh='flat if isinstance(depth, np.ndarray): assert (len(depth.shape) == 1), 'depth is not a vector' - self.gtype = GridCode.CurvilinearZGrid + self.gtype = GridType.CurvilinearZGrid self.depth = np.zeros(1, dtype=np.float32) if depth is None else depth if not self.depth.flags['C_CONTIGUOUS']: self.depth = np.array(self.depth, order='C') @@ -525,7 +525,7 @@ def __init__(self, lon, lat, depth, time=None, time_origin=None, mesh='flat'): super().__init__(lon, lat, time, time_origin, mesh) assert (isinstance(depth, np.ndarray) and len(depth.shape) in [3, 4]), 'depth is not a 4D numpy array' - self.gtype = GridCode.CurvilinearSGrid + self.gtype = GridType.CurvilinearSGrid self.depth = depth # should be a C-contiguous array of floats if not self.depth.flags['C_CONTIGUOUS']: self.depth = np.array(self.depth, order='C') diff --git a/parcels/include/index_search.h b/parcels/include/index_search.h index 39b454cfaa..bb42b2db60 100644 --- a/parcels/include/index_search.h +++ b/parcels/include/index_search.h @@ -57,7 +57,7 @@ typedef enum typedef enum { RECTILINEAR_Z_GRID=0, RECTILINEAR_S_GRID=1, CURVILINEAR_Z_GRID=2, CURVILINEAR_S_GRID=3 - } GridCode; + } GridType; // equal/closeness comparison that is equal to numpy (double) static inline bool is_close_dbl(double a, double b) { @@ -206,7 +206,7 @@ static inline void reconnect_bnd_indices(int *xi, int *yi, int xdim, int ydim, i } -static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridCode gcode, +static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridType gtype, int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta, int ti, double time, double t0, double t1, int interp_method, int gridindexingtype) @@ -284,7 +284,7 @@ static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y, StatusCode status; if (zdim > 1){ - switch(gcode){ + switch(gtype){ case RECTILINEAR_Z_GRID: status = search_indices_vertical_z(z, zdim, zvals, zi, zeta, gridindexingtype); break; @@ -316,7 +316,7 @@ static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y, } -static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridCode gcode, +static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridType gtype, int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta, int ti, double time, double t0, double t1, int interp_method, int gridindexingtype) @@ -428,7 +428,7 @@ static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y, StatusCode status; if (zdim > 1){ - switch(gcode){ + switch(gtype){ case CURVILINEAR_Z_GRID: status = search_indices_vertical_z(z, zdim, zvals, zi, zeta, gridindexingtype); break; @@ -457,18 +457,18 @@ static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y, * */ static inline StatusCode search_indices(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta, - GridCode gcode, int ti, double time, double t0, double t1, int interp_method, + GridType gtype, int ti, double time, double t0, double t1, int interp_method, int gridindexingtype) { - switch(gcode){ + switch(gtype){ case RECTILINEAR_Z_GRID: case RECTILINEAR_S_GRID: - return search_indices_rectilinear(x, y, z, grid, gcode, xi, yi, zi, xsi, eta, zeta, + return search_indices_rectilinear(x, y, z, grid, gtype, xi, yi, zi, xsi, eta, zeta, ti, time, t0, t1, interp_method, gridindexingtype); break; case CURVILINEAR_Z_GRID: case CURVILINEAR_S_GRID: - return search_indices_curvilinear(x, y, z, grid, gcode, xi, yi, zi, xsi, eta, zeta, + return search_indices_curvilinear(x, y, z, grid, gtype, xi, yi, zi, xsi, eta, zeta, ti, time, t0, t1, interp_method, gridindexingtype); break; default: diff --git a/parcels/include/parcels.h b/parcels/include/parcels.h index 61775a8192..e4d8431e31 100644 --- a/parcels/include/parcels.h +++ b/parcels/include/parcels.h @@ -435,7 +435,7 @@ static inline StatusCode getCell3D(CField *f, int xi, int yi, int zi, int ti, fl /* Linear interpolation along the time axis */ static inline StatusCode temporal_interpolation_structured_grid(type_coord x, type_coord y, type_coord z, double time, CField *f, - GridCode gcode, int *xi, int *yi, int *zi, int *ti, + GridType gtype, int *xi, int *yi, int *zi, int *ti, float *value, int interp_method, int gridindexingtype) { StatusCode status; @@ -466,7 +466,7 @@ static inline StatusCode temporal_interpolation_structured_grid(type_coord x, ty double tsrch = (tii == 2) ? time : t0; status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], - &xsi, &eta, &zeta, gcode, ti[igrid], + &xsi, &eta, &zeta, gtype, ti[igrid], tsrch, t0, t1, interp_method, gridindexingtype); CHECKSTATUS(status); @@ -564,7 +564,7 @@ static double dist(double lon1, double lon2, double lat1, double lat2, int spher /* Linear interpolation routine for 2D C grid */ static inline StatusCode spatial_interpolation_UV_c_grid(double xsi, double eta, int xi, int yi, CStructuredGrid *grid, - GridCode gcode, float dataU[2][2], float dataV[2][2], float *u, float *v) + GridType gtype, float dataU[2][2], float dataV[2][2], float *u, float *v) { /* Cast data array into data[lat][lon] as per NEMO convention */ int xdim = grid->xdim; @@ -572,7 +572,7 @@ static inline StatusCode spatial_interpolation_UV_c_grid(double xsi, double eta, double xgrid_loc[4]; double ygrid_loc[4]; int iN; - if( (gcode == RECTILINEAR_Z_GRID) || (gcode == RECTILINEAR_S_GRID) ){ + if( (gtype == RECTILINEAR_Z_GRID) || (gtype == RECTILINEAR_S_GRID) ){ float *xgrid = grid->lon; float *ygrid = grid->lat; for (iN=0; iN < 4; ++iN){ @@ -643,7 +643,7 @@ static inline StatusCode spatial_interpolation_UV_c_grid(double xsi, double eta, static inline StatusCode temporal_interpolationUV_c_grid(type_coord x, type_coord y, type_coord z, double time, CField *U, CField *V, - GridCode gcode, int *xi, int *yi, int *zi, int *ti, + GridType gtype, int *xi, int *yi, int *zi, int *ti, float *u, float *v, int gridindexingtype) { StatusCode status; @@ -663,7 +663,7 @@ static inline StatusCode temporal_interpolationUV_c_grid(type_coord x, type_coor float u0, u1, v0, v1; double t0 = grid->time[ti[igrid]]; double t1 = grid->time[ti[igrid]+1]; /* Identify grid cell to sample through local linear search */ - status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gcode, ti[igrid], time, t0, t1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); + status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gtype, ti[igrid], time, t0, t1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); if (grid->zdim==1){ float data2D_U[2][2][2], data2D_V[2][2][2]; if (gridindexingtype == NEMO) { @@ -674,8 +674,8 @@ static inline StatusCode temporal_interpolationUV_c_grid(type_coord x, type_coor status = getCell2D(U, xi[igrid], yi[igrid]-1, ti[igrid], data2D_U, 0); CHECKSTATUS(status); status = getCell2D(V, xi[igrid]-1, yi[igrid], ti[igrid], data2D_V, 0); CHECKSTATUS(status); } - status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gcode, data2D_U[0], data2D_V[0], &u0, &v0); CHECKSTATUS(status); - status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gcode, data2D_U[1], data2D_V[1], &u1, &v1); CHECKSTATUS(status); + status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gtype, data2D_U[0], data2D_V[0], &u0, &v0); CHECKSTATUS(status); + status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gtype, data2D_U[1], data2D_V[1], &u1, &v1); CHECKSTATUS(status); } else { float data3D_U[2][2][2][2], data3D_V[2][2][2][2]; @@ -687,15 +687,15 @@ static inline StatusCode temporal_interpolationUV_c_grid(type_coord x, type_coor status = getCell3D(U, xi[igrid], yi[igrid]-1, zi[igrid], ti[igrid], data3D_U, 0); CHECKSTATUS(status); status = getCell3D(V, xi[igrid]-1, yi[igrid], zi[igrid], ti[igrid], data3D_V, 0); CHECKSTATUS(status); } - status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gcode, data3D_U[0][0], data3D_V[0][0], &u0, &v0); CHECKSTATUS(status); - status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gcode, data3D_U[1][0], data3D_V[1][0], &u1, &v1); CHECKSTATUS(status); + status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gtype, data3D_U[0][0], data3D_V[0][0], &u0, &v0); CHECKSTATUS(status); + status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gtype, data3D_U[1][0], data3D_V[1][0], &u1, &v1); CHECKSTATUS(status); } *u = u0 + (u1 - u0) * (float)((time - t0) / (t1 - t0)); *v = v0 + (v1 - v0) * (float)((time - t0) / (t1 - t0)); return SUCCESS; } else { double t0 = grid->time[ti[igrid]]; - status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gcode, ti[igrid], t0, t0, t0+1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); + status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gtype, ti[igrid], t0, t0, t0+1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); if (grid->zdim==1){ float data2D_U[2][2][2], data2D_V[2][2][2]; if (gridindexingtype == NEMO) { @@ -706,7 +706,7 @@ static inline StatusCode temporal_interpolationUV_c_grid(type_coord x, type_coor status = getCell2D(U, xi[igrid], yi[igrid]-1, ti[igrid], data2D_U, 1); CHECKSTATUS(status); status = getCell2D(V, xi[igrid]-1, yi[igrid], ti[igrid], data2D_V, 1); CHECKSTATUS(status); } - status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gcode, data2D_U[0], data2D_V[0], u, v); CHECKSTATUS(status); + status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gtype, data2D_U[0], data2D_V[0], u, v); CHECKSTATUS(status); } else{ float data3D_U[2][2][2][2], data3D_V[2][2][2][2]; @@ -718,7 +718,7 @@ static inline StatusCode temporal_interpolationUV_c_grid(type_coord x, type_coor status = getCell3D(U, xi[igrid], yi[igrid]-1, zi[igrid], ti[igrid], data3D_U, 1); CHECKSTATUS(status); status = getCell3D(V, xi[igrid]-1, yi[igrid], zi[igrid], ti[igrid], data3D_V, 1); CHECKSTATUS(status); } - status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gcode, data3D_U[0][0], data3D_V[0][0], u, v); CHECKSTATUS(status); + status = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gtype, data3D_U[0][0], data3D_V[0][0], u, v); CHECKSTATUS(status); } return SUCCESS; } @@ -726,7 +726,7 @@ static inline StatusCode temporal_interpolationUV_c_grid(type_coord x, type_coor /* Quadratic interpolation routine for 3D C grid */ static inline StatusCode spatial_interpolation_UVW_c_grid(double xsi, double eta, double zet, int xi, int yi, int zi, int ti, CStructuredGrid *grid, - GridCode gcode, float dataU[2][2][2], float dataV[2][2][2], float dataW[2][2][2], float *u, float *v, float *w, int gridindexingtype) + GridType gtype, float dataU[2][2][2], float dataV[2][2][2], float dataW[2][2][2], float *u, float *v, float *w, int gridindexingtype) { /* Cast data array into data[lat][lon] as per NEMO convention */ int xdim = grid->xdim; @@ -736,7 +736,7 @@ static inline StatusCode spatial_interpolation_UVW_c_grid(double xsi, double eta float xgrid_loc[4]; float ygrid_loc[4]; int iN; - if( gcode == RECTILINEAR_S_GRID ){ + if( gtype == RECTILINEAR_S_GRID ){ float *xgrid = grid->lon; float *ygrid = grid->lat; for (iN=0; iN < 4; ++iN){ @@ -861,7 +861,7 @@ static inline StatusCode spatial_interpolation_UVW_c_grid(double xsi, double eta } static inline StatusCode temporal_interpolationUVW_c_grid(type_coord x, type_coord y, type_coord z, double time, CField *U, CField *V, CField *W, - GridCode gcode, int *xi, int *yi, int *zi, int *ti, + GridType gtype, int *xi, int *yi, int *zi, int *ti, float *u, float *v, float *w, int gridindexingtype) { StatusCode status; @@ -884,15 +884,15 @@ static inline StatusCode temporal_interpolationUVW_c_grid(type_coord x, type_coo float u0, u1, v0, v1, w0, w1; double t0 = grid->time[ti[igrid]]; double t1 = grid->time[ti[igrid]+1]; /* Identify grid cell to sample through local linear search */ - status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zet, gcode, ti[igrid], time, t0, t1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); + status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zet, gtype, ti[igrid], time, t0, t1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); status = getCell3D(U, xi[igrid], yi[igrid], zi[igrid], ti[igrid], data3D_U, 0); CHECKSTATUS(status); status = getCell3D(V, xi[igrid], yi[igrid], zi[igrid], ti[igrid], data3D_V, 0); CHECKSTATUS(status); status = getCell3D(W, xi[igrid], yi[igrid], zi[igrid], ti[igrid], data3D_W, 0); CHECKSTATUS(status); if (grid->zdim==1){ return ERROR; } else { - status = spatial_interpolation_UVW_c_grid(xsi, eta, zet, xi[igrid], yi[igrid], zi[igrid], ti[igrid], grid, gcode, data3D_U[0], data3D_V[0], data3D_W[0], &u0, &v0, &w0, gridindexingtype); CHECKSTATUS(status); - status = spatial_interpolation_UVW_c_grid(xsi, eta, zet, xi[igrid], yi[igrid], zi[igrid], ti[igrid]+1, grid, gcode, data3D_U[1], data3D_V[1], data3D_W[1], &u1, &v1, &w1, gridindexingtype); CHECKSTATUS(status); + status = spatial_interpolation_UVW_c_grid(xsi, eta, zet, xi[igrid], yi[igrid], zi[igrid], ti[igrid], grid, gtype, data3D_U[0], data3D_V[0], data3D_W[0], &u0, &v0, &w0, gridindexingtype); CHECKSTATUS(status); + status = spatial_interpolation_UVW_c_grid(xsi, eta, zet, xi[igrid], yi[igrid], zi[igrid], ti[igrid]+1, grid, gtype, data3D_U[1], data3D_V[1], data3D_W[1], &u1, &v1, &w1, gridindexingtype); CHECKSTATUS(status); } *u = u0 + (u1 - u0) * (float)((time - t0) / (t1 - t0)); *v = v0 + (v1 - v0) * (float)((time - t0) / (t1 - t0)); @@ -900,7 +900,7 @@ static inline StatusCode temporal_interpolationUVW_c_grid(type_coord x, type_coo return SUCCESS; } else { double t0 = grid->time[ti[igrid]]; - status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zet, gcode, ti[igrid], t0, t0, t0+1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); + status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zet, gtype, ti[igrid], t0, t0, t0+1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); status = getCell3D(U, xi[igrid], yi[igrid], zi[igrid], ti[igrid], data3D_U, 1); CHECKSTATUS(status); status = getCell3D(V, xi[igrid], yi[igrid], zi[igrid], ti[igrid], data3D_V, 1); CHECKSTATUS(status); status = getCell3D(W, xi[igrid], yi[igrid], zi[igrid], ti[igrid], data3D_W, 1); CHECKSTATUS(status); @@ -908,7 +908,7 @@ static inline StatusCode temporal_interpolationUVW_c_grid(type_coord x, type_coo return ERROR; } else{ - status = spatial_interpolation_UVW_c_grid(xsi, eta, zet, xi[igrid], yi[igrid], zi[igrid], ti[igrid], grid, gcode, data3D_U[0], data3D_V[0], data3D_W[0], u, v, w, gridindexingtype); CHECKSTATUS(status); + status = spatial_interpolation_UVW_c_grid(xsi, eta, zet, xi[igrid], yi[igrid], zi[igrid], ti[igrid], grid, gtype, data3D_U[0], data3D_V[0], data3D_W[0], u, v, w, gridindexingtype); CHECKSTATUS(status); } return SUCCESS; } @@ -1081,7 +1081,7 @@ static inline StatusCode calculate_slip_conditions_3D(double xsi, double eta, do } static inline StatusCode temporal_interpolation_slip(type_coord x, type_coord y, type_coord z, double time, CField *U, CField *V, CField *W, - GridCode gcode, int *xi, int *yi, int *zi, int *ti, + GridType gtype, int *xi, int *yi, int *zi, int *ti, float *u, float *v, float *w, int interp_method, int gridindexingtype, int withW) { StatusCode status; @@ -1100,7 +1100,7 @@ static inline StatusCode temporal_interpolation_slip(type_coord x, type_coord y, float u0, u1, v0, v1, w0, w1; double t0 = grid->time[ti[igrid]]; double t1 = grid->time[ti[igrid]+1]; /* Identify grid cell to sample through local linear search */ - status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gcode, ti[igrid], time, t0, t1, interp_method, gridindexingtype); CHECKSTATUS(status); + status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gtype, ti[igrid], time, t0, t1, interp_method, gridindexingtype); CHECKSTATUS(status); if (grid->zdim==1){ float data2D_U[2][2][2], data2D_V[2][2][2], data2D_W[2][2][2]; status = getCell2D(U, xi[igrid], yi[igrid], ti[igrid], data2D_U, 0); CHECKSTATUS(status); @@ -1143,7 +1143,7 @@ static inline StatusCode temporal_interpolation_slip(type_coord x, type_coord y, } else { double t0 = grid->time[ti[igrid]]; - status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gcode, ti[igrid], t0, t0, t0+1, interp_method, gridindexingtype); CHECKSTATUS(status); + status = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gtype, ti[igrid], t0, t0, t0+1, interp_method, gridindexingtype); CHECKSTATUS(status); if (grid->zdim==1){ float data2D_U[2][2][2], data2D_V[2][2][2], data2D_W[2][2][2]; status = getCell2D(U, xi[igrid], yi[igrid], ti[igrid], data2D_U, 1); CHECKSTATUS(status); @@ -1178,10 +1178,10 @@ static inline StatusCode temporal_interpolation(type_coord x, type_coord y, type float *value, int interp_method, int gridindexingtype) { CGrid *_grid = f->grid; - GridCode gcode = _grid->gtype; + GridType gtype = _grid->gtype; - if (gcode == RECTILINEAR_Z_GRID || gcode == RECTILINEAR_S_GRID || gcode == CURVILINEAR_Z_GRID || gcode == CURVILINEAR_S_GRID) - return temporal_interpolation_structured_grid(x, y, z, time, f, gcode, xi, yi, zi, ti, value, interp_method, gridindexingtype); + if (gtype == RECTILINEAR_Z_GRID || gtype == RECTILINEAR_S_GRID || gtype == CURVILINEAR_Z_GRID || gtype == CURVILINEAR_S_GRID) + return temporal_interpolation_structured_grid(x, y, z, time, f, gtype, xi, yi, zi, ti, value, interp_method, gridindexingtype); else{ printf("Only RECTILINEAR_Z_GRID, RECTILINEAR_S_GRID, CURVILINEAR_Z_GRID and CURVILINEAR_S_GRID grids are currently implemented\n"); return ERROR; @@ -1196,15 +1196,15 @@ static inline StatusCode temporal_interpolationUV(type_coord x, type_coord y, ty StatusCode status; if (interp_method == CGRID_VELOCITY){ CGrid *_grid = U->grid; - GridCode gcode = _grid->gtype; - status = temporal_interpolationUV_c_grid(x, y, z, time, U, V, gcode, xi, yi, zi, ti, valueU, valueV, gridindexingtype); CHECKSTATUS(status); + GridType gtype = _grid->gtype; + status = temporal_interpolationUV_c_grid(x, y, z, time, U, V, gtype, xi, yi, zi, ti, valueU, valueV, gridindexingtype); CHECKSTATUS(status); return SUCCESS; } else if ((interp_method == PARTIALSLIP) || (interp_method == FREESLIP)){ CGrid *_grid = U->grid; CField *W = U; - GridCode gcode = _grid->gtype; + GridType gtype = _grid->gtype; int withW = 0; - status = temporal_interpolation_slip(x, y, z, time, U, V, W, gcode, xi, yi, zi, ti, valueU, valueV, 0, interp_method, gridindexingtype, withW); CHECKSTATUS(status); + status = temporal_interpolation_slip(x, y, z, time, U, V, W, gtype, xi, yi, zi, ti, valueU, valueV, 0, interp_method, gridindexingtype, withW); CHECKSTATUS(status); return SUCCESS; } else { status = temporal_interpolation(x, y, z, time, U, xi, yi, zi, ti, valueU, interp_method, gridindexingtype); CHECKSTATUS(status); @@ -1221,16 +1221,16 @@ static inline StatusCode temporal_interpolationUVW(type_coord x, type_coord y, t StatusCode status; if (interp_method == CGRID_VELOCITY){ CGrid *_grid = U->grid; - GridCode gcode = _grid->gtype; - if (gcode == RECTILINEAR_S_GRID || gcode == CURVILINEAR_S_GRID){ - status = temporal_interpolationUVW_c_grid(x, y, z, time, U, V, W, gcode, xi, yi, zi, ti, valueU, valueV, valueW, gridindexingtype); CHECKSTATUS(status); + GridType gtype = _grid->gtype; + if (gtype == RECTILINEAR_S_GRID || gtype == CURVILINEAR_S_GRID){ + status = temporal_interpolationUVW_c_grid(x, y, z, time, U, V, W, gtype, xi, yi, zi, ti, valueU, valueV, valueW, gridindexingtype); CHECKSTATUS(status); return SUCCESS; } } else if ((interp_method == PARTIALSLIP) || (interp_method == FREESLIP)){ CGrid *_grid = U->grid; - GridCode gcode = _grid->gtype; + GridType gtype = _grid->gtype; int withW = 1; - status = temporal_interpolation_slip(x, y, z, time, U, V, W, gcode, xi, yi, zi, ti, valueU, valueV, valueW, interp_method, gridindexingtype, withW); CHECKSTATUS(status); + status = temporal_interpolation_slip(x, y, z, time, U, V, W, gtype, xi, yi, zi, ti, valueU, valueV, valueW, interp_method, gridindexingtype, withW); CHECKSTATUS(status); return SUCCESS; } status = temporal_interpolationUV(x, y, z, time, U, V, xi, yi, zi, ti, valueU, valueV, interp_method, gridindexingtype); CHECKSTATUS(status); diff --git a/parcels/kernel.py b/parcels/kernel.py index 1c9758440a..867251a68c 100644 --- a/parcels/kernel.py +++ b/parcels/kernel.py @@ -30,7 +30,7 @@ ) from parcels.compilation.codegenerator import KernelGenerator, LoopGenerator from parcels.field import Field, NestedField, VectorField -from parcels.grid import GridCode +from parcels.grid import GridType from parcels.tools.global_statics import get_cache_dir from parcels.tools.loggers import logger from parcels.tools.statuscodes import ( @@ -323,7 +323,7 @@ def check_fieldsets_in_kernels(self, pyfunc): raise NotImplementedError('Analytical Advection only works in Scipy mode') if self._fieldset.U.interp_method != 'cgrid_velocity': raise NotImplementedError('Analytical Advection only works with C-grids') - if self._fieldset.U.grid.gtype not in [GridCode.CurvilinearZGrid, GridCode.RectilinearZGrid]: + if self._fieldset.U.grid.gtype not in [GridType.CurvilinearZGrid, GridType.RectilinearZGrid]: raise NotImplementedError('Analytical Advection only works with Z-grids in the vertical') elif pyfunc is AdvectionRK45: if not hasattr(self.fieldset, 'RK45_tol'): diff --git a/parcels/particleset.py b/parcels/particleset.py index 0628452b60..09010f6507 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -23,7 +23,7 @@ from parcels.application_kernels.advection import AdvectionRK4 from parcels.compilation.codecompiler import GNUCompiler from parcels.field import NestedField -from parcels.grid import CurvilinearGrid, GridCode +from parcels.grid import CurvilinearGrid, GridType from parcels.interaction.interactionkernel import InteractionKernel from parcels.interaction.neighborsearch import ( BruteFlatNeighborSearch, @@ -569,7 +569,7 @@ def monte_carlo_sample(cls, start_field, size, mode='monte_carlo'): j, i = np.unravel_index(inds, p_interior.shape) grid = start_field.grid lon, lat = ([], []) - if grid.gtype in [GridCode.RectilinearZGrid, GridCode.RectilinearSGrid]: + if grid.gtype in [GridType.RectilinearZGrid, GridType.RectilinearSGrid]: lon = grid.lon[i] + xsi * (grid.lon[i + 1] - grid.lon[i]) lat = grid.lat[j] + eta * (grid.lat[j + 1] - grid.lat[j]) else: diff --git a/tests/test_data/create_testfields.py b/tests/test_data/create_testfields.py index f52c8d773f..bbc259d779 100644 --- a/tests/test_data/create_testfields.py +++ b/tests/test_data/create_testfields.py @@ -2,7 +2,7 @@ import numpy as np -from parcels import FieldSet, GridCode +from parcels import FieldSet, GridType try: from pympler import asizeof @@ -87,12 +87,12 @@ def write_simple_2Dt(field, filename, varname=None): varname = field.name # Create DataArray objects for file I/O - if field.grid.gtype == GridCode.RectilinearZGrid: + if field.grid.gtype == GridType.RectilinearZGrid: nav_lon = xr.DataArray(field.grid.lon + np.zeros((field.grid.ydim, field.grid.xdim), dtype=np.float32), coords=[('y', field.grid.lat), ('x', field.grid.lon)]) nav_lat = xr.DataArray(field.grid.lat.reshape(field.grid.ydim, 1) + np.zeros(field.grid.xdim, dtype=np.float32), coords=[('y', field.grid.lat), ('x', field.grid.lon)]) - elif field.grid.gtype == GridCode.CurvilinearZGrid: + elif field.grid.gtype == GridType.CurvilinearZGrid: nav_lon = xr.DataArray(field.grid.lon, coords=[('y', range(field.grid.ydim)), ('x', range(field.grid.xdim))]) nav_lat = xr.DataArray(field.grid.lat, coords=[('y', range(field.grid.ydim)), ('x', range(field.grid.xdim))]) else: From 38df5388d6fcd2a5da3e5e5bd4526f4738851703 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 26 Jul 2024 16:26:55 +0200 Subject: [PATCH 09/11] Update parcels/grid.py Great idea! Co-authored-by: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> --- parcels/grid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/grid.py b/parcels/grid.py index 5693a09264..8c949e9ab1 100644 --- a/parcels/grid.py +++ b/parcels/grid.py @@ -7,7 +7,7 @@ from parcels.tools.converters import TimeConverter from parcels.tools.loggers import logger -__all__ = ['GridType', 'RectilinearZGrid', 'RectilinearSGrid', 'CurvilinearZGrid', 'CurvilinearSGrid', 'CGrid', 'Grid'] +__all__ = ['GridType', 'GridCode', 'RectilinearZGrid', 'RectilinearSGrid', 'CurvilinearZGrid', 'CurvilinearSGrid', 'CGrid', 'Grid'] class GridType(IntEnum): From 1937f7951a8af7852885bd57dce34a65c34e03b4 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 26 Jul 2024 16:28:10 +0200 Subject: [PATCH 10/11] Adding backward compatibility of GridCode (for now) --- parcels/grid.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/parcels/grid.py b/parcels/grid.py index 8c949e9ab1..d983313a58 100644 --- a/parcels/grid.py +++ b/parcels/grid.py @@ -17,6 +17,11 @@ class GridType(IntEnum): CurvilinearSGrid = 3 +# GridCode has been renamed to GridType for consistency. +# TODO: Remove alias in Parcels v4 +GridCode = GridType + + class CGrid(Structure): _fields_ = [('gtype', c_int), ('grid', c_void_p)] From c66c28551e05e12704e9fae02de1e564490cb57e Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 29 Jul 2024 08:18:05 +0200 Subject: [PATCH 11/11] 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 22e05f8dc9..cf19dd10cc 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)