Nesting - particles stop moving when leaving the fine resolution grid

[cgrimaldi]

Question

Question

Hi
I am trying to use the nesting function in ocean parcels, but can't seem to get it to work. I do not get any error, however when I look at my particle tracks output, it stops moving when it leaves my finer grid resolution model. I am assuming I am doing something wrong in the nesting side of things, but have followed the tutorial pretty closely. See code below. Any help is appreciated thanks :)

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import pandas
from datetime import timedelta as delta
from datetime import datetime as datetime

from parcels import AdvectionRK4, Field, FieldSet, JITParticle, NestedField, ParticleSet

def DeleteParticle(particle, fieldset, time):
    particle.delete()


# LOAD COARSER RESOLUTION GRID - ROMS
fname = 'ROMS.nc'
filenames = {'U': fname, 'V': fname}
variables = {'U': 'u_eastward', 'V': 'v_northward'}
dimensions = {'U': { 'lon': 'lon_rho_utm','lat': 'lat_rho_utm','time': 'ocean_time'},
              'V': { 'lon': 'lon_rho_utm','lat': 'lat_rho_utm','time': 'ocean_time'}}
fieldset_Large = FieldSet.from_netcdf(filenames, variables, dimensions, mesh="flat",allow_time_extrapolation = True)

# PLOT COARSER RESOLUTION GRID - ROMS
fieldset_Large.computeTimeChunk()
plt.pcolormesh(fieldset_Large.U.grid.lon, fieldset_Large.U.grid.lat, fieldset_Large.U.data[0,:,:])
plt.xlabel("Northings[m]")
plt.ylabel("Eastings[m]")
plt.colorbar()
plt.show()

# READ FINER RESOLUTION GRID - DFM
fname='DFM.nc'
filenames = {'U': fname, 'V': fname}
variables = {'U': 'U', 'V': 'V'}
dimensions = {'U': {'lon': 'x','lat': 'y','time': 'time'},'V': {'lon': 'x','lat': 'y','time': 'time'}}
fieldset_Small = FieldSet.from_netcdf(filenames, variables, dimensions, mesh="flat",allow_time_extrapolation = True)

# PLOT FINER RESOLUTION GRID - DFM
fieldset_Small.computeTimeChunk()
plt.pcolormesh(fieldset_Small.U.grid.lon, fieldset_Small.U.grid.lat, fieldset_Small.U.data[1,:, :])
plt.xlabel("Northings[m]")
plt.ylabel("Eastings[m]")
plt.colorbar()
plt.show()

# CREATE NESTED FIELD
U = NestedField('U', [fieldset_Small.U,fieldset_Large.U])
V = NestedField('V', [fieldset_Small.V,fieldset_Large.V])
fieldset_nest = FieldSet(U, V)

# RUN PTRACKING
pset = ParticleSet(fieldset_nest, pclass=JITParticle,lon=[200000], lat=[7590000])

output_file = pset.ParticleFile(
    name="NestedFieldParticle_test.zarr", outputdt=delta(days=1))

pset.execute(AdvectionRK4, runtime=delta(days=30), dt=delta(days=1), output_file=output_file)

[VeckoTheGecko]

Hi Camille!

There doesn't seem to be anything wrong with the code that you are working with. Would you be able to attach the plotted trajectories as well?

Cheers,
Nick

[erikvansebille]

Oh sorry, my mistake. You should not use the np.sqrt(), as numpy functions can't be used in Kernels. I edited the proposed kernel above

[cgrimaldi]

Oh right - easy fix thanks @erikvansebille. From what I understand, here you are just creating a kernel to check when the particle is in the small domain, and if it is, you pop out an out of bond error, but how do I actually make it use the larger grid velocities when it leaves the small grid?

I was trying to do something like that, based on my code higher up, does it make sense?

def CheckInSmallDomain(particle, fieldset, time):
    u,v = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
    if (u**2+v**2) < 1e-12:
        (u,v) = fieldset_Large.UV[particle]
    else:
        (u,v) = fieldset_Small.UV[particle]

% RUN PTRACKING
pset = ParticleSet(fieldset_Small, pclass=JITParticle,lon=[200000], lat=[7630000])
output_file = pset.ParticleFile(name="NestedFieldParticle_test.zarr", outputdt=delta(days=1))
pset.execute([AdvectionRK4,CheckInSmallDomain], runtime=delta(days=30), dt=delta(days=1), output_file=output_file)

But getting this error: NotImplementedError: Cannot convert 'fieldset_Large' used in kernel to C-code

[erikvansebille]

Well, that is the point of the NestedFields; under the hood Parcels, will go to the next Field in the list of Fields when an OutofBounds error occurs; so no need to do anything with fieldset_large and fieldset_small etc. I hope that with just my proposed CheckInSmallDomain() kernel above, your nesting will work as you implement it at the top

[cgrimaldi]

I was going from the suggestion that @michaeldenes made about going around the nested field.
When I use your suggestion, it just raises an out of bond error when the particle leaves the smaller domain..

def CheckInSmallDomain(particle, fieldset, time):
u,v = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
if (u2+v2) < 1e-12:
particle.state = StatusCode.ErrorOutOfBounds

[cgrimaldi]

@VeckoTheGecko it seems like if I manually modify the NaN with a random missing value in the netcdf file, the nested field works. I will let a couple of simulations run and report back tomorrow. Hopefully that works, thanks for the NaN hint!
Camille