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

Brendan Griffen brendan.f.griffen at gmail.com
Mon Jun 9 10:48:00 PDT 2014


Hi Nathan, I'm not sure if you missed it but earlier Matt asked if there
was some copying which shouldn't be happening?

Brendan


On Mon, Jun 9, 2014 at 1:22 PM, Nathan Goldbaum <nathan12343 at gmail.com>
wrote:

>
>
>
> 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
>>
>>
>
> _______________________________________________
> 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/3127120c/attachment.htm>


More information about the yt-users mailing list