Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactoring of CField struct to include gridindexingtype and interp_method #1673

Closed
wants to merge 2 commits into from
Closed

Refactoring of CField struct to include gridindexingtype and interp_method #1673

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1dccf0d
Using gridindexingtype from CField struct in C code
erikvansebille Aug 29, 2024
0e065a5
Also making interp_method a CField attribute
erikvansebille Aug 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 17 additions & 3 deletions parcels/field.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()}, {self.gridindexingtype.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):
Expand Down Expand Up @@ -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)),
]
Expand All @@ -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,
}
Comment on lines 1479 to +1491
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is still a bit ugly, because it mirrors the information from

Parcels/parcels/include/index_search.h

Lines 22 to 30 in 997e3f8

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;
typedef enum
{
NEMO = 0, MITGCM = 1, MOM5 = 2, POP = 3
} GridIndexingType;

It would be cleaner to either

  1. Use the string here in the CField, instead of an integer. I tried using c_char_p, but that didn't work (yet)
  2. Somehow import the typedefs from index_search.h (above) here
    In any case, I think we want to avoid happing two places where these mappings between strings and integers occur...


cstruct = CField(
self.grid.xdim,
Expand All @@ -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),
)
Expand Down Expand Up @@ -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()}, {U.gridindexingtype.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()}, {U.gridindexingtype.upper()})"
+ f" &{varU}, &{varV})"
)
return ccode_str

Expand Down
2 changes: 1 addition & 1 deletion parcels/include/index_search.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand Down
51 changes: 30 additions & 21 deletions parcels/include/parcels.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -436,10 +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, int gridindexingtype)
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 */
Expand Down Expand Up @@ -644,10 +646,12 @@ 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,
GridType gtype, int *xi, int *yi, int *zi, int *ti,
float *u, float *v, int gridindexingtype)
float *u, float *v)
{
StatusCode status;
CStructuredGrid *grid = U->grid->grid;
GridIndexingType gridindexingtype = U->gridindexingtype;

int igrid = U->igrid;

/* Find time index for temporal interpolation */
Expand Down Expand Up @@ -726,7 +730,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,
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)
GridType gtype, float dataU[2][2][2], float dataV[2][2][2], float dataW[2][2][2], float *u, float *v, float *w)
{
/* Cast data array into data[lat][lon] as per NEMO convention */
int xdim = grid->xdim;
Expand Down Expand Up @@ -862,10 +866,11 @@ 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,
GridType gtype, int *xi, int *yi, int *zi, int *ti,
float *u, float *v, float *w, int gridindexingtype)
float *u, float *v, float *w)
{
StatusCode status;
CStructuredGrid *grid = U->grid->grid;
GridIndexingType gridindexingtype = U->gridindexingtype;
int igrid = U->igrid;

/* Find time index for temporal interpolation */
Expand All @@ -891,8 +896,8 @@ static inline StatusCode temporal_interpolationUVW_c_grid(type_coord x, type_coo
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, 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);
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); 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); CHECKSTATUS(status);
}
*u = u0 + (u1 - u0) * (float)((time - t0) / (t1 - t0));
*v = v0 + (v1 - v0) * (float)((time - t0) / (t1 - t0));
Expand All @@ -908,7 +913,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, gtype, 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); CHECKSTATUS(status);
}
return SUCCESS;
}
Expand Down Expand Up @@ -1082,10 +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 gridindexingtype, 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 */
Expand Down Expand Up @@ -1175,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, int gridindexingtype)
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, gridindexingtype);
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;
Expand All @@ -1191,52 +1198,54 @@ 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, int gridindexingtype)
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;
status = temporal_interpolationUV_c_grid(x, y, z, time, U, V, gtype, xi, yi, zi, ti, valueU, valueV, gridindexingtype); CHECKSTATUS(status);
status = temporal_interpolationUV_c_grid(x, y, z, time, U, V, gtype, xi, yi, zi, ti, valueU, valueV); CHECKSTATUS(status);
return SUCCESS;
} else if ((interp_method == PARTIALSLIP) || (interp_method == FREESLIP)){
CGrid *_grid = U->grid;
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, gridindexingtype, 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, gridindexingtype); CHECKSTATUS(status);
status = temporal_interpolation(x, y, z, time, V, xi, yi, zi, ti, valueV, interp_method, gridindexingtype); 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, int gridindexingtype)
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;
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);
status = temporal_interpolationUVW_c_grid(x, y, z, time, U, V, W, gtype, xi, yi, zi, ti, valueU, valueV, valueW); CHECKSTATUS(status);
return SUCCESS;
}
} else if ((interp_method == PARTIALSLIP) || (interp_method == FREESLIP)){
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, gridindexingtype, 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, gridindexingtype); 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, gridindexingtype); CHECKSTATUS(status);
status = temporal_interpolation(x, y, z, time, W, xi, yi, zi, ti, valueW); CHECKSTATUS(status);
return SUCCESS;
}

Expand Down
Loading