Mailman 3 May 2018 - yt-users

Re: Newbie question: ParticlePlot(...)
by Nick Gnedin 16 May '18

16 May '18

Nathan, Done (I tared 4 files together): /data/gnedin/REI/TMP >tar cvf gnedin_hop.tar hop.* hop.0.h5 hop.1.h5 hop.2.h5 hop.3.h5 /data/gnedin/REI/TMP >yt upload gnedin_hop.tar Stream body handler: received 9687040 bytes http://use.yt/upload/318d837f The files are available on BlueWaters at: /mnt/c/scratch/sciteam/ngnedin/PERM/B40/A.DMO/a=1.0008/hop/hop.*.h5 Thanks a lot. Nick On 05/16/2018 09:38 AM, Nathan Goldbaum wrote: > Hi, > > Is there any chance you can share the hdf5 file? I'm not able to > reproduce this warning using a test dataset I have so it might have > something to do with the data. > > You can share a dataset with "yt upload path/to/data.h5". That will > upload to a server at NCSA and then print a link out to share here once > its done. > > If the data cannot be shared publicly you can send me an e-mail off-list > with e.g. a dropbox link. > > -Nathan > > On Wed, May 16, 2018 at 10:32 AM, Nick Gnedin <ngnedin(a)gmail.com > <mailto:ngnedin@gmail.com>> wrote: > > > I am trying to make a plot of all halos found by HOP. I copied your > example from http://yt-project.org/doc/cookbook/simple_plots.html > <http://yt-project.org/doc/cookbook/simple_plots.html> > > If I just plot partile positions: > > --------------------- script ------------------- > import yt > import numpy as np > > # load the dataset > ds = yt.load('../B40/A.DMO/a=1.0008/hop/hop.0.h5') > > # create our plot > p = yt.ParticlePlot(ds, 'particle_position_x', 'particle_position_y') > > m = ds.all_data()[('all', 'particle_mass')] > print(np.amin(m),np.amax(m)) > > # pick some appropriate units > p.set_axes_unit('Mpc') > p.set_unit('particle_mass', 'Msun') > > #save result > p.save() > ----------------------------------------------- > > then everything works fine. However, if I use mass weighting: > > p = yt.ParticlePlot(ds, 'particle_position_x', > 'particle_position_y', 'particle_mass') > > I get the following error: > > ... > yt : [INFO ] 2018-05-16 09:27:39,004 Allocating for 1.997e+05 > particles (index particle type 'all') > yt : [INFO ] 2018-05-16 09:27:39,381 Identified 5.972e+04 octs > yt : [INFO ] 2018-05-16 09:27:40,021 xlim = 0.000000 > 199229564535610884216586240.000000 > yt : [INFO ] 2018-05-16 09:27:40,021 ylim = 0.000000 > 199229564535610884216586240.000000 > yt : [INFO ] 2018-05-16 09:27:40,023 xlim = 0.000000 > 199229564535610884216586240.000000 > yt : [INFO ] 2018-05-16 09:27:40,023 ylim = 0.000000 > 199229564535610884216586240.000000 > yt : [INFO ] 2018-05-16 09:27:40,042 Splatting (('all', > 'particle_mass')) onto a 800 by 800 mesh > 4.1036023233207635e+41 g 1.921697915571756e+47 g > yt : [INFO ] 2018-05-16 09:27:44,084 Saving plot > hop.0.h5_Particle_z_particle_mass.png > /u/sciteam/ngnedin/.local/lib/python3.5/site-packages/matplotlib/colors.py:1020: > RuntimeWarning: invalid value encountered in less_equal > mask |= resdat <= 0 > > Notice, that there are no nonpositive masses in the catalog > print(np.amin(m),np.amax(m)): > 4.1036023233207635e+41 g 1.921697915571756e+47 g > > > > > _______________________________________________ > yt-users mailing list -- yt-users(a)python.org > <mailto:yt-users@python.org> > To unsubscribe send an email to yt-users-leave(a)python.org > <mailto:yt-users-leave@python.org> > >

2 1

Newbie question: ParticlePlot(...)
by Nick Gnedin 16 May '18

16 May '18

