Skip to content

Commit 29a9fbb

Browse files
erikvansebilleerikvansebille
committed
Merge branch 'master' into warn_particle_times_outside_fieldset_time_bounds
2 parents 8b6f0b0 + 70e26eb commit 29a9fbb

File tree

15 files changed

+470
-469
lines changed

15 files changed

+470
-469
lines changed

‎parcels/application_kernels/advection.py‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -184,8 +184,8 @@ def AdvectionAnalytical(particle, fieldset, time):
184184
time_i = np.linspace(0, fieldset.U.grid.time[ti + 1] - fieldset.U.grid.time[ti], I_s)
185185
ds_t = min(ds_t, time_i[np.where(time - fieldset.U.grid.time[ti] < time_i)[0][0]])
186186

187-
xsi, eta, zeta, xi, yi, zi = fieldset.U._search_indices(
188-
particle.lon, particle.lat, particle.depth, particle=particle
187+
zeta, eta, xsi, zi, yi, xi = fieldset.U._search_indices(
188+
-1, particle.depth, particle.lat, particle.lon, particle=particle
189189
)
190190
if withW:
191191
if abs(xsi - 1) < tol:
@@ -232,14 +232,14 @@ def AdvectionAnalytical(particle, fieldset, time):
232232
else:
233233
dz = 1.0
234234

235-
c1 = fieldset.UV.dist(px[0], px[1], py[0], py[1], grid.mesh, np.dot(i_u.phi2D_lin(xsi, 0.0), py))
236-
c2 = fieldset.UV.dist(px[1], px[2], py[1], py[2], grid.mesh, np.dot(i_u.phi2D_lin(1.0, eta), py))
237-
c3 = fieldset.UV.dist(px[2], px[3], py[2], py[3], grid.mesh, np.dot(i_u.phi2D_lin(xsi, 1.0), py))
238-
c4 = fieldset.UV.dist(px[3], px[0], py[3], py[0], grid.mesh, np.dot(i_u.phi2D_lin(0.0, eta), py))
235+
c1 = fieldset.UV.dist(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py))
236+
c2 = fieldset.UV.dist(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
237+
c3 = fieldset.UV.dist(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
238+
c4 = fieldset.UV.dist(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
239239
rad = np.pi / 180.0
240240
deg2m = 1852 * 60.0
241241
meshJac = (deg2m * deg2m * math.cos(rad * particle.lat)) if grid.mesh == "spherical" else 1
242-
dxdy = fieldset.UV.jacobian(xsi, eta, px, py) * meshJac
242+
dxdy = fieldset.UV.jacobian(py, px, eta, xsi) * meshJac
243243

244244
if withW:
245245
U0 = direction * fieldset.U.data[ti, zi + 1, yi + 1, xi] * c4 * dz

‎parcels/compilation/codegenerator.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -834,7 +834,7 @@ def visit_FieldEvalNode(self, node):
834834
statements_croco = [
835835
c.Statement(f"float cs_w[] = {*Cs_w, }".replace("(", "{").replace(")", "}")),
836836
c.Statement(
837-
f"{node.var} = croco_from_z_to_sigma(U, H, Zeta, {args[3]}, {args[2]}, {args[1]}, time, &particles->xi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->ti[pnum*ngrid], hc, &cs_w)"
837+
f"{node.var} = croco_from_z_to_sigma(time, {args[1]}, {args[2]}, {args[3]}, U, H, Zeta, &particles->ti[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->xi[pnum*ngrid], hc, &cs_w)"
838838
),
839839
]
840840
args = (args[0], node.var, args[2], args[3])
@@ -863,7 +863,7 @@ def visit_VectorFieldEvalNode(self, node):
863863
statements_croco = [
864864
c.Statement(f"float cs_w[] = {*Cs_w, }".replace("(", "{").replace(")", "}")),
865865
c.Statement(
866-
f"{node.var4} = croco_from_z_to_sigma(U, H, Zeta, {args[3]}, {args[2]}, {args[1]}, time, &particles->xi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->ti[pnum*ngrid], hc, &cs_w)"
866+
f"{node.var4} = croco_from_z_to_sigma(time, {args[1]}, {args[2]}, {args[3]}, U, H, Zeta, &particles->ti[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->xi[pnum*ngrid], hc, &cs_w)"
867867
),
868868
]
869869
args = (args[0], node.var4, args[2], args[3])

‎parcels/field.py‎

Lines changed: 97 additions & 91 deletions
Large diffs are not rendered by default.

‎parcels/fieldfilebuffer.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -92,7 +92,7 @@ def parse_name(self, name):
9292
return name
9393

9494
@property
95-
def lonlat(self):
95+
def latlon(self):
9696
lon = self.dataset[self.dimensions["lon"]]
9797
lat = self.dataset[self.dimensions["lat"]]
9898
if self.nolonlatindices and self.gridindexingtype not in ["croco"]:
@@ -141,7 +141,7 @@ def lonlat(self):
141141
if rectilinear:
142142
lon_subset = lon_subset[0, :]
143143
lat_subset = lat_subset[:, 0]
144-
return lon_subset, lat_subset
144+
return lat_subset, lon_subset
145145

146146
@property
147147
def depth(self):

‎parcels/include/index_search.h‎

Lines changed: 38 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -127,9 +127,10 @@ static inline StatusCode search_indices_vertical_z(type_coord z, int zdim, float
127127
return SUCCESS;
128128
}
129129

130-
static inline StatusCode search_indices_vertical_s(type_coord z, int xdim, int ydim, int zdim, float *zvals,
131-
int xi, int yi, int *zi, double xsi, double eta, double *zeta,
132-
int z4d, int ti, int tdim, double time, double t0, double t1, int interp_method)
130+
static inline StatusCode search_indices_vertical_s(double time, type_coord z,
131+
int tdim, int zdim, int ydim, int xdim, float *zvals,
132+
int ti, int *zi, int yi, int xi, double *zeta, double eta, double xsi,
133+
int z4d, double t0, double t1, int interp_method)
133134
{
134135
if (interp_method == BGRID_VELOCITY || interp_method == BGRID_W_VELOCITY || interp_method == BGRID_TRACER){
135136
xsi = 1;
@@ -184,7 +185,7 @@ static inline StatusCode search_indices_vertical_s(type_coord z, int xdim, int y
184185
return SUCCESS;
185186
}
186187

187-
static inline void reconnect_bnd_indices(int *xi, int *yi, int xdim, int ydim, int onlyX, int sphere_mesh)
188+
static inline void reconnect_bnd_indices(int *yi, int *xi, int ydim, int xdim, int onlyX, int sphere_mesh)
188189
{
189190
if (*xi < 0){
190191
if (sphere_mesh)
@@ -211,10 +212,11 @@ static inline void reconnect_bnd_indices(int *xi, int *yi, int xdim, int ydim, i
211212
}
212213

213214

214-
static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridType gtype,
215-
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
216-
int ti, double time, double t0, double t1, int interp_method,
217-
int gridindexingtype)
215+
static inline StatusCode search_indices_rectilinear(double time, type_coord z, type_coord y, type_coord x,
216+
CStructuredGrid *grid, GridType gtype,
217+
int ti, int *zi, int *yi, int *xi, double *zeta, double *eta, double *xsi,
218+
double t0, double t1, int interp_method,
219+
int gridindexingtype)
218220
{
219221
int xdim = grid->xdim;
220222
int ydim = grid->ydim;
@@ -261,7 +263,7 @@ static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y,
261263
++(*xi);
262264
else if (xvalsi > x)
263265
--(*xi);
264-
reconnect_bnd_indices(xi, yi, xdim, ydim, 1, 1);
266+
reconnect_bnd_indices(yi, xi, ydim, xdim, 1, 1);
265267
xvalsi = xvals[*xi];
266268
if (xvalsi < x - 225) xvalsi += 360;
267269
if (xvalsi > x + 225) xvalsi -= 360;
@@ -294,9 +296,9 @@ static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y,
294296
status = search_indices_vertical_z(z, zdim, zvals, zi, zeta, gridindexingtype);
295297
break;
296298
case RECTILINEAR_S_GRID:
297-
status = search_indices_vertical_s(z, xdim, ydim, zdim, zvals,
298-
*xi, *yi, zi, *xsi, *eta, zeta,
299-
z4d, ti, tdim, time, t0, t1, interp_method);
299+
status = search_indices_vertical_s(time, z, tdim, zdim, ydim, xdim, zvals,
300+
ti, zi, *yi, *xi, zeta, *eta, *xsi,
301+
z4d, t0, t1, interp_method);
300302
break;
301303
default:
302304
status = ERRORINTERPOLATION;
@@ -321,10 +323,12 @@ static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y,
321323
}
322324

323325

324-
static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridType gtype,
325-
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
326-
int ti, double time, double t0, double t1, int interp_method,
327-
int gridindexingtype)
326+
static inline StatusCode search_indices_curvilinear(double time, type_coord z, type_coord y, type_coord x,
327+
CStructuredGrid *grid, GridType gtype,
328+
int ti, int *zi, int *yi, int *xi,
329+
double *zeta, double *eta, double *xsi,
330+
double t0, double t1, int interp_method,
331+
int gridindexingtype)
328332
{
329333
int xi_old = *xi;
330334
int yi_old = *yi;
@@ -407,23 +411,23 @@ static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y,
407411
(*yi)--;
408412
if (*eta > 1+tol)
409413
(*yi)++;
410-
reconnect_bnd_indices(xi, yi, xdim, ydim, 0, sphere_mesh);
414+
reconnect_bnd_indices(yi, xi, ydim, xdim, 0, sphere_mesh);
411415
it++;
412416
if ( it > maxIterSearch){
413-
printf("Correct cell not found for (%f, %f) after %d iterations\n", x, y, maxIterSearch);
417+
printf("Correct cell not found for (lat, lon) = (%f, %f) after %d iterations\n", y, x, maxIterSearch);
414418
printf("Debug info: old particle indices: (yi, xi) %d %d\n", yi_old, xi_old);
415419
printf(" new particle indices: (yi, xi) %d %d\n", *yi, *xi);
416420
printf(" Mesh 2d shape: %d %d\n", ydim, xdim);
417-
printf(" Relative particle position: (xsi, eta) %1.16e %1.16e\n", *xsi, *eta);
421+
printf(" Relative particle position: (eta, xsi) %1.16e %1.16e\n", *eta, *xsi);
418422
return ERROROUTOFBOUNDS;
419423
}
420424
}
421425
if ( (*xsi != *xsi) || (*eta != *eta) ){ // check if nan
422-
printf("Correct cell not found for (%f, %f))\n", x, y);
426+
printf("Correct cell not found for (lat, lon) = (%f, %f))\n", y, x);
423427
printf("Debug info: old particle indices: (yi, xi) %d %d\n", yi_old, xi_old);
424428
printf(" new particle indices: (yi, xi) %d %d\n", *yi, *xi);
425429
printf(" Mesh 2d shape: %d %d\n", ydim, xdim);
426-
printf(" Relative particle position: (xsi, eta) %1.16e %1.16e\n", *xsi, *eta);
430+
printf(" Relative particle position: (eta, xsi) %1.16e %1.16e\n", *eta, *xsi);
427431
return ERROROUTOFBOUNDS;
428432
}
429433
if (*xsi < 0) *xsi = 0;
@@ -438,9 +442,9 @@ static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y,
438442
status = search_indices_vertical_z(z, zdim, zvals, zi, zeta, gridindexingtype);
439443
break;
440444
case CURVILINEAR_S_GRID:
441-
status = search_indices_vertical_s(z, xdim, ydim, zdim, zvals,
442-
*xi, *yi, zi, *xsi, *eta, zeta,
443-
z4d, ti, tdim, time, t0, t1, interp_method);
445+
status = search_indices_vertical_s(time, z, tdim, ydim, xdim, zdim, zvals,
446+
ti, zi, *yi, *xi, zeta, *eta, *xsi,
447+
z4d, t0, t1, interp_method);
444448
break;
445449
default:
446450
status = ERRORINTERPOLATION;
@@ -460,21 +464,23 @@ static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y,
460464
/* Local linear search to update grid index
461465
* params ti, sizeT, time. t0, t1 are only used for 4D S grids
462466
* */
463-
static inline StatusCode search_indices(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid,
464-
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
465-
GridType gtype, int ti, double time, double t0, double t1, int interp_method,
466-
int gridindexingtype)
467+
static inline StatusCode search_indices(double time, type_coord z, type_coord y, type_coord x,
468+
CStructuredGrid *grid,
469+
int ti, int *zi, int *yi, int *xi,
470+
double *zeta, double *eta, double *xsi,
471+
GridType gtype, double t0, double t1, int interp_method,
472+
int gridindexingtype)
467473
{
468474
switch(gtype){
469475
case RECTILINEAR_Z_GRID:
470476
case RECTILINEAR_S_GRID:
471-
return search_indices_rectilinear(x, y, z, grid, gtype, xi, yi, zi, xsi, eta, zeta,
472-
ti, time, t0, t1, interp_method, gridindexingtype);
477+
return search_indices_rectilinear(time, z, y, x, grid, gtype, ti, zi, yi, xi, zeta, eta, xsi,
478+
t0, t1, interp_method, gridindexingtype);
473479
break;
474480
case CURVILINEAR_Z_GRID:
475481
case CURVILINEAR_S_GRID:
476-
return search_indices_curvilinear(x, y, z, grid, gtype, xi, yi, zi, xsi, eta, zeta,
477-
ti, time, t0, t1, interp_method, gridindexingtype);
482+
return search_indices_curvilinear(time, z, y, x, grid, gtype, ti, zi, yi, xi, zeta, eta, xsi,
483+
t0, t1, interp_method, gridindexingtype);
478484
break;
479485
default:
480486
printf("Only RECTILINEAR_Z_GRID, RECTILINEAR_S_GRID, CURVILINEAR_Z_GRID and CURVILINEAR_S_GRID grids are currently implemented\n");

‎parcels/include/interpolation_utils.h‎

Lines changed: 42 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ typedef enum
88
} Orientation;
99

1010

11-
static inline void phi2D_lin(double xsi, double eta, double *phi)
11+
static inline void phi2D_lin(double eta, double xsi, double *phi)
1212
{
1313
phi[0] = (1-xsi) * (1-eta);
1414
phi[1] = xsi * (1-eta);
@@ -25,40 +25,40 @@ static inline void phi1D_quad(double xsi, double *phi)
2525
}
2626

2727

28-
static inline void dphidxsi3D_lin(double xsi, double eta, double zet, double *dphidxsi, double *dphideta, double *dphidzet)
28+
static inline void dphidxsi3D_lin(double zeta, double eta, double xsi, double *dphidzeta, double *dphideta, double *dphidxsi)
2929
{
30-
dphidxsi[0] = - (1-eta) * (1-zet);
31-
dphidxsi[1] = (1-eta) * (1-zet);
32-
dphidxsi[2] = ( eta) * (1-zet);
33-
dphidxsi[3] = - ( eta) * (1-zet);
34-
dphidxsi[4] = - (1-eta) * ( zet);
35-
dphidxsi[5] = (1-eta) * ( zet);
36-
dphidxsi[6] = ( eta) * ( zet);
37-
dphidxsi[7] = - ( eta) * ( zet);
38-
39-
dphideta[0] = - (1-xsi) * (1-zet);
40-
dphideta[1] = - ( xsi) * (1-zet);
41-
dphideta[2] = ( xsi) * (1-zet);
42-
dphideta[3] = (1-xsi) * (1-zet);
43-
dphideta[4] = - (1-xsi) * ( zet);
44-
dphideta[5] = - ( xsi) * ( zet);
45-
dphideta[6] = ( xsi) * ( zet);
46-
dphideta[7] = (1-xsi) * ( zet);
47-
48-
dphidzet[0] = - (1-xsi) * (1-eta);
49-
dphidzet[1] = - ( xsi) * (1-eta);
50-
dphidzet[2] = - ( xsi) * ( eta);
51-
dphidzet[3] = - (1-xsi) * ( eta);
52-
dphidzet[4] = (1-xsi) * (1-eta);
53-
dphidzet[5] = ( xsi) * (1-eta);
54-
dphidzet[6] = ( xsi) * ( eta);
55-
dphidzet[7] = (1-xsi) * ( eta);
30+
dphidxsi[0] = - (1-eta) * (1-zeta);
31+
dphidxsi[1] = (1-eta) * (1-zeta);
32+
dphidxsi[2] = ( eta) * (1-zeta);
33+
dphidxsi[3] = - ( eta) * (1-zeta);
34+
dphidxsi[4] = - (1-eta) * ( zeta);
35+
dphidxsi[5] = (1-eta) * ( zeta);
36+
dphidxsi[6] = ( eta) * ( zeta);
37+
dphidxsi[7] = - ( eta) * ( zeta);
38+
39+
dphideta[0] = - (1-xsi) * (1-zeta);
40+
dphideta[1] = - ( xsi) * (1-zeta);
41+
dphideta[2] = ( xsi) * (1-zeta);
42+
dphideta[3] = (1-xsi) * (1-zeta);
43+
dphideta[4] = - (1-xsi) * ( zeta);
44+
dphideta[5] = - ( xsi) * ( zeta);
45+
dphideta[6] = ( xsi) * ( zeta);
46+
dphideta[7] = (1-xsi) * ( zeta);
47+
48+
dphidzeta[0] = - (1-xsi) * (1-eta);
49+
dphidzeta[1] = - ( xsi) * (1-eta);
50+
dphidzeta[2] = - ( xsi) * ( eta);
51+
dphidzeta[3] = - (1-xsi) * ( eta);
52+
dphidzeta[4] = (1-xsi) * (1-eta);
53+
dphidzeta[5] = ( xsi) * (1-eta);
54+
dphidzeta[6] = ( xsi) * ( eta);
55+
dphidzeta[7] = (1-xsi) * ( eta);
5656
}
5757

58-
static inline void dxdxsi3D_lin(double *px, double *py, double *pz, double xsi, double eta, double zet, double *jacM, int sphere_mesh)
58+
static inline void dxdxsi3D_lin(double *pz, double *py, double *px, double zeta, double eta, double xsi, double *jacM, int sphere_mesh)
5959
{
60-
double dphidxsi[8], dphideta[8], dphidzet[8];
61-
dphidxsi3D_lin(xsi, eta, zet, dphidxsi, dphideta, dphidzet);
60+
double dphidxsi[8], dphideta[8], dphidzeta[8];
61+
dphidxsi3D_lin(zeta, eta, xsi, dphidzeta, dphideta, dphidxsi);
6262

6363
int i;
6464
for(i=0; i<9; ++i)
@@ -76,20 +76,22 @@ static inline void dxdxsi3D_lin(double *px, double *py, double *pz, double xsi,
7676
for(i=0; i<8; ++i){
7777
jacM[3*0+0] += px[i] * dphidxsi[i] * jac_lon; // dxdxsi
7878
jacM[3*0+1] += px[i] * dphideta[i] * jac_lon; // dxdeta
79-
jacM[3*0+2] += px[i] * dphidzet[i] * jac_lon; // dxdzet
79+
jacM[3*0+2] += px[i] * dphidzeta[i] * jac_lon; // dxdzeta
8080
jacM[3*1+0] += py[i] * dphidxsi[i] * jac_lat; // dydxsi
8181
jacM[3*1+1] += py[i] * dphideta[i] * jac_lat; // dydeta
82-
jacM[3*1+2] += py[i] * dphidzet[i] * jac_lat; // dydzet
82+
jacM[3*1+2] += py[i] * dphidzeta[i] * jac_lat; // dydzeta
8383
jacM[3*2+0] += pz[i] * dphidxsi[i]; // dzdxsi
8484
jacM[3*2+1] += pz[i] * dphideta[i]; // dzdeta
85-
jacM[3*2+2] += pz[i] * dphidzet[i]; // dzdzet
85+
jacM[3*2+2] += pz[i] * dphidzeta[i]; // dzdzeta
8686
}
8787
}
8888

89-
static inline double jacobian3D_lin_face(double *px, double *py, double *pz, double xsi, double eta, double zet, Orientation orientation, int sphere_mesh)
89+
static inline double jacobian3D_lin_face(double *pz, double *py, double *px,
90+
double zeta, double eta, double xsi,
91+
Orientation orientation, int sphere_mesh)
9092
{
9193
double jacM[9];
92-
dxdxsi3D_lin(px, py, pz, xsi, eta, zet, jacM, sphere_mesh);
94+
dxdxsi3D_lin(pz, py, px, zeta, eta, xsi, jacM, sphere_mesh);
9395

9496
double j[3];
9597

@@ -112,10 +114,12 @@ static inline double jacobian3D_lin_face(double *px, double *py, double *pz, dou
112114
return sqrt(j[0]*j[0]+j[1]*j[1]+j[2]*j[2]);
113115
}
114116

115-
static inline double jacobian3D_lin(double *px, double *py, double *pz, double xsi, double eta, double zet, int sphere_mesh)
117+
static inline double jacobian3D_lin(double *pz, double *py, double *px,
118+
double zeta, double eta, double xsi,
119+
int sphere_mesh)
116120
{
117121
double jacM[9];
118-
dxdxsi3D_lin(px, py, pz, xsi, eta, zet, jacM, sphere_mesh);
122+
dxdxsi3D_lin(pz, py, px, zeta, eta, xsi, jacM, sphere_mesh);
119123

120124
double jac = jacM[3*0+0] * (jacM[3*1+1]*jacM[3*2+2] - jacM[3*2+1]*jacM[3*1+2])
121125
- jacM[3*0+1] * (jacM[3*1+0]*jacM[3*2+2] - jacM[3*2+0]*jacM[3*1+2])

0 commit comments

Comments
 (0)