Particle Generators
@ John ZuHone | Friday, Jan 4, 2013 | 5 minute read | Update at Friday, Dec 4, 2020

This post shows how to generate particle fields from pre-defined particle lists, lattice distributions, and distributions based on density fields.

Particle Generators

Notebook Download

Generating particle initial conditions is now possible in yt. The following shows how to generate particle fields from pre-defined particle lists, lattice distributions, and distributions based on density fields.

First, we define a gridded density field where the particle density field has been “cloud-in-cell” (CIC) interpolated to the grid, and define a function that assigns a set of particle indices based on a number of particles and a starting index. This is for a case where we want to add particles to an already existing set but make sure they have uniqune indices.

In[1]:

    from yt.mods import *
    from yt.utilities.particle_generator import *
    import yt.utilities.initial_conditions as ic
    import yt.utilities.flagging_methods as fm
    from yt.frontends.stream.api import refine_amr
    from yt.utilities.lib import CICDeposit_3
    
    def _pgdensity(field, data):
        blank = np.zeros(data.ActiveDimensions, dtype='float32')
        if data.NumberOfParticles == 0: return blank
        CICDeposit_3(data["particle_position_x"].astype(np.float64),
                     data["particle_position_y"].astype(np.float64),
                     data["particle_position_z"].astype(np.float64),
                     data["particle_gas_density"].astype(np.float32),
                     np.int64(data.NumberOfParticles),
                     blank, np.array(data.LeftEdge).astype(np.float64),
                     np.array(data.ActiveDimensions).astype(np.int32),
                     np.float64(data['dx']))
        return blank
    add_field("particle_density_cic", function=_pgdensity,
              validators=[ValidateGridType()], 
              display_name=r"$\mathrm{Particle}\/\mathrm{Density}$")
    
    def add_indices(npart, start_num) :
        return np.arange((npart)) + start_num

Next, we’ll set up a uniform grid with some random density data:

In[2]:

    domain_dims = (128, 128, 128)
    dens = 0.1*np.random.random(domain_dims)
    fields = {"Density": dens}
    ug = load_uniform_grid(fields, domain_dims, 1.0)

As a first example, we’ll generate particle fields from pre-existing NumPy arrays. First, we define a list of particle field names, and then assign random positions to the particles in one corner of the grid. We then call FromListParticleGenerator, which generates the particles. assign_indices assigns the indices (using numpy.arange by default). apply_to_stream applies the particle fields to the grid.

In[3]:

    num_particles1 = 10000
    field_list = ["particle_position_x","particle_position_y",
                  "particle_position_z","particle_gas_density"]
    x = np.random.uniform(low=0.0, high=0.5, size=num_particles1) # random positions
    y = np.random.uniform(low=0.0, high=0.5, size=num_particles1) # random positions
    z = np.random.uniform(low=0.0, high=0.5, size=num_particles1) # random positions
    pdata = {'particle_position_x':x,
             'particle_position_y':y,
             'particle_position_z':z}
    particles1 = FromListParticleGenerator(ug, num_particles1, pdata)
    particles1.assign_indices()
    particles1.apply_to_stream()
    yt : [INFO     ] 2013-01-01 21:24:32,484 Adding Density to list of fields
    yt : [INFO     ] 2013-01-01 21:24:32,486 Adding particle_position_z to list of fields
    yt : [INFO     ] 2013-01-01 21:24:32,487 Adding particle_index to list of fields
    yt : [INFO     ] 2013-01-01 21:24:32,487 Adding particle_position_x to list of fields
    yt : [INFO     ] 2013-01-01 21:24:32,488 Adding particle_position_y to list of fields

Now that the particles are part of the parameter file, they may be manipulated and plotted:

In[4]:

    slc = SlicePlot(ug, 2, ["Density"], center=ug.domain_center)
    slc.set_cmap("Density","spring")
    slc.annotate_particles(0.2, p_size=10.0) # Display all particles within a thick slab 0.2 times the domain width
    slc.show()

Now let’s try adding a particle distribution in a lattice-shaped spatial arrangement. Let’s choose ten particles on a side, and place them in a small region away from the random particles. We’ll use the special add_indices function we defined earlier to assign indices that are all different from the ones the already existing particles have.

In[5]:

    pdims = np.array([10,10,10]) # number of particles on a side in each dimension
    ple = np.array([0.6,0.6,0.6]) # left edge of particle positions
    pre = np.array([0.9,0.9,0.9]) # right edge of particle positions
    particles2 = LatticeParticleGenerator(ug, pdims, ple, pre, field_list)
    particles2.assign_indices(function=add_indices, npart=np.product(pdims),
                              start_num=num_particles1)
    particles2.apply_to_stream()
    yt : [INFO     ] 2013-01-01 21:24:33,957 Adding particle_gas_density to list of fields

