Correct cell not found at North Pole with a curvilinear grid

[noemieplanat]

Question

Question

Hi all,
I am trying to run parcels in the Arctic with a curvilinear Z grid. This grid (see photos) has the advantage that is has no "pole".

• I was able to modify the reconnect_bnd_indices function in index_search.h so that the initialization of the particles works (fail vs sucess images below). It's not a "beautiful" way to do it but it takes into account the curvature of the grid along y direction.
• it runs smoothly, except at the pole. I keep having particles that show "correct cell not found"

 Correct cell not found for (lat, lon) = (89.895232, -107.906885) after 1000 iterations
Debug info: old particle indices: (yi, xi) 1243 855
            new particle indices: (yi, xi) 1244 855
            Mesh 2d shape:  1801 1580
            Relative particle position:  (eta, xsi) 1.0321382495380942e+00 5.6039516532941140e-02

I have tried to run the same in python (from e.g. indices 1243, 855 and target lat-lon = (89.895232, -107.906885)) by copying in the search_indices_curvilinear function in python. Surprisingly, it works in python, which make me think the grid that I give it (ds_mesh.nav_lat, ds_mesh.nav_lon) is not exactly the one Parcels uses in that function.

So two questions :

• how is the grid handled, is it re-ordered ? Is anything special happening at the pole ?
• Is it normal that the seach_indices_curvilinear function is failing at the pole ?

Any idea would be helpful, thanks !
Noemie

Fail and sucess initialization of a particle at location 68N, -165W :


Grid in CREG :

Supporting code/error messages

Here is the modified reconnect_bnd_indices function that I use.

static inline void reconnect_bnd_indices(int *yi, int *xi, int ydim, int xdim, int onlyX, int sphere_mesh)
{
  if ((*xi < 0) || (*xi > xdim -2)){
    if (sphere_mesh){
      if (*yi > 1245 - ydim){
        (*yi) = 1245 + (1245 - *yi);
        (*xi) = 0;
      }
      else{
        if (*xi <0){(*xi) = xdim -2;}
        else (*xi) = 0;
      }
  }
    else
      (*xi) = 0;
  }
  if (*xi > xdim-2){
    if (sphere_mesh)
      (*xi) = 0;
    else
      (*xi) = xdim-2;
  }
  if (onlyX == 0){
    if (*yi < 0){
      (*yi) = 0;
    }
    if (*yi > ydim-2){
      (*yi) = ydim-2;
      if (sphere_mesh)
        (*xi) = xdim - (*xi);
    }
  }
}

Parcel's script in python :

from parcels import FieldSet, ParticleSet, Variable, JITParticle, AdvectionRK4_3D 
from operator import attrgetter
import numpy as np
from Config_file_glouton_adv import Config
from datetime import timedelta as delta


max_aging = 1; year = 2001; m0 = 1; cs = 512; chunk_part = 200; chunk_obs = 100; chunk_vert = 25

class ocean_particle(JITParticle):
    #add some variables
    age = Variable('age',dtype=np.float32, initial=0.)
    stuck = Variable('stuck', dtype=np.int32, initial=0.)
    prev_lon = Variable('prev_lon', dtype=np.float32, to_write=False,
                        initial=attrgetter('lon'))  # the previous longitude
    prev_lat = Variable('prev_lat', dtype=np.float32, to_write=False,
                        initial=attrgetter('lat'))  # the previous latitude.
    uvel = Variable('uvel', dtype=np.float32)
    vvel = Variable('vvel', dtype=np.float32)
    wvel = Variable('wvel', dtype=np.float32)

