Error accessing a user-defined field with an arbitrary_grid

[ngoldbaum]

Code for reproduction

import yt

def _tracerf(field, data):
    return data['Metal_Density']/data['gas', 'density']

ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")

ds.add_field(("gas", "tracerf"), function=_tracerf, units="dimensionless",
             take_log=False)

ad = ds.all_data()
galgas = ds.arbitrary_grid([0.4, 0.4, 0.4], [0.99, 0.99, 0.99],
                           dims=[512, 512, 512])
tracerp = galgas['tracerf']

The above script produces the following traceback:

Traceback (most recent call last):
  File "yt_arbgridtest.py", line 16, in <module>
    tracerp = galgas['tracerf']
  File "/Users/goldbaum/Documents/yt-git-fixes/yt/data_objects/data_containers.py", line 281, in __getitem__
    self.get_data(f)
  File "/Users/goldbaum/Documents/yt-git-fixes/yt/data_objects/construction_data_containers.py", line 641, in get_data
    if len(fill) > 0: self._fill_fields(fill)
  File "/Users/goldbaum/Documents/yt-git-fixes/yt/data_objects/construction_data_containers.py", line 851, in _fill_fields
    assert(len(fields) == 1)
AssertionError

[matthewturk]

It works for me if I change Metal_Density to density, so I think this might be an issue with the field definition having an error in it.

[jzuhone]

@MatthewTurk Isn't her issue for user derived fields though? That wouldn't catch it, then.

[matthewturk]

@jzuhone I mean in the field definition -- I still access the derived field.

[matthewturk]

Oddly, I am able to do galgas["Metal_Density"]/galgas["gas", "density"] so it's not as clear cut as I thought with the error. It's something more subtle, like a dependency not getting picked up during field dependency computation.

[stonnes]

Also, if you call for the user-defined field with a cut_region it works. So it is purely when you call for the user-defined field for the arbitrary_grid. Maybe this was already clear, but just in case!

[jzuhone]

FWIW, I tried your script by making a different user-defined field without Metal_Density, and that did work.

[stonnes]

Good point, I can also find user-defined fields that work. This one, with no Metal_Density call, also works for a cut region and not for the arbitrary grid--this is actually the first fail that clued me in to the problem (you can either take the absolute value as I do here or not divide and get the value with dimensions g/cm/s**2 and both fail):

def _boundval(field, data):
    Mdisk = 1.15e11 * Msun
    Mbulge = 1e10 * Msun
    a = 3.5e-3 * Mpc
    b = 7.0e-4 * Mpc
    rbulge = 6.0e-4 * Mpc
    rdm = 2.3e-2 * Mpc
    rhodm = 3.8e-25 * g/cm**3
    GravConst = 6.67e-8 *cm**3/g/s**2

    a = ds.arr(a,'Mpc').to('cm')
    b = ds.arr(b,'Mpc').to('cm')
    rbulge = ds.arr(rbulge,'Mpc').to('cm')
    rdm = ds.arr(rdm,'Mpc').to('cm')
    Mdisk = ds.arr(Mdisk,'Msun').to('g')
    Mbulge = ds.arr(Mbulge,'Msun').to('g')

    AMx = 0.0
    AMy = 0.342
    AMz = 0.93969

    xc = ds.arr(0.5,'code_length') #* code_length
    yc = ds.arr(0.5,'code_length') #* code_length
    zc = ds.arr(0.5,'code_length') #* code_length

    xd = data["x"].in_units('code_length') - xc
    yd = data["y"].in_units('code_length') - yc
    zd = data["z"].in_units('code_length') - zc

    xd = ds.arr(xd,'code_length').to('cm')
    yd = ds.arr(yd,'code_length').to('cm')
    zd = ds.arr(zd,'code_length').to('cm')

    zheight = AMx*xd + AMy*yd + AMz*zd
    xcyl = xd - zheight*AMx
    ycyl = yd - zheight*AMy
    zcyl = zd - zheight*AMz
    rsphere = np.sqrt(xd*xd + yd*yd + zd*zd)
    rcyl = np.sqrt(xcyl*xcyl+ycyl*ycyl+zcyl*zcyl)

    Epotbulge = -Mbulge/(rsphere + rbulge)
    Epotdisk = -Mdisk/(pow(rcyl*rcyl + pow(a + pow(zheight*zheight + b*b,0.5),2.0),0.5))
    Epotdm = -1.0*np.pi*rhodm*rdm*rdm*(np.pi - 2.0*(1.0+rdm/rsphere)* np.arctan(rsphere/rdm) + 2.0*(1.0+rdm/rsphere)*np.log(1.0+rsphere/rdm) - (1.0-rdm/rsphere)*np.log(1.0+(pow(rsphere/rdm,2))))
    Epottot = Epotbulge + Epotdisk + Epotdm
    return (GravConst*Epottot*data["gas","density"].in_cgs() + data["gas","thermal_energy"]*data["gas","density"].in_cgs() + data["gas","kinetic_energy"].in_cgs())/np.abs(GravConst*Epottot*data["gas","density"].in_cgs() + data["gas","thermal_energy"]*data["gas","density"].in_cgs() + data["gas","kinetic_energy"].in_cgs())

[ngoldbaum]

@stonnes I opened a PR, see #1530.