I am trying to make a plot of all halos found by HOP. I copied your example from http://yt-project.org/doc/cookbook/simple_plots.html If I just plot partile positions: --------------------- script ------------------- import yt import numpy as np # load the dataset ds = yt.load('../B40/A.DMO/a=1.0008/hop/hop.0.h5') # create our plot p = yt.ParticlePlot(ds, 'particle_position_x', 'particle_position_y') m = ds.all_data()[('all', 'particle_mass')] print(np.amin(m),np.amax(m)) # pick some appropriate units p.set_axes_unit('Mpc') p.set_unit('particle_mass', 'Msun') #save result p.save() ----------------------------------------------- then everything works fine. However, if I use mass weighting: p = yt.ParticlePlot(ds, 'particle_position_x', 'particle_position_y', 'particle_mass') I get the following error: ... yt : [INFO ] 2018-05-16 09:27:39,004 Allocating for 1.997e+05 particles (index particle type 'all') yt : [INFO ] 2018-05-16 09:27:39,381 Identified 5.972e+04 octs yt : [INFO ] 2018-05-16 09:27:40,021 xlim = 0.000000 199229564535610884216586240.000000 yt : [INFO ] 2018-05-16 09:27:40,021 ylim = 0.000000 199229564535610884216586240.000000 yt : [INFO ] 2018-05-16 09:27:40,023 xlim = 0.000000 199229564535610884216586240.000000 yt : [INFO ] 2018-05-16 09:27:40,023 ylim = 0.000000 199229564535610884216586240.000000 yt : [INFO ] 2018-05-16 09:27:40,042 Splatting (('all', 'particle_mass')) onto a 800 by 800 mesh 4.1036023233207635e+41 g 1.921697915571756e+47 g yt : [INFO ] 2018-05-16 09:27:44,084 Saving plot hop.0.h5_Particle_z_particle_mass.png /u/sciteam/ngnedin/.local/lib/python3.5/site-packages/matplotlib/colors.py:1020: RuntimeWarning: invalid value encountered in less_equal mask |= resdat <= 0 Notice, that there are no nonpositive masses in the catalog print(np.amin(m),np.amax(m)): 4.1036023233207635e+41 g 1.921697915571756e+47 g

2 1

reading in Fenics on point tetrahedral ouputs to yt
by Steven Armour 15 May '18

15 May '18

FENICS is one of the most popular FEM solvers and considered one of the best FEM solvers for Python as well as being a sibling NumFocues project. The outputs of fenics is the dolfin-XML or .pvd/.vtu . And thus Fenics has traditionally relied on VTK to view its outputs. I would like to, however, move from the nightmare install that is VTK to using yt to view the output. The issue is that fenics only supports tetrahedral meshes with on mesh point data. Whereas from reading http://yt-project.org/doc/examining/loading_data.html sections: Generic Array Data, Semi-Structured Grid Data as well as http://yt-project.org/doc/visualizing/unstructured_mesh_rendering.html that yt does or does not load tetrahedral data and does not load on mesh point data but instead cell-centered data. I have tried using the meshio library to convert from .vtu to .e (exodus II) with no success. There is also an example notebook based on one of the examples from the fenics tutorials here https://github.com/GProtoZeroW/IssueNotebooks/blob/master/FENICS%20YT.ipynb So what would be ideal is to get a loader method that reads in a fenics object to yt to bypass any writing to a file so that one can get visulisations like what has been done with the MOOSE ex (http://yt-project.org/doc/visualizing/unstructured_mesh_rendering.html) but with FENICS

2 1

yt.analysis_modules.halo_analysis.enzofof_merger_tree module
by lianxiaoli87＠sjtu.edu.cn 14 May '18

14 May '18

Dear all: Recently, I try to get the merger-tree for my Enzo simulation. The halo catalogs can be obtained by HaloCatalog or FOFHaloFinder correctly. But it will be wrong when I use these halo lists and catalogs to get their merger tree. Here is my code and the Error-message: Code: ts = yt.load("DD????/DD????") for ds in ts: halo_list = FOFHaloFinder(ds) i = int(ds.basename[2:]) halo_list.write_out("FOF/groups_%05i.txt" % i) halo_list.write_particle_lists("FOF/particles_%05i" % i) mt = EnzoFOFMergerTree(external_FOF=False) for i in range(20): mt.build_tree(i) mt.save_halo_evolution('halos.h5') for i in range(20): plot_halo_evolution('halos.h5', i) Error Message: TypeError: '<' not supported between instances of 'str' and 'int'. This error is from the commander of mt.build_tree(i), for the other commanders can work well. I try to figure it out and check the yt code, but it still show me this wrong message. Hopefully for your help! Xiaoli Lian

2 9

Re: Incorrect particle_identifier in HOP for 2048^3 runs
by Nick Gnedin 10 May '18

10 May '18

