[yt-users] cloud in cell mesh construction for particle data
Brendan Griffen
brendan.f.griffen at gmail.com
Sun Jun 8 12:26:51 PDT 2014
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/8170d289/attachment.html>
More information about the yt-users
mailing list