Out of boundaries error

[mainarabgouveia]

Hi,

How do I make my released particles that leave the domain or arrive on the land don't cause my simulations to stop?

I'm getting this error:

INFO: Compiled ArrayJITParticleAdvectionRK45AdvectionDiffusionM1 ==> /tmp/parcels-1000/lib77962480406540cb0477e74d24ef8364_0.so
98% (14938200.0 of 15120000.0) |####### | Elapsed Time: 0:15:43 ETA: 0:00:11

OutOfBoundsError Traceback (most recent call last)
/tmp/ipykernel_32755/3033041448.py in
1 # Now run the remaining ParticleSet (using less accurate Euler Forward for speedup)
2
----> 3 pset.execute(pset.Kernel(AdvectionRK45) + AdvectionDiffusionM1, # the kernel (which defines how particles move)
4 runtime=timedelta(days=175), # the total length of the run
5 dt=-timedelta(minutes=5), # the timestep of the kernel (negative timestep for backwards run)

~/anaconda3/envs/py3_parcels/lib/python3.9/site-packages/parcels/particleset/baseparticleset.py in execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, moviedt, recovery, output_file, movie_background_field, verbose_progress, postIterationCallbacks, callbackdt)
456 # If we don't perform interaction, only execute the normal kernel efficiently.
457 if self.interaction_kernel is None:
--> 458 self.kernel.execute(self, endtime=next_time, dt=dt, recovery=recovery, output_file=output_file,
459 execute_once=execute_once)
460 # Interaction: interleave the interaction and non-interaction kernel for each time step.

~/anaconda3/envs/py3_parcels/lib/python3.9/site-packages/parcels/kernel/kernelsoa.py in execute(self, pset, endtime, dt, recovery, output_file, execute_once)
225 recovery_kernel = recovery_map[p.state]
226 p.set_state(StateCode.Success)
--> 227 recovery_kernel(p, self.fieldset, p.time)
228 if p.isComputed():
229 p.reset_state()

~/anaconda3/envs/py3_parcels/lib/python3.9/site-packages/parcels/tools/statuscodes.py in recovery_kernel_out_of_bounds(particle, fieldset, time)
194 # TODO: JIT does not yet provide the context that created
195 # the exception. We need to pass that info back from C.
--> 196 raise OutOfBoundsError(particle, fieldset)
197 else:
198 error = particle.exception

OutOfBoundsError: 0
Particle P[4](lon=17.006130, lat=-29.650196, depth=0.000000, time=36208800.000000)
Time: 2021-02-23T02:30:00.000000000, timestep dt: -3600.000000
Out-of-bounds sampling by particle at (17.006130, -29.650196, 0.000000)

[erikvansebille]

Dear @mainarabgouveia. You will have to include a recovery kernel to your execution:

pset.execute(..., recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})

And then create a Recovery Kernel to delete the particle:

def DeleteParticle(particle, fieldset, time):
    particle.delete()

See also https://github.com/OceanParcels/parcels/blob/master/parcels/examples/tutorial_Agulhasparticles.ipynb for a full example

[mollyjames2]

Hi Erik, I have included the DeleteErrorParticle kernel you mentioned above and included it in my set execute

    particle.age += particle.dt 

def Sample_land(particle, fieldset, time):
    particle.on_land = fieldset.landmask[time, particle.depth, particle.lat, particle.lon]

def Particle_death(particle, fieldset, time):
    runtime = 2419200 # 28 days (in seconds)
    if particle.age > runtime:
        particle.delete()
        
def DeleteParticle(particle, fieldset, time):
    particle.delete()

def Distance(particle, fieldset, time):
    lat_dist = (particle.lat - particle.prev_lat) * 1.11e2 # latitudinal distance
    lon_dist = (particle.lon - particle.prev_lon) * 1.11e2 * math.cos(particle.lat * math.pi/180) # longitudinal distance
    particle.distance += math.sqrt(math.pow(lon_dist, 2) + math.pow(lat_dist, 2)) # calculate Euclidean distance
    particle.prev_lon = particle.lon # Set the values for the next time step   
    particle.prev_lat = particle.lat


def DeleteErrorParticle(particle, fieldset, time):
    if particle.state == StatusCode.ErrorOutOfBounds:
        particle.delete()

#Execute an advection kernel on the given fieldset

output_file = pset.ParticleFile(name="/users/modellers/moja/projects/SANH/NEMO-PARCELS/NEMO_parcels_SANH.zarr", outputdt=timedelta(hours=1)) # the file name and the time step of the outputs

# Cutom kernels need to be cast into a kernel object using pset.Kernel before they can function 
kernels = pset.Kernel(AdvectionRK4)+ pset.Kernel(SampleAge)+ pset.Kernel(Sample_land)+pset.Kernel(Distance)+pset.Kernel(DeleteErrorParticle)+pset.Kernel(DiffusionUniformKh)

