diff --git a/parcels/field.py b/parcels/field.py index ece0f9c53..7af24caf1 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 cc9bcaf2d..d983313a5 100644 --- a/parcels/grid.py +++ b/parcels/grid.py @@ -7,16 +7,21 @@ from parcels.tools.converters import TimeConverter from parcels.tools.loggers import logger -__all__ = ['GridCode', 'RectilinearZGrid', 'RectilinearSGrid', 'CurvilinearZGrid', 'CurvilinearSGrid', 'CGrid', 'Grid'] +__all__ = ['GridType', 'GridCode', 'RectilinearZGrid', 'RectilinearSGrid', 'CurvilinearZGrid', 'CurvilinearSGrid', 'CGrid', 'Grid'] -class GridCode(IntEnum): +class GridType(IntEnum): RectilinearZGrid = 0 RectilinearSGrid = 1 CurvilinearZGrid = 2 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)] @@ -326,7 +331,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 +376,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 +486,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 +530,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 39b454cfa..bb42b2db6 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 61775a819..e4d8431e3 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 1c9758440..867251a68 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 0628452b6..09010f650 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 f52c8d773..bbc259d77 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: