Problem with zonal periodic BC and 3D advection - OutOfBoundsError

[eavellashaw]

Hi everyone!

I am currently struggling with the zonal periodic BC for my domain (see figure attached).

My domain goes from -280 to 80 deg E roughly, so I had to add a halo west and a halo_east to my fieldset:

fieldset.add_constant("halo_west", np.nanmin(fieldset.U.grid.lon))
fieldset.add_constant("halo_east", np.nanmax(fieldset.U.grid.lon))

fieldset.add_periodic_halo(zonal=True)

As suggested in #1333 checks are done at the different longitudes needed for the interpolation:

def Modified_AdvectionRK4_3D(particle, fieldset, time):
    # make sure we don't exceed the surface and the boundaries
    depmin = 5.033549785614014

    if particle.depth <depmin:
        particle.stuck_surf += 1
        particle.delete()

    if particle.lon < fieldset.halo_west:
        particle.lon += fieldset.halo_east - fieldset.halo_west
    elif particle.lon > fieldset.halo_east:
        particle.lon -= fieldset.halo_east - fieldset.halo_west
    
    (u1, v1, w1) = fieldset.UVW[time, particle.depth, particle.lat, particle.lon]

    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.stuck_surf += 1
        particle.delete()
    
    if lon1 < fieldset.halo_west:
        lon1 += fieldset.halo_east - fieldset.halo_west
    elif lon1 > fieldset.halo_east:
        lon1 -= fieldset.halo_east - fieldset.halo_west
  
    (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.stuck_surf += 1
        particle.delete()
    
    if lon2 < fieldset.halo_west:
        lon2 += fieldset.halo_east - fieldset.halo_west
    elif lon2 > fieldset.halo_east:
        lon2 -= fieldset.halo_east - fieldset.halo_west

    (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.stuck_surf += 1
        particle.delete()
    
    if lon3 < fieldset.halo_west:
        lon3 += fieldset.halo_east - fieldset.halo_west
    elif lon3 > fieldset.halo_east:
        lon3 -= fieldset.halo_east - fieldset.halo_west

    (u4, v4, w4) = fieldset.UVW[time + particle.dt, dep3, lat3, lon3]

    particle.lon += (u1 + 2*u2 + 2*u3 + u4) / 6. * particle.dt
    particle.lat += (v1 + 2*v2 + 2*v3 + v4) / 6. * particle.dt
    particle.depth += (w1 + 2*w2 + 2*w3 + w4) / 6. * particle.dt

    if particle.depth <depmin:
        particle.stuck_surf += 1
        particle.delete()
    
    if particle.lon < fieldset.halo_west:
        particle.lon += fieldset.halo_east - fieldset.halo_west
    elif particle.lon > fieldset.halo_east:
        particle.lon -= fieldset.halo_east - fieldset.halo_west

Still, I got the following error triggered:

Traceback (most recent call last):
  File "/aos/home/eashaw/Desktop/OP_Run_Scripts/Particles_Adv_annual.py", line 279, in <module>
    pset.execute(kernels,
  File "/aos/home/eashaw/miniconda3/envs/py3_parcels/lib/python3.11/site-packages/parcels/particleset/baseparticleset.py", line 532, in execute
    self.kernel.execute(self, endtime=next_time, dt=dt, recovery=recovery, output_file=output_file,
  File "/aos/home/eashaw/miniconda3/envs/py3_parcels/lib/python3.11/site-packages/parcels/kernel/kernelsoa.py", line 213, in execute
    recovery_kernel(p, self.fieldset, p.time)
  File "/aos/home/eashaw/miniconda3/envs/py3_parcels/lib/python3.11/site-packages/parcels/tools/statuscodes.py", line 182, in recovery_kernel_out_of_bounds
    raise OutOfBoundsError(particle, fieldset)
parcels.tools.statuscodes.OutOfBoundsError: 0
Particle P[208](lon=79.597166, lat=83.043988, depth=302.538240, stuck_land=0.000000, stuck_surf=0.000000, uvel=0.000001, vvel=0.000000, wvel=-0.000001, time=6688200.000000)
Time: 0191-03-18 20:05:00
timestep dt: 600.0
Out-of-bounds sampling by particle at (79.59716616254354, 83.04398831132214, 302.53824025989746)

Do you have an idea about what is wrong in my code?

Thanks in advance!

Esteban Avella Shaw

[michaeldenes]

Hi Esteban,

At a preliminary glance your code looks okay to me.

Could you share your fieldset creation code? As a first step, could you share what the following:
np.nanmin(fieldset.U.grid.lon))
np.nanmax(fieldset.U.grid.lon))
np.nanmin(fieldset.U.grid.lat))
np.nanmax(fieldset.U.grid.lat))
produce?
Are you providing any indices to your latitude/longitude dimensions that would cause you not to load the entire grid?