pset.execute(kernels,              # the kernel (which defines how particles move)
            runtime=timedelta(days=14),    # the total length of the run
            dt=timedelta(minutes=1),      # the timestep of the kernel (negative timestep for backwards run)
            output_file=output_file)

However, the model is still terminating

---------------------------------------------------------------------------
FieldOutOfBoundError                      Traceback (most recent call last)
Cell In[33], line 15
     12 # Cutom kernels need to be cast into a kernel object using pset.Kernel before they can function 
     13 kernels = pset.Kernel(AdvectionRK4)+ pset.Kernel(SampleAge)+ pset.Kernel(Sample_land)+pset.Kernel(Distance)+pset.Kernel(DeleteErrorParticle)+pset.Kernel(DiffusionUniformKh)
---> 15 pset.execute(kernels,              # the kernel (which defines how particles move)
     16             runtime=timedelta(days=14),    # the total length of the run
     17             dt=timedelta(minutes=1),      # the timestep of the kernel (negative timestep for backwards run)
     18             output_file=output_file)

File ~/miniforge3/envs/parcels/lib/python3.11/site-packages/parcels/particleset.py:948, in ParticleSet.execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, output_file, verbose_progress, postIterationCallbacks, callbackdt)
    946 # If we don't perform interaction, only execute the normal kernel efficiently.
    947 if self.interaction_kernel is None:
--> 948     res = self.kernel.execute(self, endtime=next_time, dt=dt)
    949     if res == StatusCode.StopAllExecution:
    950         return StatusCode.StopAllExecution

File ~/miniforge3/envs/parcels/lib/python3.11/site-packages/parcels/kernel.py:585, in Kernel.execute(self, pset, endtime, dt)
    583     raise TimeExtrapolationError(p.time)
    584 elif p.state == StatusCode.ErrorOutOfBounds:
--> 585     raise FieldOutOfBoundError(p.lon, p.lat, p.depth)
    586 elif p.state == StatusCode.ErrorThroughSurface:
    587     raise FieldOutOfBoundSurfaceError(p.lon, p.lat, p.depth)

FieldOutOfBoundError: Field sampled out-of-bound, at (81.19326152956675, 5.092092458112055, 0.0)

Do you have any insight on why this may be happening?
I am running parcels version 3.0.0

Thanks in advance

[erikvansebille]

It's really important that the DeleteErrorParticle kernel is the last in you list of kernels. So make it

kernels = [AdvectionRK4, DiffusionUniformKh, SampleAge, Sample_land, Distance, DeleteErrorParticle]

See also the last paragraph of https://docs.oceanparcels.org/en/latest/examples/tutorial_kernelloop.html#Working-with-Status-Codes

[nntaganou]

Hi Erik,

I am getting an out-of-bound error while using the AdvectionAnalytical advection kernel. Even though I did use the above kernel that deletes the particles at the boundaries, the code still gives me the out-of-bound error for U velocity, while the particle state is still not over 40.

This is the code:

variables = {'U': 'u',
             'V': 'v'}   
dimensions = {'lat': 'Latitude',
              'lon': 'Longitude',
              'time': 'MT'}
 
fieldset = FieldSet.from_xarray_dataset(ds_all_data,variables,dimensions,allow_time_extrapolation=True)
fieldset.mindepth = fieldset.U.depth[0]
fieldset.mindepth = fieldset.U.depth[0]
fieldset.U.interp_method = "cgrid_velocity"
fieldset.V.interp_method = "cgrid_velocity"   

x = lat_value
y = ptcs_lon_all
lons, lats = np.meshgrid(y,x)
psetAA = ParticleSet(fieldset=fieldset, 
                    pclass=ScipyParticle, 
                    lon=lons, 
                    lat=lats,
                    time=0) 

folder_output = '/unity/g2/nntaganou/parcels_gom/'
name = 'test_release_yc_1993_3h_inter_step_daily_data'
file_output = join(f"{folder_output}", f"{name}.zarr")
output = psetAA.ParticleFile(name=file_output, outputdt=timedelta(hours=3))


def CheckOutOfBounds(particle, fieldset, time):
    if particle.state == StatusCode.ErrorOutOfBounds:
        particle.delete()
    
start_date = datetime(1993,1,2,0)
end_date = datetime(1994,1,1,0)
run_time = end_date-start_date
kernels = [AdvectionAnalytical,CheckOutOfBounds]

psetAA.execute(kernels,
            runtime=run_time, 
            dt=timedelta(hours=3), 
            output_file=output)

netcdf_file = join(f'{folder_output}', f'{name}.nc')
ds= xr.open_zarr(f"{file_output}")

This is the error:

---------------------------------------------------------------------------
FieldOutOfBoundError                      Traceback (most recent call last)
Cell In[13], line 1
----> 1 psetAA.execute(AdvectionAnalytical,
      2             runtime=run_time, 
      3             #endtime=end_date)
      4             dt=timedelta(hours=3), 
      5             output_file=output)

File ~/mambaforge/envs/testenv/lib/python3.11/site-packages/parcels/particleset.py:951, in ParticleSet.execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, output_file, verbose_progress, postIterationCallbacks, callbackdt)
    949 # If we don't perform interaction, only execute the normal kernel efficiently.
    950 if self.interaction_kernel is None:
--> 951     res = self.kernel.execute(self, endtime=next_time, dt=dt)
    952     if res == StatusCode.StopAllExecution:
    953         return StatusCode.StopAllExecution

File ~/mambaforge/envs/testenv/lib/python3.11/site-packages/parcels/kernel.py:564, in Kernel.execute(self, pset, endtime, dt)
    562     self.execute_jit(pset, endtime, dt)
    563 else:
--> 564     self.execute_python(pset, endtime, dt)
    566 # Remove all particles that signalled deletion
    567 self.remove_deleted(pset)

File ~/mambaforge/envs/testenv/lib/python3.11/site-packages/parcels/kernel.py:543, in Kernel.execute_python(self, pset, endtime, dt)
    540     self.scipy_positionupdate_kernels_added = True
    542 for p in pset:
--> 543     self.evaluate_particle(p, endtime)
    544     if p.state == StatusCode.StopAllExecution:
    545         return StatusCode.StopAllExecution

File ~/mambaforge/envs/testenv/lib/python3.11/site-packages/parcels/kernel.py:628, in Kernel.evaluate_particle(self, p, endtime)
    626 if abs(endtime - p.time_nextloop) < abs(p.dt)-1e-6:
    627     p.dt = abs(endtime - p.time_nextloop) * sign_dt
--> 628 res = self._pyfunc(p, self._fieldset, p.time_nextloop)
    630 if res is None:
    631     if sign_dt*p.time < sign_dt*endtime and p.state == StatusCode.Success:

File <ast>:24, in SetcoordsAdvectionAnalyticalUpdatecoords(particle, fieldset, time)

File ~/mambaforge/envs/testenv/lib/python3.11/site-packages/parcels/field.py:998, in Field.search_indices(self, x, y, z, ti, time, particle, search2D)
    996 def search_indices(self, x, y, z, ti=-1, time=-1, particle=None, search2D=False):
    997     if self.grid.gtype in [GridCode.RectilinearSGrid, GridCode.RectilinearZGrid]:
--> 998         return self.search_indices_rectilinear(x, y, z, ti, time, particle=particle, search2D=search2D)
    999     else:
   1000         return self.search_indices_curvilinear(x, y, z, ti, time, particle=particle, search2D=search2D)

File ~/mambaforge/envs/testenv/lib/python3.11/site-packages/parcels/field.py:822, in Field.search_indices_rectilinear(self, x, y, z, ti, time, particle, search2D)
    820         raise FieldOutOfBoundError(x, y, z, field=self)
    821 if grid.ydim > 1 and (y < grid.lonlat_minmax[2] or y > grid.lonlat_minmax[3]):
--> 822     raise FieldOutOfBoundError(x, y, z, field=self)
    824 if grid.xdim > 1:
    825     if grid.mesh != 'spherical':

FieldOutOfBoundError: U sampled out-of-bound, at (-79.0750963318386, 31.960647583007816, 0.0)

The (x, y) at (-79.0750963318386, 31.960647583007816) are the upper bounds of my domain.

Do you have any ideas as to why I am getting this error? This is the first time I am using analytical advection to advect particles and I am not sure if I am missing something.

Thank you in advance!

[erikvansebille]

I think in this case, the error is in the AnalyticalAdvection scheme. This was coded up as a quick showcase that it would in principle be possible to implement such a scheme within parcels, but wasn't well-tested. I.e. no proper error-handling. You've run into the limits of its capabilities, I fear...

However, my next big project the coming weeks/months is to 'properly' implement AnalyticalAdvection into Parcels, also in JIT-mode. I will then make sure that these out-of-bounds errors are properly treated, too.

So I suggest you either wait a few more weeks and/or use RungeKutta for the timing being; or write your own removal kernel to remove particles before they even get close to the boundary (if that is even an option for your use-case)

[nntaganou]

Thank you for letting me know!

I might be able to implement a removal kernel after the parcels exit the Straits of Florida -my region is the Gulf of Mexico, so something like that could work.

In your future project, is the fully time-dependent scheme going to be implemented or is it going to be a further development of the stepwise scheme?

Thank you!

[mainarabgouveia]

Tks for your help and time.