def Modified_AdvectionRK4_3D(particle, fieldset, time):
    # make sure we don't exceed the surface and the boundaries 
    #depmin = 1# 0.50576001 #premier point de deptht
    #maxlat = 89.95
    #if particle.depth <depmin:
    #   particle.delete()
    if particle.lat <66:
       particle.delete()

    lon0 = particle.lon
    if lon0>180:
        lon0 -=360
    lat0 = particle.lat
    depth0 = particle.depth
    (u1, v1, w1) = fieldset.UVW[time, depth0, lat0, lon0]

    lon1 = particle.lon + u1*.5*particle.dt
    lat1 = particle.lat + v1*.5*particle.dt
    dep1 = particle.depth + w1*.5*particle.dt
    #if dep1 <depmin :
    #    particle.delete() #delete if surface
    if lat1<66:
        particle.delete()
    if lon1>180:
        lon1 -=360
    (u2, v2, w2) = fieldset.UVW[time + .5 * particle.dt, dep1, lat1, lon1]

    lon2 = particle.lon + u2*.5*particle.dt
    lat2 = particle.lat + v2*.5*particle.dt
    dep2 = particle.depth + w2*.5*particle.dt
    #if dep2 < depmin:
    #    particle.delete()
    if lat2<66:
        particle.delete()
    if lon2>180:
        lon2-=360
    (u3, v3, w3) = fieldset.UVW[time + .5 * particle.dt, dep2, lat2, lon2]
    lon3 = particle.lon + u3*particle.dt
    lat3 = particle.lat + v3*particle.dt
    dep3 = particle.depth + w3*particle.dt
    #if dep3 < depmin:
    #    particle.delete()
    if lat3<66:
        particle.delete()
    if lon3 > 180:
        lon3 -=360
    (u4, v4, w4) = fieldset.UVW[time + particle.dt, dep3, lat3, lon3]

    particle_dlon = (u1 + 2.*u2 + 2.*u3 + u4) / 6. * particle.dt
    particle_dlat = (v1 + 2.*v2 + 2.*v3 + v4) / 6. * particle.dt
    particle_ddepth = (w1 + 2.*w2 + 2.*w3 + w4) / 6. * particle.dt
   
    if particle.lat+particle_dlat<66:
         particle.delete()

    # if particle.depth <depmin:
    #     particle.delete()
    # if particle.lat>maxlat:
    #     particle.delete()

def Recorvery_CheckOutOfBounds_all(particle, fieldset, time):
    if particle.state == StatusCode.ErrorOutOfBounds: # type: ignore
        particle.delete()

def Recorvery_CheckOutOfBounds_Bering(particle, fieldset, time):
    # remove particles advected southward through BS
    if particle.state == StatusCode.ErrorOutOfBounds and particle.lat<66: # type: ignore
        particle.delete()

def Recorvery_KeepInOcean(particle, fieldset, time):
    if particle.state == StatusCode.ErrorThroughSurface: # type: ignore
        particle_ddepth = 2
        particle.state = StatusCode.Success # type: ignore


mesh_file_h = Config['domain_path']
nmb_days = int((max_aging)*365)
print('nmb_days : ', nmb_days)

listeU = [Config['datapath']+'CREG12.L75-REF12_y'+str(y)+'m{0:02}'.format(m)+'.1d_gridU.nc' for y in range(year, year+int(max_aging+1)+1)  for m in range(1, 13) if (y-year)*365+(m-m0)*30  <= nmb_days]
listeV = [Config['datapath']+'CREG12.L75-REF12_y'+str(y)+'m{0:02}'.format(m)+'.1d_gridV.nc' for y in range(year, year+int(max_aging+1)+1)  for m in range(1, 13) if (y-year)*365+(m-m0)*30 <= nmb_days]
listeW = [Config['datapath']+'CREG12.L75-REF12_y'+str(y)+'m{0:02}'.format(m)+'.1d_gridW.nc' for y in range(year, year+int(max_aging+1)+1)  for m in range(1, 13) if (y-year)*365+(m-m0)*30 <= nmb_days]


# dictionaries
filenames = {'U': {'lon': mesh_file_h, 'lat': mesh_file_h, 'depth':mesh_file_h, 'data': listeU},
             'V': {'lon': mesh_file_h, 'lat': mesh_file_h, 'depth':mesh_file_h, 'data': listeV},
             'W': {'lon': mesh_file_h, 'lat': mesh_file_h, 'depth':mesh_file_h, 'data': listeW},
             }
variables = {'U': 'vozocrtx',
             'V': 'vomecrty', 
             'W': 'vovecrtz',
             }
dimensions = {'U': {'lon': 'glamf', 'lat': 'gphif', 'depth':'nav_lev', 'time':'time_counter'},
              'V': {'lon': 'glamf', 'lat': 'gphif', 'depth':'nav_lev', 'time':'time_counter'},
              'W': {'lon': 'glamf', 'lat': 'gphif','depth':'nav_lev',  'time':'time_counter'}}

