Sinking kernel error

[BillyBoy996]

Question

Question

I seem to have an issue when trying to apply a constant sinking rate to particles so they act as propagules. I am using MITgcm model outputs. I tried to base myself off the ARGO float tutorial. However, I keep getting the following error message when trying to run the model.
This error only appears if I am running the kernel. I have tried positioning the kernels in different orders and also changing the simulation timestep but to no avail.

Supporting code/error messages

File /Applications/miniconda3/envs/parcels/lib/python3.13/site-packages/parcels/particleset.py:1163, in ParticleSet.execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, output_file, verbose_progress, postIterationCallbacks, callbackdt, delete_cfiles)
   [1160](https://file+.vscode-resource.vscode-cdn.net/Applications/miniconda3/envs/parcels/lib/python3.13/site-packages/parcels/particleset.py:1160) while (time < endtime and dt > 0) or (time > endtime and dt < 0):
   [1161](https://file+.vscode-resource.vscode-cdn.net/Applications/miniconda3/envs/parcels/lib/python3.13/site-packages/parcels/particleset.py:1161)     # Check if we can fast-forward to the next time needed for the particles
   [1162](https://file+.vscode-resource.vscode-cdn.net/Applications/miniconda3/envs/parcels/lib/python3.13/site-packages/parcels/particleset.py:1162)     if dt > 0:
-> [1163](https://file+.vscode-resource.vscode-cdn.net/Applications/miniconda3/envs/parcels/lib/python3.13/site-packages/parcels/particleset.py:1163)         skip_kernel = True if min(self.time) > (time + dt) else False
   [1164](https://file+.vscode-resource.vscode-cdn.net/Applications/miniconda3/envs/parcels/lib/python3.13/site-packages/parcels/particleset.py:1164)     else:
   [1165](https://file+.vscode-resource.vscode-cdn.net/Applications/miniconda3/envs/parcels/lib/python3.13/site-packages/parcels/particleset.py:1165)         skip_kernel = True if max(self.time) < (time + dt) else False

ValueError: min() iterable argument is empty

Paste your code within this block

def Sinking(particle, fieldset, time):
    """Sink particles at a constant rate"""
    w = 0.01
    # Update depth
    particle_ddepth += w * particle.dt

N = len(lon_start)   
times_np = np.repeat(np.datetime64(sim_start_date), N)
pset = parcels.ParticleSet(fieldset=fieldset, pclass=InvParticle, lon=lon_start, lat=lat_start, depth=start_depths, time=times_np, repeatdt=timedelta(days=seeding_frequency))
pfile = pset.ParticleFile(name=output_name, outputdt=timedelta(hours=output_frequency), chunks=(int(1e4), 1)) #save particles position at output freq
kernels = [Sinking, Surface, parcels.AdvectionRK4_3D, Distance, Ages, Settlement, CheckOutOfBounds]
pset.execute(kernels, runtime=timedelta(days=length_simulation), dt=timedelta(hours=sim_timestep), output_file=pfile)

[VeckoTheGecko]

I think the problem is in this line:

parcels.ParticleSet(fieldset=fieldset, pclass=InvParticle, lon=lon_start, lat=lat_start, depth=start_depths, time=times_np, repeatdt=timedelta(days=seeding_frequency))

Can you provide the whole script? (particularly I want to know how times_np is defined)

---

tip: you can use markdown syntax for quoting code. i.e.,

```python
print("Surround your code in triple backticks and specify the language")
```

which produces this when rendered

print("Surround your code in triple backticks and specify the language")

[BillyBoy996]

Ok sure I'll provide the full code below. The times_np is a datetime array defined by the simulation start date and any delay I might add on.

# ---------------------------
# Main loop: create fieldset and run parcels
# ---------------------------
for period in simulation_periods:
    year = period["year"]
    day_delay = period["day_delay"]
    ufiles = sorted(glob(f"{data_dir}/RGEMS3_{year}_v5.nc"))
    if len(ufiles) == 0:
        warnings.warn(f"No RGEMS3 files found for year {year} in {data_dir}. Skipping.")
        continue

    filenames = {
        "U": {"lon": mesh_mask, "lat": mesh_mask, "depth": mesh_mask, "data": ufiles},
        "V": {"lon": mesh_mask, "lat": mesh_mask, "depth": mesh_mask, "data": ufiles},
        "W": {"lon": mesh_mask, "lat": mesh_mask, "depth": mesh_mask, "data": ufiles},
    }
    c_grid_dimensions = {
        "lon": "XG",
        "lat": "YG",
        "depth": "Zl",
        "time": "time",
    }
    dimensions = {"U": c_grid_dimensions, "V": c_grid_dimensions, "W": c_grid_dimensions}
    # depth_indices = np.where(depths >= max_depth)[0] 
    # indices = {'depth': depth_indices.tolist()}

    fieldset = parcels.FieldSet.from_mitgcm(filenames, variables, dimensions, allow_time_extrapolation=True)

    ## Fieldset
    ## Update FieldSet (reload with new year’s data)
    fieldset.add_field(parcels.Field('HR_bathymetry', bathymetry_data, lon=bat_lon, lat=bat_lat, interp_method='linear'))
    fieldset.add_field(parcels.Field('bathymetry', depth_array.values, lon=depth_array.XC.values, lat=depth_array.YC.values, interp_method='linear'))
    fieldset.add_field(parcels.Field('HR_landmask', landmask, lon=bat_lon, lat=bat_lat, interp_method='nearest'))
    fieldset.add_field(parcels.Field('HR_stranding_mask', combined_mask, lon=bat_lon, lat=bat_lat, interp_method='nearest'))
    fieldset.add_field(parcels.Field('GMR_mask', gmr_mask, lon=bat_lon, lat=bat_lat, interp_method='linear'))
    fieldset.add_constant('advection_duration', advection_duration)
    fieldset.add_constant('max_depth', abs(max_depth))
    fieldset.add_constant('sink_velocity', 0.01)  # sinking velocity in m/s
    fieldset.add_field(parcels.Field('distance', distance_to_coastline,lon=bat_lon,lat=bat_lat, interp_method='linear'))
    fieldset.add_field(parcels.Field('nearest_edge_lat', nearest_edge_lat, lon=bat_lon, lat=bat_lat, interp_method='nearest'))
    fieldset.add_field(parcels.Field('nearest_edge_lon', nearest_edge_lon, lon=bat_lon, lat=bat_lat, interp_method='nearest'))

    # Update simulation start and end dates
    sim_start_date = datetime(year, 1, 1, 3, 0, 0) + timedelta(days=day_delay) # Compute simulation start and end dates
    sim_end_date = sim_start_date + timedelta(days=length_simulation) 
    start_str = sim_start_date.strftime("%d%m%y") # Format date range for output file name
    end_str = sim_end_date.strftime("%d%m%y")
    bat_grid_name = 'GEBCO' if use_GEBCO else 'ETOPO'
    Settle_depth_str = f"{abs(stranding_depth)}m" # Format depth string (remove minus sign)
    output_name = f"Inv2_InsideGMR_{year}_{Settle_depth_str}_{start_str}-{end_str}_{seeding_frequency}dayfreq_{bat_grid_name}.zarr"

    N = len(deploy_lons)   
    times_np = np.repeat(np.datetime64(sim_start_date), N)
    pset = parcels.ParticleSet(fieldset=fieldset, pclass=InvParticle, lon=deploy_lons, lat=deploy_lats, depth=start_depths, time=times_np, repeatdt=timedelta(days=seeding_frequency))
    pfile = pset.ParticleFile(name=output_name, outputdt=timedelta(hours=output_frequency), chunks=(int(1e4), 1)) #save particles position at output freq
    kernels = [Surface, parcels.AdvectionRK4_3D, Distance, Ages, Settlement, CheckOutOfBounds]
    print(f"Running simulation for {year} starting day {day_delay} with {N} particles...")
    pset.execute(kernels, runtime=timedelta(days=length_simulation), dt=timedelta(hours=sim_timestep), output_file=pfile)
    print(f"Simulation complete: {output_name}")

[nilodna]

Hi @BillyBoy996, not sure if this will fix, but can you try to convert your times_np to a list of seconds instead of a list of datetime64?

Just to test, do the following to release particles at the first available time on your hydrodynamic field.:

times_np = np.repeat([0], N)

If this test works, then probably the issue is the way you're defining your times_np, as @VeckoTheGecko pointed out initially.