[eavellashaw]

Hi @michaeldenes, thanks a lot for your answer! Here is my fieldset creation code:

times = []
for f in files:
    ds = xr.open_dataset(f)
    times.append(ds.time.values)

file_names ={'U': {'lon': mesh_file_h, 'lat': mesh_file_h, 'depth': mesh_file_v, 'data': files},
	        'V': {'lon': mesh_file_h, 'lat': mesh_file_h, 'depth': mesh_file_v, 'data': files},
	        'W': {'lon': mesh_file_h, 'lat': mesh_file_h, 'depth': mesh_file_v, 'data': files}
            #  'T': {'lon': mesh_file_h, 'lat': mesh_file_h, 'depth': mesh_file_v, 'data': files},
            #  'S': {'lon': mesh_file_h, 'lat': mesh_file_h, 'depth': mesh_file_v, 'data': files}
            }
variables = {'U':'u_te',
	         'V':'v_tn',
	         'W':'wt'
            #  'T':'temp',
            #  'S':'salt'
             }
dimensions = {'U': {'lon': 'geolon_t', 'lat':'geolat_t', 'depth': 'st_ocean', 'time': 'time'},
                  'V': {'lon': 'geolon_t', 'lat':'geolat_t', 'depth': 'st_ocean', 'time': 'time'},
                  'W': {'lon': 'geolon_t', 'lat':'geolat_t', 'depth': 'st_ocean', 'time': 'time'}
                #   'T': {'lon': 'geolon_t', 'lat':'geolat_t', 'depth': 'st_ocean', 'time': 'time'},
                #   'S': {'lon': 'geolon_t', 'lat':'geolat_t', 'depth': 'st_ocean', 'time': 'time'}
                }

fieldset = FieldSet.from_c_grid_dataset(file_names,variables,dimensions, timestamps=times, allow_time_extrapolation=True) # only got one timestep

and here are the numbers you requested:

np.nanmin(fieldset.U.grid.lon) = -279.99948
np.nanmax(fieldset.U.grid.lon) = 79.99948
np.nanmin(fieldset.U.grid.lat) = -80.02101
np.nanmax(fieldset.U.grid.lat) = 89.94043

As far as I am concerned I am not providing any indices that would cause the entire grid to not be loaded. Also maybe worth to mention, I am working on a tripolar curvilinear grid (see the screenshots attached).

[michaeldenes]

The particle seems to be inside the domain. While it seems to be coming up close to the RHS boundary (lon=79.597166), it should be in the top part of the grid(with the two poles). I'm not all that sure about how things work above 60°N in these c-grids, my focus was always in the Southern hemisphere!

Parcels should be able to handle c-grids without the need for a periodic BC kernel or halos. As a quick check, have you tried running your code with no halos and with the simple RK4 advection (without the periodic BC checks)?

[eavellashaw]

The thing is that the velocity fields have been interpolated to be as if they were on a C-grid but are coming from outputs on a B-grid. Maybe that could also be what is causing the problem as there may have been an error during the interpolation. But in my opinion that should not cause too much of a problem.

I will try running it without the BC checks and halos for the interpolated velocities, as well as with the BC checks and halos but with the B-grid outputs, I will keep you posted!

[michaeldenes]

You can also try (in your current configuration) sampling the velocity field manually:
fieldset.UV[6688200, 302.53824025989746, 83.04398831132214, 79.59716616254354]
This will sample the velocity field at the current particle location. You can perturb the lat/lon values a bit to see if you can reproduce the out of bounds error. If you get valid values, and no out of bounds error, then something strange is going on!