css = {"time": ("time_counter", 1), "lat": ("gphif", 300), "lon": ("glamf", cs), "depth":('nav_lev', chunk_vert)}
fset = FieldSet.from_nemo(filenames, variables, dimensions, allow_time_extrapolation =False, transpose=False, chunksize=css)#field_chunksize=False)
fset.add_constant("Max_aging", max_aging)

fset.add_constant("Max_aging", max_aging)
print('Field set - Done')


pset = ParticleSet(fieldset=fset,   # the fields on which the particles are advected
                   pclass=ocean_particle,  # the type of particles (JITParticle or ScipyParticle)
                   lon=[-105.761814],
                   lat=[89.891613], 
                   depth = 30,
                   time=[0], 
                   lonlatdepth_dtype=np.float64)

print('Particle set - Done')
kernels = pset.Kernel(Modified_AdvectionRK4_3D)
kernels = kernels +pset.Kernel(Recorvery_KeepInOcean) +pset.Kernel(Recorvery_CheckOutOfBounds_Bering)+pset.Kernel(Recorvery_CheckOutOfBounds_all)
print('Kernel - Done')
chunks_particles = chunk_part
chunks_obs = chunk_obs
output_file = pset.ParticleFile('/storage/nplanat/toto8', outputdt=delta(hours = 24), chunks=(chunks_particles, chunks_obs))

pset.execute(kernels, runtime=delta(days=3), dt=delta(minutes=1440), output_file=output_file, verbose_progress=True)

[erikvansebille]

Thanks for reaching out, @noemieplanat. Quite strange indeed that the search_indices_curvilinear() function does work, but pset.execute() not. Can you try running the simulation in Scipy-mode, i.e. by changing your ocean_particle class to

class ocean_particle(ScipyParticle)

[erikvansebille]

OK, but does that FieldOutOfBoundError happens immediately? Or only later in the execution? Is there an output file? What does the trajectory look like until then? Is it different from the JITParticle trajectory?

[noemieplanat]

I initialize the particle very close to the pole. It crashes after a couple of iterations - at the same moment when the JITParticle fails to find the correct cell.
Is there an artificial U-V file that I could design to help find the bug ? Poleward, Eastward ? Along x,y directions of the model ?

[erikvansebille]

Hmm, so the bug sits deeper than I hoped...

Our development efforts now are focussing strongly on the new version 4, which is a complete overhaul of the Field class and will include unstructured grid-support. I expect this to be more robust for complex grids like yours. But it does mean it's hard for us to find time to dig into this issue in the v3-codebase. So I'm afraid we can't really help you here now; you can either wait until the first releases of v4, or interpolate your data to a regular grid?

[noemieplanat]

Hi Erik,
Thanks for the reply. Unfortunately it is too much data for me to interpolate it. When do you expect V4 to be out ? Do you have somewhere some info on how the grid is used in the search_indices_curvilinear() function ? I can try to dig a bit into it to see if I can find the bug or a way to overpass it. I've tried printf the grid corners when entering the function, and they do not equal the ones that I have in my mesh grid file, which makes me wonder how the grid is used ... ? Do you know where I can find the info ?
Alternatively I will delete all my particles north of 87deg N.

[erikvansebille]

Hi @noemieplanat, we expect the first 'alpha'-versions of v4 to be out in September.

There are two versions of the search_indices_curvilinear() function:

1. One is written in C and is used for JITParticles:

   Parcels/parcels/include/index_search.h

   Lines 326 to 462 in f3a987e

   	static inline StatusCode search_indices_curvilinear(double time, type_coord z, type_coord y, type_coord x,
   	CStructuredGrid *grid, GridType gtype,
   	int ti, int *zi, int *yi, int *xi,
   	double *zeta, double *eta, double *xsi,
   	double t0, double t1, int interp_method,
   	int gridindexingtype)
   	{
   	int xi_old = *xi;
   	int yi_old = *yi;
   	int xdim = grid->xdim;
   	int ydim = grid->ydim;
   	int zdim = grid->zdim;
   	int tdim = grid->tdim;
   	float *xvals = grid->lon;
   	float *yvals = grid->lat;
   	float *zvals = grid->depth;
   	float *xy_minmax = grid->lonlat_minmax;
   	int sphere_mesh = grid->sphere_mesh;
   	int zonal_periodic = grid->zonal_periodic;
   	int z4d = grid->z4d;
   	// NEMO convention
   	float (* xgrid)[xdim] = (float (*)[xdim]) xvals;
   	float (* ygrid)[xdim] = (float (*)[xdim]) yvals;
   	if (zonal_periodic == 0){
   	if ((x < xy_minmax[0]) || (x > xy_minmax[1])){
   	if (xgrid[0][0] < xgrid[0][xdim-1]) {return ERROROUTOFBOUNDS;}
   	else if (x < xgrid[0][0] && x > xgrid[0][xdim-1]) {return ERROROUTOFBOUNDS;}
   	}
   	}
   	if ((y < xy_minmax[2]) || (y > xy_minmax[3]))
   	return ERROROUTOFBOUNDS;
   	double a[4], b[4];
   	*xsi = *eta = -1;
   	int maxIterSearch = 1e6, it = 0;
   	double tol = 1e-10;
   	while ( (*xsi < -tol) || (*xsi > 1+tol) || (*eta < -tol) || (*eta > 1+tol) ){
   	double xgrid_loc[4] = {xgrid[*yi][*xi], xgrid[*yi][*xi+1], xgrid[*yi+1][*xi+1], xgrid[*yi+1][*xi]};
   	if (sphere_mesh){ //we are on the sphere
   	int i4;
   	if (xgrid_loc[0] < x - 225) xgrid_loc[0] += 360;
   	if (xgrid_loc[0] > x + 225) xgrid_loc[0] -= 360;
   	for (i4 = 1; i4 < 4; ++i4){
   	if (xgrid_loc[i4] < xgrid_loc[0] - 180) xgrid_loc[i4] += 360;
   	if (xgrid_loc[i4] > xgrid_loc[0] + 180) xgrid_loc[i4] -= 360;
   	}
   	}
   	double ygrid_loc[4] = {ygrid[*yi][*xi], ygrid[*yi][*xi+1], ygrid[*yi+1][*xi+1], ygrid[*yi+1][*xi]};
   	a[0] = xgrid_loc[0];
   	a[1] = -xgrid_loc[0] + xgrid_loc[1];
   	a[2] = -xgrid_loc[0] + xgrid_loc[3];
   	a[3] = xgrid_loc[0] - xgrid_loc[1] + xgrid_loc[2] - xgrid_loc[3];
   	b[0] = ygrid_loc[0];
   	b[1] = -ygrid_loc[0] + ygrid_loc[1];
   	b[2] = -ygrid_loc[0] + ygrid_loc[3];
   	b[3] = ygrid_loc[0] - ygrid_loc[1] + ygrid_loc[2] - ygrid_loc[3];
   	double aa = a[3]*b[2] - a[2]*b[3];
   	double bb = a[3]*b[0] - a[0]*b[3] + a[1]*b[2] - a[2]*b[1] + x*b[3] - y*a[3];
   	double cc = a[1]*b[0] - a[0]*b[1] + x*b[1] - y*a[1];
   	if (fabs(aa) < 1e-12) // Rectilinear cell, or quasi
   	*eta = -cc / bb;
   	else{
   	double det = sqrt(bb*bb-4*aa*cc);
   	if (det == det) // so, if det is nan we keep the xsi, eta from previous iter
   	*eta = (-bb+det)/(2*aa);
   	}
   	if ( fabs(a[1]+a[3]*(*eta)) < 1e-12 ) // this happens when recti cell rotated of 90deg
   	*xsi = ( (y-ygrid_loc[0]) / (ygrid_loc[1]-ygrid_loc[0]) +
   	(y-ygrid_loc[3]) / (ygrid_loc[2]-ygrid_loc[3]) ) * .5;
   	else
   	*xsi = (x-a[0]-a[2]* (*eta)) / (a[1]+a[3]* (*eta));
   	if ( (*xsi < 0) && (*eta < 0) && (*xi == 0) && (*yi == 0) )
   	return ERROROUTOFBOUNDS;
   	if ( (*xsi > 1) && (*eta > 1) && (*xi == xdim-1) && (*yi == ydim-1) )
   	return ERROROUTOFBOUNDS;
   	if (*xsi < -tol)
   	(*xi)--;
   	if (*xsi > 1+tol)
   	(*xi)++;
   	if (*eta < -tol)
   	(*yi)--;
   	if (*eta > 1+tol)
   	(*yi)++;
   	reconnect_bnd_indices(yi, xi, ydim, xdim, 0, sphere_mesh);
   	it++;
   	if ( it > maxIterSearch){
   	printf("Correct cell not found for (lat, lon) = (%f, %f) after %d iterations\n", y, x, maxIterSearch);
   	printf("Debug info: old particle indices: (yi, xi) %d %d\n", yi_old, xi_old);
   	printf(" new particle indices: (yi, xi) %d %d\n", *yi, *xi);
   	printf(" Mesh 2d shape: %d %d\n", ydim, xdim);
   	printf(" Relative particle position: (eta, xsi) %1.16e %1.16e\n", *eta, *xsi);
   	return ERROROUTOFBOUNDS;
   	}
   	}
   	if ( (*xsi != *xsi) || (*eta != *eta) ){ // check if nan
   	printf("Correct cell not found for (lat, lon) = (%f, %f))\n", y, x);
   	printf("Debug info: old particle indices: (yi, xi) %d %d\n", yi_old, xi_old);
   	printf(" new particle indices: (yi, xi) %d %d\n", *yi, *xi);
   	printf(" Mesh 2d shape: %d %d\n", ydim, xdim);
   	printf(" Relative particle position: (eta, xsi) %1.16e %1.16e\n", *eta, *xsi);
   	return ERROROUTOFBOUNDS;
   	}
   	if (*xsi < 0) *xsi = 0;
   	if (*xsi > 1) *xsi = 1;
   	if (*eta < 0) *eta = 0;
   	if (*eta > 1) *eta = 1;
   	StatusCode status;
   	if (zdim > 1){
   	switch(gtype){
   	case CURVILINEAR_Z_GRID:
   	status = search_indices_vertical_z(z, zdim, zvals, zi, zeta, gridindexingtype);
   	break;
   	case CURVILINEAR_S_GRID:
   	status = search_indices_vertical_s(time, z, tdim, ydim, xdim, zdim, zvals,
   	ti, zi, *yi, *xi, zeta, *eta, *xsi,
   	z4d, t0, t1, interp_method);
   	break;
   	default:
   	status = ERRORINTERPOLATION;
   	}
   	CHECKSTATUS(status);
   	}
   	else
   	*zeta = 0;
   	if ( (*xsi < 0) || (*xsi > 1) ) return ERRORINTERPOLATION;
   	if ( (*eta < 0) || (*eta > 1) ) return ERRORINTERPOLATION;
   	if ( (*zeta < 0) || (*zeta > 1) ) return ERRORINTERPOLATION;
   	return SUCCESS;
   	}

