Cycle_age updating faster than time step when trying to replicate the ArgoVerticalMovement example

[hrsdawson]

Question

Question

I'm trying to simulate some virtual Argo Floats and have starting by replicating the ArgoVerticalMovement example. However I'm running into some problems:

1. The cycle_age seems to be updated multiple times within a time step during the descent and ascent phases.

2. The float depth doesn't get updated to the drift depth at the end of Phase 0 when it enters the if particle.depth + particle_ddepth >= drift depth: loop. I can only get the float to my desired drifting depth if I set the depth at the start of the Phase 1 loop. However the float seems to skip a depth when I do this. In my example the float moves 90m in one time step (15 mins). So the float should descend from 181 m to 271 m to 361 m at which point is should be reset to 300 m (the drift depth). However in the output the float moves straight from 181 m to 300 m. Does this mean that Parcels is not executing any depth changes at the end of the Phase 1 loop and if so, why?

I've only recently updated to Parcels version 3 (it's been a while since I've run simulations!) so it's possible I've missed some key updates/understanding of what's changed. This code was run using version 3.1.0. Any pointers are much appreciated, thank you.

Supporting code/error messages

from parcels import FieldSet, ParticleSet, Variable, JITParticle, ScipyParticle, AdvectionRK4

output_directory = '/scratch/e14/hd4873/particle_tracking/virtual_floats/logistically_feasible/'
mesh_mask = output_directory + '/ocean.nc'

# Define the files containing the U and V files (2D advection).
data_path = '/g/data/cj50/access-om2/raw-output/access-om2-01/01deg_jra55v140_iaf/'
ufiles = sorted(glob(data_path + f'output*/ocean/ocean-3d-u-1-daily-mean-ym_*.nc'))
vfiles = sorted(glob(data_path + f'output*/ocean/ocean-3d-v-1-daily-mean-ym_*.nc'))

filenames = {
    'U': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': mesh_mask, 'data': ufiles},
    'V': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': mesh_mask, 'data': vfiles},
    }
 
variables = {'U': 'u', 'V': 'v'}

dimensions = {'U': {'lon': 'xu_ocean', 'lat': 'yu_ocean', 'depth': 'sw_ocean', 'time': 'time'},
              'V': {'lon': 'xu_ocean', 'lat': 'yu_ocean', 'depth': 'sw_ocean', 'time': 'time'},
             }
 
cs = {"U": {"lon": ("xu_ocean", 400), "lat": ("yu_ocean", 300), "depth": ("st_ocean", 2), "time": ("time", 1)}, 
      "V": {"lon": ("xu_ocean", 400), "lat": ("yu_ocean", 300), "depth": ("st_ocean", 2), "time": ("time", 1)},
     }

fieldset = FieldSet.from_mom5(filenames, variables, dimensions, 
                              mesh = 'spherical', 
                              chunksize = cs, 
                              tracer_interp_method = 'bgrid_tracer',
                             )

# Set a minimum depth == the uppermost layer in the hydrodynamic data
fieldset.mindepth = fieldset.U.depth[0]
print(fieldset.mindepth) 
# this is 1.08 m

# Define a new Particle type including extra Variables
ArgoParticle = parcels.JITParticle.add_variables(
    [
        parcels.Variable("cycle_phase", dtype=np.int32, initial=0.0),
        parcels.Variable("cycle_age", dtype=np.float32, initial=0.0),
        parcels.Variable("drift_age", dtype=np.float32, initial=0.0),
    ]
)

