Hi Nathan,

Thanks for your prompt response.

Earlier I've included the ghost zones, but the script gives me "unexpected keyword argument" error. And I just looked into the source code, and found the keyword is "ghost_zones" instead of "ghost_zone" appearing in the docs example. After this correction, the script works.

Thanks again!

All best,
Suoqing


On Jan 23, 2013, at 9:39 PM, Nathan Goldbaum <nathan12343@gmail.com> wrote:

Hi Suoqing,

Others on the list know more about this than me, but I want to note that you do need to include the ghost zones.  Whenever yt needs to difference a field, it will generate ghost zones by doing a cascading interpolation.  The ghost zones do not need to be written to disk.

There might be small differences between the ghost zones yt generates and the ones that were actually used during the calculation, so it's possible there will be some small artifacts in the final result.  Usually the ghost zones generated by yt will give excellent results.

-Nathan

On Jan 23, 2013, at 6:31 PM, Ji Suoqing wrote:

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
_______________________________________________
yt-users mailing list
yt-users@lists.spacepope.org
http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org

_______________________________________________
yt-users mailing list
yt-users@lists.spacepope.org
http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org