Also making interp_method a CField attribute

OceanParcels · Aug 29, 2024 · 0e065a5 · 0e065a5
1 parent 1dccf0d
commit 0e065a5
Show file tree

Hide file tree

Showing 3 changed files with 35 additions and 17 deletions.
diff --git 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
@@ -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
@@ -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,21 +1212,22 @@ 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;
   }
 }
 
 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;
 }