Separate fieldsets working properly, but loop kernels failing systematicaly

[AlxLhrNc]

Question

Question

Following the parcels tutorial tutorial_kernelloop.ipynb, I tried to simulate 14 drifters being impacted by currents and wind, accounting for a 3% wind impact. While the current OR wind fieldsets work great independently under the native AdvectionRK4 kernel, using both together with either custom kernels (WindInteraction or AdvectionRK4_wind) always result in the particles leaving the fieldset in ~4 dt every time.

Any idea why? Thanks in advance

Code

import os
import parcels
from glob import glob

def WindInteraction(particle, fieldset, time):
    '''
    doi for wind_factor: https://doi.org/10.1175/JPO-D-20-0275.1
    '''
    wind_factor = 0.03 # adjust to control the influence of the wind 0.03
    particle_dlon += (
        fieldset.UWind[time, particle.depth, particle.lat, particle.lon] * wind_factor * particle.dt
    )
    particle_dlat += (
        fieldset.VWind[time, particle.depth, particle.lat, particle.lon] * wind_factor * particle.dt
    )

def AdvectionRK4_wind(particle, fieldset, time):  # pragma: no cover
    """Advection of particles using fourth-order Runge-Kutta integration."""
    (u1, v1) = fieldset.UWind[particle], fieldset.VWind[particle]
    lon1, lat1 = (particle.lon + u1 * 0.5 * particle.dt, particle.lat + v1 * 0.5 * particle.dt)
    (u2, v2) = (fieldset.UWind[time + 0.5 * particle.dt, particle.depth, lat1, lon1, particle],
                fieldset.VWind[time + 0.5 * particle.dt, particle.depth, lat1, lon1, particle])
    lon2, lat2 = (particle.lon + u2 * 0.5 * particle.dt, particle.lat + v2 * 0.5 * particle.dt)
    (u3, v3) = (fieldset.UWind[time + 0.5 * particle.dt, particle.depth, lat2, lon2, particle],
                fieldset.VWind[time + 0.5 * particle.dt, particle.depth, lat2, lon2, particle])
    lon3, lat3 = (particle.lon + u3 * particle.dt, particle.lat + v3 * particle.dt)
    (u4, v4) = (fieldset.UWind[time + particle.dt, particle.depth, lat3, lon3, particle],
                fieldset.VWind[time + particle.dt, particle.depth, lat3, lon3, particle])
    particle_dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * particle.dt  # noqa
    particle_dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * particle.dt  # noqa

# delete-out-of-bound-partices kernel
def DeleteParticle(particle, fieldset, time):
    if particle.state >= 40: #== parcels.StatusCode.ErrorOutOfBounds:
        particle.delete()

# currents from https://doi.org/10.48670/moi-00021
nc_files = sorted(glob(os.path.join(PATH_DATA, "NetCDF", "Multi", "*.nc")))
fieldset = parcels.FieldSet.from_netcdf(nc_files,
                                                variables={"U": "uo",
                                                           "V" : "vo",
                                                           },
                                                dimensions={"time" : "time",
                                                            "lat" : "latitude",
                                                            "lon" : "longitude"
                                                            },
                                                )
# wind from https://doi.org/10.48670/moi-00185
nc_files = sorted(glob(os.path.join(PATH_DATA, "NetCDF", "Wind", "*.nc")))
fieldset_wind = parcels.FieldSet.from_netcdf(nc_files,
                                                variables={"UWind": "eastward_wind",
                                                           "VWind" : "northward_wind"},
                                                dimensions={"time" : "time",
                                                            "lat" : "latitude",
                                                            "lon" : "longitude"},
                                                )
# Add wind to fieldset
fieldset.add_field(fieldset_wind.UWind, name="UWind")
fieldset.add_field(fieldset_wind.VWind, name="VWind")

# Create the particle set
print(f"Initiating ParticleSet ...")
pset = parcels.ParticleSet.from_list(
            fieldset=fieldset,  # the fields on which the particles are advected
            pclass=parcels.JITParticle,  # the type of particles (JITParticle or ScipyParticle)
            lon=initial_locations["Longitude"],  # a vector of release longitudes
            lat=initial_locations["Latitude"],  # a vector of release latitudes
            time=initial_locations["DeviceDateTime"],  # a vector of release times
            id=initial_locations["id"],  # a vector of release times
        )
# Generating the output file in zarr format
output_file = pset.ParticleFile(
            name= os.path.join(PATH_DATA, "dual_projection14.zarr"),  # the file name
            outputdt=pd.Timedelta(hours=1),  # the time step of the outputs
        )

 # Run the advection
