Particles out of bounds in idealized curvilinear domain

[wingtorres]

I am trying to use Parcels for offline Lagrangian particle integration in fields derived from ROMS (coupled to SWAN via COAWST) simulations of wave-driven flow on an idealized island. The annulus grid is curvilinear, laterally periodic, and irregularly-spaced (shown below) on an f-plane at latitude -30°.

The velocity field in polar coordinates after 14 days of constant radially symmetric wave-forcing on the outer boundary looks like

(sorry for the oddly formatted figure - I'm still getting familiar with hvplot)

Anyway, for some reason Parcels has problems integrating particle tracks on some parts of the domain, but not others. Here are trajectories for single particles seeded in each reef pass jet integrated for 12 hours just using the surface (u,v) fields (particle scatter color is mapped to time).

Based on print output from my recovery kernel, Parcels seems to think the particles seeded in most of the jets are not in bounds from the get go. One thing that has me a bit worried is that when I plot a field using fieldset.V.show(), only one section of the domain is displayed (although particle integration seems to work for some locations outside of this slice and not others), even though fieldset.V.lon and fieldset.V.lat clearly have the correct shape and values for the full annulus.

Below is the code I use to load the model output, create the particleset+fieldset, and do the integration calling "mpirun -np x python driftersIntegrate.py". I've also provided a Dropbox link to the ROMS netcdf file in case that is helpful. ncks was used to extract just the top sigma level and every 6th time slice to reduce file size. I haven't tried to implement the periodic boundary condition yet because I want to sort this issue out first, but hopefully that will come next. Thank you very much in advance - Parcels is an awesome code!

https://www.dropbox.com/s/hmzgyo7d4mzi5lk/anntest_strided.nc?dl=0

Import Modules

import os
import sys
import xarray as xr
import netCDF4 as nc
import numpy as np
import cmocean
from datetime import timedelta as delta
from parcels import AdvectionRK4, ErrorCode, FieldSet, JITParticle, ParticleFile, ParticleSet, Variable, plotting

Load data

home = os.environ['HOME']
filename = home + '/Desktop/anntest_strided2D.nc'
Ds = xr.open_dataset(filename)
lon = Ds['lon_psi'].values
lat = Ds['lat_psi'].values
angle_psi = 0.25*( Ds.angle.isel(eta_rho = slice(0, -1), xi_rho = slice(0, -1) ) \
         + Ds.angle.isel(eta_rho = slice(1, None), xi_rho = slice(1, None) ) \
         + Ds.angle.isel(eta_rho = slice(1, None), xi_rho = slice(0, -1) ) \
         + Ds.angle.isel(eta_rho = slice(0,-1), xi_rho = slice(1,None)) ).rename({'xi_rho': 'xi_psi', 'eta_rho': 'eta_psi'}) #to psi points

Interpolate Arakawa C-grid velocity components to psi/f points, calculate Lagrangian velocity, then rotate to eastward/northward orientation

#velocity to psi points
u_eulerian = 0.5*( Ds['u'].isel(eta_u = slice(0,-1)) + Ds['u'].isel(eta_u = slice(1,None)) ).rename({'xi_u': 'xi_psi', 'eta_u': 'eta_psi'}) 
v_eulerian = 0.5*( Ds['v'].isel(xi_v = slice(0,-1)) + Ds['v'].isel(xi_v = slice(1,None)) ).rename({'xi_v': 'xi_psi', 'eta_v': 'eta_psi'}) 
u_stokes = 0.5*( Ds['u_stokes'].isel(eta_u = slice(0,-1)) + Ds['u_stokes'].isel(eta_u = slice(1,None)) ).rename({'xi_u': 'xi_psi', 'eta_u': 'eta_psi'}) 
v_stokes = 0.5*( Ds['v_stokes'].isel(xi_v = slice(0,-1)) + Ds['v_stokes'].isel(xi_v = slice(1,None)) ).rename({'xi_v': 'xi_psi', 'eta_v': 'eta_psi'}) 
u = u_eulerian + u_stokes
v = v_eulerian + v_stokes
uveitheta = (u + 1j*v)*np.exp(1j*angle_psi) #rotate to eastward and northward
u_eastward  = np.real(uveitheta.isel(s_rho = -1)) #top slice
v_northward = np.imag(uveitheta.isel(s_rho = -1)) #top slice

Initialize fieldset, particleset, and recovery+advection kernel

#set fieldset
data = {'U': u_eastward.values, 'V': v_northward.values} #2D velocity
dimensions = {'lon': lon, 'lat': lat, 'time': t }
fieldset = FieldSet.from_data(data, dimensions, transpose = False, mesh = 'spherical', time_periodic = t[-1] - t[0])
fieldset.V.show()

#recovery kernel
def DeleteParticle(particle, fieldset, time):
    print("Deleting particle")
    particle.delete()

recovery = {ErrorCode.ErrorOutOfBounds: DeleteParticle,
            ErrorCode.ErrorThroughSurface: DeleteParticle}

#initialize particleset
theta = np.linspace(0, 2*np.pi, 12, endpoint = False) + np.pi/12
lonp = .1*np.cos(theta)
latp = .1*np.sin(theta)
pset = ParticleSet.from_list(fieldset = fieldset, pclass = JITParticle, lon = lonp.tolist(), lat = latp.tolist() )
kernels = AdvectionRK4
output_file = pset.ParticleFile(name= "PeriodicParticle", outputdt = delta(seconds = 300) )

Integrate for 12 hours and save output

pset.execute( kernels, runtime = delta(hours = 12.0), dt = delta(seconds = 60), output_file = output_file, recovery = recovery )
output_file.export()
output_file.close()

[erikvansebille]

Hi @wingtorres, thanks for reporting this Issue. I am travelling at the moment with poor internet, so can't download your data file to test myself. I will have a try later next week. However, it may very well be that your model setup is too complicated for Parcels, and that making it work on such a torus requires significant changes under the hood. So don't get your hopes too high that we can make this work. But I'll have a look next week

[wingtorres]

Hey @erikvansebille thanks so much for being willing to look into it! Maybe there are circumpolar applications of Parcels we could reference?

[erikvansebille]

Dear @wingtorres. I now had a look at your data file, and as I feared the topology of this simulation is really quite a challenge to Parcels. When we developed the curvilinear capabilities, we just didn't expect these extremely curved grids. The point is that Parcels expects monotonically increasing lot and lat grids, which is not the case in your simulation

If you really want to keep using Parcels, the easiest solution probably is to interpolate your fields onto a regular grid; although regridding is never ideal because it introduces errors in the interpolation. Do you have another Lagrangian code that does work on this topology?

[wingtorres]

Hi @erikvansebille thank you for investigating, and it is totally understandable that Parcels was not designed with this kind of case in mind. Might another solution be to just re-generate the grid as if it were at a pole? That way lon/lat would be monotically increasing and I could then conceivably use zonally periodic BC's. The central reference point is arbitrary because it is an idealized domain, and I am manually setting it to an f-plane anyway.

ROMS' online particle integration does work with this domain, but enabling floats grievously slowed down the computation so I ended up not going that route. Below is a snapshot of 5000 online particle tracks initially seeded evenly throughout the lagoon after three days time (Coriolis parameter set to 0).

[erikvansebille]

Wow, this figure looks amazingly beautiful! Would be great to check if Parcels can produce similar figures, and also be efficient. Indeed, moving the pole to the centre of the torus could be a very good strategy. Let us know if it works!

[wingtorres]

hey @erikvansebille good news - moving the grid to the North pole worked! I'm not sure if ROMS coordinates are referenced on a geoid or spherical earth, but a spherical transformation seemed to work fine at the very least for visualization purposes.

r = 6378e3 #radius of earth at equator. does ROMS account for geoid?
z = (r**2 - x**2 - y**2)**0.5
lon, lat = (180/np.pi)*np.arctan2(y,x) , 90-(180/np.pi)*np.arccos( z/r )

then I followed the Periodic Boundary tutorial to create the zonal halo and periodic BC kernel.

Regarding the figure... I made that figure in Blender by reading in particle tracks through its Python API and animating them as curves with halo/particle system effects, but I've found that the Python package Datashader can produce really nice graphics as well with much less of a learning curve . Datashader seems really well equipped to handle tons of data points even with minimal computing resources (i.e. just a laptop), so it may be of interest to you.

Here are the results for 2^14 particles seeded initially in the backreef+lagoon on the annulus/torus problem (this time w/ coriolis turned on) and integrated for 48 hours. The color here should be proportional to the density of pathlines over the integration period. I think this nicely illustrates how some particles that leave via the wave-driven jets come back over the reef crest and recirculate, but many continue on and join neighboring jets instead. This was done just with a static slice of output for testing purposes though - I suspect the full offline integration will look more chaotic.

import os 
import xarray as xr
import numpy as np
import hvplot.xarray
import hvplot.pandas
import datashader as ds
import datashader.transfer_functions as tf
import holoviews as hv
import cmocean
hv.extension('bokeh')
renderer = hv.renderer('bokeh')

particlename = "~/Dropbox/Oceanography/Projects/annulus/floats/annFloats.nc"
Dp = xr.open_dataset(particlename, decode_times = False)
df = Dp.to_dataframe() 
df.loc[[0], ['lon','lat']] = np.nan #add intermittent row with nans (first observation) to separate lines
df = df.sort_values(['trajectory', 'obs']) #organize by trajectory then order by time

cvs = ds.Canvas(plot_width = 1200, plot_height = 1200)
agg = cvs.line(dp, 'lon', y = 'lat', agg = ds.count())
img = tf.shade(agg, cmap = cmocean.cm.ice, how = "eq_hist")
plot = tf.set_background(tf.dynspread(img, threshold = 0.25, max_px = 4), 'black')
plot

[erikvansebille]

Thanks @wingtorres, this is really great. Beautiful visualization! I’m ping-ing @CKehl, who is also very interested in visualizations. Great that you’ve been able to make this work with a relatively simple coordinate transformation

[CKehl]

Hi - yeah, the renderings are aesthetically quite pleasing! I also know the bokeh rendering engine - there are actually quite some InfoVis frameworks build on top of it. Thanks for the forward of 'Datashader' - that seems to deserve some more attention indeed. From what I see, the coriolis force is the driving factor in the simulation, with the toroidal drift of the particles. The darkerning mid-bad around the torus looks a bit 'artificial' - would be interesting to see what happened with the particles there, or if that's just a visualisation artefact. Thanks for sharing your results, and giving some inspiration to those advanced plotting libraries!

[wingtorres]

@CKehl thank you for your comment! The darkening mid band I think you are referring is inside of the wave break where the reef is very shallow and flows are fast so the circulation isn't as much affected by Coriolis and the particles initially seeded there are quickly flushed out through the jets. So I believe the visualization is reflecting the physics accurately, but an animation would likely help to verify this.

As an aside I'm really impressed by the efficiency of both Parcels and Datashader - both the integration and visualization were done in <2 min compute time on just my personal laptop.