Matt, Thanks. A simple workaround is not to use the 'particle_identifier' field at all (assuming the halo finder works correctly). Perhaps, it can be deleted or flagged as "do not use". n On 5/10/2018 2:53 PM, Matthew Turk wrote: > Hi Nick, you're right. It looks like there's a cast-to-float that > shouldn't be happening. I'm going to try to track this down, but I think > I've heard of other folks addressing this in the past in other frontends, > so I'll start by asking around. > On Thu, May 10, 2018 at 2:33 PM Nick Gnedin <gnedin(a)fnal.gov> wrote: > > >> I am running HOP halo finder for several large simulations like this: >> import yt >> from yt.analysis_modules.halo_analysis.api import HaloCatalog >> yt.enable_parallelism() >> path = "/mnt/c/scratch/sciteam/ngnedin/PERM/B40/B" >> d = yt.load(path+"/rei40/rei40.art") >> hc = > > HaloCatalog(data_ds=d,finder_method='hop',output_dir=path+"/hop",finder_kwargs={"dm_only":False,"ptype":"N-BODY"}) >> hc.create() > >> For 1024^3 and smaller simulation the halo catalog is created correctly, >> but for 2048^3 runs the 'particle_identifier' field is saved in the out >> pout incorrectly. For example: > >> # 1024^3 >> >>> h = yt.load("../B40/A/a=0.1600/hop/hop.0.h5") >> >>> d = h.all_data()[('halos','particle_identifier')] >> >>> d.sort() >> >>> print(len(d),d[-1]) >> 293749 293748.0 dimensionless >> >>> print(d[0:10]) >> [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.] dimensionless > >> # 2048^3 >> >>> h = yt.load("../B80/A/a=0.1606/hop/hop.0.h5") >> >>> d = h.all_data()[('halos','particle_identifier')] >> >>> d.sort() >> >>> print(len(d),d[-1]) >> 1974910 43782.0 dimensionless >> >>> print(d[0:10]) >> [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] dimensionless > >> It looks like the identifier gets wrapped somehow, even if the number of >> halos in the 2048^3 run is less than 2M. > >> Nick Gnedin >> _______________________________________________ >> yt-users mailing list -- yt-users(a)python.org >> To unsubscribe send an email to yt-users-leave(a)python.org > > --- > This email has been checked for viruses by AVG. > http://www.avg.com >

1 0

Incorrect particle_identifier in HOP for 2048^3 runs
by Nick Gnedin 10 May '18

10 May '18

I am running HOP halo finder for several large simulations like this: import yt from yt.analysis_modules.halo_analysis.api import HaloCatalog yt.enable_parallelism() path = "/mnt/c/scratch/sciteam/ngnedin/PERM/B40/B" d = yt.load(path+"/rei40/rei40.art") hc = HaloCatalog(data_ds=d,finder_method='hop',output_dir=path+"/hop",finder_kwargs={"dm_only":False,"ptype":"N-BODY"}) hc.create() For 1024^3 and smaller simulation the halo catalog is created correctly, but for 2048^3 runs the 'particle_identifier' field is saved in the out pout incorrectly. For example: # 1024^3 >>> h = yt.load("../B40/A/a=0.1600/hop/hop.0.h5") >>> d = h.all_data()[('halos','particle_identifier')] >>> d.sort() >>> print(len(d),d[-1]) 293749 293748.0 dimensionless >>> print(d[0:10]) [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.] dimensionless # 2048^3 >>> h = yt.load("../B80/A/a=0.1606/hop/hop.0.h5") >>> d = h.all_data()[('halos','particle_identifier')] >>> d.sort() >>> print(len(d),d[-1]) 1974910 43782.0 dimensionless >>> print(d[0:10]) [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] dimensionless It looks like the identifier gets wrapped somehow, even if the number of halos in the 2048^3 run is less than 2M. Nick Gnedin

2 1

How does YT create the gas velocity field from particle data (e.g., EAGLE, or SPH in general)?
by Stephanie Ho 10 May '18

10 May '18

Hi, I am trying to understand how YT creates the field: ['gas', 'velocity'] (or ['deposit', 'PartType0_smoothed_velocity_x'] or y/z) from EAGLE, or generally from SPH datasets. I am tracing the code to yt/fields/field_info_container.py, yt/fields/particle_fields.py, which then they would somehow link to yt/geometry/particle_deposit.pyx and yt/geometry/particle_smooth.pyx. However, I am still not sure where and how the velocity field is created. The main reason for doing so, is to understand what does the output gas velocity represent, and how is it calculated, e.g., is it weighted by the all gas density (or mass), is it weighted by HI density (or mass) etc. It would be a great help if anyone can point to the relevant codes for the calculation. Thanks, Stephanie

4 4

TypeError: only length-1 arrays can be converted to Python scalars
by Guido granda muñoz 02 May '18

02 May '18

