ParticleSet Empty at Beginning When Using Custom Kernel

[jamespolly]

Hello! The Kernel loop example shows how Parcels v3.0 allows for vectorized summing from multiple fields.

I've created a custom kernel that is slightly modified from the pre-built AdvectionRK4. Only the variable names have been changed.

The problem I'm experiencing has been recreated in the following script. Everything above idemo = 0 is from the above linked example and works. Everything below is my reproducible example resulting in the error: WARNING: ParticleSet is empty on writing as array at time 0

What I want to happen is to bring in multiple fields from multiple NetCDFs, add them to the first fieldset, and have different kernels operate on these fields (via VectorField objects).

What am I doing wrong to have no particles at the beginning of the simulation in my broken example below idemo = 0?

Here is a reproducible example below:

from datetime import timedelta

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

from parcels import (
    AdvectionRK4,
    Field,
    FieldSet,
    JITParticle,
    ParticleSet,
    download_example_dataset,
    StatusCode,
)
from parcels.field import VectorField


def delete_particle(
    particle: JITParticle, fieldset: FieldSet, time: timedelta
) -> None:  # pragma: no cover
    """
    Delete particles when they drift out of defined vector field domain.

    !! Note that including a return statement here will result in compilation errors
    for parcels v3.0.0 !!

    :param particle: particle to be deleted.
    :param fieldset: ``FieldSet`` instance that contains the ``particle``.
    :param time: time delta between steps by the ``particle``.
    """
    if particle.state == StatusCode.ErrorOutOfBounds:
        particle.delete()


# Load the GlobCurrent data in the Agulhas region from the example_data
example_dataset_folder = download_example_dataset("GlobCurrent_example_data")
filenames = {
    "U": f"{example_dataset_folder}/20*.nc",
    "V": f"{example_dataset_folder}/20*.nc",
}
variables = {
    "U": "eastward_eulerian_current_velocity",
    "V": "northward_eulerian_current_velocity",
}
dimensions = {"lat": "lat", "lon": "lon", "time": "time"}
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)
# uppermost layer in the hydrodynamic data
fieldset.mindepth = fieldset.U.depth[0]

lons = 26.0 * np.ones(10)
lats = np.linspace(-37.5, -34.5, 10)

idemo = 0  # set this to 1 to run the example in the docs:
# https://docs.oceanparcels.org/en/latest/examples/tutorial_kernelloop.html#Example-use
if idemo == 1:
    pset = ParticleSet(fieldset, pclass=JITParticle, lon=lons, lat=lats)
    output_file = pset.ParticleFile(
        name="advection_then_wind.zarr", outputdt=timedelta(hours=6)
    )
    pset.execute(
        [AdvectionRK4, delete_particle],
        runtime=timedelta(days=5),
        dt=timedelta(hours=1),
        output_file=output_file,
    )
else:
    filenames_wind = {
        "UWind": f"{example_dataset_folder}/20*.nc",
        "VWind": f"{example_dataset_folder}/20*.nc",
    }
    variables_wind = {
        "UWind": "eastward_eulerian_current_velocity",
        "VWind": "northward_eulerian_current_velocity",
    }
    fieldset_wind = FieldSet.from_netcdf(filenames_wind, variables_wind, dimensions)

    fieldset.add_field(fieldset_wind.UWind, name="UWind")
    fieldset.add_field(fieldset_wind.VWind, name="VWind")
    wind_scale_factor = 0.5
    fieldset.UWind.set_scaling_factor(wind_scale_factor)
    fieldset.VWind.set_scaling_factor(wind_scale_factor)
    fieldset.add_vector_field(VectorField("UVWind", fieldset.UWind, fieldset.VWind))

    def WindAdvectionRK4(particle, fieldset, time):
        """Advection of particles using fourth-order Runge-Kutta integration.

        Function needs to be converted to Kernel object before execution.
        Based on pre-built Kernel `AdvectionRK4`.
        Requires fieldset with `UWind` and `VWind` fields.
        """
        (uw1, vw1) = fieldset.UVWind[particle]
        lonw1, latw1 = (
            particle.lon + uw1 * 0.5 * particle.dt,
            particle.lat + vw1 * 0.5 * particle.dt,
        )
        (uw2, vw2) = fieldset.UVWind[
            time + 0.5 * particle.dt, particle.depth, latw1, lonw1, particle
        ]
        lonw2, latw2 = (
            particle.lon + uw2 * 0.5 * particle.dt,
            particle.lat + vw2 * 0.5 * particle.dt,
        )
        (uw3, vw3) = fieldset.UVWind[
            time + 0.5 * particle.dt, particle.depth, latw2, lonw2, particle
        ]
        lonw3, latw3 = (
            particle.lon + uw3 * particle.dt,
            particle.lat + vw3 * particle.dt,
        )
        (uw4, vw4) = fieldset.UVWind[
            time + particle.dt, particle.depth, latw3, lonw3, particle
        ]
        particle_dlon += (uw1 + 2 * uw2 + 2 * uw3 + uw4) / 6.0 * particle.dt  # noqa
        particle_dlat += (vw1 + 2 * vw2 + 2 * vw3 + vw4) / 6.0 * particle.dt  # noqa

    pset = ParticleSet(fieldset, pclass=JITParticle, lon=lons, lat=lats)
    output_file = pset.ParticleFile(
        name="advection_then_wind_reproducible.zarr", outputdt=timedelta(hours=6)
    )
    pset.execute(
        [WindAdvectionRK4, delete_particle],
        runtime=timedelta(days=5),
        dt=timedelta(hours=1),
        output_file=output_file,
    )

