[yt-dev] Clump finding in 2D

Thu Jul 4 10:02:26 PDT 2013

Hi John,

On Wed, Jul 3, 2013 at 6:39 PM, John Forbes <jforbes at ucolick.org> wrote:
> Hi folks,
>
> I'm attempting to run yt's clump finding algorithm on a 2D FLASH dataset.
> Out of the box, following the basic example in the docs, I get the session
> at the bottom of this email, including an error mentioning an inability to
> get certain data from ghost zones.

Ah ... I think this will actually be somewhat problematic.  I am
thinking about the places this might be a problem, and I think the
biggest will be that ghost zone generation in 2D may not be possible.
The problem is likely the periodicity, which could wrap around
unsuccessfully.  I believe that in principle periodicity in the 2D
clump finding wouldn't be an issue if ghost zone generation functioned
properly for 2D datasets.

>
> I am able to get around this by commenting out the following lines in
> yt/analysis_modules/level_sets/contour_finder.py
>
> 105     for gi,grid in enumerate(grids):
> 106         pbar.update(gi)
> 107         cg = grid.retrieve_ghost_zones(1, "tempContours",
> smoothed=False)
> 108         grid_set.update(set(cg._grids))
> 109         fd = cg["tempContours"].astype('int64')
> 110         boundary_tree = amr_utils.construct_boundary_relationships(fd)
> 111         tree.update(((a, b) for a, b in boundary_tree))
>
> which I assume simply allows the code to ignore the ghost zones (which -
> wild speculation - maybe yt expects to exist in the 3rd dimension but
> don't?). My questions are:
> 1) Is this likely to have ill effects? (with this change find_clumps()
> appears to work, but I haven't inspected the results closely enough to know
> if they're reasonable).

The lines that are commented out above, unfortunately, also include
all of the joining across grid datasets.  So for instance, you'll be
able to identify clumps in a single patch, but not any that cross
patches.  I think the general solution here is to figure out how to
make ghost zone generation in 2D work properly.  I may have time to
work on this next week -- is there a chance that you have a small
dataset that I could test it out on?  And, could you file a ticket?

https://bitbucket.org/yt_analysis/yt/issues/new

-Matt

> 2) Is there a different recommended fix?
>
> Thank you,
> John
>
>
>
>
>
>
> In [1]: from yt.analysis_modules.level_sets.api import *
>
> In [2]: field="Density"
>
> In [3]: step=2.0
>
> In [4]: data_source=pf.h.all_data()
> yt : [WARNING  ] 2013-07-03 15:22:06,608 Dimensionality less than 3:
> reverting to overlap_proj
>
> In [5]: cmin = 10**np.floor(np.log10(data_source[field]).min())
> yt : [INFO     ] 2013-07-03 15:22:24,519 Getting field dens from 7872
>
> In [6]: cmax = 10**np.floor(np.log10(data_source[field]).max()+1)
>
> In [7]: function='self.data["Density"].size>20'
>
> In [8]: master_clump = Clump(data_source,None,field,function=function)
>
> In [9]: find_clumps(master_clump,cmin,cmax,step)
> Finding clumps: min: 1.000000e-28, max: 1.000000e-24, step: 2.000000
> yt : [ERROR    ] 2013-07-03 15:23:23,911 Covering problem: 190 cells are
> uncovered
> ---------------------------------------------------------------------------
> KeyError                                  Traceback (most recent call last)
> /home/jforbes/yt-x86_64/src/yt-hg/yt/mods.pyc in <module>()
> ----> 1 find_clumps(master_clump,cmin,cmax,step)
>
> /home/jforbes/yt-x86_64/src/yt-hg/yt/analysis_modules/level_sets/clump_handling.pyc
> in find_clumps(clump, min_val, max_val, d_clump)
>     177     print "Finding clumps: min: %e, max: %e, step: %f" % (min_val,
> max_val, d_clump)
>     178     if min_val >= max_val: return
> --> 179     clump.find_children(min_val)
>     180
>     181     if (len(clump.children) == 1):
>
> /home/jforbes/yt-x86_64/src/yt-hg/yt/analysis_modules/level_sets/clump_handling.pyc
> in find_children(self, min_val, max_val)
>     118         if max_val is None: max_val = self.max_val
>     119         contour_info = identify_contours(self.data, self.field,
> min_val, max_val,
> --> 120                                          self.cached_fields)
>     121         for cid in contour_info:
>     122             new_clump = self.data.extract_region(contour_info[cid])
>
> /home/jforbes/yt-x86_64/src/yt-hg/yt/analysis_modules/level_sets/contour_finder.py
> in identify_contours(data_source, field, min_val, max_val, cached_fields)
>     105     for gi,grid in enumerate(grids):
>     106         pbar.update(gi)
> --> 107         cg = grid.retrieve_ghost_zones(1, "tempContours",
> smoothed=False)
>     108         grid_set.update(set(cg._grids))
>     109         fd = cg["tempContours"].astype('int64')
>
> /home/jforbes/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.pyc in
> retrieve_ghost_zones(self, n_zones, fields, all_levels, smoothed)
>     477                 level, new_left_edge, **kwargs)
>     478         else:
> --> 479             cube = self.hierarchy.covering_grid(level,
> new_left_edge, **kwargs)
>     480
>     481         return cube
>
> /home/jforbes/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in
> __init__(self, level, left_edge, dims, fields, pf, num_ghost_zones,
> use_pbar, **kwargs)
>    3712         self.domain_width = np.rint((self.pf.domain_right_edge -
>    3713
> self.pf.domain_left_edge)/self.dds).astype('int64')
> -> 3714         self._refresh_data()
>    3715
>    3716     def _get_list_of_grids(self, buffer = 0.0):
>
> /home/jforbes/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in
> _refresh_data(self)
>    3729
>    3730     def _refresh_data(self):
> -> 3731         AMR3DData._refresh_data(self)
>    3732         self['dx'] = self.dds[0] * np.ones(self.ActiveDimensions,
> dtype='float64')
>    3733         self['dy'] = self.dds[1] * np.ones(self.ActiveDimensions,
> dtype='float64')
>
> /home/jforbes/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in
> _refresh_data(self)
>     318         """
>     319         self.clear_data()
> --> 320         self.get_data()
>     321
>     322     def keys(self):
>
> /home/jforbes/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in
> get_data(self, fields)
>    3773             n_bad = np.where(self[obtain_fields[0]]==-999)[0].size
>    3774             mylog.error("Covering problem: %s cells are uncovered",
> n_bad)
> -> 3775             raise KeyError(n_bad)
>    3776
>    3777     def _generate_field(self, field):
>
> KeyError: 190
>
>
> _______________________________________________
> yt-dev mailing list
> yt-dev at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
>