[yt-users] cloud in cell mesh construction for particle data

Nathan Goldbaum nathan12343 at gmail.com
Mon Jun 9 10:22:31 PDT 2014


On Mon, Jun 9, 2014 at 12:32 AM, Brendan Griffen <
brendan.f.griffen at gmail.com> wrote:

> Thanks Nathan. Comparing the two though:
>
> 26 1536.664 MiB 1130.480 MiB       uniform_array = grid_object['deposit', 'all_cic']
>
> 24 1285.703 MiB  890.539 MiB       cic_density = ad["deposit", "all_cic"]
>
> Isn't it the uniform grid which uses the most memory? I'm running the profiling now and will get back to you. I hope in spite of the fact that it does crash it will still give me some useful output.
>
>
That's true - I was thinking it might be more memory efficient since there
is no need to construct the covering grid in addition to the octree.


> Brendan
>
>
>
> On Mon, Jun 9, 2014 at 2:05 AM, Nathan Goldbaum <nathan12343 at gmail.com>
> wrote:
>
>> Hey Brendan,
>>
>> Could you try running your script using the memory_profiler module on
>> pypi?
>>
>> Here's an example script that uses the memory profiler:
>> http://paste.yt-project.org/show/4748/
>>
>> and the output for that script: http://paste.yt-project.org/show/4749/
>>
>> For what it's worth, it does indeed look like matt's suggestion to use
>> load_uniform_grid is an option, and might be more memory efficient in the
>> end since you will go directly to the uniform grid you want without
>> creating an octree.  Here's an example:
>> http://paste.yt-project.org/show/4750/
>>
>> Here's the memory usage information for that example:
>> http://paste.yt-project.org/show/4751/
>>
>> I used a 256^3 uniform grid with normally distributed random data - I'm
>> not sure whether it will also be more memory efficient in your case.
>>
>>
>>
>> On Sun, Jun 8, 2014 at 9:32 PM, Brendan Griffen <
>> brendan.f.griffen at gmail.com> wrote:
>>
>>> This is the full error if it helps at all? It is indeed, loading in all
>>> of the quantities.
>>>
>>> Loading particles...
>>>  --> Loading particle type: 1
>>> yt : [INFO     ] 2014-06-08 15:12:05,540 Parameters: current_time
>>>        = 0.0
>>> yt : [INFO     ] 2014-06-08 15:12:05,540 Parameters: domain_dimensions
>>>       = [2 2 2]
>>> yt : [INFO     ] 2014-06-08 15:12:05,541 Parameters: domain_left_edge
>>>        = [ 0.  0.  0.]
>>> yt : [INFO     ] 2014-06-08 15:12:05,542 Parameters: domain_right_edge
>>>       = [ 100.  100.  100.]
>>> yt : [INFO     ] 2014-06-08 15:12:05,542 Parameters:
>>> cosmological_simulation   = 0.0
>>> yt : [INFO     ] 2014-06-08 15:12:05,548 Allocating for 1.074e+09
>>> particles
>>> yt : [INFO     ] 2014-06-08 15:16:16,195 Identified 7.584e+07 octs
>>> yt : [INFO     ] 2014-06-08 15:16:16,299 Loading field plugins.
>>> yt : [INFO     ] 2014-06-08 15:16:16,299 Loaded angular_momentum (8 new
>>> fields)
>>> yt : [INFO     ] 2014-06-08 15:16:16,299 Loaded astro (14 new fields)
>>> yt : [INFO     ] 2014-06-08 15:16:16,300 Loaded cosmology (20 new fields)
>>> yt : [INFO     ] 2014-06-08 15:16:16,300 Loaded fluid (56 new fields)
>>> yt : [INFO     ] 2014-06-08 15:16:16,301 Loaded fluid_vector (88 new
>>> fields)
>>> yt : [INFO     ] 2014-06-08 15:16:16,301 Loaded geometric (103 new
>>> fields)
>>> yt : [INFO     ] 2014-06-08 15:16:16,301 Loaded local (103 new fields)
>>> yt : [INFO     ] 2014-06-08 15:16:16,302 Loaded magnetic_field (109 new
>>> fields)
>>> yt : [INFO     ] 2014-06-08 15:16:16,302 Loaded species (109 new fields)
>>>
>>> ---------------------------------------------------------------------------
>>> MemoryError                               Traceback (most recent call
>>> last)
>>>  /nfs/blank/h4231/bgriffen/data/lib/yt-x86_64/lib/python2.7/site-packages/IPython/utils/py3compat.pyc
>>> in execfile(fname, *where)
>>>     202             else:
>>>     203                 filename = fname
>>> --> 204             __builtin__.execfile(filename, *where)
>>>
>>> /nfs/blank/h4231/bgriffen/work/projects/caterpillar/c2ray/cic/ytcic.py
>>> in <module>()
>>>      99         slc.set_figure_size(4)
>>>     100         slc.save()
>>> --> 101
>>>     102     for ndim in dimlist:
>>>     103         print
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>> in __getitem__(self, key)
>>>     218                 return self.field_data[f]
>>>     219             else:
>>> --> 220                 self.get_data(f)
>>>     221         # fi.units is the unit expression string. We depend on
>>> the registry
>>>     222         # hanging off the dataset to define this unit object.
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>> in get_data(self, fields)
>>>     627
>>>     628         fields_to_generate += gen_fluids + gen_particles
>>> --> 629         self._generate_fields(fields_to_generate)
>>>     630
>>>     631     def _generate_fields(self, fields_to_generate):
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>> in _generate_fields(self, fields_to_generate)
>>>     644                 fi = self.pf._get_field_info(*field)
>>>     645                 try:
>>> --> 646                     fd = self._generate_field(field)
>>>     647                     if type(fd) == np.ndarray:
>>>     648                         fd = self.pf.arr(fd, fi.units)
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>> in _generate_field(self, field)
>>>     255                 tr = self._generate_particle_field(field)
>>>     256             else:
>>> --> 257                 tr = self._generate_fluid_field(field)
>>>     258             if tr is None:
>>>     259                 raise YTCouldNotGenerateField(field, self.pf)
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>> in _generate_fluid_field(self, field)
>>>     273             finfo.check_available(gen_obj)
>>>     274         except NeedsGridType as ngt_exception:
>>> --> 275             rv = self._generate_spatial_fluid(field,
>>> ngt_exception.ghost_zones)
>>>     276         else:
>>>     277             rv = finfo(gen_obj)
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>> in _generate_spatial_fluid(self, field, ngz)
>>>     289                     o = self._current_chunk.objs[0]
>>>     290                     with o._activate_cache():
>>> --> 291                         ind += o.select(self.selector,
>>> self[field], rv, ind)
>>>     292         else:
>>>     293             chunks = self.index._chunk(self, "spatial", ngz =
>>> ngz)
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>> in __getitem__(self, key)
>>>     218                 return self.field_data[f]
>>>     219             else:
>>> --> 220                 self.get_data(f)
>>>     221         # fi.units is the unit expression string. We depend on
>>> the registry
>>>     222         # hanging off the dataset to define this unit object.
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>> in get_data(self, fields)
>>>     627
>>>     628         fields_to_generate += gen_fluids + gen_particles
>>> --> 629         self._generate_fields(fields_to_generate)
>>>     630
>>>     631     def _generate_fields(self, fields_to_generate):
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>> in _generate_fields(self, fields_to_generate)
>>>     644                 fi = self.pf._get_field_info(*field)
>>>     645                 try:
>>> --> 646                     fd = self._generate_field(field)
>>>     647                     if type(fd) == np.ndarray:
>>>     648                         fd = self.pf.arr(fd, fi.units)
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>> in _generate_field(self, field)
>>>     255                 tr = self._generate_particle_field(field)
>>>     256             else:
>>> --> 257                 tr = self._generate_fluid_field(field)
>>>     258             if tr is None:
>>>     259                 raise YTCouldNotGenerateField(field, self.pf)
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>> in _generate_fluid_field(self, field)
>>>     275             rv = self._generate_spatial_fluid(field,
>>> ngt_exception.ghost_zones)
>>>     276         else:
>>> --> 277             rv = finfo(gen_obj)
>>>     278         return rv
>>>     279
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/fields/derived_field.pyc
>>> in __call__(self, data)
>>>     176                 "for %s" % (self.name,))
>>>     177         with self.unit_registry(data):
>>> --> 178             dd = self._function(self, data)
>>>     179         for field_name in data.keys():
>>>     180             if field_name not in original_fields:
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/fields/particle_fields.pyc
>>> in particle_cic(field, data)
>>>     113     def particle_cic(field, data):
>>>     114         pos = data[ptype, coord_name]
>>> --> 115         d = data.deposit(pos, [data[ptype, mass_name]], method =
>>> "cic")
>>>     116         d = data.apply_units(d, data[ptype, mass_name].units)
>>>     117         d /= data["index", "cell_volume"]
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/octree_subset.pyc
>>> in deposit(self, positions, fields, method)
>>>     167         mylog.debug("Depositing %s (%s^3) particles into %s
>>> Octs",
>>>     168             positions.shape[0], positions.shape[0]**0.3333333,
>>> nvals[-1])
>>> --> 169         pos =
>>> np.asarray(positions.convert_to_units("code_length"),
>>>     170                          dtype="float64")
>>>     171         # We should not need the following if we know in advance
>>> all our fields
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/units/yt_array.pyc in
>>> convert_to_units(self, units)
>>>     366
>>>     367         self.units = new_units
>>> --> 368         self *= conversion_factor
>>>     369         return self
>>>     370
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/units/yt_array.pyc in
>>> __imul__(self, other)
>>>     667         """ See __mul__. """
>>>     668         oth = sanitize_units_mul(self, other)
>>> --> 669         return np.multiply(self, oth, out=self)
>>>     670
>>>     671     def __div__(self, right_object):
>>>
>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/units/yt_array.pyc in
>>> __array_wrap__(self, out_arr, context)
>>>     966                 # casting to YTArray avoids creating a
>>> YTQuantity with size > 1
>>>     967                 return YTArray(np.array(out_arr, unit))
>>> --> 968             return ret_class(np.array(out_arr), unit)
>>>     969
>>>     970
>>>
>>> MemoryError:
>>>
>>>
>>> On Sun, Jun 8, 2014 at 10:50 PM, Matthew Turk <matthewturk at gmail.com>
>>> wrote:
>>>
>>>> Hi Brendan,
>>>>
>>>> On Sun, Jun 8, 2014 at 9:21 PM, Brendan Griffen
>>>> <brendan.f.griffen at gmail.com> wrote:
>>>> > Hi Matt,
>>>> >
>>>> > Thanks for your detailed email. Forgive my naivety but why do you
>>>> need the
>>>> > oct-tree in the first place? I have a my own fortran code for
>>>> constructing a
>>>> > cloud in cell mesh and it uses very little overhead (just the n^3
>>>> grid and
>>>> > the particle data itself). I then calculate the dx,dy,dzs to the
>>>> nearest 8
>>>> > grid points and distribute accordingly in a omp loop which is done in
>>>> a
>>>> > fraction of a second. Does the situation with yt come about (oct tree
>>>> etc.)
>>>> > necessarily because of the way it handles particle data? Is it
>>>> essentially
>>>> > used to map the particles to domains in the grid or something?
>>>>
>>>> That's not naive at all.  There are two reasons --
>>>>
>>>> 1) The octree is used for indexing for neighbor lookups and
>>>> early-termination of region selection for particles
>>>> 2) The octree is used to estimate the "required resolution" for any
>>>> operation that requires a space-filling value.  (i.e., any time that a
>>>> particle becomes a volume.)
>>>>
>>>> Projections in yt are adaptive, in that they project down to the
>>>> finest appropriate resolution.  There's also the "arbitrary_grid"
>>>> operation, which does precisely what you're describing, but as it
>>>> stands right now the octree gets constructed at time of instantiation
>>>> of the indexing system.  Thinking it over, you may be able to avoid
>>>> that completely by not using load_particles and instead using
>>>> load_uniform_grid and supplying your desired dimensions.  The field
>>>> names should be the same.
>>>>
>>>> >
>>>> > The machine I max memory on has 128GB and the snapshots are using
>>>> 1024^3
>>>> > particles. Do you have any idea of how much memory the oct-tree uses
>>>> as a
>>>> > function of particle/grid number? I am going to try on a 256GB machine
>>>> > (though this is a bit of a hassle). I'll see how I go.
>>>>
>>>> I am disappointed that it's blowing out your RAM.  This week I will
>>>> try to get some memory profiling done.  Could you file a bug to this
>>>> effect, which will help me track it?  Peak memory usage during
>>>> indexing should only be 64 bits * Nparticles, unless you're using
>>>> load_particles, in which case all the fields will *also* have to be in
>>>> memory.  It's about 8 gigabytes per field.  So, I think there's
>>>> something going wrong.
>>>>
>>>> >
>>>> > Thanks.
>>>> >
>>>> > Brendan
>>>> >
>>>> >
>>>> > On Sun, Jun 8, 2014 at 6:25 PM, Matthew Turk <matthewturk at gmail.com>
>>>> wrote:
>>>> >>
>>>> >> Hi all,
>>>> >>
>>>> >> I feel like I owe a brief explanation of why things are tricky right
>>>> >> now, what we're planning on doing, and how we're experimenting and
>>>> >> developing.
>>>> >>
>>>> >> Presently, the particle geometry handlers build a single mesh from
>>>> all
>>>> >> particles in the dataset, along with a coarse bitmask that correlates
>>>> >> files to regions in the domain.  This requires the allocation of a
>>>> >> single int64 array of size Nparticles, which is sorted in place and
>>>> >> then fed into an octree construction algorithm that then spits back
>>>> >> out the mesh.  Each octree component contains 3 64-bit integers and
>>>> >> eitehr a void pointer or a pointer to eight other octs.  Increasing
>>>> >> n_ref decreases the number of octs in this mesh; when smoothing
>>>> >> operaitons are conducted, a second "index" mesh is created for
>>>> looking
>>>> >> up particles near mesh points.  Mesh points are used for adaptive
>>>> >> resolution smoothing and other "deposit particles on the grid
>>>> somehow"
>>>> >> operations (including SPH kernel).
>>>> >>
>>>> >> Anyway, because right now it requires a global mesh to be
>>>> constructed,
>>>> >> this is expensive and requires holding a 64-bit integer in memory for
>>>> >> each particle.  I think if you're loading the particles in
>>>> differently
>>>> >> there is some additional overhead as well, but I'm still a bit
>>>> >> surprised you OOM on a 1024^3 dataset.
>>>> >>
>>>> >> In general, we don't *need* this global mesh; is can be constructed
>>>> as
>>>> >> required, which would speed up both the initial index phase as well
>>>> as
>>>> >> the final meshing process.  I got about 50% of the way to
>>>> implementing
>>>> >> this last fall, but because of various concerns and deadlines I
>>>> >> haven't finished it.  I intend to get back to it probably in July,
>>>> >> right after we put out a 3.0, so that we can have it in time for 3.1.
>>>> >> In principle this will make the particle codes much more similar to
>>>> >> ARTIO, in that the mesh will be constructed only as required and
>>>> >> discarded when no longer required, which will make them much more
>>>> >> memory efficient.
>>>> >>
>>>> >> But, getting a single mesh for extremely large data is a very high
>>>> >> priority; right now for the 10240^3 run we've been loading up
>>>> >> individual sub-chunks, which I want to stop doing.
>>>> >>
>>>> >> From the technical perspective, these are the things that need to
>>>> >> happen on the yt side for particle datasets to move to this "lazy"
>>>> >> mode of loading; most of this is based on things learned from 2HOT
>>>> and
>>>> >> ARTIO, and will involve converting to a forest-of-octrees.
>>>> >>
>>>> >>  * Split into spatially-organized subchunks of ParticleOctreeSubset
>>>> >> objects, such that these map 1:Nfiles, and that can be constructed on
>>>> >> the fly.
>>>> >>  * Construct a dual-mesh of the bitmask "ParticleRegion" object that
>>>> >> will help with identifying neighbors to a given oct cell, so that if
>>>> >> we're inside one octree we know which neighbor octrees to grab if we
>>>> >> need particles for smoothing things (fast boundary particle
>>>> >> identification is later down the road)
>>>> >>  * Parallel sort of particles, or using the parallel ring function;
>>>> >> may not be necessary after all
>>>> >>
>>>> >> All of this is doable, and I'd be happy to work with people if they'd
>>>> >> like to take a shot at implementing it, but I've mostly put it on my
>>>> >> list for post-3.0.
>>>> >>
>>>> >> -Matt
>>>> >>
>>>> >> On Sun, Jun 8, 2014 at 2:43 PM, Nathan Goldbaum <
>>>> nathan12343 at gmail.com>
>>>> >> wrote:
>>>> >> >
>>>> >> >
>>>> >> >
>>>> >> > On Sun, Jun 8, 2014 at 12:27 PM, Brendan Griffen
>>>> >> > <brendan.f.griffen at gmail.com> wrote:
>>>> >> >>
>>>> >> >> Also, how do I construct just a zero filled yt array with
>>>> dimensions
>>>> >> >> (ndim,ndim,ndim)? Thanks
>>>> >> >
>>>> >> >
>>>> >> >
>>>> >> > from yt import YTArray
>>>> >> > from numpy import np
>>>> >> >
>>>> >> > arr = YTArray(np.zeros([ndim, ndim, ndim]),
>>>> input_units=units_string)
>>>> >> >
>>>> >> > or alternatively:
>>>> >> >
>>>> >> > from yt.units import kiloparsec
>>>> >> >
>>>> >> > arr = kiloparsec*np.zeros([ndim, ndim, ndim])
>>>> >> >
>>>> >> > it doesn't have to be kiloparsec - you can compose the units you
>>>> want
>>>> >> > out of
>>>> >> > any of the unit symbols that live in yt.units.
>>>> >> >
>>>> >> > See this page for a ton more detail about yt's new unit system:
>>>> >> > http://yt-project.org/docs/dev-3.0/analyzing/units/index.html
>>>> >> >
>>>> >> >>
>>>> >> >>
>>>> >> >> Brendan
>>>> >> >>
>>>> >> >>
>>>> >> >> On Sun, Jun 8, 2014 at 3:26 PM, Brendan Griffen
>>>> >> >> <brendan.f.griffen at gmail.com> wrote:
>>>> >> >>>
>>>> >> >>> Hi,
>>>> >> >>>
>>>> >> >>> Since I get memory errors. Could I not just read in the blocks
>>>> of the
>>>> >> >>> output individually then basically stack the mesh each time.
>>>> That way
>>>> >> >>> not
>>>> >> >>> every single particle of the snapshot has to be loaded at the
>>>> same
>>>> >> >>> time.
>>>> >> >>> Would that just be a case of doing
>>>> >> >>>
>>>> >> >>> level = int(math.log(ndim,2))
>>>> >> >>> cg = ds.covering_grid(level=level,
>>>> >> >>> left_edge=[0,0,0],dims=[ndim,ndim,ndim])
>>>> >> >>> arr = cg['deposit', 'all_density']
>>>> >> >>> arrall += arr
>>>> >> >>>
>>>> >> >>> in a loop over each HDF5 block?
>>>> >> >
>>>> >> >
>>>> >> > It's likely that the memory use is dominated by the octree rather
>>>> than
>>>> >> > the
>>>> >> > covering grid.  This is with 1024^3 particles, correct?
>>>> >> >
>>>> >> > You can probably significantly reduce the memory used by the
>>>> octree by
>>>> >> > increasing n_ref in the call to load_particles.
>>>> >> >
>>>> >> > See this page for more detail about load_particles:
>>>> >> >
>>>> >> >
>>>> http://yt-project.org/docs/dev-3.0/examining/loading_data.html#generic-particle-data
>>>> >> >
>>>> >> > Larger n_ref means fewer octree cells (lower resolution), but it
>>>> also
>>>> >> > means
>>>> >> > lower poisson noise and lower memory use.
>>>> >> >
>>>> >> > Alternatively, as Matt suggested, you could break your 1024^3
>>>> ensemble
>>>> >> > of
>>>> >> > particles up into chunks, loop over the chunk, creating a particle
>>>> >> > octree
>>>> >> > and then a covering grid for each subset of the particles.  Your
>>>> final
>>>> >> > covering grid is just the sub of the covering grids for each
>>>> subset of
>>>> >> > particles.
>>>> >> >
>>>> >> >>>
>>>> >> >>>
>>>> >> >>> Thanks.
>>>> >> >>> Brendan
>>>> >> >>>
>>>> >> >>>
>>>> >> >>>
>>>> >> >>>
>>>> >> >>> On Fri, Jun 6, 2014 at 7:26 PM, Matthew Turk <
>>>> matthewturk at gmail.com>
>>>> >> >>> wrote:
>>>> >> >>>>
>>>> >> >>>>
>>>> >> >>>> On Jun 6, 2014 4:54 PM, "Brendan Griffen"
>>>> >> >>>> <brendan.f.griffen at gmail.com>
>>>> >> >>>> wrote:
>>>> >> >>>> >
>>>> >> >>>> > OK great. It is very low resolution but it worked. Thanks for
>>>> all
>>>> >> >>>> > your
>>>> >> >>>> > help. My higher resolution run 1024^3 in 100 Mpc seems to
>>>> crash on
>>>> >> >>>> > 128GB
>>>> >> >>>> > memory machine. I might have to look elsewhere.
>>>> >> >>>> >
>>>> >> >>>>
>>>> >> >>>> If you are looking for  low resolution extraction you can tune
>>>> the
>>>> >> >>>> memory usage by changing the parameter n_ref to something
>>>> higher.
>>>> >> >>>>
>>>> >> >>>> Supporting extremely large datasets in a single mesh is on the
>>>> >> >>>> roadmap
>>>> >> >>>> for the late summer or fall, after a 3.0 release goes out. For
>>>> now
>>>> >> >>>> you can
>>>> >> >>>> also extract before you load in; this is sort of how we are
>>>> >> >>>> supporting an
>>>> >> >>>> INCITE project with very large particle counts.
>>>> >> >>>>
>>>> >> >>>>
>>>> >> >>>> > Also, I normally use Canopy distribution but I just use an
>>>> alias to
>>>> >> >>>> > loadyt which erases my PYTHONPATH and I can't access scipy
>>>> and a
>>>> >> >>>> > few other
>>>> >> >>>> > libraries any more. What is the best practice here? Should I
>>>> just
>>>> >> >>>> > manually
>>>> >> >>>> > export PYTHONPATH and point to the libraries need in canopy
>>>> or can
>>>> >> >>>> > they play
>>>> >> >>>> > nice together?
>>>> >> >>>> >
>>>> >> >>>> > Thanks.
>>>> >> >>>> >
>>>> >> >>>> > BG
>>>> >> >>>> >
>>>> >> >>>> >
>>>> >> >>>> > On Fri, Jun 6, 2014 at 2:54 PM, Nathan Goldbaum
>>>> >> >>>> > <nathan12343 at gmail.com> wrote:
>>>> >> >>>> >>
>>>> >> >>>> >>
>>>> >> >>>> >>
>>>> >> >>>> >>
>>>> >> >>>> >> On Fri, Jun 6, 2014 at 11:48 AM, Brendan Griffen
>>>> >> >>>> >> <brendan.f.griffen at gmail.com> wrote:
>>>> >> >>>> >>>
>>>> >> >>>> >>> OK great. Thanks. I just wanted a homogeneous mesh. 512^3
>>>> with no
>>>> >> >>>> >>> nesting of any kind. Though when I plot the image it looks
>>>> like
>>>> >> >>>> >>> it is
>>>> >> >>>> >>> assigning particles incorrectly (low resolution on the
>>>> outside).
>>>> >> >>>> >>> This is
>>>> >> >>>> >>> just a test image.
>>>> >> >>>> >>>
>>>> >> >>>> >>
>>>> >> >>>> >> The SlicePlot is visualizing the octree so there is less
>>>> >> >>>> >> resolution
>>>> >> >>>> >> where there are fewer particles. If you want to visualize the
>>>> >> >>>> >> covering grid
>>>> >> >>>> >> you're going to need to visualize that separately.
>>>> >> >>>> >>
>>>> >> >>>> >>>
>>>> >> >>>> >>> ds = yt.load_particles(data, length_unit=3.08e24,
>>>> >> >>>> >>> mass_unit=1.9891e33,bbox=bbox)
>>>> >> >>>> >>>
>>>> >> >>>> >>> ad = ds.all_data()
>>>> >> >>>> >>> print ad['deposit', 'all_cic']
>>>> >> >>>> >>> slc = yt.SlicePlot(ds, 2, ('deposit', 'all_cic'))
>>>> >> >>>> >>> slc.set_figure_size(4)
>>>> >> >>>> >>> cg = ds.covering_grid(level=9,
>>>> >> >>>> >>> left_edge=[0,0,0],dims=[512,512,512])
>>>> >> >>>> >>>
>>>> >> >>>> >>
>>>> >> >>>> >> To actually produce the uniform resolution ndarray, you're
>>>> going
>>>> >> >>>> >> to
>>>> >> >>>> >> need to do something like:
>>>> >> >>>> >>
>>>> >> >>>> >> array = cg[('deposit', 'all_cic')]
>>>> >> >>>> >>
>>>> >> >>>> >> array will then be a 3D array you can do whatever you want
>>>> with.
>>>> >> >>>> >> By
>>>> >> >>>> >> default it has units, but to strip them off you'll just need
>>>> to
>>>> >> >>>> >> cast to
>>>> >> >>>> >> ndarray:
>>>> >> >>>> >>
>>>> >> >>>> >> array_without_units = array.v
>>>> >> >>>> >>
>>>> >> >>>> >>
>>>> >> >>>> >>>
>>>> >> >>>> >>> Also, is there a way to load multiple particle types?
>>>> >> >>>> >>>
>>>> >> >>>> >>> Do I just need to stack the particles into the array here?
>>>> >> >>>> >>>
>>>> >> >>>> >>> data = {'particle_position_x': pos[:,0],
>>>> >> >>>> >>>         'particle_position_y': pos[:,1],
>>>> >> >>>> >>>         'particle_position_z': pos[:,2],
>>>> >> >>>> >>>         'particle_mass': np.array([mpart]*npart)}
>>>> >> >>>> >>>
>>>> >> >>>> >>> Then feed it in as usual?
>>>> >> >>>> >>
>>>> >> >>>> >>
>>>> >> >>>> >> That's right, although if the particle masses are different
>>>> for
>>>> >> >>>> >> the
>>>> >> >>>> >> different particle types that code snippet will need to be
>>>> >> >>>> >> generalized to
>>>> >> >>>> >> handle that.
>>>> >> >>>> >>
>>>> >> >>>> >> I think in principle it should be possible to make
>>>> load_particles
>>>> >> >>>> >> handle different particle types just like an SPH dataset that
>>>> >> >>>> >> contains
>>>> >> >>>> >> multiple particle types, but right now that hasn't been
>>>> >> >>>> >> implemented yet.
>>>> >> >>>> >>
>>>> >> >>>> >>>
>>>> >> >>>> >>>
>>>> >> >>>> >>> Brendan
>>>> >> >>>> >>>
>>>> >> >>>> >>>
>>>> >> >>>> >>> On Thu, Jun 5, 2014 at 9:44 PM, Nathan Goldbaum
>>>> >> >>>> >>> <nathan12343 at gmail.com> wrote:
>>>> >> >>>> >>>>
>>>> >> >>>> >>>> That's right, you can set that via the bbox keyword
>>>> parameter
>>>> >> >>>> >>>> for
>>>> >> >>>> >>>> load_particles.  I'd urge you to take a look at the
>>>> docstrings
>>>> >> >>>> >>>> and source
>>>> >> >>>> >>>> code for load_particles.
>>>> >> >>>> >>>>
>>>> >> >>>> >>>>
>>>> >> >>>> >>>> On Thu, Jun 5, 2014 at 6:34 PM, Brendan Griffen
>>>> >> >>>> >>>> <brendan.f.griffen at gmail.com> wrote:
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>> Thanks very much Nathan. I tried to load in my own data
>>>> but I
>>>> >> >>>> >>>>> think there are too many particles or I have to
>>>> specifically
>>>> >> >>>> >>>>> set the domain
>>>> >> >>>> >>>>> size.
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>> In this area:
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>> data = {'particle_position_x': pos[:,0],
>>>> >> >>>> >>>>>         'particle_position_y': pos[:,1],
>>>> >> >>>> >>>>>         'particle_position_z': pos[:,2],
>>>> >> >>>> >>>>>         'particle_mass': np.array([mpart]*npart)}
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>> ds = yt.load_particles(data, length_unit=3.08e24,
>>>> >> >>>> >>>>> mass_unit=1.9891e36)
>>>> >> >>>> >>>>> ad = ds.all_data()
>>>> >> >>>> >>>>> print ad['deposit', 'all_cic']
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>> In [3]: run ytcic.py
>>>> >> >>>> >>>>> yt : [INFO     ] 2014-06-05 21:29:06,183 Parameters:
>>>> >> >>>> >>>>> current_time
>>>> >> >>>> >>>>> = 0.0
>>>> >> >>>> >>>>> yt : [INFO     ] 2014-06-05 21:29:06,183 Parameters:
>>>> >> >>>> >>>>> domain_dimensions         = [2 2 2]
>>>> >> >>>> >>>>> yt : [INFO     ] 2014-06-05 21:29:06,184 Parameters:
>>>> >> >>>> >>>>> domain_left_edge          = [ 0.  0.  0.]
>>>> >> >>>> >>>>> yt : [INFO     ] 2014-06-05 21:29:06,185 Parameters:
>>>> >> >>>> >>>>> domain_right_edge         = [ 1.  1.  1.]
>>>> >> >>>> >>>>> yt : [INFO     ] 2014-06-05 21:29:06,185 Parameters:
>>>> >> >>>> >>>>> cosmological_simulation   = 0.0
>>>> >> >>>> >>>>> yt : [INFO     ] 2014-06-05 21:29:06,188 Allocating for
>>>> >> >>>> >>>>> 1.342e+08
>>>> >> >>>> >>>>> particles
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> ---------------------------------------------------------------------------
>>>> >> >>>> >>>>> YTDomainOverflow                          Traceback (most
>>>> >> >>>> >>>>> recent
>>>> >> >>>> >>>>> call last)
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /nfs/blank/h4231/bgriffen/data/lib/yt-x86_64/lib/python2.7/site-packages/IPython/utils/py3compat.pyc
>>>> >> >>>> >>>>> in execfile(fname, *where)
>>>> >> >>>> >>>>>     202             else:
>>>> >> >>>> >>>>>     203                 filename = fname
>>>> >> >>>> >>>>> --> 204             __builtin__.execfile(filename, *where)
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /nfs/blank/h4231/bgriffen/work/projects/caterpillar/c2ray/cic/ytcic.py in
>>>> >> >>>> >>>>> <module>()
>>>> >> >>>> >>>>>      52
>>>> >> >>>> >>>>>      53 ad = ds.all_data()
>>>> >> >>>> >>>>> ---> 54 print ad['deposit', 'all_cic']
>>>> >> >>>> >>>>>      55 slc = yt.SlicePlot(ds, 2, ('deposit', 'all_cic'))
>>>> >> >>>> >>>>>      56 slc.set_figure_size(4)
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>>> >> >>>> >>>>> in __getitem__(self, key)
>>>> >> >>>> >>>>>     205         Returns a single field.  Will add if
>>>> necessary.
>>>> >> >>>> >>>>>     206         """
>>>> >> >>>> >>>>> --> 207         f = self._determine_fields([key])[0]
>>>> >> >>>> >>>>>     208         if f not in self.field_data and key not in
>>>> >> >>>> >>>>> self.field_data:
>>>> >> >>>> >>>>>     209             if f in self._container_fields:
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc
>>>> >> >>>> >>>>> in _determine_fields(self, fields)
>>>> >> >>>> >>>>>     453                     raise
>>>> YTFieldNotParseable(field)
>>>> >> >>>> >>>>>     454                 ftype, fname = field
>>>> >> >>>> >>>>> --> 455                 finfo =
>>>> self.pf._get_field_info(ftype,
>>>> >> >>>> >>>>> fname)
>>>> >> >>>> >>>>>     456             else:
>>>> >> >>>> >>>>>     457                 fname = field
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/static_output.pyc
>>>> >> >>>> >>>>> in _get_field_info(self, ftype, fname)
>>>> >> >>>> >>>>>     445     _last_finfo = None
>>>> >> >>>> >>>>>     446     def _get_field_info(self, ftype, fname =
>>>> None):
>>>> >> >>>> >>>>> --> 447         self.index
>>>> >> >>>> >>>>>     448         if fname is None:
>>>> >> >>>> >>>>>     449             ftype, fname = "unknown", ftype
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/data_objects/static_output.pyc
>>>> >> >>>> >>>>> in index(self)
>>>> >> >>>> >>>>>     277                 raise RuntimeError("You should not
>>>> >> >>>> >>>>> instantiate Dataset.")
>>>> >> >>>> >>>>>     278             self._instantiated_index =
>>>> >> >>>> >>>>> self._index_class(
>>>> >> >>>> >>>>> --> 279                 self,
>>>> dataset_type=self.dataset_type)
>>>> >> >>>> >>>>>     280             # Now we do things that we need an
>>>> >> >>>> >>>>> instantiated index for
>>>> >> >>>> >>>>>     281             # ...first off, we create our
>>>> field_info
>>>> >> >>>> >>>>> now.
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/frontends/stream/data_structures.pyc
>>>> >> >>>> >>>>> in __init__(self, pf, dataset_type)
>>>> >> >>>> >>>>>     942     def __init__(self, pf, dataset_type = None):
>>>> >> >>>> >>>>>     943         self.stream_handler = pf.stream_handler
>>>> >> >>>> >>>>> --> 944         super(StreamParticleIndex,
>>>> self).__init__(pf,
>>>> >> >>>> >>>>> dataset_type)
>>>> >> >>>> >>>>>     945
>>>> >> >>>> >>>>>     946     def _setup_data_io(self):
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/geometry/particle_geometry_handler.pyc
>>>> >> >>>> >>>>> in __init__(self, pf, dataset_type)
>>>> >> >>>> >>>>>      48         self.directory =
>>>> >> >>>> >>>>> os.path.dirname(self.index_filename)
>>>> >> >>>> >>>>>      49         self.float_type = np.float64
>>>> >> >>>> >>>>> ---> 50         super(ParticleIndex, self).__init__(pf,
>>>> >> >>>> >>>>> dataset_type)
>>>> >> >>>> >>>>>      51
>>>> >> >>>> >>>>>      52     def _setup_geometry(self):
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/geometry/geometry_handler.pyc
>>>> >> >>>> >>>>> in __init__(self, pf, dataset_type)
>>>> >> >>>> >>>>>      54
>>>> >> >>>> >>>>>      55         mylog.debug("Setting up domain geometry.")
>>>> >> >>>> >>>>> ---> 56         self._setup_geometry()
>>>> >> >>>> >>>>>      57
>>>> >> >>>> >>>>>      58         mylog.debug("Initializing data grid data
>>>> IO")
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/geometry/particle_geometry_handler.pyc
>>>> >> >>>> >>>>> in _setup_geometry(self)
>>>> >> >>>> >>>>>      52     def _setup_geometry(self):
>>>> >> >>>> >>>>>      53         mylog.debug("Initializing Particle
>>>> Geometry
>>>> >> >>>> >>>>> Handler.")
>>>> >> >>>> >>>>> ---> 54         self._initialize_particle_handler()
>>>> >> >>>> >>>>>      55
>>>> >> >>>> >>>>>      56
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/geometry/particle_geometry_handler.pyc
>>>> >> >>>> >>>>> in _initialize_particle_handler(self)
>>>> >> >>>> >>>>>      87                 pf.domain_left_edge,
>>>> >> >>>> >>>>> pf.domain_right_edge,
>>>> >> >>>> >>>>>      88                 [N, N, N], len(self.data_files))
>>>> >> >>>> >>>>> ---> 89         self._initialize_indices()
>>>> >> >>>> >>>>>      90         self.oct_handler.finalize()
>>>> >> >>>> >>>>>      91         self.max_level =
>>>> self.oct_handler.max_level
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/geometry/particle_geometry_handler.pyc
>>>> >> >>>> >>>>> in _initialize_indices(self)
>>>> >> >>>> >>>>>     109             npart =
>>>> >> >>>> >>>>> sum(data_file.total_particles.values())
>>>> >> >>>> >>>>>     110             morton[ind:ind + npart] = \
>>>> >> >>>> >>>>> --> 111
>>>> self.io._initialize_index(data_file,
>>>> >> >>>> >>>>> self.regions)
>>>> >> >>>> >>>>>     112             ind += npart
>>>> >> >>>> >>>>>     113         morton.sort()
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> /bigbang/data/bgriffen/lib/yt-x86_64/src/yt-hg/yt/frontends/stream/io.pyc in
>>>> >> >>>> >>>>> _initialize_index(self, data_file, regions)
>>>> >> >>>> >>>>>     144                 raise
>>>> YTDomainOverflow(pos.min(axis=0),
>>>> >> >>>> >>>>> pos.max(axis=0),
>>>> >> >>>> >>>>>     145
>>>> >> >>>> >>>>> data_file.pf.domain_left_edge,
>>>> >> >>>> >>>>> --> 146
>>>> >> >>>> >>>>> data_file.pf.domain_right_edge)
>>>> >> >>>> >>>>>     147             regions.add_data_file(pos,
>>>> >> >>>> >>>>> data_file.file_id)
>>>> >> >>>> >>>>>     148             morton.append(compute_morton(
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>> YTDomainOverflow: Particle bounds [ 0.  0.  0.] and [
>>>> >> >>>> >>>>> 99.99999237
>>>> >> >>>> >>>>> 99.99999237  99.99999237] exceed domain bounds [ 0.  0.
>>>>  0.]
>>>> >> >>>> >>>>> code_length and
>>>> >> >>>> >>>>> [ 1.  1.  1.] code_length
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>> On Thu, Jun 5, 2014 at 8:22 PM, Nathan Goldbaum
>>>> >> >>>> >>>>> <nathan12343 at gmail.com> wrote:
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>> Here's a worked out example that does what you're
>>>> looking for
>>>> >> >>>> >>>>>> using a fake 1 million particle dataset:
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>>
>>>> http://nbviewer.ipython.org/gist/ngoldbaum/546d37869aafe71cfe38
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>> In this notebook I make use of two key yt features:
>>>> >> >>>> >>>>>> `load_particles`, and `covering_grid`.
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>> load_particles creates a "stream" dataset based on
>>>> in-memory
>>>> >> >>>> >>>>>> data
>>>> >> >>>> >>>>>> fed in as a numpy array.  This dataset acts just like an
>>>> >> >>>> >>>>>> on-disk simulation
>>>> >> >>>> >>>>>> dataset, but doesn't come with the baggage of needing to
>>>> write
>>>> >> >>>> >>>>>> a custom
>>>> >> >>>> >>>>>> frontend to read a specific data format off disk.
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>> covering_grid is a way to generate uniform resolution
>>>> data
>>>> >> >>>> >>>>>> from
>>>> >> >>>> >>>>>> an AMR dataset. It acts like a python dictionary where
>>>> the
>>>> >> >>>> >>>>>> keys are field
>>>> >> >>>> >>>>>> names and returns 3D numpy arrays of whatever uniform
>>>> >> >>>> >>>>>> resolution you specify
>>>> >> >>>> >>>>>> when you create the covering_grid.
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>> Note that if you're using load_particles all of your data
>>>> >> >>>> >>>>>> needs
>>>> >> >>>> >>>>>> to live in memory.  If your data is too big for that
>>>> you'll
>>>> >> >>>> >>>>>> need to write a
>>>> >> >>>> >>>>>> frontend for your data format or use a memmap to an
>>>> on-disk
>>>> >> >>>> >>>>>> file somehow.
>>>> >> >>>> >>>>>> I'm not an expert on that but others on the list should
>>>> be
>>>> >> >>>> >>>>>> able to help out.
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>> Hope that gets you well on your way :)
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>> -Nathan
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>> On Thu, Jun 5, 2014 at 5:04 PM, Desika Narayanan
>>>> >> >>>> >>>>>> <dnarayan at haverford.edu> wrote:
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>> Hey Brendan,
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>> A couple of extra tools you might find helpful in
>>>> conjunction
>>>> >> >>>> >>>>>>> with Nathan's example of depositing the particles onto
>>>> an
>>>> >> >>>> >>>>>>> octree are at:
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>> http://paste.yt-project.org/show/4737/
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>> Where I load a gadget snapshot, and then recover the
>>>> >> >>>> >>>>>>> coordinates
>>>> >> >>>> >>>>>>> and width of each cell.
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>> In response to your last question - the particles are
>>>> >> >>>> >>>>>>> deposited
>>>> >> >>>> >>>>>>> into an octree grid (so, you'll see that the cell sizes
>>>> >> >>>> >>>>>>> aren't all the same
>>>> >> >>>> >>>>>>> size).   I don't know if depositing onto a regular
>>>> NxNxN mesh
>>>> >> >>>> >>>>>>> is possible,
>>>> >> >>>> >>>>>>> though would be interested to hear if so.
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>> -d
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>> On Thu, Jun 5, 2014 at 7:58 PM, Brendan Griffen
>>>> >> >>>> >>>>>>> <brendan.f.griffen at gmail.com> wrote:
>>>> >> >>>> >>>>>>>>
>>>> >> >>>> >>>>>>>> Thanks. I'll get the "bleeding edge" version first
>>>> then try
>>>> >> >>>> >>>>>>>> your suggestions. Though I want to return the NxNxN
>>>> array
>>>> >> >>>> >>>>>>>> and be able to
>>>> >> >>>> >>>>>>>> write this mesh to a file. It is *only* using the cic
>>>> part
>>>> >> >>>> >>>>>>>> of yt and it
>>>> >> >>>> >>>>>>>> should return the mesh to be written? Just wanted to
>>>> >> >>>> >>>>>>>> clarify?
>>>> >> >>>> >>>>>>>>
>>>> >> >>>> >>>>>>>> Thanks.
>>>> >> >>>> >>>>>>>> Brendan
>>>> >> >>>> >>>>>>>>
>>>> >> >>>> >>>>>>>>
>>>> >> >>>> >>>>>>>> On Thu, Jun 5, 2014 at 6:49 PM, Nathan Goldbaum
>>>> >> >>>> >>>>>>>> <nathan12343 at gmail.com> wrote:
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>> On Thu, Jun 5, 2014 at 3:36 PM, John ZuHone
>>>> >> >>>> >>>>>>>>> <jzuhone at gmail.com> wrote:
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>> Hi Brendan,
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>> Which version of yt are you using?
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>> If you're using 3.0, this is actually fairly easy.
>>>> If you
>>>> >> >>>> >>>>>>>>>> look in yt.fields.particle_fields.py, around line
>>>> 85, you
>>>> >> >>>> >>>>>>>>>> can see how this
>>>> >> >>>> >>>>>>>>>> is done for the "particle_density" and
>>>> "particle_mass"
>>>> >> >>>> >>>>>>>>>> fields. Basically you
>>>> >> >>>> >>>>>>>>>> can call a "deposit" method which takes the particle
>>>> field
>>>> >> >>>> >>>>>>>>>> quantity you want
>>>> >> >>>> >>>>>>>>>> deposited and deposits it into cells. The underlying
>>>> >> >>>> >>>>>>>>>> calculation is done
>>>> >> >>>> >>>>>>>>>> using Cython, so it's fast.
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>> And you shouldn't ever actually need to call these
>>>> >> >>>> >>>>>>>>> "deposit"
>>>> >> >>>> >>>>>>>>> functions, since "deposit" is exposed as a field type
>>>> for
>>>> >> >>>> >>>>>>>>> all datasets that
>>>> >> >>>> >>>>>>>>> contain particles.
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>> Here is a notebook that does this for Enzo AMR data:
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>
>>>> http://nbviewer.ipython.org/gist/ngoldbaum/5e19e4e6cc2bf330149c
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>> This dataset contains about a million particles and
>>>> >> >>>> >>>>>>>>> generates
>>>> >> >>>> >>>>>>>>> a CIC deposition for the whole domain in about 6
>>>> seconds
>>>> >> >>>> >>>>>>>>> from a cold start.
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>> If you're using 2.x, then you can do the same thing,
>>>> but
>>>> >> >>>> >>>>>>>>>> it's
>>>> >> >>>> >>>>>>>>>> not as straightforward. You can see how this works in
>>>> >> >>>> >>>>>>>>>> yt.data_objects.universal_fields.py, around line
>>>> 986,
>>>> >> >>>> >>>>>>>>>> where the
>>>> >> >>>> >>>>>>>>>> "particle_density" field is defined. Basically, it
>>>> calls
>>>> >> >>>> >>>>>>>>>> CICDeposit_3, which
>>>> >> >>>> >>>>>>>>>> is in yt.utilities.lib.CICDeposit.pyx.
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>> Let me know if you need any more clarification.
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>> Best,
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>> John Z
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>> On Jun 5, 2014, at 6:07 PM, Brendan Griffen
>>>> >> >>>> >>>>>>>>>> <brendan.f.griffen at gmail.com> wrote:
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>> > Hi,
>>>> >> >>>> >>>>>>>>>> >
>>>> >> >>>> >>>>>>>>>> > I was wondering if there were any Cython routines
>>>> within
>>>> >> >>>> >>>>>>>>>> > yt
>>>> >> >>>> >>>>>>>>>> > which takes particle data and converts it into a
>>>> >> >>>> >>>>>>>>>> > cloud-in-cell based mesh
>>>> >> >>>> >>>>>>>>>> > which can be written to a file of my choosing.
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>> What sort of mesh were you looking for?  yt will
>>>> internally
>>>> >> >>>> >>>>>>>>> construct an octree if it is fed particle data.  I'm
>>>> not
>>>> >> >>>> >>>>>>>>> sure whether this
>>>> >> >>>> >>>>>>>>> octree can be saved to disk for later analysis.
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>> It's also possible to create a uniform resolution
>>>> covering
>>>> >> >>>> >>>>>>>>> grid containing field data for a deposited quantity,
>>>> which
>>>> >> >>>> >>>>>>>>> can be quite
>>>> >> >>>> >>>>>>>>> easily saved to disk in a number of ways.
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>> I heard a while ago there was some such
>>>> functionality but
>>>> >> >>>> >>>>>>>>>> it
>>>> >> >>>> >>>>>>>>>> could be too far down the yt rabbit hole to be used
>>>> as a
>>>> >> >>>> >>>>>>>>>> standalone? Is this
>>>> >> >>>> >>>>>>>>>> true? I have my own Python code for doing it but it
>>>> just
>>>> >> >>>> >>>>>>>>>> isn't fast enough
>>>> >> >>>> >>>>>>>>>> and thought I'd ask the yt community if there were
>>>> any
>>>> >> >>>> >>>>>>>>>> wrapper tools
>>>> >> >>>> >>>>>>>>>> available to boost the speed.
>>>> >> >>>> >>>>>>>>>> >
>>>> >> >>>> >>>>>>>>>> > Thanks.
>>>> >> >>>> >>>>>>>>>> > Brendan
>>>> >> >>>> >>>>>>>>>> > _______________________________________________
>>>> >> >>>> >>>>>>>>>> > yt-users mailing list
>>>> >> >>>> >>>>>>>>>> > yt-users at lists.spacepope.org
>>>> >> >>>> >>>>>>>>>> >
>>>> >> >>>> >>>>>>>>>> >
>>>> >> >>>> >>>>>>>>>> >
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>> _______________________________________________
>>>> >> >>>> >>>>>>>>>> yt-users mailing list
>>>> >> >>>> >>>>>>>>>> yt-users at lists.spacepope.org
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>>
>>>> >> >>>> >>>>>>>>>>
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>> _______________________________________________
>>>> >> >>>> >>>>>>>>> yt-users mailing list
>>>> >> >>>> >>>>>>>>> yt-users at lists.spacepope.org
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>>
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>> >>>>>>>>>
>>>> >> >>>> >>>>>>>>
>>>> >> >>>> >>>>>>>>
>>>> >> >>>> >>>>>>>> _______________________________________________
>>>> >> >>>> >>>>>>>> yt-users mailing list
>>>> >> >>>> >>>>>>>> yt-users at lists.spacepope.org
>>>> >> >>>> >>>>>>>>
>>>> >> >>>> >>>>>>>>
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>> >>>>>>>>
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>> _______________________________________________
>>>> >> >>>> >>>>>>> yt-users mailing list
>>>> >> >>>> >>>>>>> yt-users at lists.spacepope.org
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>>
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>> >>>>>>>
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>> _______________________________________________
>>>> >> >>>> >>>>>> yt-users mailing list
>>>> >> >>>> >>>>>> yt-users at lists.spacepope.org
>>>> >> >>>> >>>>>>
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>> >>>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>> _______________________________________________
>>>> >> >>>> >>>>> yt-users mailing list
>>>> >> >>>> >>>>> yt-users at lists.spacepope.org
>>>> >> >>>> >>>>>
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>> >>>>>
>>>> >> >>>> >>>>
>>>> >> >>>> >>>>
>>>> >> >>>> >>>> _______________________________________________
>>>> >> >>>> >>>> yt-users mailing list
>>>> >> >>>> >>>> yt-users at lists.spacepope.org
>>>> >> >>>> >>>>
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>> >>>>
>>>> >> >>>> >>>
>>>> >> >>>> >>>
>>>> >> >>>> >>> _______________________________________________
>>>> >> >>>> >>> yt-users mailing list
>>>> >> >>>> >>> yt-users at lists.spacepope.org
>>>> >> >>>> >>>
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>> >>>
>>>> >> >>>> >>
>>>> >> >>>> >>
>>>> >> >>>> >> _______________________________________________
>>>> >> >>>> >> yt-users mailing list
>>>> >> >>>> >> yt-users at lists.spacepope.org
>>>> >> >>>> >>
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>> >>
>>>> >> >>>> >
>>>> >> >>>> >
>>>> >> >>>> > _______________________________________________
>>>> >> >>>> > yt-users mailing list
>>>> >> >>>> > yt-users at lists.spacepope.org
>>>> >> >>>> >
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>> >
>>>> >> >>>>
>>>> >> >>>>
>>>> >> >>>> _______________________________________________
>>>> >> >>>> yt-users mailing list
>>>> >> >>>> yt-users at lists.spacepope.org
>>>> >> >>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>>>
>>>> >> >>>
>>>> >> >>
>>>> >> >>
>>>> >> >> _______________________________________________
>>>> >> >> yt-users mailing list
>>>> >> >> yt-users at lists.spacepope.org
>>>> >> >> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >>
>>>> >> >
>>>> >> >
>>>> >> > _______________________________________________
>>>> >> > yt-users mailing list
>>>> >> > yt-users at lists.spacepope.org
>>>> >> > http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >> >
>>>> >> _______________________________________________
>>>> >> yt-users mailing list
>>>> >> yt-users at lists.spacepope.org
>>>> >> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >
>>>> >
>>>> >
>>>> > _______________________________________________
>>>> > yt-users mailing list
>>>> > yt-users at lists.spacepope.org
>>>> > http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>> >
>>>> _______________________________________________
>>>> yt-users mailing list
>>>> yt-users at lists.spacepope.org
>>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>>
>>>
>>>
>>> _______________________________________________
>>> yt-users mailing list
>>> yt-users at lists.spacepope.org
>>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>>
>>>
>>
>> _______________________________________________
>> yt-users mailing list
>> yt-users at lists.spacepope.org
>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>
>>
>
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20140609/f271acb3/attachment.html>


More information about the yt-users mailing list