Ocean Glider Behavior (Active Particle) in ROMS output

[ChadDGibson]

Hello Community,

Any help with this error would be greatly appreciated...Anyone experienced similar results?
I am encountering a Value Error while running a virtual glider (active ScipyParticle) in a ROMS output with AdvectionRK4_3D.
I have made the glider behavior work in the Argo Tutorial flow field, and I can run a passive particle in the ROMS output. However, when I use AdvectionRK4_3D and GliderMovement Kernels together, I get the following error:

ValueError                                Traceback (most recent call last)
Cell In[38], [line 137](vscode-notebook-cell:?execution_count=38&line=137)
[130](vscode-notebook-cell:?execution_count=38&line=130) output_file = pset.ParticleFile(
[131](vscode-notebook-cell:?execution_count=38&line=131)     name="glider_float",
[132](vscode-notebook-cell:?execution_count=38&line=132)     outputdt=timedelta(minutes=10))#,
[133](vscode-notebook-cell:?execution_count=38&line=133)   #  chunks=(1, 500),  # setting to write in chunks of 500 observations ***Messes Up Color Plot*** Deleted
[134](vscode-notebook-cell:?execution_count=38&line=134) #)
[135](vscode-notebook-cell:?execution_count=38&line=135) 
[136](vscode-notebook-cell:?execution_count=38&line=136) # Now execute the kernels for 30 days, saving data every 30 minutes
--> [137](vscode-notebook-cell:?execution_count=38&line=137) pset.execute(
[138](vscode-notebook-cell:?execution_count=38&line=138)     kernels,
[139](vscode-notebook-cell:?execution_count=38&line=139)     runtime=timedelta(hours=12),
[140](vscode-notebook-cell:?execution_count=38&line=140)     dt=timedelta(minutes=5),
[141](vscode-notebook-cell:?execution_count=38&line=141)     output_file=output_file,
[142](vscode-notebook-cell:?execution_count=38&line=142) )

