Accessing particle_density in a Castro pure nbody simulation
Hey,
I'm trying to access die derived field "particle_density" via the following script:
from yt.mods import *
pf = load('plt00000') my_sphere = pf.h.sphere([0,0,0], 1000.0/pf['kpc']) print my_sphere["particle_density"]
Result: http://pastebin.com/Hukh7ahB
The data is from a pure nbody Castro run and contains no "Density"-Field in the output.
Like this:
gurke: ~/castro/runs/test_my_ic > yt stats plt00000 yt : [WARNING ] 2011-06-10 10:29:56,591 Setting 1.0 in code units to be 1.0 cm yt : [WARNING ] 2011-06-10 10:29:56,591 No time units. Setting 1.0 = 1 second. yt : [INFO ] 2011-06-10 10:29:56,592 Parameters: current_time = 0.0 yt : [INFO ] 2011-06-10 10:29:56,592 Parameters: domain_dimensions = [64, 64, 64] yt : [INFO ] 2011-06-10 10:29:56,592 Parameters: domain_left_edge = [ 0. 0. 0.] yt : [INFO ] 2011-06-10 10:29:56,592 Parameters: domain_right_edge = [ 90. 90. 90.] yt : [INFO ] 2011-06-10 10:29:56,593 Parameters: cosmological_simulation = 0.0 yt : [INFO ] 2011-06-10 10:29:56,598 Adding particle_count to list of fields yt : [INFO ] 2011-06-10 10:29:56,598 Adding particle_mass_density to list of fields Warning: invalid value encountered in sqrt Warning: divide by zero encountered in divide Warning: divide by zero encountered in divide Warning: invalid value encountered in sqrt level # grids # cells --------------------------- 0 8 262144 ---------------------------- 8 262144
t = 0.00000000e+00 = 0.00000000e+00 s = 0.00000000e+00 years
Smallest Cell: Width: 4.557e-25 mpc Width: 4.557e-22 kpc Width: 4.557e-19 pc Width: 9.402e-14 au Width: 2.022e-11 rsun Width: 8.736e-06 miles Width: 1.562e-02 unitary Width: 1.406e+00 1 Width: 1.406e+00 cm Traceback (most recent call last): File "/home/sklemer/yt/yt-x86_64/bin/yt", line 9, in <module> load_entry_point('yt==2.1dev', 'console_scripts', 'yt')() File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/command_line.py", line 1061, in run_main sys.exit(YT.main()) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/cmdln.py", line 257, in main return self.cmd(args) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/cmdln.py", line 280, in cmd retval = self.onecmd(argv) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/cmdln.py", line 412, in onecmd return self._dispatch_cmd(handler, argv) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/cmdln.py", line 1100, in _dispatch_cmd return handler(argv[0], opts, *args) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/command_line.py", line 224, in arg_iterate func(self, subcmd, opts, arg) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/command_line.py", line 535, in do_stats v, c = pf.h.find_max("Density") File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/object_finding_mixin.py", line 55, in find_max mg, mc, mv, pos = self.find_max_cell_location(field, finest_levels) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/object_finding_mixin.py", line 68, in find_max_cell_location source.quantities["MaxLocation"]( field, lazy_reader=True) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.py", line 91, in __call__ return self._call_func_lazy(args, kwargs) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.py", line 98, in _call_func_lazy rv = self.func(GridChildMaskWrapper(g, self._data_source), *args, **kwargs) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.py", line 603, in _MaxLocation if data[field].size > 0: File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.py", line 60, in __getitem__ data = self.data_source._get_data_from_grid(self.grid, item) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.py", line 73, in save_state tr = func(self, grid, field) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.py", line 2313, in _get_data_from_grid if grid[field].size == 1: # dx, dy, dz, cellvolume File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 147, in __getitem__ self.get_data(key) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 190, in get_data self._generate_field(field) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 135, in _generate_field self[field] = self.pf.field_info[field](self) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/field_info_container.py", line 309, in __call__ dd = self._function(self, data) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/frontends/castro/fields.py", line 84, in <lambda> add_field(theirs, function=lambda a, b: b[mine], take_log=True) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 147, in __getitem__ self.get_data(key) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 190, in get_data self._generate_field(field) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 125, in _generate_field self.pf.field_info[field].check_available(self) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/field_info_container.py", line 276, in check_available validator(data) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/field_info_container.py", line 374, in __call__ raise NeedsDataField(doesnt_have) yt.data_objects.field_info_container.NeedsDataField: (['density'])
Or here like this:
In [27]: pf.keys() Out[27]: ['miles', 'rsun', 'mpc', 'cm', 'unitary', '1', 'pc', 'kpc', 'au', '1', 'days', 'years', 'DualEnergyFormalism', 'cosmology_calculator', 'TopGridDimensions', 'CosmologyHubbleConstantNow', 'TopGridRank', 'CosmologyOmegaLambdaNow', 'Time', 'HydroMethod', 'CosmologyOmegaMatterNow', 'EOSType', 'particle_velocity_z', 'particle_position_z', 'particle_position_x', 'particle_position_y', 'particle_mass', 'particle_velocity_x', 'particle_velocity_y']
bye /Steffen
PS: Sorry for not using the yt-pastebin script but the lodge_it script seems to not work via the http proxy from my mashine.
Hi Steffen,
This is, funnily enough, two problems.
The one for yt stats is actually just telling you that it wanted to do something with the Density, but you have an N-body only calculation, so it can't. I've made this conditional.
The second issue is that of units. Right now the Nyx units are a bit uncertain and unknown to yt, so I have erred on the side of failing rather than returning bad data. I think that you should be able to project "particle_mass_density" without trouble. Today and early next week, in a weird coincidence, was when we were supposed to work on units for this frontend; I'll keep you in the loop on that. The ticket I was going to keep track of this with is here:
http://hg.enzotools.org/yt/issue/281/nyx-units-are-incorrect
If you wouldn't mind, it'd be very helpful if you could fill in a bit about what you know -- from a pragmatic point of view -- about the units in that ticket.
Thanks,
Matt
On Fri, Jun 10, 2011 at 1:32 AM, Steffen Klemer moh@gmx.org wrote:
Hey,
I'm trying to access die derived field "particle_density" via the following script:
from yt.mods import *
pf = load('plt00000') my_sphere = pf.h.sphere([0,0,0], 1000.0/pf['kpc']) print my_sphere["particle_density"]
Result: http://pastebin.com/Hukh7ahB
The data is from a pure nbody Castro run and contains no "Density"-Field in the output.
Like this:
gurke: ~/castro/runs/test_my_ic > yt stats plt00000 yt : [WARNING ] 2011-06-10 10:29:56,591 Setting 1.0 in code units to be 1.0 cm yt : [WARNING ] 2011-06-10 10:29:56,591 No time units. Setting 1.0 = 1 second. yt : [INFO ] 2011-06-10 10:29:56,592 Parameters: current_time = 0.0 yt : [INFO ] 2011-06-10 10:29:56,592 Parameters: domain_dimensions = [64, 64, 64] yt : [INFO ] 2011-06-10 10:29:56,592 Parameters: domain_left_edge = [ 0. 0. 0.] yt : [INFO ] 2011-06-10 10:29:56,592 Parameters: domain_right_edge = [ 90. 90. 90.] yt : [INFO ] 2011-06-10 10:29:56,593 Parameters: cosmological_simulation = 0.0 yt : [INFO ] 2011-06-10 10:29:56,598 Adding particle_count to list of fields yt : [INFO ] 2011-06-10 10:29:56,598 Adding particle_mass_density to list of fields Warning: invalid value encountered in sqrt Warning: divide by zero encountered in divide Warning: divide by zero encountered in divide Warning: invalid value encountered in sqrt level # grids # cells
0 8 262144
8 262144
t = 0.00000000e+00 = 0.00000000e+00 s = 0.00000000e+00 years
Smallest Cell: Width: 4.557e-25 mpc Width: 4.557e-22 kpc Width: 4.557e-19 pc Width: 9.402e-14 au Width: 2.022e-11 rsun Width: 8.736e-06 miles Width: 1.562e-02 unitary Width: 1.406e+00 1 Width: 1.406e+00 cm Traceback (most recent call last): File "/home/sklemer/yt/yt-x86_64/bin/yt", line 9, in <module> load_entry_point('yt==2.1dev', 'console_scripts', 'yt')() File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/command_line.py", line 1061, in run_main sys.exit(YT.main()) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/cmdln.py", line 257, in main return self.cmd(args) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/cmdln.py", line 280, in cmd retval = self.onecmd(argv) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/cmdln.py", line 412, in onecmd return self._dispatch_cmd(handler, argv) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/cmdln.py", line 1100, in _dispatch_cmd return handler(argv[0], opts, *args) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/command_line.py", line 224, in arg_iterate func(self, subcmd, opts, arg) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/utilities/command_line.py", line 535, in do_stats v, c = pf.h.find_max("Density") File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/object_finding_mixin.py", line 55, in find_max mg, mc, mv, pos = self.find_max_cell_location(field, finest_levels) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/object_finding_mixin.py", line 68, in find_max_cell_location source.quantities["MaxLocation"]( field, lazy_reader=True) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.py", line 91, in __call__ return self._call_func_lazy(args, kwargs) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.py", line 98, in _call_func_lazy rv = self.func(GridChildMaskWrapper(g, self._data_source), *args, **kwargs) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.py", line 603, in _MaxLocation if data[field].size > 0: File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.py", line 60, in __getitem__ data = self.data_source._get_data_from_grid(self.grid, item) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.py", line 73, in save_state tr = func(self, grid, field) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.py", line 2313, in _get_data_from_grid if grid[field].size == 1: # dx, dy, dz, cellvolume File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 147, in __getitem__ self.get_data(key) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 190, in get_data self._generate_field(field) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 135, in _generate_field self[field] = self.pf.field_info[field](self) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/field_info_container.py", line 309, in __call__ dd = self._function(self, data) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/frontends/castro/fields.py", line 84, in <lambda> add_field(theirs, function=lambda a, b: b[mine], take_log=True) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 147, in __getitem__ self.get_data(key) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 190, in get_data self._generate_field(field) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.py", line 125, in _generate_field self.pf.field_info[field].check_available(self) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/field_info_container.py", line 276, in check_available validator(data) File "/home/sklemer/yt/yt-x86_64/src/yt-hg/yt/data_objects/field_info_container.py", line 374, in __call__ raise NeedsDataField(doesnt_have) yt.data_objects.field_info_container.NeedsDataField: (['density'])
Or here like this:
In [27]: pf.keys() Out[27]: ['miles', 'rsun', 'mpc', 'cm', 'unitary', '1', 'pc', 'kpc', 'au', '1', 'days', 'years', 'DualEnergyFormalism', 'cosmology_calculator', 'TopGridDimensions', 'CosmologyHubbleConstantNow', 'TopGridRank', 'CosmologyOmegaLambdaNow', 'Time', 'HydroMethod', 'CosmologyOmegaMatterNow', 'EOSType', 'particle_velocity_z', 'particle_position_z', 'particle_position_x', 'particle_position_y', 'particle_mass', 'particle_velocity_x', 'particle_velocity_y']
bye /Steffen
PS: Sorry for not using the yt-pastebin script but the lodge_it script seems to not work via the http proxy from my mashine.
-- () ascii ribbon campaign - against html e-mail /\ www.asciiribbon.org - against proprietary attachments _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Hey Matt,
The second issue is that of units.
Aaaahhhh... ...found this as well after checking this a bit more. Sometimes it's a bit hard to not get confused by the yt codebase :-) I couldn't imagine that yt tries to convert the result into the right units... ...but it's a nice feature!! :-)
http://hg.enzotools.org/yt/issue/281/nyx-units-are-incorrect
If you wouldn't mind, it'd be very helpful if you could fill in a bit about what you know -- from a pragmatic point of view -- about the units in that ticket.
For the Particle-Plotfiles it is:
x,y,z in MPc velocities in km/s/a (where a is the comoving a) masses in M_{\sun}
bye /Steffen
Hi Steffen,
On Tue, Jun 14, 2011 at 7:15 AM, Steffen Klemer moh@gmx.org wrote:
Hey Matt,
The second issue is that of units.
Aaaahhhh... ...found this as well after checking this a bit more. Sometimes it's a bit hard to not get confused by the yt codebase :-)
We have tried to separate it into fundamental submodules, to ease this transition. Like many scientific codes it has grown organically, but we have actively tried to organize it in a logical way. That being said, we would be eager to hear your suggestions for improving the organization, so please feel free to either issue pull requests from your fork on BitBucket or fill out a ticket.
I couldn't imagine that yt tries to convert the result into the right units... ...but it's a nice feature!! :-)
Thanks. One of the principal ideas behind yt has been to address physically-relevant quantities, rather than simulation quantities. (Not that the latter isn't available, however.)
http://hg.enzotools.org/yt/issue/281/nyx-units-are-incorrect
If you wouldn't mind, it'd be very helpful if you could fill in a bit about what you know -- from a pragmatic point of view -- about the units in that ticket.
For the Particle-Plotfiles it is:
x,y,z in MPc velocities in km/s/a (where a is the comoving a) masses in M_{\sun}
Excellent! I have spoken to Ann about this, and she's preparing a more comprehensive solution. For now I will put this in. For clarification, to translate from the units in the plotfile to cgs, it would be multiplication by 10^4 * (current scale factor), right?
Thanks very much,
Matt
bye /Steffen
-- () ascii ribbon campaign - against html e-mail /\ www.asciiribbon.org - against proprietary attachments _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Am Tue, 14.06.2011 um 13:57 schrieb Matthew Turk matthewturk@gmail.com:
velocities in km/s/a (where a is the comoving a)
Excellent! I have spoken to Ann about this, and she's preparing a more comprehensive solution. For now I will put this in. For clarification, to translate from the units in the plotfile to cgs, it would be multiplication by 10^4 * (current scale factor), right?
I think it's 10^5 from km to cm but otherwise yes:
10^5 * (comoving_a)
thanks for your help and looking forward to a generic nyx unit scheme :-) /Steffen
Hi Steffen,
On Wed, Jun 15, 2011 at 12:40 AM, Steffen Klemer moh@gmx.org wrote:
Am Tue, 14.06.2011 um 13:57 schrieb Matthew Turk matthewturk@gmail.com:
velocities in km/s/a (where a is the comoving a)
Excellent! I have spoken to Ann about this, and she's preparing a more comprehensive solution. For now I will put this in. For clarification, to translate from the units in the plotfile to cgs, it would be multiplication by 10^4 * (current scale factor), right?
I think it's 10^5 from km to cm but otherwise yes:
10^5 * (comoving_a)
thanks for your help and looking forward to a generic nyx unit scheme :-) /Steffen
Hehe, yes, it is 10^5 ... thanks for being kind about my typo. :)
I have made the change to convert these units, but it looks to me like -- my own fault here -- the unit system is composed of a number of copy/pasted items that I wrote quite a while ago. I'll take a closer look at these and fix them as we go along.
Best,
Matt
-- () ascii ribbon campaign - against html e-mail /\ www.asciiribbon.org - against proprietary attachments _______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
participants (2)
-
Matthew Turk -
Steffen Klemer