print(f"Advecting the particles ...")
pset.execute(
            [parcels.AdvectionRK4, AdvectionRK4_wind, DeleteParticle],  # the kernel (which defines how particles move)
            runtime=pd.Timedelta(days=RUNTIME),  # the total length of the run
            dt=pd.Timedelta(minutes=15),  # the timestep of the kernel, try to match the frequency of the drifters
            output_file=output_file,
        )

Output

Currents with AdvectionRK4

Winds with AdvectionRK4

Currents+Wind with AdvectionRK4+AdvectionRK4_wind

The dual kernels always project the particles North-East, which do not match the expected West direction. The only difference between wind and current fieldsets is the resolution (wind: 0.125°, currents: 0.083°). Could this be the problem? If it is, any possibilities of interpolating either current or wind to the other resolution?

Thanks again!

[erikvansebille]

Thanks for using Parcels, @AlxLhrNc. The issue that you're describing very likely is related to the under-the-hood unit conversion from m/s to degree/s. This works automatically for fields with names U and V, but you will need to tell Parcels explicitly what the units are if fields have other names (like UWind)

See the Unit Conversion tutorial for more information. Hope this helps!

[AlxLhrNc]

Thank you for the quick reply. However, it doesn't seem to be unit related, as wind and currents work as intended independently as stated in my post, and both share the same unit (m/s) as per their respective doi (currents, winds). My only guess would be a mistake in the kernels, either AdvectionRK4_wind or WindInteraction, but I can't figure out where.

If it helps, they both result in the exact same figure, with particles leaving after exactly 4 dt. This '4 dt exit' also happens (but westward) when I replace the wind_factor variable by a negative value (wind_factor = -0.03), see image there after.

[AlxLhrNc]

Many apologies, I realized I misunderstood your message after reading into the tutorial a bit more. However, I ran into an error that would be worth investigating.

Code from tutorial

fieldset_wind = parcels.FieldSet.from_netcdf(nc_files,
                                                variables={"U": "eastward_wind",
                                                           "V" : "northward_wind"},
                                                dimensions={"time" : "time",
                                                            "lat" : "latitude",
                                                            "lon" : "longitude"},
                                                )
    # Add wind to fieldset
    fieldset.add_field(fieldset_wind.U, name="UWind")
    fieldset.add_field(fieldset_wind.V, name="VWind")
    pset = parcels.ParticleSet.from_list(
            fieldset=fieldset,  # the fields on which the particles are advected
            pclass=parcels.JITParticle,  # the type of particles (JITParticle or ScipyParticle)
            lon=initial_locations["Longitude"],  # a vector of release longitudes
            lat=initial_locations["Latitude"],  # a vector of release latitudes
            time=initial_locations["DeviceDateTime"],  # a vector of release times
            id=initial_locations["id"],  # a vector of release times
        )

Error raised with my data

Exception has occurred: AssertionError
Field U.name (UWind) is not consistent
  File "C:\Users\ALN\Documents\Jobs\Oceanly\Scripts\drifters9_parcels_currents_winds.py", line 153, in <module>
    fieldset=fieldset,  # the fields on which the particles are advected
            pclass=parcels.JITParticle,  # the type of particles (JITParticle or ScipyParticle)
            lon=initial_locations["Longitude"],  # a vector of release longitudes
            lat=initial_locations["Latitude"],  # a vector of release latitudes
            time=initial_locations["DeviceDateTime"],  # a vector of release times
            id=initial_locations["id"],  # a vector of release times
        )
AssertionError: Field U.name (UWind) is not consisten

Fix

After checking, it seems that in my case, fieldset.add_field(fieldset_wind.U, name="UWind") is unable to rename the added field. I managed to fix it using the following code:

    fieldset_wind = parcels.FieldSet.from_netcdf(nc_files,
                                                variables={"UWind": "eastward_wind",
                                                           "VWind" : "northward_wind"},
                                                dimensions={"time" : "time",
                                                            "lat" : "latitude",
                                                            "lon" : "longitude"},
                                                )
    # Add wind to fieldset
    fieldset.add_field(fieldset_wind.UWind, name="UWind")
    fieldset.UWind.units = fieldset.U.units
    fieldset.add_field(fieldset_wind.VWind, name="VWind")
    fieldset.VWind.units = fieldset.V.units

It might be a very case-by-case error, but I find it strange.

Thanks a lot for your help.