diff --git a/parcels/field.py b/parcels/field.py index 210d32a0d..a1d38b199 100644 --- a/parcels/field.py +++ b/parcels/field.py @@ -1371,7 +1371,7 @@ def eval(self, time, z, y, x, particle=None, applyConversion=True): def ccode_eval(self, var, t, z, y, x): self._check_velocitysampling() - ccode_str = f"temporal_interpolation({x}, {y}, {z}, {t}, {self.ccode_name}, &particles->xi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->ti[pnum*ngrid], &{var}, {self.interp_method.upper()})" + ccode_str = f"temporal_interpolation({x}, {y}, {z}, {t}, {self.ccode_name}, &particles->xi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->ti[pnum*ngrid], &{var})" return ccode_str def ccode_convert(self, _, z, y, x): @@ -1456,6 +1456,7 @@ class CField(Structure): ("allow_time_extrapolation", c_int), ("time_periodic", c_int), ("gridindexingtype", c_int), + ("interp_method", c_int), ("data_chunks", POINTER(POINTER(POINTER(c_float)))), ("grid", POINTER(CGrid)), ] @@ -1476,6 +1477,18 @@ class CField(Structure): self.c_data_chunks[i] = None gridindexingtype_mapping = {"nemo": 0, "mitgcm": 1, "mom5": 2, "pop": 3} + interp_method_mapping = { + "LINEAR": 0, + "NEAREST": 1, + "CGRID_VELOCITY": 2, + "CGRID_TRACER": 3, + "BGRID_VELOCITY": 4, + "BGRID_W_VELOCITY": 5, + "BGRID_TRACER": 6, + "LINEAR_INVDIST_LAND_TRACER": 7, + "PARTIALSLIP": 8, + "FREESLIP": 9, + } cstruct = CField( self.grid.xdim, @@ -1486,6 +1499,7 @@ class CField(Structure): allow_time_extrapolation, time_periodic, gridindexingtype_mapping[self.gridindexingtype], + interp_method_mapping[self.interp_method.upper()], (POINTER(POINTER(c_float)) * len(self.c_data_chunks))(*self.c_data_chunks), pointer(self.grid.ctypes_struct), ) @@ -2223,13 +2237,13 @@ def ccode_eval(self, varU, varV, varW, U, V, W, t, z, y, x): ccode_str = ( f"temporal_interpolationUVW({x}, {y}, {z}, {t}, {U.ccode_name}, {V.ccode_name}, {W.ccode_name}, " + "&particles->xi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->ti[pnum*ngrid]," - + f"&{varU}, &{varV}, &{varW}, {U.interp_method.upper()})" + + f"&{varU}, &{varV}, &{varW})" ) else: ccode_str = ( f"temporal_interpolationUV({x}, {y}, {z}, {t}, {U.ccode_name}, {V.ccode_name}, " + "&particles->xi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->ti[pnum*ngrid]," - + f" &{varU}, &{varV}, {U.interp_method.upper()})" + + f" &{varU}, &{varV})" ) return ccode_str diff --git a/parcels/include/index_search.h b/parcels/include/index_search.h index bb42b2db6..ceaada507 100644 --- a/parcels/include/index_search.h +++ b/parcels/include/index_search.h @@ -22,7 +22,7 @@ typedef float type_coord; typedef enum { LINEAR=0, NEAREST=1, CGRID_VELOCITY=2, CGRID_TRACER=3, BGRID_VELOCITY=4, BGRID_W_VELOCITY=5, BGRID_TRACER=6, LINEAR_INVDIST_LAND_TRACER=7, PARTIALSLIP=8, FREESLIP=9 - } InterpCode; + } InterpMethod; typedef enum { diff --git a/parcels/include/parcels.h b/parcels/include/parcels.h index 3ec15943c..d8e9bb020 100644 --- a/parcels/include/parcels.h +++ b/parcels/include/parcels.h @@ -18,7 +18,7 @@ extern "C" { typedef struct { - int xdim, ydim, zdim, tdim, igrid, allow_time_extrapolation, time_periodic, gridindexingtype; + int xdim, ydim, zdim, tdim, igrid, allow_time_extrapolation, time_periodic, gridindexingtype, interp_method; float ****data_chunks; CGrid *grid; } CField; @@ -436,11 +436,12 @@ 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, GridType gtype, int *xi, int *yi, int *zi, int *ti, - float *value, int interp_method) + float *value) { StatusCode status; CStructuredGrid *grid = f->grid->grid; GridIndexingType gridindexingtype = f->gridindexingtype; + InterpMethod interp_method = f->interp_method; int igrid = f->igrid; /* Find time index for temporal interpolation */ @@ -1086,11 +1087,12 @@ 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, GridType gtype, int *xi, int *yi, int *zi, int *ti, - float *u, float *v, float *w, int interp_method, int withW) + float *u, float *v, float *w, int withW) { StatusCode status; CStructuredGrid *grid = U->grid->grid; GridIndexingType gridindexingtype = U->gridindexingtype; + InterpMethod interp_method = U->interp_method; int igrid = U->igrid; /* Find time index for temporal interpolation */ @@ -1180,13 +1182,13 @@ static inline StatusCode temporal_interpolation_slip(type_coord x, type_coord y, static inline StatusCode temporal_interpolation(type_coord x, type_coord y, type_coord z, double time, CField *f, int *xi, int *yi, int *zi, int *ti, - float *value, int interp_method) + float *value) { CGrid *_grid = f->grid; GridType gtype = _grid->gtype; 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); + return temporal_interpolation_structured_grid(x, y, z, time, f, gtype, xi, yi, zi, ti, value); else{ printf("Only RECTILINEAR_Z_GRID, RECTILINEAR_S_GRID, CURVILINEAR_Z_GRID and CURVILINEAR_S_GRID grids are currently implemented\n"); return ERROR; @@ -1196,9 +1198,10 @@ static inline StatusCode temporal_interpolation(type_coord x, type_coord y, type static inline StatusCode temporal_interpolationUV(type_coord x, type_coord y, type_coord z, double time, CField *U, CField *V, int *xi, int *yi, int *zi, int *ti, - float *valueU, float *valueV, int interp_method) + float *valueU, float *valueV) { StatusCode status; + InterpMethod interp_method = U->interp_method; if (interp_method == CGRID_VELOCITY){ CGrid *_grid = U->grid; GridType gtype = _grid->gtype; @@ -1209,11 +1212,11 @@ static inline StatusCode temporal_interpolationUV(type_coord x, type_coord y, ty CField *W = U; GridType gtype = _grid->gtype; int withW = 0; - status = temporal_interpolation_slip(x, y, z, time, U, V, W, gtype, xi, yi, zi, ti, valueU, valueV, 0, interp_method, withW); CHECKSTATUS(status); + status = temporal_interpolation_slip(x, y, z, time, U, V, W, gtype, xi, yi, zi, ti, valueU, valueV, 0, withW); CHECKSTATUS(status); return SUCCESS; } else { - status = temporal_interpolation(x, y, z, time, U, xi, yi, zi, ti, valueU, interp_method); CHECKSTATUS(status); - status = temporal_interpolation(x, y, z, time, V, xi, yi, zi, ti, valueV, interp_method); CHECKSTATUS(status); + status = temporal_interpolation(x, y, z, time, U, xi, yi, zi, ti, valueU); CHECKSTATUS(status); + status = temporal_interpolation(x, y, z, time, V, xi, yi, zi, ti, valueV); CHECKSTATUS(status); return SUCCESS; } } @@ -1221,9 +1224,10 @@ static inline StatusCode temporal_interpolationUV(type_coord x, type_coord y, ty static inline StatusCode temporal_interpolationUVW(type_coord x, type_coord y, type_coord z, double time, CField *U, CField *V, CField *W, int *xi, int *yi, int *zi, int *ti, - float *valueU, float *valueV, float *valueW, int interp_method) + float *valueU, float *valueV, float *valueW) { StatusCode status; + InterpMethod interp_method = U->interp_method; if (interp_method == CGRID_VELOCITY){ CGrid *_grid = U->grid; GridType gtype = _grid->gtype; @@ -1235,13 +1239,13 @@ static inline StatusCode temporal_interpolationUVW(type_coord x, type_coord y, t CGrid *_grid = U->grid; GridType gtype = _grid->gtype; int withW = 1; - status = temporal_interpolation_slip(x, y, z, time, U, V, W, gtype, xi, yi, zi, ti, valueU, valueV, valueW, interp_method, withW); CHECKSTATUS(status); + status = temporal_interpolation_slip(x, y, z, time, U, V, W, gtype, xi, yi, zi, ti, valueU, valueV, valueW, withW); CHECKSTATUS(status); return SUCCESS; } - status = temporal_interpolationUV(x, y, z, time, U, V, xi, yi, zi, ti, valueU, valueV, interp_method); CHECKSTATUS(status); + status = temporal_interpolationUV(x, y, z, time, U, V, xi, yi, zi, ti, valueU, valueV); CHECKSTATUS(status); if (interp_method == BGRID_VELOCITY) interp_method = BGRID_W_VELOCITY; - status = temporal_interpolation(x, y, z, time, W, xi, yi, zi, ti, valueW, interp_method); CHECKSTATUS(status); + status = temporal_interpolation(x, y, z, time, W, xi, yi, zi, ti, valueW); CHECKSTATUS(status); return SUCCESS; }