[yt-users] Error slicing 2D spherical FLASH data

Matthew Turk matthewturk at gmail.com
Tue Jun 24 06:22:23 PDT 2014


Hi all,

Here's a pull request that should address the issue.

https://bitbucket.org/yt_analysis/yt/pull-request/974/fixing-non-cartesian-flash-bounding-box

Please leave a comment if it looks good to you.

-Matt

On Mon, Jun 23, 2014 at 9:08 PM, Matthew Turk <matthewturk at gmail.com> wrote:
> On Mon, Jun 23, 2014 at 7:59 PM, Luke Zoltan Kelley
> <lkelley at cfa.harvard.edu> wrote:
>> I'm very confused.  When I load this file (that I uploaded earlier), I get
>>
>> In [2]: pf = load('starwind_hdf5_plt_cnt_0000')
>> yt : [WARNING  ] 2014-06-23 20:55:08,536 integer runtime parameter checkpointfilenumber overwrites a simulation scalar of the same name
>> yt : [INFO     ] 2014-06-23 20:55:08,556 Parameters: current_time              = 0.0
>> yt : [INFO     ] 2014-06-23 20:55:08,556 Parameters: domain_dimensions         = [ 8 16  1]
>> yt : [INFO     ] 2014-06-23 20:55:08,557 Parameters: domain_left_edge          = [ 0.  0.  0.]
>> yt : [INFO     ] 2014-06-23 20:55:08,557 Parameters: domain_right_edge         = [  4.00000000e+11   3.14159265e+00   1.74532925e-02]
>> yt : [INFO     ] 2014-06-23 20:55:08,557 Parameters: cosmological_simulation   = 0.0
>>
>> so it seems to think it is 2D (which is *correct*), and the axis-0 (rad) and axis-1 (theta) bounds seem correct.  The phi bound is pi/180 for some reason (but this can be overwritten in the FLASH parameter file).  I tried manually resetting the pf.domain_right_edge to be [4e11, pi, pi] but I still get the exact same error --- I don't understand Matt's fix.
>>
>
> Ah, I misread your earlier email and thought it was supposed to be 1D.
>
> Now, that being said, I think yt is actually doing what John suggested
> after I made one minor change.  The thing I'm a bit confused by is how
> it's setting the left/right edges of the grids inside yt.  What seems
> to be right for Cartesian isn't right for non-Cartesian.
>
> John, the tricky logic seems to all be inside where we use ND and rdx
> to lock against integer multiples.  I'm going to try again to take a
> look at this tomorrow (have to go tonight) but basically everything
> seems to be in place, there's just this confusion with blocks, etc,
> and the domain_dimensions not adding up to the domain_width.
>
> -Matt
>
>>
>>
>> On Jun 23, 2014, at 6:38 PM, John ZuHone <jzuhone at gmail.com> wrote:
>>
>>> AFAIK this is partially due to undesirable behavior in the way that FLASH writes bounding boxes in spherical coordinates when one or more dimension is suppressed. Buried deep within the code somewhere (which I cannot remember but I have run into before) there is logic that sets the bounding box along the suppressed dimensions to pi/nblock[xyz], where nblock[xyz] is the corresponding value of the number of blocks on a side for that dimension, even though the value of nblock[xyz] is completely irrelevant.
>>>
>>> My suggestion would be that within yt we should re-set the bounding boxes in 1D/2D datasets to be 0 -> pi or 0 -> 2*pi for theta and phi. This will also help us get cell_volume right.
>>>
>>> Are you sure that Luke's simulation is 1D? It can't be if the blocks are 8x8x1.
>>>
>>> On Jun 23, 2014, at 6:17 PM, Matthew Turk <matthewturk at gmail.com> wrote:
>>>
>>>> The NaNs are on purpose, to blank out the regions that hadn't been
>>>> seen by the pixelizer.
>>>>
>>>> I grabbed Luke's data and played around with it, and I think I figured
>>>> out the problem.  The dimensionality shows up as 2, although Luke ran
>>>> a 1D simulation.  The blocks are all 8x8x1, too.
>>>>
>>>> Now, if I set the right edge of the grids in axis one to be pi, I get
>>>> a plot like the one attached.
>>>>
>>>> If we can figure out in what situations dimensionality will be
>>>> reported as 2, but the left/right edges arent' set correctly, we can
>>>> get this working out of the box.  Ideas, Luke, John, or anyone else?
>>>>
>>>> On Mon, Jun 23, 2014 at 4:46 PM, Nathan Goldbaum <nathan12343 at gmail.com> wrote:
>>>>> The issue seems to be that the pixelizer is returning NaNs:
>>>>>
>>>>> ipdb> up
>>>>>>
>>>>>> /Users/goldbaum/Documents/yt-hg/yt/visualization/base_plot_types.py(116)_init_image()
>>>>>   115                                       aspect=aspect, vmax=self.zmax,
>>>>> cmap=cmap)
>>>>> --> 116         self.cb = self.figure.colorbar(self.image, self.cax)
>>>>>   117
>>>>>
>>>>> ipdb> l
>>>>>   111             norm = matplotlib.colors.Normalize()
>>>>>   112         extent = [float(e) for e in extent]
>>>>>   113         self.image = self.axes.imshow(data.to_ndarray(),
>>>>> origin='lower',
>>>>>   114                                       extent=extent, norm=norm,
>>>>> vmin=self.zmin,
>>>>>   115                                       aspect=aspect, vmax=self.zmax,
>>>>> cmap=cmap)
>>>>> --> 116         self.cb = self.figure.colorbar(self.image, self.cax)
>>>>>   117
>>>>>   118     def _repr_png_(self):
>>>>>   119         canvas = FigureCanvasAgg(self.figure)
>>>>>   120         f = StringIO()
>>>>>   121         canvas.print_figure(f)
>>>>>
>>>>> ipdb> print data
>>>>> [[ nan  nan  nan ...,  nan  nan  nan]
>>>>> [ nan  nan  nan ...,  nan  nan  nan]
>>>>> [ nan  nan  nan ...,  nan  nan  nan]
>>>>> ...,
>>>>> [ nan  nan  nan ...,  nan  nan  nan]
>>>>> [ nan  nan  nan ...,  nan  nan  nan]
>>>>> [ nan  nan  nan ...,  nan  nan  nan]] g/cm**3
>>>>> ipdb> print data.max()
>>>>> nan g/cm**3
>>>>> ipdb> print data.min()
>>>>> nan g/cm**3
>>>>> ipdb> print np.any(np.isfinite(data))
>>>>> False
>>>>>
>>>>> Unfortunately I'm not sure why this is happening.  FWIW, Matt knows a lot
>>>>> more than I do about how the pixelizer works for non-cartesian geometries.
>>>>>
>>>>>
>>>>> On Mon, Jun 23, 2014 at 2:15 PM, Luke Zoltan Kelley
>>>>> <lkelley at cfa.harvard.edu> wrote:
>>>>>>
>>>>>> Sure, thanks Nathan!
>>>>>> I've uploaded the file here
>>>>>>
>>>>>>
>>>>>> https://drive.google.com/file/d/0B4FXIXq2e3y4eVpHMHNPMVJCLVU/edit?usp=sharing
>>>>>>
>>>>>> Luke
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> On Jun 23, 2014, at 5:08 PM, Nathan Goldbaum <nathan12343 at gmail.com>
>>>>>> wrote:
>>>>>>
>>>>>> For some reason it's falling over when it's trying to construct the
>>>>>> colorbar for the plot.  I'm not sure why and would need to look closer at
>>>>>> the data file.  If you'd be willing to share it with me I'm happy to do
>>>>>> that.
>>>>>>
>>>>>> If you want to debug this on your own my first suggestion would be to run
>>>>>> the script under gdb and try to figure out why the colorbar construction is
>>>>>> going wonky.
>>>>>>
>>>>>> -Nathan
>>>>>>
>>>>>>
>>>>>> On Mon, Jun 23, 2014 at 2:02 PM, Luke Zoltan Kelley
>>>>>> <lkelley at cfa.harvard.edu> wrote:
>>>>>>>
>>>>>>> They seem fine:
>>>>>>>
>>>>>>> In [33]: dens = pf.all_data()['density']
>>>>>>>
>>>>>>> In [34]: print dens.min()
>>>>>>> 9.99999968266e-20 g/cm**3
>>>>>>>
>>>>>>> In [35]: print dens.max()
>>>>>>> 157.337768555 g/cm**3
>>>>>>>
>>>>>>> These are the values I expected.
>>>>>>>
>>>>>>>
>>>>>>> On Jun 23, 2014, at 4:47 PM, Nathan Goldbaum <nathan12343 at gmail.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>> Hey Luke,
>>>>>>>
>>>>>>> What do the raw density values in your data look like?
>>>>>>>
>>>>>>> dens = ds.all_data()['density']
>>>>>>> print dens.max()
>>>>>>> print dens.min()
>>>>>>>
>>>>>>> -Nathan
>>>>>>>
>>>>>>>
>>>>>>> On Mon, Jun 23, 2014 at 1:37 PM, Luke Zoltan Kelley
>>>>>>> <lkelley at cfa.harvard.edu> wrote:
>>>>>>>>
>>>>>>>> Hello again yt-users,
>>>>>>>>
>>>>>>>> I'm trying to plot density from a FLASH simulation in 2D spherical.  I
>>>>>>>> assumed that a 'SlicePlot' would be the natural choice, slicing along axis
>>>>>>>> '2'/'z'/'phi'.
>>>>>>>>
>>>>>>>> pf = load('starwind_hdf5_plt_cnt_0000')
>>>>>>>> slc = SlicePlot(pf, 2, "density")
>>>>>>>>
>>>>>>>> However, I get the error:
>>>>>>>>
>>>>>>>> yt : [INFO     ] 2014-06-23 16:31:42,631 xlim = -400000000000.000000
>>>>>>>> 400000000000.000000
>>>>>>>> yt : [INFO     ] 2014-06-23 16:31:42,632 ylim = -400000000000.000000
>>>>>>>> 400000000000.000000
>>>>>>>> yt : [INFO     ] 2014-06-23 16:31:42,632 Making a fixed resolution
>>>>>>>> buffer of (('gas', 'density')) 800 by 800
>>>>>>>> yt : [INFO     ] 2014-06-23 16:31:42,698 xlim = -400000000000.000000
>>>>>>>> 400000000000.000000
>>>>>>>> yt : [INFO     ] 2014-06-23 16:31:42,698 ylim = -400000000000.000000
>>>>>>>> 400000000000.000000
>>>>>>>> yt : [INFO     ] 2014-06-23 16:31:42,699 Making a fixed resolution
>>>>>>>> buffer of (('gas', 'density')) 800 by 800
>>>>>>>> yt : [INFO     ] 2014-06-23 16:31:42,703 Making a fixed resolution
>>>>>>>> buffer of (('gas', 'density')) 800 by 800
>>>>>>>> ERROR: ValueError: The truth value of an array with more than one
>>>>>>>> element is ambiguous. Use a.any() or a.all() [numpy.ma.core]
>>>>>>>> astropy: [ERROR    ] 2014-06-23 16:31:42,864 ValueError: The truth value
>>>>>>>> of an array with more than one element is ambiguous. Use a.any() or a.all()
>>>>>>>>
>>>>>>>> ---------------------------------------------------------------------------
>>>>>>>> ValueError                                Traceback (most recent call
>>>>>>>> last)
>>>>>>>> <ipython-input-26-62c5de6a0872> in <module>()
>>>>>>>> ----> 1 slc = SlicePlot(pf, 2, "density")
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Applications/yt/yt-x86_64/src/yt-hg/yt/visualization/plot_window.pyc
>>>>>>>> in SlicePlot(pf, normal, fields, axis, *args, **kwargs)
>>>>>>>>  1872             del kwargs['north_vector']
>>>>>>>>  1873
>>>>>>>> -> 1874         return AxisAlignedSlicePlot(pf, normal, fields, *args,
>>>>>>>> **kwargs)
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Applications/yt/yt-x86_64/src/yt-hg/yt/visualization/plot_window.pyc
>>>>>>>> in __init__(self, pf, axis, fields, center, width, axes_unit, origin,
>>>>>>>> fontsize, field_parameters, window_size, aspect)
>>>>>>>>  1093         if axes_unit is None:
>>>>>>>>  1094             axes_unit = get_axes_unit(width, pf)
>>>>>>>> -> 1095         self.set_axes_unit(axes_unit)
>>>>>>>>  1096
>>>>>>>>  1097 class ProjectionPlot(PWViewerMPL):
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Applications/yt/yt-x86_64/src/yt-hg/yt/visualization/plot_container.pyc
>>>>>>>> in newfunc(*args, **kwargs)
>>>>>>>>    67         rv = f(*args, **kwargs)
>>>>>>>>    68         args[0]._plot_valid = False
>>>>>>>> ---> 69         args[0]._setup_plots()
>>>>>>>>    70         return rv
>>>>>>>>    71     return newfunc
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Applications/yt/yt-x86_64/src/yt-hg/yt/visualization/plot_window.pyc
>>>>>>>> in _setup_plots(self)
>>>>>>>>   836                 self._colormaps[f], extent, zlim,
>>>>>>>>   837                 self.figure_size, fp.get_size(),
>>>>>>>> --> 838                 self.aspect, fig, axes, cax)
>>>>>>>>   839
>>>>>>>>   840             axes_unit_labels = ['', '']
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Applications/yt/yt-x86_64/src/yt-hg/yt/visualization/plot_window.pyc
>>>>>>>> in __init__(self, data, cbname, cmap, extent, zlim, figure_size, fontsize,
>>>>>>>> unit_aspect, figure, axes, cax)
>>>>>>>>  1703             size, axrect, caxrect, zlim, figure, axes, cax)
>>>>>>>>  1704
>>>>>>>> -> 1705         self._init_image(data, cbname, cmap, extent,
>>>>>>>> unit_aspect)
>>>>>>>>  1706
>>>>>>>>  1707         self.image.axes.ticklabel_format(scilimits=(-2, 3))
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Applications/yt/yt-x86_64/src/yt-hg/yt/visualization/base_plot_types.pyc
>>>>>>>> in _init_image(self, data, cbnorm, cmap, extent, aspect)
>>>>>>>>   114                                       extent=extent, norm=norm,
>>>>>>>> vmin=self.zmin,
>>>>>>>>   115                                       aspect=aspect,
>>>>>>>> vmax=self.zmax, cmap=cmap)
>>>>>>>> --> 116         self.cb = self.figure.colorbar(self.image, self.cax)
>>>>>>>>   117
>>>>>>>>   118     def _repr_png_(self):
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/figure.pyc
>>>>>>>> in colorbar(self, mappable, cax, ax, use_gridspec, **kw)
>>>>>>>>  1449                 cax, kw = cbar.make_axes(ax, **kw)
>>>>>>>>  1450         cax.hold(True)
>>>>>>>> -> 1451         cb = cbar.colorbar_factory(cax, mappable, **kw)
>>>>>>>>  1452
>>>>>>>>  1453         self.sca(current_ax)
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/colorbar.pyc
>>>>>>>> in colorbar_factory(cax, mappable, **kwargs)
>>>>>>>>  1272         cb = ColorbarPatch(cax, mappable, **kwargs)
>>>>>>>>  1273     else:
>>>>>>>> -> 1274         cb = Colorbar(cax, mappable, **kwargs)
>>>>>>>>  1275
>>>>>>>>  1276     mappable.callbacksSM.connect('changed',
>>>>>>>> cb.on_mappable_changed)
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/colorbar.pyc
>>>>>>>> in __init__(self, ax, mappable, **kw)
>>>>>>>>   875                 kw['alpha'] = mappable.get_alpha()
>>>>>>>>   876
>>>>>>>> --> 877             ColorbarBase.__init__(self, ax, **kw)
>>>>>>>>   878
>>>>>>>>   879     def on_mappable_changed(self, mappable):
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/colorbar.pyc
>>>>>>>> in __init__(self, ax, cmap, norm, alpha, values, boundaries, orientation,
>>>>>>>> ticklocation, extend, spacing, ticks, format, drawedges, filled, extendfrac,
>>>>>>>> extendrect, label)
>>>>>>>>   315         # The rest is in a method so we can recalculate when
>>>>>>>> clim changes.
>>>>>>>>   316         self.config_axis()
>>>>>>>> --> 317         self.draw_all()
>>>>>>>>   318
>>>>>>>>   319     def _extend_lower(self):
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/colorbar.pyc
>>>>>>>> in draw_all(self)
>>>>>>>>   336         and do all the drawing.
>>>>>>>>   337         '''
>>>>>>>> --> 338         self._process_values()
>>>>>>>>   339         self._find_range()
>>>>>>>>   340         X, Y = self._mesh()
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/colorbar.pyc
>>>>>>>> in _process_values(self, b)
>>>>>>>>   650                 self.norm.vmin = 0
>>>>>>>>   651                 self.norm.vmax = 1
>>>>>>>> --> 652             b = self.norm.inverse(self._uniform_y(self.cmap.N +
>>>>>>>> 1))
>>>>>>>>   653             if self._extend_lower():
>>>>>>>>   654                 b[0] = b[0] - 1
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/colors.pyc
>>>>>>>> in inverse(self, value)
>>>>>>>>  1001         if cbook.iterable(value):
>>>>>>>>  1002             val = ma.asarray(value)
>>>>>>>> -> 1003             return vmin * ma.power((vmax / vmin), val)
>>>>>>>>  1004         else:
>>>>>>>>  1005             return vmin * pow((vmax / vmin), value)
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/ma/core.pyc
>>>>>>>> in __mul__(self, other)
>>>>>>>>  3705     def __mul__(self, other):
>>>>>>>>  3706         "Multiply other by self, and return a new masked array."
>>>>>>>> -> 3707         return multiply(self, other)
>>>>>>>>  3708     #
>>>>>>>>  3709     def __rmul__(self, other):
>>>>>>>>
>>>>>>>>
>>>>>>>> /Users/lzkelley/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/ma/core.pyc
>>>>>>>> in __call__(self, a, b, *args, **kwargs)
>>>>>>>>   936         # Case 1. : scalar
>>>>>>>>   937         if not result.ndim:
>>>>>>>> --> 938             if m:
>>>>>>>>   939                 return masked
>>>>>>>>   940             return result
>>>>>>>>
>>>>>>>> ValueError: The truth value of an array with more than one element is
>>>>>>>> ambiguous. Use a.any() or a.all()
>>>>>>>>
>>>>>>>>
>>>>>>>> yt instinfo
>>>>>>>> Version = 3.0-dev
>>>>>>>> Changeset = 90a724768665
>>>>>>>>
>>>>>>>>
>>>>>>>> Any help would be greatly appreciated,
>>>>>>>> Thanks so much!
>>>>>>>> Luke
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> yt-users mailing list
>>>>>>>> yt-users at lists.spacepope.org
>>>>>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> yt-users mailing list
>>>>>>> yt-users at lists.spacepope.org
>>>>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> yt-users mailing list
>>>>>>> yt-users at lists.spacepope.org
>>>>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> yt-users mailing list
>>>>>> yt-users at lists.spacepope.org
>>>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>>>>
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> yt-users mailing list
>>>>>> yt-users at lists.spacepope.org
>>>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>>>>
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> yt-users mailing list
>>>>> yt-users at lists.spacepope.org
>>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>>>
>>>> <starwind_hdf5_plt_cnt_0000_Slice_phi_dens.png>_______________________________________________
>>>> yt-users mailing list
>>>> yt-users at lists.spacepope.org
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>
>>> _______________________________________________
>>> yt-users mailing list
>>> yt-users at lists.spacepope.org
>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>
>> _______________________________________________
>> yt-users mailing list
>> yt-users at lists.spacepope.org
>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org



More information about the yt-users mailing list