Hello yt-users, I'm trying to read 2D flash data using a disk data object. The code is the following: import yt import numpy as np import matplotlib.pyplot as plt def angular_velocity(x,y,vx,vy): r=np.sqrt(x**2+y**2) theta=np.arctan2(y/x) w=(-vx*np.sin(theta)+vy*np.cos(theta))/r return w data_f="/home/guido/FLASH4.5/object/qsub/" data_name="rotator_mhd_2d_hdf5_chk_0015" property="density" field="dynamical_time" ds=yt.load(data_f+data_name) disk=ds.disk(ds.domain_center,[0.0,0.0,1.0],(ds.domain_ center/2.0,'cm'),(1,'cm')) time=ds.current_time.in_units("yr") print("field list \n:",ds.field_list) rp1=yt.create_profile(disk,'radius','temp') And the error I got is: yt : [INFO ] 2018-04-25 19:30:40,053 Particle file found: rotator_mhd_2d_hdf5_chk_0015 yt : [INFO ] 2018-04-25 19:30:40,079 Parameters: current_time = 0.75008895406 yt : [INFO ] 2018-04-25 19:30:40,079 Parameters: domain_dimensions = [8 8 1] yt : [INFO ] 2018-04-25 19:30:40,080 Parameters: domain_left_edge = [ 0. 0. 0.] yt : [INFO ] 2018-04-25 19:30:40,080 Parameters: domain_right_edge = [ 1. 1. 1.] yt : [INFO ] 2018-04-25 19:30:40,080 Parameters: cosmological_simulation = 0.0 Traceback (most recent call last): File "plot_ang.py", line 28, in <module> rp1=yt.create_profile(disk,'radius','temp') File "/home/guido/.conda/envs/kuro/lib/python3.5/site-packages/yt/data_objects/profiles.py", line 1034, in create_profile for f, l in zip(bin_fields, logs)] File "/home/guido/.conda/envs/kuro/lib/python3.5/site-packages/yt/data_objects/profiles.py", line 1034, in <listcomp> for f, l in zip(bin_fields, logs)] File "/home/guido/.conda/envs/kuro/lib/python3.5/site-packages/yt/data_objects/derived_quantities.py", line 511, in __call__ rv = super(Extrema, self).__call__(fields, non_zero) File "/home/guido/.conda/envs/kuro/lib/python3.5/site-packages/yt/data_objects/derived_quantities.py", line 66, in __call__ for sto, ds in parallel_objects(chunks, -1, storage = storage): File "/home/guido/.conda/envs/kuro/lib/python3.5/site-packages/yt/utilities/parallel_tools/parallel_analysis_interface.py", line 507, in parallel_objects for result_id, obj in oiter: File "/home/guido/.conda/envs/kuro/lib/python3.5/site-packages/yt/data_objects/data_containers.py", line 1223, in chunks self.get_data() # Ensure we have built ourselves File "/home/guido/.conda/envs/kuro/lib/python3.5/site-packages/yt/data_objects/data_containers.py", line 1266, in get_data self.index._identify_base_chunk(self) File "/home/guido/.conda/envs/kuro/lib/python3.5/site-packages/yt/geometry/grid_geometry_handler.py", line 293, in _identify_base_chunk gi = dobj.selector.select_grids(self.grid_left_edge, File "/home/guido/.conda/envs/kuro/lib/python3.5/site-packages/yt/data_objects/data_containers.py", line 1218, in selector self._selector = sclass(self) File "yt/geometry/selection_routines.pyx", line 1006, in yt.geometry.selection_routines.DiskSelector.__init__ TypeError: only length-1 arrays can be converted to Python scalars I am using : yt ,3.4.1 ,py35_0 And the output of the statistics of the input file used is: $ yt stats ~/FLASH4.5/wind_tunnel/qsub_folder/windtunnel_4lev_hdf5_chk_0000 yt : [INFO ] 2018-04-25 19:33:56,124 Particle file found: windtunnel_4lev_hdf5_chk_0000 yt : [INFO ] 2018-04-25 19:33:56,149 Parameters: current_time = 0.0 yt : [INFO ] 2018-04-25 19:33:56,150 Parameters: domain_dimensions = [120 40 1] yt : [INFO ] 2018-04-25 19:33:56,150 Parameters: domain_left_edge = [ 0. 0. 0.] yt : [INFO ] 2018-04-25 19:33:56,150 Parameters: domain_right_edge = [ 3. 1. 1.] yt : [INFO ] 2018-04-25 19:33:56,150 Parameters: cosmological_simulation = 0.0 level # grids # cells # cells^3 ---------------------------------------------- 0 63 4032 16 ---------------------------------------------- 63 4032 t = 0.00000000e+00 = 0.00000000e+00 s = 0.00000000e+00 years Smallest Cell: Width: 8.102e-27 Mpc Width: 8.102e-21 pc Width: 1.671e-15 AU Width: 2.500e-02 cm Could you please help me to find out what I am doing wrong, I'm new to this yt version. Kind Regards, -- Guido

3 3