[yt-users] cloud in cell mesh construction for particle data
Brendan Griffen
brendan.f.griffen at gmail.com
Sun Jun 8 12:27:42 PDT 2014
Also, how do I construct just a zero filled yt array with dimensions
(ndim,ndim,ndim)? Thanks
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?
>
> 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
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20140608/8b4b497a/attachment.htm>
More information about the yt-users
mailing list