yt-project / yt

Main yt repository
http://yt-project.org
Other
468 stars 277 forks source link

doc: example Particle Trajectories analysis fails #4590

Open Xarthisius opened 1 year ago

Xarthisius commented 1 year ago

Bug report

Bug summary

See: https://yt-project.org/doc/analyzing/particle_trajectories.html

Code for reproduction

import glob
from os.path import join

import yt
from yt.config import ytcfg

path = ytcfg.get("yt", "test_data_dir")
my_fns = glob.glob(join(path, "Orbit", "orbit_hdf5_chk_00[0-9][0-9]"))
my_fns.sort()
fields = ["particle_velocity_x", "particle_velocity_y", "particle_velocity_z"]
ds = yt.load(my_fns[0])
dd = ds.all_data()
indices = dd["all", "particle_index"].astype("int")
ts = yt.DatasetSeries(my_fns)
trajs = ts.particle_trajectories(indices, fields=fields, suppress_logging=True)

Actual outcome

Generating [particle_velocity_x, particle_velocity_y, particle_velocity_z] fields in trajectories:   0%| | 0/31 [00:00<Traceback (most recent call last):
  File "/home/xarth/yt/traj/ala.py", line 19, in <module>
    trajs = ts.particle_trajectories(indices, fields=fields, suppress_logging=True)
  File "/home/xarth/codes/xarthisius/yt/yt/data_objects/time_series.py", line 434, in particle_trajectories
    return ParticleTrajectories(
  File "/home/xarth/codes/xarthisius/yt/yt/data_objects/particle_trajectories.py", line 147, in __init__
    self._get_data(fields)
  File "/home/xarth/codes/xarthisius/yt/yt/data_objects/particle_trajectories.py", line 278, in _get_data
    cube = grid.retrieve_ghost_zones(1, grid_fields)
  File "/home/xarth/codes/xarthisius/yt/yt/data_objects/index_subobjects/grid_patch.py", line 265, in retrieve_ghost_zones
    cube = self.ds.covering_grid(
  File "/home/xarth/codes/xarthisius/yt/yt/data_objects/construction_data_containers.py", line 723, in __init__
    self.get_data(fields)
  File "/home/xarth/codes/xarthisius/yt/yt/data_objects/construction_data_containers.py", line 910, in get_data
    self[a] = self._data_source[f]
  File "/home/xarth/codes/xarthisius/yt/yt/data_objects/data_containers.py", line 235, in __getitem__
    self.get_data(f)
  File "/home/xarth/codes/xarthisius/yt/yt/data_objects/selection_objects/data_selection_objects.py", line 140, in get_data
    self.index._identify_base_chunk(self)
  File "/home/xarth/codes/xarthisius/yt/yt/geometry/grid_geometry_handler.py", line 333, in _identify_base_chunk
    gi = dobj.selector.select_grids(
  File "/home/xarth/codes/xarthisius/yt/yt/data_objects/selection_objects/data_selection_objects.py", line 90, in selector
    self._selector = sclass(self)
  File "yt/geometry/_selection_routines/region_selector.pxi", line 60, in yt.geometry.selection_routines.RegionSelector.__init__
RuntimeError: yt attempted to read outside the boundaries of a non-periodic domain along dimension 0.
Region left edge = -0.0625 code_length, Region right edge = 0.5625 code_length
Dataset left edge = 0.0 code_length, Dataset right edge = 1.0 code_length

This commonly happens when trying to compute ghost cells up to the domain boundary. Two possible solutions are to select a smaller region that does not border domain edge (see https://yt-project.org/docs/analyzing/objects.html?highlight=region)
or override the periodicity with
ds.force_periodicity()

Expected outcome

No error

mtryan83 commented 10 months ago

Per @neutrinoceros, I'm adding some of the info from #4765.

Briefly, it looks like the problem occurs when yt tries to add ("all","particle_position_x") as a grid field in the _get_data() function. I don't know if that's expected or not. Note that ("io","particle_position_x") works just fine, presumably because it's not a grid field (i.e. passes the if ftype in self.data_series[0].particle_types: statement in line 237). I don't know the FLASH format, so I don't know if "all" and "io" are equivalent. Alternatively, it might be a problem in how the grid is loaded?

So from what I saw, it looked like there were three possible sources for the bug:

  1. The "all" particle type isn't included in the FLASH dataset list of particle types and/or isn't mapped to "io" correctly
  2. if ftype in self.data_series[0].particle_types: isn't a robust way to check if a field is a grid field
  3. There's a problem with how the grid is loaded, possibly related to #4137?

Since it seems to load just fine when using "io", I would guess 1. or 2. would be the issue, but I don't know if grid information is relevant. That seems like a question for @jzuhone .

Also, I'd expect the test test_particle_trajectories.py::test_orbit_traj to fail based on this, but wasn't able to test.

mtryan83 commented 10 months ago

Just FYI, it looks like test_orbit_traj and test_etc_traj aren't even running. For a little bit more, see #4768.

I think the test problem can be resolved simply by changing the require_ds decorator to require_file (a la #4692) which results in both tests failing on my machine. The test_orbit_traj is failing for the reason described here. The other test failure might be related, it appears to occur in the same grid_field section, but a slightly different part:

# This will fail for non-grid index objects
for grid in particle_grids:
    cube = grid.retrieve_ghost_zones(1, grid_fields) # test_orbit_traj fails here
    for field in grid_fields:
        CICSample_3(  # test_etc_traj fails here
            x,
            y,
            z,
            pfield[field],
            self.num_indices,
            cube[fds[field]],
            np.array(grid.LeftEdge, dtype="float64"),
            np.array(grid.ActiveDimensions, dtype="int32"),
            grid.dds[0],
        )

with the error File "yt/utilities/lib/particle_mesh_operations.pyx", line 214, in yt.utilities.lib.particle_mesh_operations.CICSample_3 ValueError: Buffer has wrong number of dimensions (expected 3, got 1)

Based on the code comment above, I'm guessing these fields aren't supposed to be considered grid fields, and the bug is either 1 or 2 from my previous comment