[erikvansebille]

Hi @jamespolly, thanks for using parcels!

It's hard to be sure without having access to the hydrodynamic data, but one thing we've seen before is that particles get deleted if the longitude definitions between the fields don't align. For example when the longitudes of your ocean flow fields are defined between 0 and 360 degrees, and for the wind fields between -180 and 180 degrees. Particles could then get an out-of-bounds in one field, but not the other.

I assume that if you don't add delete_particle to your kernel list, you will get an OutOfBounds error? The details of that error (which field triggered it where) might help you get to the bottom of this

[jamespolly]

Thanks so much for your response @erikvansebille!

A quick response to answer the points you raised prior to me digging deeper:

1. The hydrodynamic data I'm using in the above example is coming from the kernel loop example (i.e. download_example_dataset("GlobCurrent_example_data")). I've tried to use that example from the docs to recreate my end goal which is summing fields created via fieldset.from_netcdf().
2. Removing the delete_particle() method from the kernel list is a great idea! The resulting error is FieldOutOfBoundError: Field sampled out-of-bound, at (26.0, -37.5, 0.0). Note that in the example above only one set of points is defined, and everything works for idemo = 1.
3. The custom kernel I define has two differences from the pre-built AdvectionRK4 kernel. First, the VectorField name is changed from fieldset.UV to fieldset.UVWind. Second, I have renamed variables (e.g., u1 becomes uw1to avoid the warning when having multiple kernels with the same variable names. So, to check the FieldOutOfBoundError I check the latitude and longitude attributes of these vector fields (only U shown, V is the same):

In [14]: fieldset.UV.U.lon
Out[14]: 
array([14.875, 15.125, 15.375, 15.625, 15.875, 16.125, 16.375, 16.625,
       16.875, 17.125, 17.375, 17.625, 17.875, 18.125, 18.375, 18.625,
       18.875, 19.125, 19.375, 19.625, 19.875, 20.125, 20.375, 20.625,
       20.875, 21.125, 21.375, 21.625, 21.875, 22.125, 22.375, 22.625,
       22.875, 23.125, 23.375, 23.625, 23.875, 24.125, 24.375, 24.625,
       24.875, 25.125, 25.375, 25.625, 25.875, 26.125, 26.375, 26.625,
       26.875, 27.125, 27.375, 27.625, 27.875, 28.125, 28.375, 28.625,
       28.875, 29.125, 29.375, 29.625, 29.875, 30.125, 30.375, 30.625,
       30.875, 31.125, 31.375, 31.625, 31.875, 32.125, 32.375, 32.625,
       32.875, 33.125, 33.375, 33.625, 33.875, 34.125, 34.375, 34.625,
       34.875], dtype=float32)

In [15]: fieldset.UVWind.U.lon
Out[15]: 
array([14.875, 15.125, 15.375, 15.625, 15.875, 16.125, 16.375, 16.625,
       16.875, 17.125, 17.375, 17.625, 17.875, 18.125, 18.375, 18.625,
       18.875, 19.125, 19.375, 19.625, 19.875, 20.125, 20.375, 20.625,
       20.875, 21.125, 21.375, 21.625, 21.875, 22.125, 22.375, 22.625,
       22.875, 23.125, 23.375, 23.625, 23.875, 24.125, 24.375, 24.625,
       24.875, 25.125, 25.375, 25.625, 25.875, 26.125, 26.375, 26.625,
       26.875, 27.125, 27.375, 27.625, 27.875, 28.125, 28.375, 28.625,
       28.875, 29.125, 29.375, 29.625, 29.875, 30.125, 30.375, 30.625,
       30.875, 31.125, 31.375, 31.625, 31.875, 32.125, 32.375, 32.625,
       32.875, 33.125, 33.375, 33.625, 33.875, 34.125, 34.375, 34.625,
       34.875], dtype=float32)

In [16]: fieldset.UV.U.lat
Out[16]: 
array([-40.125, -39.875, -39.625, -39.375, -39.125, -38.875, -38.625,
       -38.375, -38.125, -37.875, -37.625, -37.375, -37.125, -36.875,
       -36.625, -36.375, -36.125, -35.875, -35.625, -35.375, -35.125,
       -34.875, -34.625, -34.375, -34.125, -33.875, -33.625, -33.375,
       -33.125, -32.875, -32.625, -32.375, -32.125, -31.875, -31.625,
       -31.375, -31.125, -30.875, -30.625, -30.375, -30.125],
      dtype=float32)

In [17]: fieldset.UVWind.U.lat
Out[17]: 
array([-40.125, -39.875, -39.625, -39.375, -39.125, -38.875, -38.625,
       -38.375, -38.125, -37.875, -37.625, -37.375, -37.125, -36.875,
       -36.625, -36.375, -36.125, -35.875, -35.625, -35.375, -35.125,
       -34.875, -34.625, -34.375, -34.125, -33.875, -33.625, -33.375,
       -33.125, -32.875, -32.625, -32.375, -32.125, -31.875, -31.625,
       -31.375, -31.125, -30.875, -30.625, -30.375, -30.125],
      dtype=float32)

So the error indicating that the field is sampled out of bounds at (26.0, -37.5, 0.0) is odd given that this point seems to be within the bounds of the domain.

[erikvansebille]

Ah, in that case what's probably happening is that the unitconverters are not properly set on the 'wind'-fields.

If you're on a spherical mesh (i.e. in lon/lat), parcels automatically converts velocities in m/s to degrees/m using these so-called unitconverters.

It does so based on the Field -names, where U and V (and als Kh_zonal and Kh_meridional) are 'special', see the list here.

These unitconverters are probably not set for your UWind and VWind fields, meaning that the velocities (in m/s) are interpreted as in degrees/s by the advection kernel.

Your particles hence immediately fly off the grid (with a dt of 1 hour, they could move 3600 degrees in one time step!)

So the solution is (likely) to add these unitconverters

from parcels.tools.converters import GeographicPolar

fieldset.Uwind.units = GeographicPolar()
fieldset.Vwind.units = Geographic()

See for much more details the Spherical grids and UnitConversion tutorial.

[jamespolly]

This was indeed the problem @erikvansebille, thank you so much!

As I was investigating, and as you pointed out, the computation of the particle_dlat and particle_dlon were really big values...

I created two fieldsets using .from_netcdf() but with different keys for the filenames and variables dictionaries provided as arguments to the .from_netcdf() method. One set of dictionaries used 'U' and 'V' for the keys, whereas the other dictionary used 'UWind' and 'VWind' for the keys (as I have in the opening post).

I could see that fieldset.U.data was the same as fieldset.UWind.data and several other checks of this kind on .data didn't reveal any discrepancies.

However, when I would test some lines from the kernel, differences started to appear.

pset = ParticleSet(fieldset, pclass=JITParticle, lon=lons, lat=lats)
output_file = pset.ParticleFile(
    name="advection_then_wind_reproducible.zarr", outputdt=timedelta(hours=6)
)

fieldset.computeTimeChunk()
particle = pset[0]
(u1, v1) = fieldset.UV[particle]
print(u1, v1)

Would not agree with the output of:

pset_wind = ParticleSet(fieldset_wind, pclass=JITParticle, lon=lons, lat=lats)
output_file_wind = pset.ParticleFile(
    name="advection_then_wind_reproducible_wind.zarr", outputdt=timedelta(hours=6)
)

fieldset_wind.computeTimeChunk()
particle_wind = pset_wind[0]
(uw1, vw1) = fieldset_wind.UVWind[particle_wind]
print(uw1, vw1)

Despite fieldset.UV.U.data being equal to fieldset.UVWind.U.data.

Just adding this extra info here in case it is useful to anyone down the line. Thank you again @erikvansebille for your help!

[erikvansebille]

Great to hear this fixed your issue.

Now that you looked so much at this: Do you have any idea how we could make this issue of unit converters for user-defined Fields either clearer in the documentation, or to somehow infer their intended unit converter in the code? I see this going wrong more often... Suggestions very welcome!

[jamespolly]

@erikvansebille thank you for asking!

I have a couple thoughts:

1. The quickest thing to do would be to make a notebook or docs page similar to the kernel loop example but specifically for the case of using multiple fields originating from the .from_netcdf() method. This would need to go through the steps we've identified in this discussion. 😄
2. When a user tries to name their own fields, as I was doing above (e.g., "UWind"), I found that there were some failures due to an assertion in Fieldset.check_complete() which wants to see the name and the type (e.g. U or V) be the same. Without understanding fully why this check is there, I can imagine there being some similar checks on these quantities, and on the grids/converters that are associated with these fields that might generate warnings with instructions when grids/converters are not specified.

I'd love to contribute here and maybe this is a good place to start. I have a couple other ideas that I'll try to record here later.

[erikvansebille]

Thanks @jamespolly, these seem like good ideas. Your first option of a tutorial could be appended/part of the unit converter tutorial then?

And yes, your second option of throwing a warning in fieldset.check_complete() could also be interesting to explore, especially if we can write some rules that make it likely the user supplies a velocity field. E.g.:

1. The mesh is set to spherical (otherwise no need for UnitConverter, so this is a necessary condition)
2. The unitconverter is not explicitly set
3. The field-name starts with U or V; or the field name contains a u or v?
4. All values in the data are between say -2 and 2 (m/s); although note that wind speeds may be higher(?)
5. ...