def ArgoVerticalMovement(particle, fieldset, time):
    driftdepth = 300       # drifting depth in m
    vertical_speed = 0.1   # sink and rise rate in m/s 
    cycletime = 4 * 86400  # Just using 4 days for now
    drifttime = cycletime - (5*60*60)  # given the shallow drifting depth, just allowing for 5 hour rise/rall rate for the purposes of this test

    # Phase 0: sinking to drift depth with vertical_speed
    if particle.cycle_phase == 0:
        particle_ddepth += vertical_speed * particle.dt
        if particle.depth + particle_ddepth >= driftdepth:
            # Note: This code did not work, so updating particle depth directly in next loop
            particle_ddepth = driftdepth - particle.depth
            particle.cycle_phase = 1

    # Phase 1: Drifting at depth for drifttime seconds
    elif particle.cycle_phase == 1:
        # updates to depth at end of cycle_phase 0 don't seem to be working, so updating directly here instead
        particle.depth = driftdepth 
        particle_ddepth = 0
        particle.drift_age += particle.dt
        if particle.drift_age >= drifttime:
            particle.drift_age = 0  # reset drift_age for next cycle
            particle.cycle_phase = 3

    # Skipping Phase 2 (second descent phase) for testing purposes

    # Phase 3: Rising with vertical_speed until at surface
    elif particle.cycle_phase == 3:
        particle_ddepth -= vertical_speed * particle.dt
        if particle.depth + particle_ddepth <= fieldset.mindepth:
            particle_ddepth = fieldset.mindepth - particle.depth 
            particle.cycle_phase = 4

    # Phase 4: Transmitting at surface until cycletime is reached
    elif particle.cycle_phase == 4:
        # updates to depth at end of cycle_phase 3 don't seem to be working, so updating directly here instead
        particle.depth = fieldset.mindepth 
        particle_ddepth = 0       
        if particle.cycle_age > cycletime:
            particle.cycle_phase = 0
            particle.cycle_age = 0

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


lons = np.array(77.)
lats = np.array(-68.24674827)
depths = np.array(1.0825615)

pset = ParticleSet.from_list(fieldset=fieldset,   
                             pclass=ArgoParticle,  
                             lon=lons,          
                             lat=lats,           
                             depth=depths,       
                            ) 



output_filename = 'TestArgoParticle_MinimalExample.zarr'

# Save to file at same timestep as advection kernels
output_file = pset.ParticleFile(
    name=output_directory + output_filename, 
    outputdt=timedelta(minutes=15), 
    chunks=(1, 500), 
)

kernels = [ArgoVerticalMovement, AdvectionRK4]

pset.execute(
    kernels,
    runtime=timedelta(days=7),    # short run time for example
    dt=timedelta(minutes=15),     
    output_file=output_file,      
)

Then, reading in the output:

ds = xr.open_zarr(output_directory + "TestArgoParticle_MinimalExample.zarr")


# Problem 1: During the descent phase the cycle_age does not increment at the same rate as the timestep: 
# Cycle age:

# Cycle_age difference between timesteps:

# Problem 2: When descending during the first cycle the floats moves from 181 m to 300 m. There should be another recorded depth at 271 m before reaching the 300 m drifting depth. The same thing happens when the float ascends (moves from 121 m straight to 1 m)

ds.z.values

[erikvansebille]

Hi @hrsdawson, thanks for reaching out with this question!

It could be that the version of the Argo kernel that you're using is a bit outdated. I remember that we had to make some updates to the Kernel when we were developing our Virtual Ship Classroom.

Could you check the code at https://github.com/OceanParcels/virtualship/blob/main/src/virtualship/instruments/argo_float.py and see if that does work?

If so, we may need to update the example here in the Parcels repo too

[hrsdawson]