2. One is written in python and is used for ScipyParticles and when you directly call the function:

   Parcels/parcels/_index_search.py

   Lines 228 to 317 in f3a987e

   	def _search_indices_curvilinear(field: Field, time, z, y, x, ti=-1, particle=None, search2D=False):
   	if particle:
   	xi = particle.xi[field.igrid]
   	yi = particle.yi[field.igrid]
   	else:
   	xi = int(field.grid.xdim / 2) - 1
   	yi = int(field.grid.ydim / 2) - 1
   	xsi = eta = -1.0
   	grid = field.grid
   	invA = np.array([[1, 0, 0, 0], [-1, 1, 0, 0], [-1, 0, 0, 1], [1, -1, 1, -1]])
   	maxIterSearch = 1e6
   	it = 0
   	tol = 1.0e-10
   	if not grid.zonal_periodic:
   	if x < grid.lonlat_minmax[0] or x > grid.lonlat_minmax[1]:
   	if grid.lon[0, 0] < grid.lon[0, -1]:
   	_raise_field_out_of_bound_error(z, y, x)
   	elif x < grid.lon[0, 0] and x > grid.lon[0, -1]: # This prevents from crashing in [160, -160]
   	_raise_field_out_of_bound_error(z, y, x)
   	if y < grid.lonlat_minmax[2] or y > grid.lonlat_minmax[3]:
   	_raise_field_out_of_bound_error(z, y, x)
   	while xsi < -tol or xsi > 1 + tol or eta < -tol or eta > 1 + tol:
   	px = np.array([grid.lon[yi, xi], grid.lon[yi, xi + 1], grid.lon[yi + 1, xi + 1], grid.lon[yi + 1, xi]])
   	if grid.mesh == "spherical":
   	px[0] = px[0] + 360 if px[0] < x - 225 else px[0]
   	px[0] = px[0] - 360 if px[0] > x + 225 else px[0]
   	px[1:] = np.where(px[1:] - px[0] > 180, px[1:] - 360, px[1:])
   	px[1:] = np.where(-px[1:] + px[0] > 180, px[1:] + 360, px[1:])
   	py = np.array([grid.lat[yi, xi], grid.lat[yi, xi + 1], grid.lat[yi + 1, xi + 1], grid.lat[yi + 1, xi]])
   	a = np.dot(invA, px)
   	b = np.dot(invA, py)
   	aa = a[3] * b[2] - a[2] * b[3]
   	bb = a[3] * b[0] - a[0] * b[3] + a[1] * b[2] - a[2] * b[1] + x * b[3] - y * a[3]
   	cc = a[1] * b[0] - a[0] * b[1] + x * b[1] - y * a[1]
   	if abs(aa) < 1e-12: # Rectilinear cell, or quasi
   	eta = -cc / bb
   	else:
   	det2 = bb * bb - 4 * aa * cc
   	if det2 > 0: # so, if det is nan we keep the xsi, eta from previous iter
   	det = np.sqrt(det2)
   	eta = (-bb + det) / (2 * aa)
   	if abs(a[1] + a[3] * eta) < 1e-12: # this happens when recti cell rotated of 90deg
   	xsi = ((y - py[0]) / (py[1] - py[0]) + (y - py[3]) / (py[2] - py[3])) * 0.5
   	else:
   	xsi = (x - a[0] - a[2] * eta) / (a[1] + a[3] * eta)
   	if xsi < 0 and eta < 0 and xi == 0 and yi == 0:
   	_raise_field_out_of_bound_error(0, y, x)
   	if xsi > 1 and eta > 1 and xi == grid.xdim - 1 and yi == grid.ydim - 1:
   	_raise_field_out_of_bound_error(0, y, x)
   	if xsi < -tol:
   	xi -= 1
   	elif xsi > 1 + tol:
   	xi += 1
   	if eta < -tol:
   	yi -= 1
   	elif eta > 1 + tol:
   	yi += 1
   	(yi, xi) = _reconnect_bnd_indices(yi, xi, grid.ydim, grid.xdim, grid.mesh)
   	it += 1
   	if it > maxIterSearch:
   	print(f"Correct cell not found after {maxIterSearch} iterations")
   	_raise_field_out_of_bound_error(0, y, x)
   	xsi = max(0.0, xsi)
   	eta = max(0.0, eta)
   	xsi = min(1.0, xsi)
   	eta = min(1.0, eta)
   	if grid.zdim > 1 and not search2D:
   	if grid._gtype == GridType.CurvilinearZGrid:
   	try:
   	(zi, zeta) = search_indices_vertical_z(field.grid, field.gridindexingtype, z)
   	except FieldOutOfBoundError:
   	_raise_field_out_of_bound_error(z, y, x)
   	elif grid._gtype == GridType.CurvilinearSGrid:
   	(zi, zeta) = search_indices_vertical_s(field.grid, field.interp_method, time, z, y, x, ti, yi, xi, eta, xsi)
   	else:
   	zi = -1
   	zeta = 0
   	if not ((0 <= xsi <= 1) and (0 <= eta <= 1) and (0 <= zeta <= 1)):
   	_raise_field_sampling_error(z, y, x)
   	if particle:
   	particle.xi[field.igrid] = xi
   	particle.yi[field.igrid] = yi
   	particle.zi[field.igrid] = zi
   	return (zeta, eta, xsi, zi, yi, xi)

Good luck digging!