OceanParcels / Parcels

Main code for Parcels (Probably A Really Computationally Efficient Lagrangian Simulator)
https://www.oceanparcels.org
MIT License
294 stars 135 forks source link

BrownianMotion2D kernel issue 'AttributeError: 'NoneType' object has no attribute 'flags'' #798

Closed HaydenSchilling closed 4 years ago

HaydenSchilling commented 4 years ago

Hi, I'm pretty new to Python and I'm trying to run some simulations but am having trouble getting the BrownianMotion2D kernel to work. The simulations run fine without this kernel but with this kernal it gives the following error:


  File "C:\Users\hayde\Crab-Dispersal\Crab_Dispersal.py", line 164, in <module>
    endtime=end_time)

  File "C:\Users\hayde\.conda\envs\py3_parcels\lib\site-packages\parcels\particleset.py", line 527, in execute
    self.kernel.execute(self, endtime=time, dt=dt, recovery=recovery, output_file=output_file)

  File "C:\Users\hayde\.conda\envs\py3_parcels\lib\site-packages\parcels\kernel.py", line 341, in execute
    self.execute_jit(pset, endtime, dt)

  File "C:\Users\hayde\.conda\envs\py3_parcels\lib\site-packages\parcels\kernel.py", line 250, in execute_jit
    fargs = [byref(f.ctypes_struct) for f in self.field_args.values()]

  File "C:\Users\hayde\.conda\envs\py3_parcels\lib\site-packages\parcels\kernel.py", line 250, in <listcomp>
    fargs = [byref(f.ctypes_struct) for f in self.field_args.values()]

  File "C:\Users\hayde\.conda\envs\py3_parcels\lib\site-packages\parcels\field.py", line 964, in ctypes_struct
    if not self.data_chunks[i].flags.c_contiguous:

AttributeError: 'NoneType' object has no attribute 'flags'

I believe I'm using a 3D NEMO style grid. (model detailed here: https://doi.org/10.1029/2019JC015227 ). My input parameters are:

# Set diffusion constants.
Kh_zonal = 100
Kh_meridional = 100

fieldset = FieldSet.from_nemo(filenames, variables, dimensions, indices, allow_time_extrapolation=True)#, transpose=True)

# Create field of Kh_zonal and Kh_meridional, using same grid as U
#[time, depth, particle.lon, particle.lat] # Think this order is correct for here
size4D = (30,30,fieldset.U.grid.ydim, fieldset.U.grid.xdim)
fieldset.add_field(Field('Kh_zonal', Kh_zonal*np.ones(size4D), grid=fieldset.temp.grid))
fieldset.add_field(Field('Kh_meridional', Kh_meridional*np.ones(size4D), grid=fieldset.temp.grid))

Thanks for any help,

erikvansebille commented 4 years ago

Thanks @HaydenSchilling, for reporting this Issue. The problem is that you're using the Grid of fieldset.temp, which is based on a netcdf file to create 'simple' numpy Field. That does not work, because of all kinds of reasons with what we call 'deferred loading' in Parcels. I have just pushed an update that will trigger a clearer error message in the future: #797.

In the meantime, you could fix your code by using:

fieldset.add_field(Field('Kh_zonal', Kh_zonal*np.ones(size4D), lon=fieldset.temp.grid.lon, 
                                          lat=fieldset.temp.grid.lat, depth=fieldset.temp.grid.depth,
                                          time=fieldset.temp.grid.time))

Also note that in general you only need a 2D array (lon,lat) for these diffusion fields, because Parcels can work without for depth and time dimensions. So

size2D = (fieldset.U.grid.ydim, fieldset.U.grid.xdim)
fieldset.add_field(Field('Kh_zonal', Kh_zonal*np.ones(size2D), lon=fieldset.temp.grid.lon, lat=fieldset.temp.grid.lat))

And once #784 is fixed, you'd could do with a simple scalar value

HaydenSchilling commented 4 years ago

Thanks so much @erikvansebille! The 2D method makes sense and works! But the first suggestion keeping the 4D size does not (FYI), it gives an error code:

ValueError: cannot reshape array of size 77072400 into shape (30,316,271)

It doesn't appear to be having issues the 4D size, with both time and depth (which I don't think should be necessary, as you pointed out). I think it's trying to put too many points into too small an array.

DEHewitt commented 4 years ago

Hi @erikvansebille, I am working with Hayden on the simulations that the above issue relates to. Thanks so much for the help so far.

A quick follow up question: would the fix suggested (i.e. the '2D method') have any implications for the units of our diffusion constants? In subsequent simulations we have had to have them extremely small (1x10^-7), compared to previous work of Haydens (~100).

Thanks for any help you can offer!

erikvansebille commented 4 years ago

Ah good point! You'll probably have to add mesh='spherical' to the creation of the new Field. So

size2D = (fieldset.U.grid.ydim, fieldset.U.grid.xdim)
fieldset.add_field(Field('Kh_zonal', Kh_zonal*np.ones(size2D), lon=fieldset.temp.grid.lon, 
                         lat=fieldset.temp.grid.lat, mesh='spherical'))

In that way, Parcels knows that the units of lon and lat are in degrees, and it will automagically convert the 'Kh' from m2/s to deg2/s. This wasn't needed before, when you used grid=fieldset.temp.grid, because that grid already had the information of a spherical grid.

Hope this helps!

DEHewitt commented 4 years ago

Thanks @erikvansebille! Seems to have solved it!

DEHewitt commented 4 years ago

Hi @erikvansebille, following on from the issue above - we are now attempting to nest models and I think I have done that correctly, however when I attempt to run the code from your last comment I get the following error: AttributeError: 'NestedField' object has no attribute 'grid'. Any help you can offer would be great!

DEHewitt commented 4 years ago

Just adding the lines which nest our models: U = NestedField('U', [fieldset_ROMS.U, fieldset_BRAN.U]) V = NestedField('V', [fieldset_ROMS.V, fieldset_BRAN.V]) fieldset = FieldSet(U, V)

erikvansebille commented 4 years ago

Hi @DEHewitt. The AttributeError is because fieldset.U is now a vector, so you should do something like

size2D_0 = (fieldset.U[0].grid.ydim, fieldset.U[0].grid.xdim)
fieldset.add_field(Field('Kh_zonal', Kh_zonal*np.ones(size2D_0), lon=fieldset.U[0].grid.lon, 
                          lat=fieldset.U[0].grid.lat, mesh='spherical'))

Note that this means that you use the same 'Kh_zonal' in the entire domain.

If you want different Kh_zonal values in the different nests, you'll have to do something like (untested)

size2D_0 = (fieldset.U[0].grid.ydim, fieldset.U[0].grid.xdim)
size2D_1 = (fieldset.U[1].grid.ydim, fieldset.U[1].grid.xdim)
fieldset.add_field([Field('Kh_zonal', Kh_zonal_0*np.ones(size2D_0), lon=fieldset.U[0].grid.lon, 
                          lat=fieldset.U[0].grid.lat, mesh='spherical'),
                    Field('Kh_zonal', Kh_zonal_1*np.ones(size2D_1), lon=fieldset.U[1].grid.lon, 
                          lat=fieldset.U[1].grid.lat, mesh='spherical')])
DEHewitt commented 4 years ago

Thanks @erikvansebille! That seems to have worked. I haven't tested the second method you suggest, but will notify here if we do and it works.