Athena_pp frontend crashes when used for unstructured meshes with refinement

[Yurlungur]

Bug report

Bug summary

The Athena++ frontend crashes when reading in unstructured meshes with refinement. The problem is on line 109 of yt/frontends/athena_pp/data_structures.py, where bc[i] is an array with no shape. I'm using yt-dev but I was able to reproduce the problem in yt 3.4.0.

Code for reproduction

To reproduce, just generate the initial date for the Fishbone Moncrief Torus as described here in the Athena++ tutorial. File has been curldropped here.

Then try to read it in with the following code (I used python2, but should be the same for python3):

units_override = {"length_unit":(1.0,"Mpc"),
                  "time_unit":(1.0,"Myr"),
                  "mass_unit":(1.0e14,"Msun")}
ds = yt.load("fm_torus.cons.00000.athdf", units_override=units_override)
yt.SlicePlot(ds,'phi',('athena_pp','dens'))

Actual outcome

The result is the following traceback.

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-16-5aaf470c0415> in <module>()
----> 1 ds.field_list

/home/jonahm/programming/yt/yt/data_objects/static_output.pyc in field_list(self)
    546     @property
    547     def field_list(self):
--> 548         return self.index.field_list
    549 
    550     def create_field_info(self):

/home/jonahm/programming/yt/yt/data_objects/static_output.pyc in index(self)
    504                 raise RuntimeError("You should not instantiate Dataset.")
    505             self._instantiated_index = self._index_class(
--> 506                 self, dataset_type=self.dataset_type)
    507             # Now we do things that we need an instantiated index for
    508             # ...first off, we create our field_info now.

/home/jonahm/programming/yt/yt/frontends/athena_pp/data_structures.py in __init__(self, ds, dataset_type)
     65     def __init__(self, ds, dataset_type = 'athena_pp'):
     66         self._handle = ds._handle
---> 67         super(AthenaPPLogarithmicIndex, self).__init__(ds, dataset_type)
     68         self.index_filename = self.dataset.filename
     69         self.directory = os.path.dirname(self.dataset.filename)

/home/jonahm/programming/yt/yt/geometry/unstructured_mesh_handler.pyc in __init__(self, ds, dataset_type)
     34         self.directory = os.path.dirname(self.index_filename)
     35         self.float_type = np.float64
---> 36         super(UnstructuredIndex, self).__init__(ds, dataset_type)
     37 
     38     def _setup_geometry(self):

/home/jonahm/programming/yt/yt/geometry/geometry_handler.pyc in __init__(self, ds, dataset_type)
     48 
     49         mylog.debug("Setting up domain geometry.")
---> 50         self._setup_geometry()
     51 
     52         mylog.debug("Initializing data grid data IO")

/home/jonahm/programming/yt/yt/geometry/unstructured_mesh_handler.pyc in _setup_geometry(self)
     38     def _setup_geometry(self):
     39         mylog.debug("Initializing Unstructured Mesh Geometry Handler.")
---> 40         self._initialize_mesh()
     41 
     42     def get_smallest_dx(self):

/home/jonahm/programming/yt/yt/frontends/athena_pp/data_structures.py in _initialize_mesh(self)
    106         pbar = get_pbar("Constructing meshes", num_meshes)
    107         for i in range(num_meshes):
--> 108             ob = bc[i][0]
    109             x = x1f[ob,:]
    110             y = x2f[ob,:]

IndexError: too many indices for array

Expected outcome

What should happen is that the code should load the dataset completely and I should get a plot.

I've traced the problem to an interaction with some lines above. The bc object is constructed above with

neigh = block_grid[ii:ii+2,jj:jj+2,kk:kk+2,levels[i]]
                if np.all(neigh > -1):
                    loc_ids = neigh.transpose().flatten()
                    bc.append(loc_ids)
                    block_list[loc_ids] = -1
                else:
                    bc.append(np.array(i))
                    block_list[i] = -1

meaning the bc object contains some arrays with no shape... i.e., len(bc[i].shape) == 0. This causes a problem in the loop below which contains multiple calls to bc[i][j] for various i and j.

I can make the code run to completion by adding the following lines in the loop that accesses bc[i][j], but I worry that this is causing yt to read in incomplete data.

if len(bc[i].shape) == 0:
    continue

I was hoping John, or whomever is currently caring for the Athena++ frontend could take a look at this.

Thanks!

Version Information

• Operating System: Linux/Debian Testing
• Python Version: 2.7 (but the problem appears in python 3 too)
• yt version: 3.5-dev but can reproduce on 3.4.0

I installed yt via git.

[jzuhone]

This has likely been resolved by #4562. Given how long it's been, I don't know if @Yurlungur wants to try again, but I'm closing this unless it needs reopening.

[neutrinoceros]

I just checked and I confirm that #4562 fixed this, thank you John !