Particle Generators

Author: John ZuHone
Published on: Jan 4, 2013, 7:05:00 AM
Permalink - Source code

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.


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 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
                 blank, np.array(data.LeftEdge).astype(np.float64),
    return blank
add_field("particle_density_cic", function=_pgdensity,

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:


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.


num_particles1 = 10000
field_list = ["particle_position_x","particle_position_y",
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,
particles1 = FromListParticleGenerator(ug, num_particles1, pdata)
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:


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

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.


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),
yt : [INFO     ] 2013-01-01 21:24:33,957 Adding particle_gas_density to list of fields

We now have both sets of particles:


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

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


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.


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.


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,
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.


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")