Hello,
I’m trying to get yt to work with Enzo datasets with ghost zones. I’m having an Enzo crash in a large simulation I’m trying to debug and it would be helpful to inspect the ghost zones in yt. I see there has been significant work to build up the framework in the Enzo frontend to allow for this but much of it is inconsistent and not working yet. Sorry for this being long and overly verbose. I want to make sure I explain things in enough detail since there are a lot of moving parts. It also helps me learn how the frontend works if I explain everything. There's a TLDR at the end. I've opened a yt issue but thought I should post a message here too.
There are 3 main problems I’ve identified. One is trivial to fix, the other 2 require more work and a well thought out plan so that current working parts of the code are not broken in the process.
In `yt/frontends/enzo/data_structures.py` there’s a class `EnzoGridGZ` derived from `EnzoGrid`. It adds a member function `retrieve_ghost_zones` which is supposed to make a covering_grid including the ghost zones. This is exactly what I need.
The first problem is the call
`self.index.covering_grid(level, new_left_edge, **kwargs)`
should be
`self.ds.covering_grid(level, new_left_edge, **kwargs)`
Also, the kwargs dictionary is not created properly. This is an easy fix.
It was
```kwargs = {'dims': self.ActiveDimensions + 2*n_zones,
'num_ghost_zones':n_zones,
'use_pbar':False}
kwargs.update(self.field_parameters)```
it should be
```kwargs = {'dims': self.ActiveDimensions + 2*n_zones,
'fields': fields,
'ds': self.ds,
'num_ghost_zones':n_zones,
'use_pbar':False,
'field_parameters': self.field_parameters}
```
Next are the real problems. Here’s the call stack when yt fails. I’ll walk through it explaining flow of control and where the problems arise. I'm using a the `<enzo_dir>/run/Cosmology/AdiabaticExpansion` simulation with a 16^3 root grid and 3 ghost zones.
```
/Users/coreybrummel-smith/simulations/testing/WriteGhostZones/ghost.py(21)<module>()
-> gz = grid_gz.retrieve_ghost_zones(3, ['density'])
yt/frontends/enzo/data_structures.py(146)retrieve_ghost_zones()
-> level, new_left_edge, **kwargs
yt/data_objects/construction_data_containers.py(555)__init__()
-> self.get_data(fields)
yt/data_objects/construction_data_containers.py(642)get_data()
-> if len(fill) > 0: self._fill_fields(fill)
yt/data_objects/construction_data_containers.py(698)_fill_fields()
-> for chunk in parallel_objects(self._data_source.chunks(fields, "io")):
yt/utilities/parallel_tools/parallel_analysis_interface.py(507)parallel_objects()
-> for result_id, obj in oiter:
yt/data_objects/data_containers.py(1276)chunks()
-> self.get_data(fields)
yt/data_objects/data_containers.py(1370)get_data()
-> (fluids, self, self._current_chunk)
yt/geometry/geometry_handler.py(245)_read_fluid_fields()
-> chunk_size
/utilities/io_handler.py(148)_read_fluid_selection()
-> ind[field] += obj.select(selector, data, rv[field], ind[field])
yt/data_objects/grid_patch.py(423)select()
-> dest[offset:offset+count, i] = source[sl][np.squeeze(mask)]
```
Things start out okay, the covering grid has the correct dimensions including ghost zones and correct edges
```
ipdb> self.ActiveDimensions
array([22, 22, 22], dtype=int32)
ipdb> self.left_edge
YTArray([-0.1875, -0.1875, -0.1875]) code_length
ipdb> self.right_edge
YTArray([1.1875, 1.1875, 1.1875]) code_length
```
Now `_fill_fields()` is called from `get_data()`. Here we do
```
output_fields = [np.zeros(self.ActiveDimensions, dtype="float64")
for field in fields]
domain_dims = self.ds.domain_dimensions.astype("int64") \
* self.ds.relative_refinement(0, self.level)
```
`output_fields[0].shape` is correct (22,22,22), but `domain_dims` is (16,16,16) because it comes from TopGridDimensions in the Enzo parameter file.
Next, `_fill_fields()` loops over chunks while calling `fill_region()`
As soon as `self._data_source.chunks()` is called, the YTRegion (self._data_source) associated with the covering grid calls `get_data()`. Since there is no _current_chunk yet, `self.index._identify_base_chunk()` gets called. By doing this, `dobj.selector.select_grids()` gets called which modifies the YTRegion’s left edge and right edge.
```
ipdb> dobj.left_edge
YTArray([0.8125, 0.8125, 0.8125]) code_length
ipdb> dobj.right_edge
YTArray([2.1875, 2.1875, 2.1875]) code_length
```
Then the size is computed
```dobj.size = self._count_selection(dobj, fast_index = fast_index)```
_count_selection()
The size ultimately comes from ActiveDimensions of the EnzoGrid in YTRegion's base chunk which is (16,16,16) not (22,22,22).
Now back in `self._chunk()` where self is the EnzoHierarchy of the YTCoveringGrid. We call `self._chunk_io()` and yield the chunk back to the loop in `YTRegion.chunks()`
Now that we've prepared the chunk we try to load in the density field with
`self.get_data(fields)` which calls `self.index._read_fluid_fields(fluids, self, self._current_chunk)` which calls `self.io._read_fluid_selection(self._chunk_io(dobj), selector, fields_to_read, chunk_size)`
here a dictionary named `rv` is created which will store the field data.
`rv[field] = np.empty(size, dtype="=f8")`
But `size` is determined by the chunk size which is 16^3 (b/c the grid's ActiveDimensions are 16^3). If it is to include the ghost zones, it should be 22^3 instead.
Now we do
```for field, obj, data in self.io_iter(chunks, fields):
...
ind[field] += obj.select(selector, data, rv[field], ind[field])
```
`self.io_iter(chunks, fields)` loops over the chunks, grids, and fields, opens the appropriate Enzo hdf5 dataset, then does
```
dims = obj.ActiveDimensions[::-1] + nodal_flag[::-1]
data = np.empty(dims, dtype=h5_dtype)
yield field, obj, self._read_obj_field(
obj, field, (fid, data))
```
`data` a 16^3 ndarray again because the grid's ActiveDimensions are 16^3.
`self._read_obj_field` reads the data from the hdf5 dataset and stores it in the `data` array and returns the transpose of this array to the io_itter function. Here's tho code.
```
dg = h5py.h5d.open(fid, b(node))
dg.read(h5py.h5s.ALL, h5py.h5s.ALL, data)
return data.T
```
Interestingly this doesn't fail even though the dataset is 22^3 not 16^3. The `h5py.h5d.read()` doc says "It is your responsibility to ensure that the memory dataspace provided is compatible with the shape of the Numpy array. ... You can easily crash Python by reading in data from too large a dataspace." Not sure why it doesn't crash since these two don't match in this case.
Anyway, here's the second key problem. Because `WriteGhostZones = 1` in the Enzo parameter file, `self` is an `IOHandlerPackedHDF5GhostZones` object and its version of `_read_obj_field()` returns the data with the ghost zones trimmed off. Here's the class definition so you can see what I mean.
```
class IOHandlerPackedHDF5GhostZones(IOHandlerPackedHDF5):
_dataset_type = "enzo_packed_3d_gz"
def __init__(self, *args, **kwargs):
super(IOHandlerPackedHDF5GhostZones, self).__init__(*args, **kwargs)
NGZ = self.ds.parameters.get("NumberOfGhostZones", 3)
self._base = (slice(NGZ, -NGZ),
slice(NGZ, -NGZ),
slice(NGZ, -NGZ))
def _read_obj_field(self, *args, **kwargs):
return super(IOHandlerPackedHDF5GhostZones, self)._read_obj_field(
*args, **kwargs)[self._base]
```
Obviously this is supposed to be for when you have a dataset with ghost zones, most of the time you don't want to include them when inspecting data with yt. For instance, since my root grid is 16^3 with 3 ghost zones, instead of returning a 22^3 data cube it would return the 16^3 datacube. But here we want to use the ghost zones and furthermore, the dimensions weren't set properly so it returns a 10^3 array.
Now we yield the 10^3 data array back to the loop in `_read_fluid_selection`. Then we do
```
ind[field] += obj.select(selector, data, rv[field], ind[field])
```
which is supposed to copy `data` into `rv`. Here we encounter another problem. Here's the code for `obj.select()`
```
def select(self, selector, source, dest, offset):
mask = self._get_selector_mask(selector)
count = self.count(selector)
if count == 0: return 0
dim = np.squeeze(self.ds.dimensionality)
nodal_flag = source.shape[:dim] - self.ActiveDimensions[:dim]
# DEBUG NOTE:
# source.shape = [10,10,10]
# self.ActiveDimensions = [16,16,16]
if sum(nodal_flag) == 0:
dest[offset:offset+count] = source[mask]
else:
slices = get_nodal_slices(source.shape, nodal_flag, dim)
for i , sl in enumerate(slices):
dest[offset:offset+count, i] = source[sl][np.squeeze(mask)]
# DEBUG NOTE:
# source size originally comes from chunk size and was then modified to not include ghost zones
return count
```
`nodal_flag` should be `[0,0,0]` for the density field but since `source.shape` is incorrect, `nodal_flag = [-6,-6,-6]`.
Thus, we enter the else: block which should not happen for the density field. Finally we execute the
`dest[offset:offset+count, i] = source[sl][np.squeeze(mask)]` command and fail because `source[sl].shape` is `(9,9,9)` and `mask.shape` is `(16,16,16)`
TLDR:
In summary, there are two main problems with trying to create a YTCoveringGrid of an Enzo grid from a dataset with ghost zones. When trying to load in the field data, yt prepares numpy arrays whose shape is determined from the EnzoGrid's ActiveDimensions which does not include the ghost zones. This problem is compounded with the second problem that yt assumes we have loaded in the the full grid including ghost zones and then trims the edges of the array to eliminate the would-be ghost cells.
I'm not sure the best way to go about solving this. I don't think it is wise or elegant to just throw if statements everywhere to check if we are are trying to create a covering grid including ghost zones.
I think there should be two ways to handle Enzo datasets with ghost zones.
1) Most common: Load in an Enzo dataset that has ghost zones but treat it like a regular dataset and exclude the ghost zones when making a CoveringGrid. This is the idea behind the `IOHandlerPackedHDF5GhostZones` class.
2) My interest: Make a covering grid of an individual Enzo grid including its ghost zones. I think this was the idea for the `EnzoGridGZ.retrieve_ghost_zones()` function.
If anyone has any ideas, I'd appreciate your input.
Thanks,
Corey
