[yt-users] Define a field of spacial derivatives for FLASH4 2D cylindrical uniform grid

[Ji Suoqing]

Hi,

I would like to compute some spacial derivatives for FLASH4 2D cylindrical uniform grid datafile with YT. 

For simplicity, I tried the example given in YT docs, although I know the derivative form is only valid for Cartesian, and this is just for code test purposes. My script are:

@derived_field(name='DivV')  # I do not use ValidateSpatial… ghost_zone… because the ghost cells are not stored in FLSAH check point file

def DivV(field, data):

        sl_left = slice(None,-2)
        sl_right = slice(2,None)
        div_fac = 2.0

        ds = div_fac * data['dx'].flat[0]
        f  = data["velx"][sl_right,1:-1]/ds
        f -= data["vely"][sl_left ,1:-1]/ds
        
        ds = div_fac * data['dy'].flat[0]
        f += data["velx"][1:-1,sl_right]/ds
        f -= data["vely"][1:-1,sl_left]/ds
        
        new_field = na.zeros(data["velx"].shape, dtype='float64')
        new_field[1:-1,1:-1] = f
        
        return new_field

Then, if I add:

pf = load('relax_hdf5_chk_0001')
pf.h.find_max('DivV')

The script works, but if I use:

dd = pf.h.all_data()
print dd['DivV'].sum()

and it will give error:

yt : [INFO     ] 2013-01-23 21:17:01,551 Getting field dx from 256
yt : [INFO     ] 2013-01-23 21:17:01,625 Getting field velx from 256
Traceback (most recent call last):
  File "test.py", line 28, in <module>
    print dd['DivV'].sum()
  File "/usr/local/yt-i386/src/yt-hg/yt/data_objects/data_containers.py", line 330, in __getitem__
    self.get_data(key)
  File "/usr/local/yt-i386/src/yt-hg/yt/data_objects/data_containers.py", line 2607, in get_data
    if self._generate_field(field):
  File "/usr/local/yt-i386/src/yt-hg/yt/data_objects/data_containers.py", line 360, in _generate_field
    self[field] = self.pf.field_info[field](self)
  File "/usr/local/yt-i386/src/yt-hg/yt/data_objects/field_info_container.py", line 378, in __call__
    dd = self._function(self, data)
  File "test.py", line 11, in DivV
    f  = data["velx"][sl_right,1:-1]/ds
IndexError: too many indices

Nor can I make slice plot for this field.

I would appreciate your suggestions for this issue.

Best wishes,
Suoqing
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20130123/422d3a1a/attachment.htm>