Hi @erikvansebille, thanks for the quick response. I've tried the code you linked and unfortunately that doesn't work for me either. It increased the cycle_age increment during the cycle 3 phase due to this code (particle.cycle_age += (particle.dt) # solve issue of not updating cycle_age during ascent. It also didn't solve the drifting depth problem. In the first cycle my float will descend only to 270m before ascending to the surface. In the second cycle it descends to 300m but then only ascends to 30m. In the cycles thereafter it descends to 210m and then ascends to 30m.

Since my problem is the cycle_age incrementing too fast, I changed the code to subtract the timestep ((particle.cycle_age -= (particle.dt)) but that led to zero increment in the cycle_age during ascent rather than the correct (or closer to correct) increment. However! Using this negative age increment partially solved the drifting depth problem I just described. The float still doesn't make it to the desired drifting depth of 300m but it does consistently descend to 270m and ascend to the surface in all sampling cycles. I have no idea why that change results in a more consistent depth cycling pattern....

Note: I had to switch the cycle 1 if statement in the virtual ship code from if particle.depth + particle_ddepth <= particle.drift_depth: to if particle.depth + particle_ddepth >= particle.drift_depth: for my particular case otherwise the float just stayed at the surface. Is this an error in the virtual ship argo float kernel, or is the depth direction configured differently for other Parcels fieldset methods? I'm using the from_mom5 method.

[erikvansebille]

Hi @hrsdawson, frustrating that you still have the issue with the depths. Could it be that the virtualship code indeed expects depths to be negative (since it's an educational tool, we want to be strict in our right-hand cartesian coordinate system). We use Grid.negate_depth for that, see e.g.

Parcels/tests/test_grids.py

Lines 117 to 124 in f3a987e

	def test_negate_depth():
	depth = np.linspace(0, 5, 10, dtype=np.float32)
	fieldset = FieldSet.from_data(
	{"U": np.zeros((10, 1, 1)), "V": np.zeros((10, 1, 1))}, {"lon": [0], "lat": [0], "depth": depth}
	)
	assert np.all(fieldset.gridset.grids[0].depth == depth)
	fieldset.U.grid.negate_depth()
	assert np.all(fieldset.gridset.grids[0].depth == -depth)

Can you try whether it works if you also negate the depth in your simulation?

[hrsdawson]

Hmm, I tried negating the depths but got this error: AttributeError: 'RectilinearZGrid' object has no attribute 'negate_depth'. I guess I can't use that with the hydrodynamic model I have?

[erikvansebille]

Hi @hrsdawson, if you don't have a negate_depth function, it means you have a slightly older version of Parcels installed. No problem, you can do the same with

fieldset.U.grid._depth = -fieldset.U.grid._depth

Or, if that doesn't work (again, depending on your installed version of Parcels)

fieldset.U.grid.depth = -fieldset.U.grid.depth

[hrsdawson]

Okay I've now tried negating the depth. That did not fix the problem. Same issues arise as when depth is positive (inconsistent cycle age updates during descent/ascent and float consistently not reaching desired drifting depth).

[erikvansebille]

OK, in that case perhaps it's easiest if I try to debug your code myself. Can you send me a minimal reproducible example, also with the hydrodynamic data you're using? Can be directly to my UU email address, if that's easier

[hrsdawson]

Thanks @erikvansebille. I just sent you some hydrodynamic files via WeTransfer, which I'd subsetted to a small geographic domain (too big to send otherwise). Intriguingly, my issues disappeared when using that subsetted field. After some digging I've discovered the problem was the chunk size I was giving Parcels for the depth dimension of the hydrodynamic fields. I was originally using a chunk size of 2 which causes the cycle age and depth problems to occur, but if I increase the chunk size to something much larger (e.g. 38 - half the length of the depth dimension) then the cycle_age and depth issues disappear.

Is this a bug, or should I have known better re the chunk sizes? I've just sent you an example script via email that illustrates the problem when using small chunk in depth, but if it's not a bug and you don't wish to dig further, then the problem is simply solved by increasing the depth chunk size!

[erikvansebille]

Thanks @hrsdawson, for this response. Quite surprising that the chunk-size settings had such an impact on your kenrel; I haven't seen that before... Glad that you now found the issue!

As we're currently focussing all our development efforts on Parcels v4, and we remove the fieldset-chunking options then in favour of native xarray-support, I don't think we have the time and capacity to delve deeper into this bug for now. But good to keep in mind that chunking can lead to unexpected kernel-behaviour....