File ~/miniforge/envs/parcels/lib/python3.12/site-packages/parcels/particleset.py:948, in ParticleSet.execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, output_file, verbose_progress, postIterationCallbacks, callbackdt)
[946](https://file+.vscode-resource.vscode-cdn.net/Users/chadman/Documents/Projects/Velo_Recon/~/miniforge/envs/parcels/lib/python3.12/site-packages/parcels/particleset.py:946) # If we don't perform interaction, only execute the normal kernel efficiently.
[947](https://file+.vscode-resource.vscode-cdn.net/Users/chadman/Documents/Projects/Velo_Recon/~/miniforge/envs/parcels/lib/python3.12/site-packages/parcels/particleset.py:947) if self.interaction_kernel is None:
--> [948](https://file+.vscode-resource.vscode-cdn.net/Users/chadman/Documents/Projects/Velo_Recon/~/miniforge/envs/parcels/lib/python3.12/site-packages/parcels/particleset.py:948)     res = self.kernel.execute(self, endtime=next_time, dt=dt)
[949](https://file+.vscode-resource.vscode-cdn.net/Users/chadman/Documents/Projects/Velo_Recon/~/miniforge/envs/parcels/lib/python3.12/site-packages/parcels/particleset.py:949)     if res == StatusCode.StopAllExecution:
[950](https://file+.vscode-resource.vscode-cdn.net/Users/chadman/Documents/Projects/Velo_Recon/~/miniforge/envs/parcels/lib/python3.12/site-packages/parcels/particleset.py:950)         return StatusCode.StopAllExecution

File ~/miniforge/envs/parcels/lib/python3.12/site-packages/parcels/kernel.py:564, in Kernel.execute(self, pset, endtime, dt)
...
[631](https://file+.vscode-resource.vscode-cdn.net/Users/chadman/Documents/Projects/Velo_Recon/~/miniforge/envs/parcels/lib/python3.12/site-packages/parcels/kernel.py:631)     if sign_dt*p.time < sign_dt*endtime and p.state == StatusCode.Success:

File <ast>:42, in SetcoordsGliderMovementAdvectionRK4_3DUpdatecoords(particle, fieldset, time)

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

MY KERNEL Definition

Define the new Kernel that mimics Argo vertical movement

def GliderMovement(particle, fieldset, time):
driftdepth = 500 # maximum depth in m
maxdepth = 1000 # maximum depth in m
vertical_speed = 0.1 # sink and rise speed in m/s

uvel = .00000025 
vvel = 0
cycletime = 1 * 86400  # total time of cycle in seconds
drifttime = 5 * 60  # time of deep drift in seconds

# glider_u = glider_U * np.sin(heading)
# glider_v = glider_U * np.cos(heading)

if particle.cycle_phase == 0:
    # Phase 0: Sinking with vertical_speed until depth is driftdepth
    particle_ddepth += vertical_speed * particle.dt
    particle_dlon += uvel * particle.dt
    particle_dlat += vvel * particle.dt
    if particle.depth >= driftdepth:
        particle.cycle_phase = 1

elif particle.cycle_phase == 1:
    # Phase 1: Drifting at depth for drifttime seconds
    # particle_ddepth += particle.dt
    particle.drift_age += particle.dt
    if particle.drift_age >= .2 * drifttime:
        particle.drift_age = 0  # reset drift_age for next cycle
        particle.cycle_phase = 3

elif particle.cycle_phase == 2:
    # Phase 2: Sinking further to maxdepth
    particle_ddepth += vertical_speed * particle.dt
    if particle.depth >= maxdepth:
        particle.cycle_phase = 3

elif particle.cycle_phase == 3:
    # Phase 3: Rising with vertical_speed until at surface
    particle_ddepth -= vertical_speed * particle.dt
    particle_dlon += uvel * particle.dt
    particle_dlat += vvel * particle.dt
    # particle.temp = fieldset.temp[time, particle.depth, particle.lat, particle.lon]  # if fieldset has temperature
    if particle.depth <= fieldset.mindepth:
        particle_ddepth = fieldset.mindepth
        # particle.temp = 0./0.  # reset temperature to NaN at end of sampling cycle
        particle.cycle_phase = 4

elif particle.cycle_phase == 4:
    # Phase 4: Transmitting at surface until cycletime is reached
    particle.drift_age += particle.dt
    if particle.drift_age >= drifttime:
        particle.drift_age = 0  # reset drift_age for next cycle
        particle.cycle_phase = 0
        # particle.cycle_age = 0

    # if particle.cycle_age > cycletime:
    #     particle.cycle_phase = 0
    #     particle.cycle_age = 0

if particle.state == StatusCode.Evaluate:
    particle.cycle_age += particle.dt  # update cycle_age

Execution: ``

from datetime import timedelta

import numpy as np

from parcels import (
    AdvectionRK4,
    AdvectionRK4_3D,
    FieldSet,
    JITParticle,
    ScipyParticle,
    ParticleSet,
    StatusCode,
    Variable,
    download_example_dataset,
    # plotTrajectoriesFiles,
)


filenames = {
    "U": {
        "lon": f"{dataset_folder}/modified_g*.nc",
        "lat": f"{dataset_folder}/modified_g*.nc",
        "depth": f"{dataset_folder}/modified_g*.nc",
        "data": f"{dataset_folder}/modified_d*.nc",
    },
    "V": {
        "lon": f"{dataset_folder}/modified_g*.nc",
        "lat": f"{dataset_folder}/modified_g*.nc",
        "depth": f"{dataset_folder}/modified_g*.nc",
        "data": f"{dataset_folder}/modified_d*.nc",
    },
    "W": {
        "lon": f"{dataset_folder}/modified_g*.nc",
        "lat": f"{dataset_folder}/modified_g*.nc",
        "depth": f"{dataset_folder}/modified_g*.nc",
        "data": f"{dataset_folder}/modified_d*.nc",
    },
    "T": {
        "lon": f"{dataset_folder}/modified_g*.nc",
        "lat": f"{dataset_folder}/modified_g*.nc",
        "depth": f"{dataset_folder}/modified_g*.nc",
        "data": f"{dataset_folder}/modified_d*.nc",
    },
}
variables = {
    "U": "u",
    "V": "v",
    "W": "w",
    "T": "temp"
    # "time": "ocean_time",
}
dimensions = {"lon": "lon_u", "lat": "lat_v", "depth": "z_rho", "time": "ocean_time"}
fieldset = FieldSet.from_c_grid_dataset(filenames, variables, dimensions)#, allow_time_extrapolation=True)
# fieldset = FieldSet.from_nemo(filenames, variables, dimensions)#, allow_time_extrapolation=True)
# uppermost layer in the hydrodynamic data
fieldset.mindepth = fieldset.U.depth[10]


# # Define a new Particle type including extra Variables
class GliderParticle(ScipyParticle):
    # Phase of cycle:
    # init_descend=0,
    # drift=1,
    # profile_descend=2,
    # profile_ascend=3,
    # transmit=4
    cycle_phase = Variable("cycle_phase", dtype=np.int32, initial=0)
    cycle_age = Variable("cycle_age", dtype=np.float32, initial=0.0)
    drift_age = Variable("drift_age", dtype=np.float32, initial=0.0)
    # if fieldset has temperature
    # Temp = Variable('T', dtype=np.float32, initial=np.nan)


# Initiate one Argo float in the Agulhas Current
# pset = ParticleSet(
#     fieldset=fieldset, pclass=ArgoParticle, lon=[128.5], lat=[20], depth=[100]
# )

# Initiate one Argo float in the Agulhas Current
pset = ParticleSet(
    fieldset=fieldset, pclass=GliderParticle, lon=[128.5], lat=[20], depth=[100]
)


# combine Argo vertical movement kernel with built-in Advection kernel
kernels = [GliderMovement, AdvectionRK4_3D] #Turned this off and it runs???
# kernels = [AdvectionRK4_3D, GliderMovement]
# kernels = [AdvectionRK4_3D, ArgoVerticalMovement]
# Create a ParticleFile object to store the output
output_file = pset.ParticleFile(
    name="glider_float",
    outputdt=timedelta(minutes=10),
    chunks=(1, 500),  # setting to write in chunks of 500 observations ***Messes Up Color Plot*** Deleted
)

# Now execute the kernels for 30 days, saving data every 30 minutes
pset.execute(
    kernels,
    runtime=timedelta(days=1),
    dt=timedelta(minutes=5),
    output_file=output_file,
)

[erikvansebille]

Hi @ChadDGibson, thanks for reporting. I have not seen this error before; not really sure where it comes from.

But before I dig deeper, can you check whether you're running with the latest version of parcels (v3.0.2)? And if not, could you upgrade to that version? We've made some changes in the kernel so it could be that this fixes your issue already...

[ChadDGibson]

Hi @erikvansebille, thank you for getting back to me so quickly.
I was running v3.0.0. Just upgraded to v3.0.2 and have the same error message.
I am also running Python 3.12.0.

[erikvansebille]

Hmm, very strange then. Your Glider Kernel is (almost?) identical to the Argo Kernel. What happens if you run the actual Argo example? Do you get the same error?

[ChadDGibson]

I have run the ArgoMovement Kernel in this flow field. Same/Similar Error. When running as a JITParticle I get a "not a scalar argument error", when run as ScipyParticle, I get the original error as shown above.
Perhaps it is in the flow field I am using, and the way datetime64 is used (1 hour time-steps).
My next step is to use the Nemo Tutorial flow field, and add the ArgoMovement/GliderMovement Kernels there.
Does that seem to be a logical next step/test?

[erikvansebille]

Yes, indeed running with the default Nemo flow field would be a next logical step. If the code works there, then I agree it's likely to be something in your FieldSet. Let me know how it goes!