Seeking advice on "ParticleSet is empty" warning during 3D simulation

[Cyuan-hub]

Question

Question

Hi, @erikvansebille. The version of Parcels I used is 3.0.2.
I am trying to run a 3D simulation using CMEMS 3D data and the AdvectionRK4_3D and CheckOutOfBounds kernel. However, at every output time step, I encounter the following warning: WARNING: ParticleSet is empty on writing as array.

I have also tried reducing the output time (e.g., changing it to 300 seconds), but the issue still persists. The warning continues to appear at all output times. Could you kindly provide guidance on how I can further debug this issue?

Thank you very much in advance.

Best regards,
Cyuan

Supporting code/error messages

# My partial code is as follows.
# fieldset
uvfiles = sorted(glob(dataUV_path + '*_uvcur.nc'))
wfiles = sorted(glob(dataW_path + '*_wcur.nc'))
mesh_mask = mesh_path + 'coordinate_lonup.nc'

filenames = {"U": {"lon": mesh_mask, "lat": mesh_mask, "depth": wfiles[0], "data": uvfiles},
             "V": {"lon": mesh_mask, "lat": mesh_mask, "depth": wfiles[0], "data": uvfiles},
             "W": {"lon": mesh_mask, "lat": mesh_mask, "depth": wfiles[0], "data": wfiles}, }

variables = {"U": "uo",  "V": "vo",   "W": "wo",}

grid_dimensions = {"lon": "longitude", "lat": "latitude", "depth": "depth", }

dimensions = {'U': grid_dimensions, 'V': grid_dimensions,'W': grid_dimensions, }

fieldset = FieldSet.from_netcdf(filenames, variables, dimensions, 
allow_time_extrapolation = True, mesh='spherical', timestamps = timestamps)

# kernel
def CheckOutOfBounds(particle, fieldset, time):  
    if particle.state == StatusCode.ErrorOutOfBounds:
        particle.delete()

# ParticleSet
part_lons = np.array(lons)
part_lats = np.array(lats)
rls_depth = np.ones(part_lons.shape)*0.5
pset = ParticleSet.from_list(fieldset = fieldset, pclass = MyParticle, lon = part_lons,  lat = part_lats, depth = rls_depth, )

# Kernels and execution
kernels = [AdvectionRK4_3D, CheckOutOfBounds] #, 

outputdt = timedelta(hours = 1)
runtime =  timedelta(days = 6) # the total length of the run 

pset.execute(kernels, runtime = runtime, dt = 300, output_file = output_file)

# Some other information:

![field_units](https://github.com/user-attachments/assets/77a8e4db-df2c-467a-9b50-e162a06c68cd)

![plot_code](https://github.com/user-attachments/assets/f694ae94-2fad-4519-b16c-f09bd6ae3509)
![plot_result](https://github.com/user-attachments/assets/0091d1e8-0fb3-40ff-abf0-fb76208f950e)

# Paste your error message within this block
INFO: Output files are stored in H:\JupyterNotebookData\SCS_case3\3D\output\3D_test2.zarr.
  0%|                                                                                     | 0/518400.0 [00:00<?, ?it/s]WARNING: ParticleSet is empty on writing as array at time 0
  0%|                                                                       | 300.0/518400.0 [00:00<04:53, 1762.26it/s]WARNING: ParticleSet is empty on writing as array at time 3600
  1%|▌                                                                    | 3900.0/518400.0 [00:00<00:30, 17002.92it/s]WARNING: ParticleSet is empty on writing as array at time 7200
  1%|▉                                                                    | 7500.0/518400.0 [00:00<00:21, 24050.12it/s]WARNING: ParticleSet is empty on writing as array at time 10800
  3%|█▉                                                                  | 14400.0/518400.0 [00:00<00:16, 30510.38it/s]WARNING: ParticleSet is empty on writing as array at time 14400
  3%|██▎                                                                 | 18000.0/518400.0 [00:00<00:15, 31327.94it/s]WARNING: ParticleSet is empty on writing as array at time 18000
  4%|██▊                                                                 | 21600.0/518400.0 [00:00<00:15, 32452.41it/s]WARNING: ParticleSet is empty on writing as array at time 21600
WARNING: ParticleSet is empty on writing as array at time 25200
  5%|███▎                                                                | 25500.0/518400.0 [00:00<00:16, 30601.35it/s]WARNING: ParticleSet is empty on writing as array at time 28800
  6%|████▎                                                               | 32400.0/518400.0 [00:01<00:13, 37374.72it/s]WARNING: ParticleSet is empty on writing as array at time 32400
WARNING: ParticleSet is empty on writing as array at time 36000
  7%|████▊                                                               | 36300.0/518400.0 [00:01<00:14, 34400.13it/s]WARNING: ParticleSet is empty on writing as array at time 39600
  8%|█████▋                                                              | 43200.0/518400.0 [00:01<00:12, 38924.56it/s]WARNING: ParticleSet is empty on writing as array at time 43200

[michaeldenes]

Hi @Cyuan-hub, likely what is happening is that all your particles are going out of bounds, and are hence deleted due to your CheckOutOfBounds kernel. It looks like there are some timesteps where your particles are advected without going out of bounds, maybe try plot the trajectories from your output file, something like

import xarray as xr
import matplotlib.pyplot as plt

# Load the ParticleFile
ds = xr.open_zarr(output_file)

plt.plot(ds['lon'].T, ds['lat'].T)
plt.show()

If you see the trajectories move to the boundaries of your domain (also check through bathymetry!) then it's simply that all of your particles are being deleted!

[Cyuan-hub]

Perhaps you can try printing the particle lon, lat, and depth in your CheckOutOfBounds kernel, which will tell you where your particle is when it throws the out of bounds error. Alternatively, remove that kernel and let the pset.execute() command fail, and see where the out of bounds error is thrown.

OK, I will try it.
When remove the CheckOutOfBounds kernel and let the pset.execute() command to fail, the result is as follows:

[michaeldenes]

Without seeing the mesh file or the data, it's hard to tell, but from the print statements and plots, I think it's likely fine. Your release longitudes look to be between your min and max fieldset longitudes, and same for latitudes. Your top depth layer is at 0.494025m, so releasing particles at 0.5m should be fine also...

When you plot the fieldset, why are you performing np.fliplr()? What does the following produce:

plt.pcolormesh(fieldset.U.grid.lon, fieldset.U.grid.lat, fieldset.U.data[0,0,:,:]
plt.show()
plt.pcolormesh(fieldset.V.grid.lon, fieldset.V.grid.lat, fieldset.V.data[0,0,:,:]
plt.show()
plt.pcolormesh(fieldset.W.grid.lon, fieldset.W.grid.lat, fieldset.W.data[0,0,:,:]
plt.show()

Just to see what the data looks like, and what the actual lon/lat limits of the data are for each field. It's possible that we've just focussed on the U grid, but something funky is going on in the V or W. Worth a shot!

[Cyuan-hub]

Hi, @michaeldenes. Thanks so much for your help with this! Your detailed feedback and suggestions have really boosted my skills and given me more confidence in getting better at physical oceanography and coding.

When I previously plotted the fieldset.U without using np.fliplr(), I noticed that the figure of the fieldset.U appeared clearly flipped horizontally (compared to the coastline), which seemed obviously incorrect. That’s why I added np.fliplr() to plot the fieldset.U.data.
The related codes and results are as follows:

It seems like something might have happened with the longitude and latitude in my U, V, and W data or coordinate data?

[michaeldenes]

@Cyuan-hub no problem at all! Yep, I think in your coordinate data you have your lons and lats flipped.
Parcels 'sees' the data like it is in the unflipped cases, so it's very likely that you have just got the order of your lons and lats wrong in the mesh_mask file you created!

[Cyuan-hub]

@Cyuan-hub no problem at all! Yep, I think in your coordinate data you have your lons and lats flipped. Parcels 'sees' the data like it is in the unflipped cases, so it's very likely that you have just got the order of your lons and lats wrong in the mesh_mask file you created!

@michaeldenes, thank you very very much!
OK. I'll recreate the mesh_mask file, and then pass it into Parcels. And I'll let you know how it goes!

[Cyuan-hub]

Hi, @michaeldenes. After correcting the longitude and latitude in the coordinate data, the issue with the warning 'ParticleSet is empty on writing as array at every output time step' has been resolved. Thanks again for your help!