One can create particle trajectories from a DatasetSeries object for a specified list of particles identified by their unique indices using the particle_trajectories method.
[1]:
%matplotlib inline
import glob
from os.path import join
import yt
from yt.config import ytcfg
path = ytcfg.get("yt", "test_data_dir")
import matplotlib.pyplot as plt
First, let’s start off with a FLASH dataset containing only two particles in a mutual circular orbit. We can get the list of filenames this way:
[2]:
my_fns = glob.glob(join(path, "Orbit", "orbit_hdf5_chk_00[0-9][0-9]"))
my_fns.sort()
And let’s define a list of fields that we want to include in the trajectories. The position fields will be included by default, so let’s just ask for the velocity fields:
[3]:
fields = ["particle_velocity_x", "particle_velocity_y", "particle_velocity_z"]
There are only two particles, but for consistency’s sake let’s grab their indices from the dataset itself:
[4]:
ds = yt.load(my_fns[0])
dd = ds.all_data()
indices = dd["all", "particle_index"].astype("int")
print(indices)
[1 2] dimensionless
which is what we expected them to be. Now we’re ready to create a DatasetSeries object and use it to create particle trajectories:
[5]:
ts = yt.DatasetSeries(my_fns)
# suppress_logging=True cuts down on a lot of noise
trajs = ts.particle_trajectories(indices, fields=fields, suppress_logging=True)
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
Cell In[5], line 3
1 ts = yt.DatasetSeries(my_fns)
2 # suppress_logging=True cuts down on a lot of noise
----> 3 trajs = ts.particle_trajectories(indices, fields=fields, suppress_logging=True)
File /tmp/yt/yt/data_objects/time_series.py:434, in DatasetSeries.particle_trajectories(self, indices, fields, suppress_logging, ptype)
384 def particle_trajectories(
385 self, indices, fields=None, suppress_logging=False, ptype=None
386 ):
387 r"""Create a collection of particle trajectories in time over a series of
388 datasets.
389
(...)
432 particle disappear.
433 """
--> 434 return ParticleTrajectories(
435 self, indices, fields=fields, suppress_logging=suppress_logging, ptype=ptype
436 )
File /tmp/yt/yt/data_objects/particle_trajectories.py:147, in ParticleTrajectories.__init__(self, outputs, indices, fields, suppress_logging, ptype)
144 self.particle_fields.append(field)
146 # Instantiate fields the caller requested
--> 147 self._get_data(fields)
File /tmp/yt/yt/data_objects/particle_trajectories.py:278, in ParticleTrajectories._get_data(self, fields)
276 # This will fail for non-grid index objects
277 for grid in particle_grids:
--> 278 cube = grid.retrieve_ghost_zones(1, grid_fields)
279 for field in grid_fields:
280 CICSample_3(
281 x,
282 y,
(...)
289 grid.dds[0],
290 )
File /tmp/yt/yt/data_objects/index_subobjects/grid_patch.py:264, in AMRGridPatch.retrieve_ghost_zones(self, n_zones, fields, all_levels, smoothed)
260 cube = self.ds.smoothed_covering_grid(
261 level, new_left_edge, field_parameters=field_parameters, **kwargs
262 )
263 else:
--> 264 cube = self.ds.covering_grid(
265 level, new_left_edge, field_parameters=field_parameters, **kwargs
266 )
267 cube._base_grid = self
268 return cube
File /tmp/yt/yt/data_objects/construction_data_containers.py:723, in YTCoveringGrid.__init__(self, level, left_edge, dims, fields, ds, num_ghost_zones, use_pbar, field_parameters, data_source)
716 self.global_startindex = (
717 np.rint((self.left_edge - self.ds.domain_left_edge) / self.dds).astype(
718 "int64"
719 )
720 + self.ds.domain_offset
721 )
722 self._setup_data_source()
--> 723 self.get_data(fields)
File /tmp/yt/yt/data_objects/construction_data_containers.py:910, in YTCoveringGrid.get_data(self, fields)
908 for a, f in sorted(alias.items()):
909 if f.sampling_type == "particle" and not is_sph_field:
--> 910 self[a] = self._data_source[f]
911 else:
912 self[a] = f(self)
File /tmp/yt/yt/data_objects/data_containers.py:235, in YTDataContainer.__getitem__(self, key)
233 return self.field_data[f]
234 else:
--> 235 self.get_data(f)
236 # fi.units is the unit expression string. We depend on the registry
237 # hanging off the dataset to define this unit object.
238 # Note that this is less succinct so that we can account for the case
239 # when there are, for example, no elements in the object.
240 try:
File /tmp/yt/yt/data_objects/selection_objects/data_selection_objects.py:139, in YTSelectionContainer.get_data(self, fields)
137 def get_data(self, fields=None):
138 if self._current_chunk is None:
--> 139 self.index._identify_base_chunk(self)
140 if fields is None:
141 return
File /tmp/yt/yt/geometry/grid_geometry_handler.py:333, in GridIndex._identify_base_chunk(self, dobj)
331 dobj._chunk_info[0] = weakref.proxy(dobj)
332 elif getattr(dobj, "_grids", None) is None:
--> 333 gi = dobj.selector.select_grids(
334 self.grid_left_edge, self.grid_right_edge, self.grid_levels
335 )
336 if any(g.filename is not None for g in self.grids[gi]):
337 _gsort = _grid_sort_mixed
File /tmp/yt/yt/data_objects/selection_objects/data_selection_objects.py:89, in YTSelectionContainer.selector(self)
85 self._selector = compose_selector(
86 self, self._data_source.selector, sclass(self)
87 )
88 else:
---> 89 self._selector = sclass(self)
90 return self._selector
File /tmp/yt/yt/geometry/_selection_routines/region_selector.pxi: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()
The ParticleTrajectories object trajs is essentially a dictionary-like container for the particle fields along the trajectory, and can be accessed as such:
[6]:
print(trajs["all", "particle_position_x"])
print(trajs["all", "particle_position_x"].shape)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 1
----> 1 print(trajs["all", "particle_position_x"])
2 print(trajs["all", "particle_position_x"].shape)
NameError: name 'trajs' is not defined
Note that each field is a 2D NumPy array with the different particle indices along the first dimension and the times along the second dimension. As such, we can access them individually by indexing the field:
[7]:
plt.figure(figsize=(6, 6))
plt.plot(trajs["all", "particle_position_x"][0], trajs["all", "particle_position_y"][0])
plt.plot(trajs["all", "particle_position_x"][1], trajs["all", "particle_position_y"][1])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 2
1 plt.figure(figsize=(6, 6))
----> 2 plt.plot(trajs["all", "particle_position_x"][0], trajs["all", "particle_position_y"][0])
3 plt.plot(trajs["all", "particle_position_x"][1], trajs["all", "particle_position_y"][1])
NameError: name 'trajs' is not defined
<Figure size 600x600 with 0 Axes>
And we can plot the velocity fields as well:
[8]:
plt.figure(figsize=(6, 6))
plt.plot(trajs["all", "particle_velocity_x"][0], trajs["all", "particle_velocity_y"][0])
plt.plot(trajs["all", "particle_velocity_x"][1], trajs["all", "particle_velocity_y"][1])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 2
1 plt.figure(figsize=(6, 6))
----> 2 plt.plot(trajs["all", "particle_velocity_x"][0], trajs["all", "particle_velocity_y"][0])
3 plt.plot(trajs["all", "particle_velocity_x"][1], trajs["all", "particle_velocity_y"][1])
NameError: name 'trajs' is not defined
<Figure size 600x600 with 0 Axes>
If we want to access the time along the trajectory, we use the key "particle_time":
[9]:
plt.figure(figsize=(6, 6))
plt.plot(trajs["particle_time"], trajs["particle_velocity_x"][1])
plt.plot(trajs["particle_time"], trajs["particle_velocity_y"][1])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 2
1 plt.figure(figsize=(6, 6))
----> 2 plt.plot(trajs["particle_time"], trajs["particle_velocity_x"][1])
3 plt.plot(trajs["particle_time"], trajs["particle_velocity_y"][1])
NameError: name 'trajs' is not defined
<Figure size 600x600 with 0 Axes>
Alternatively, if we know the particle index we’d like to examine, we can get an individual trajectory corresponding to that index:
[10]:
particle1 = trajs.trajectory_from_index(1)
plt.figure(figsize=(6, 6))
plt.plot(particle1["all", "particle_time"], particle1["all", "particle_position_x"])
plt.plot(particle1["all", "particle_time"], particle1["all", "particle_position_y"])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 1
----> 1 particle1 = trajs.trajectory_from_index(1)
2 plt.figure(figsize=(6, 6))
3 plt.plot(particle1["all", "particle_time"], particle1["all", "particle_position_x"])
NameError: name 'trajs' is not defined
Now let’s look at a more complicated (and fun!) example. We’ll use an Enzo cosmology dataset. First, we’ll find the maximum density in the domain, and obtain the indices of the particles within some radius of the center. First, let’s have a look at what we’re getting:
[11]:
ds = yt.load("enzo_tiny_cosmology/DD0046/DD0046")
slc = yt.SlicePlot(
ds,
"x",
[("gas", "density"), ("gas", "dark_matter_density")],
center="max",
width=(3.0, "Mpc"),
)
slc.show()
So far, so good–it looks like we’ve centered on a galaxy cluster. Let’s grab all of the dark matter particles within a sphere of 0.5 Mpc (identified by "particle_type == 1"):
[12]:
sp = ds.sphere("max", (0.5, "Mpc"))
indices = sp["all", "particle_index"][sp["all", "particle_type"] == 1]
Next we’ll get the list of datasets we want, and create trajectories for these particles:
[13]:
my_fns = glob.glob(join(path, "enzo_tiny_cosmology/DD*/*.hierarchy"))
my_fns.sort()
ts = yt.DatasetSeries(my_fns)
trajs = ts.particle_trajectories(indices, fields=fields, suppress_logging=True)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[13], line 4
2 my_fns.sort()
3 ts = yt.DatasetSeries(my_fns)
----> 4 trajs = ts.particle_trajectories(indices, fields=fields, suppress_logging=True)
File /tmp/yt/yt/data_objects/time_series.py:434, in DatasetSeries.particle_trajectories(self, indices, fields, suppress_logging, ptype)
384 def particle_trajectories(
385 self, indices, fields=None, suppress_logging=False, ptype=None
386 ):
387 r"""Create a collection of particle trajectories in time over a series of
388 datasets.
389
(...)
432 particle disappear.
433 """
--> 434 return ParticleTrajectories(
435 self, indices, fields=fields, suppress_logging=suppress_logging, ptype=ptype
436 )
File /tmp/yt/yt/data_objects/particle_trajectories.py:147, in ParticleTrajectories.__init__(self, outputs, indices, fields, suppress_logging, ptype)
144 self.particle_fields.append(field)
146 # Instantiate fields the caller requested
--> 147 self._get_data(fields)
File /tmp/yt/yt/data_objects/particle_trajectories.py:280, in ParticleTrajectories._get_data(self, fields)
278 cube = grid.retrieve_ghost_zones(1, grid_fields)
279 for field in grid_fields:
--> 280 CICSample_3(
281 x,
282 y,
283 z,
284 pfield[field],
285 self.num_indices,
286 cube[fds[field]],
287 np.array(grid.LeftEdge, dtype="float64"),
288 np.array(grid.ActiveDimensions, dtype="int32"),
289 grid.dds[0],
290 )
291 sto.result_id = ds.parameter_filename
292 sto.result = (self.array_indices[i], pfield)
File /tmp/yt/yt/utilities/lib/particle_mesh_operations.pyx:214, in yt.utilities.lib.particle_mesh_operations.CICSample_3()
ValueError: Buffer has wrong number of dimensions (expected 3, got 1)
Matplotlib can make 3D plots, so let’s pick three particle trajectories at random and look at them in the volume:
[14]:
fig = plt.figure(figsize=(8.0, 8.0))
ax = fig.add_subplot(111, projection="3d")
ax.plot(
trajs["all", "particle_position_x"][100],
trajs["all", "particle_position_y"][100],
trajs["all", "particle_position_z"][100],
)
ax.plot(
trajs["all", "particle_position_x"][8],
trajs["all", "particle_position_y"][8],
trajs["all", "particle_position_z"][8],
)
ax.plot(
trajs["all", "particle_position_x"][25],
trajs["all", "particle_position_y"][25],
trajs["all", "particle_position_z"][25],
)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 4
1 fig = plt.figure(figsize=(8.0, 8.0))
2 ax = fig.add_subplot(111, projection="3d")
3 ax.plot(
----> 4 trajs["all", "particle_position_x"][100],
5 trajs["all", "particle_position_y"][100],
6 trajs["all", "particle_position_z"][100],
7 )
8 ax.plot(
9 trajs["all", "particle_position_x"][8],
10 trajs["all", "particle_position_y"][8],
11 trajs["all", "particle_position_z"][8],
12 )
13 ax.plot(
14 trajs["all", "particle_position_x"][25],
15 trajs["all", "particle_position_y"][25],
16 trajs["all", "particle_position_z"][25],
17 )
NameError: name 'trajs' is not defined
It looks like these three different particles fell into the cluster along different filaments. We can also look at their x-positions only as a function of time:
[15]:
plt.figure(figsize=(6, 6))
plt.plot(trajs["all", "particle_time"], trajs["all", "particle_position_x"][100])
plt.plot(trajs["all", "particle_time"], trajs["all", "particle_position_x"][8])
plt.plot(trajs["all", "particle_time"], trajs["all", "particle_position_x"][25])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 2
1 plt.figure(figsize=(6, 6))
----> 2 plt.plot(trajs["all", "particle_time"], trajs["all", "particle_position_x"][100])
3 plt.plot(trajs["all", "particle_time"], trajs["all", "particle_position_x"][8])
4 plt.plot(trajs["all", "particle_time"], trajs["all", "particle_position_x"][25])
NameError: name 'trajs' is not defined
<Figure size 600x600 with 0 Axes>
Suppose we wanted to know the gas density along the particle trajectory, but there wasn’t a particle field corresponding to that in our dataset. Never fear! If the field exists as a grid field, yt will interpolate this field to the particle positions and add the interpolated field to the trajectory. To add such a field (or any field, including additional particle fields) we can call the add_fields method:
[16]:
trajs.add_fields([("gas", "density")])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 1
----> 1 trajs.add_fields([("gas", "density")])
NameError: name 'trajs' is not defined
We also could have included "density" in our original field list. Now, plot up the gas density for each particle as a function of time:
[17]:
plt.figure(figsize=(6, 6))
plt.plot(trajs["all", "particle_time"], trajs["gas", "density"][100])
plt.plot(trajs["all", "particle_time"], trajs["gas", "density"][8])
plt.plot(trajs["all", "particle_time"], trajs["gas", "density"][25])
plt.yscale("log")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 2
1 plt.figure(figsize=(6, 6))
----> 2 plt.plot(trajs["all", "particle_time"], trajs["gas", "density"][100])
3 plt.plot(trajs["all", "particle_time"], trajs["gas", "density"][8])
4 plt.plot(trajs["all", "particle_time"], trajs["gas", "density"][25])
NameError: name 'trajs' is not defined
<Figure size 600x600 with 0 Axes>
Finally, the particle trajectories can be written to disk. Two options are provided: ASCII text files with a column for each field and the time, and HDF5 files:
[18]:
trajs.write_out(
"halo_trajectories"
) # This will write a separate file for each trajectory
trajs.write_out_h5(
"halo_trajectories.h5"
) # This will write all trajectories to a single file
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 trajs.write_out(
2 "halo_trajectories"
3 ) # This will write a separate file for each trajectory
4 trajs.write_out_h5(
5 "halo_trajectories.h5"
6 ) # This will write all trajectories to a single file
NameError: name 'trajs' is not defined