We now have both sets of particles:

In[6]:

    slc = SlicePlot(ug, 2, ["Density"], center=ug.domain_center)
    slc.set_cmap("Density","spring")
    slc.annotate_particles(0.2, p_size=10.0)
    slc.show()

And by sorting all of the indices we can check that all of them are unique, as advertised:

In[7]:

    dd = ug.h.all_data()
    indices = np.sort(np.int32(dd["particle_index"]))
    print "All indices unique = ", np.all(np.unique(indices) == indices)
    All indices unique =  True

Now let’s get fancy. We will use the initial conditions capabilities of yt to apply a spherically symmetric density distribution based on the “beta-model” functional form, and set up a refinement method based on overdensity. Then, we will call refine_amr to apply this density distribution and refine the grid based on the overdensity over some value.

In[8]:

    fo = [ic.BetaModelSphere(1.0,0.1,0.5,[0.5,0.5,0.5],{"Density":(10.0)})] 
    rc = [fm.flagging_method_registry["overdensity"](4.0)]
    pf = refine_amr(ug, rc, fo, 3)

Now, we have an interesting density field to serve as a distribution function for particle positions. What we do next is define a spherical region over which particle positions will be generated based on the local grid density. We also will map the grid density to a particle density field using cloud-in-cell interpolation. Finally, when we apply these particles, we will set the optional argument clobber=True, which will remove the particles we already created.

In[9]:

    num_particles3 = 100000
    map_dict = {"Density": "particle_gas_density"} # key is grid field, value is particle field
    sphere = pf.h.sphere(pf.domain_center, (0.5, "unitary"))
    particles3 = WithDensityParticleGenerator(pf, sphere, num_particles3,
                                              field_list)
    particles3.assign_indices()
    particles3.map_grid_fields_to_particles(map_dict) # Map density fields to particle fields
    particles3.apply_to_stream(clobber=True) # Get rid of all pre-existing particles

Now we’ll plot up both the grid density field and the “particle_density_cic” field (defined at the top of the script), which is mapped from the particles onto the grid. We also overplot the particle positions. These should roughly correspond to the non-zero values of “particle_density_cic”, but there will be some discrepancies due to the fact that they are taken from a thick slab and only a slice of the grid-based field is shown.

In[10]:

    slc = SlicePlot(pf, 2, ["Density","particle_density_cic"], center=pf.domain_center)
    slc.set_log("Density", True)
    slc.set_log("particle_density_cic", True)
    slc.set_cmap("all", "spring")
    slc.annotate_grids()
    slc.annotate_particles(0.01,p_size=3)
    slc.show()

yt extension modules

yt has many extension packages to help you in your scientific workflow! Check these out, or create your own.

ytini

ytini is set of tools and tutorials for using yt as a tool inside the 3D visual effects software Houdini or a data pre-processor externally to Houdini.

Trident

Trident is a full-featured tool that projects arbitrary sightlines through astrophysical hydrodynamics simulations for generating mock spectral observations of the IGM and CGM.

pyXSIM

pyXSIM is a Python package for simulating X-ray observations from astrophysical sources.

ytree

Analyze merger tree data from multiple sources. It’s yt for merger trees!

yt_idv

yt_idv is a package for interactive volume rendering with yt! It provides interactive visualization using OpenGL for datasets loaded in yt. It is written to provide both scripting and interactive access.

widgyts

widgyts is a jupyter widgets extension for yt, backed by rust/webassembly to allow for browser-based, interactive exploration of data from yt.

yt_astro_analysis

yt_astro_analysis is the yt extension package for astrophysical analysis.

Make your own!!

Finally, check out our development docs on writing your own yt extensions!

Contributing to the Blog

Are you interested in contributing to the yt blog?

Check out our post on contributing to the blog for a guide!

We welcome contributions from all members of the yt community. Feel free to reach out if you need any help.

the yt data hub

The yt hub at https://girder.hub.yt/ has a ton of resources to check out, whether you have yt installed or not.

The collections host all sorts of data that can be loaded with yt. Some have been used in publications, and others are used as sample frontend data for yt. Maybe there’s data from your simulation software?

The rafts host the yt quickstart notebooks, where you can interact with yt in the browser, without needing to install it locally. Check out some of the other rafts too, like the widgyts release notebooks – a demo of the widgyts yt extension pacakge; or the notebooks from the CCA workshop – a user’s workshop on using yt.

Social Links