Issue with changing from code units to real units on Gadget-2 files
I've been attempting to look at a Gadget-2 file's temperature by converting the code units into real units, by using the unit_base argument but I haven't been able to alter the temperature in any way. I am wondering if anyone knows how to convert the gadget units of temperature? The temperature for my simulation is not individually given to each particle but is derived from the internal energy and chemical abundance. Any help would be greatly appreciated. These are my imports import yt import soxs import pyxsim import numpy as np This is the code fname = '/Users/Andrew/Downloads/output_00056/new_180' unit_base = {'UnitLength_in_cm' : 1.2e+17, 'UnitTime_in_yr' : 61.03633, 'UnitVelocity_in_cm_per_s' : 6.23e+7, 'UnitMass_in_g' : 6.97e+39, 'UnitTemp_in_K' : 3.568e+7} bbox_lim = 850 #kpc bbox = [[-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim]] ds = yt.load(fname,unit_base=unit_base,bounding_box=bbox,default_species_fields="ionized") ds.index ad= ds.all_data() This is the ds field list: ('Bndry', 'Coordinates'), ('Bndry', 'Mass'), ('Bndry', 'ParticleIDs'), ('Bndry', 'Velocities'), ('Gas', 'Coordinates'), ('Gas', 'Density'), ('Gas', 'InternalEnergy'), ('Gas', 'Mass'), ('Gas', 'ParticleIDs'), ('Gas', 'SmoothingLength'), ('Gas', 'Velocities'), ('Halo', 'Coordinates'), ('Halo', 'Mass'), ('Halo', 'ParticleIDs'), ('Halo', 'Velocities'), ('Stars', 'Coordinates'), ('Stars', 'Mass'), ('Stars', 'ParticleIDs'), ('Stars', 'Velocities'), ('all', 'Coordinates'), ('all', 'Mass'), ('all', 'ParticleIDs'), ('all', 'Velocities'), ('nbody', 'Coordinates'), ('nbody', 'Mass'), ('nbody', 'ParticleIDs'), ('nbody', 'Velocities')]
Hi! So I think in principle, what you'll want to do is compute a derived field: https://yt-project.org/doc/developing/creating_derived_fields.html I believe that it is safe to say that the input into these will be fields that are associated with units -- unless the fields are not known by yt (but the fields you list all should be) they will be assigned cgs units. So when you get them in, you can apply the appropriate formulation. If you want a bit more help, and can share the exact formulation you want to use to translate from the abundances and energies into temperature, I can go deeper into this. -Matt On Sun, Apr 10, 2022 at 3:26 PM <moonpigeona@gmail.com> wrote:
Hello! Thank you so much for your help! I am still having a slight issue though. I keep getting this error: YTNonIndexedDataContainer: The data container (<class 'yt.data_objects.index_subobjects.particle_container.ParticleContainer'>) is an unindexed type. Operations such as ires, icoords, fcoords and fwidth will not work on it. This is my code. (The way in which the temperature should be found is by the molecular weight and the internal energy) -------------------------------------------------------------------------------------------------------------------------------- Code Version 1: import yt from yt.units import dimensions UnitTemp_in_K = 3.568e+7 def my_temperature(field, data): return (data["Gas","Temperature"] * UnitTemp_in_K) yt.add_field( ("my_temperature"), function=_temperatures, sampling_type="local", units='code_mass*code_specific_energy'", force_override=True) fname = '/Users/Andrew/Downloads/output_00056/new_180' unit_base = {'UnitLength_in_cm' : 1.2e+17, 'UnitTime_in_yr' : 61.03633, 'UnitVelocity_in_cm_per_s' : 6.23e+7, 'UnitMass_in_g' : 6.97e+39} bbox_lim = 850 #kpc bbox = [[-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim]] ds = yt.load(fname,unit_base=unit_base,bounding_box=bbox,default_species_fields="ionized") ds.index ad= ds.all_data() slc = yt.SlicePlot(ds, "z", [("my_temperature")], width=(1.5,"pc")) slc.show() Then the error appears. -------------------------------------------------------------------------------------------------------------------------------- Code Version 2: import yt from yt.units import dimensions UnitTemp_in_K = 3.568e+7 def my_temperature(field, data): return (data["Gas","InternalEnergy"] * data["Gas","Mass"] * UnitTemp_in_K) yt.add_field( ("my_temperature"), function=_temperatures, sampling_type="local", units='code_mass*code_specific_energy'", force_override=True) fname = '/Users/Andrew/Downloads/output_00056/new_180' unit_base = {'UnitLength_in_cm' : 1.2e+17, 'UnitTime_in_yr' : 61.03633, 'UnitVelocity_in_cm_per_s' : 6.23e+7, 'UnitMass_in_g' : 6.97e+39} bbox_lim = 850 #kpc bbox = [[-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim]] ds = yt.load(fname,unit_base=unit_base,bounding_box=bbox,default_species_fields="ionized") ds.index ad= ds.all_data() slc = yt.SlicePlot(ds, "z", [("my_temperature")], width=(1.5,"pc")) slc.show() Then the error appears. -------------------------------------------------------------------------------------------------------------------------------- In addition, just in case this could help. The overall goal is to use the following code source_model = pyxsim.ThermalSourceModel("apec", 0.05, 11.0, 1000, Zmet=1.0, kT_min=0.025,temperature_field=None,emission_measure_field=("gas","emission_measure")) exp_time = (2000.0, "ks") # Exposure time area = (15.0, "kpc**2") # Collecting area redshift = 0.0 # because the source is local dist = (8.18, "kpc") # distance to the galactic center n_photons, n_cells = pyxsim.make_photons("gadg_photonstest",ad, redshift, area, exp_time, source_model,dist = dist) This should generate photons, however, none are being made. I know this part is pyxsim but if you have any information you think may help, that would be greatly appreciated. Thank you again!
Hi Andrew, I see a couple of issues here: 1. First, the derivation of the temperature field. In your first example, you have: UnitTemp_in_K = 3.568e+7 def my_temperature(field, data): return (data["Gas","Temperature"] * UnitTemp_in_K) yt.add_field( ("my_temperature"), function=_temperatures, sampling_type="local", units='code_mass*code_specific_energy'", force_override=True) But yt may recognize data["Gas","Temperature"] as a temperature field and return it in "K", or it may not (I’m not quite sure), and assume it is dimensionless. In either case, these are not the units of "code_mass*code_specific_energy" that you specify. The second field definition doesn’t have this problem, but I should note that pyxsim expects a temperature field to have units of K (or at least the same dimension). 2. The second potential issue, which likely causes your error, is: unit_base = {'UnitLength_in_cm' : 1.2e+17, 'UnitTime_in_yr' : 61.03633, 'UnitVelocity_in_cm_per_s' : 6.23e+7, 'UnitMass_in_g' : 6.97e+39} bbox_lim = 850 #kpc bbox = [[-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim]] Here in the unit_base, the code units of length are 1.2e17 cm, but is the bbox_lim really 850 kpc? Best, John
Here is the full traceback error: yt : [INFO yt : [INFO yt : [INFO yt : [INFO yt : [INFO ] 2022-04-16 17:25:18,009 xlim = -19.285485 19.285485 ] 2022-04-16 17:25:18,010 ylim = -19.285485 19.285485 ] 2022-04-16 17:25:18,011 xlim = -19.285485 19.285485 ] 2022-04-16 17:25:18,012 ylim = -19.285485 19.285485 ] 2022-04-16 17:25:18,014 Making a fixed resolution buffer of (('gas', 'my_temperature')) 800 by 800 3 slc.show() File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:2528, in SlicePlot(ds, normal, fields, axis, *args, **kwargs) 2522 2523 2524 2525 2526 -> 2528 mylog.warning( "Ignoring 'north_vector' keyword as it is ill-defined for " "an AxisAlignedSlicePlot object." ) del kwargs["north vector"] return (ds, normal, fields, *args, **kwargs) _ AxisAlignedSlicePlot File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:1658, in AxisAlignedSlicePlot.__init__(self, ds, axis, fields, center, width, axes_unit , origin, right_handed, fontsize, field_parameters, window_size, aspect, data_source, buff_size) 1656 1657 -> 1658 slc.get_data(fields) validate mesh_fields(slc, fields) slc, bounds, origin=origin, fontsize=fontsize, fields=fields, window_size=window_size, aspect=aspect, right_handed=right_handed, buff_size=buff_size, 253 254 --> 255 256 257 if f not in self.field_data and key not in self.field_data: if f in self._container_fields: self.field_data[f] = self.ds.arr(self._generate_container_field(f)) return self.field_data[f] else: _ PWViewerMPL ] 2022-04-16 15:46:16,784 Omega Lambda is 0.0, so we are turning off Cosmology. ] 2022-04-16 15:46:16,838 Parameters: current_time ] 2022-04-16 15:46:16,839 Parameters: domain_dimensions ] 2022-04-16 15:46:16,839 Parameters: domain_left_edge ] 2022-04-16 15:46:16,841 Parameters: domain_right_edge ] 2022-04-16 15:46:16,842 Parameters: cosmological_simulation ] 2022-04-16 15:46:16,851 Your Gadget-2 file may have extra columns or different precision! (16925336 diff => ['Bndry', 'Halo']?) ] 2022-04-16 15:46:16,857 Allocating for 2.116e+06 particles .__init__( self, 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 ) 1670 if axes_unit is None: 1671 axes_unit = get_axes_unit(width, ds) = 17.999999932944775 = [1 1 1] = [-850. -850. -850.] = [850. 850. 850.] =0 fname = '/Users/Andrew/Downloads/output_00056/new_180' unit_base = {'UnitLength_in_cm' 'UnitTime_in_yr' 'UnitVelocity_in_cm_per_s' 'UnitMass_in_g' 'UnitTemp_in_K' #3.568e+7 bbox_lim = 850 #kpc bbox = [[-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim]] : 1.2e+17, : 61.03633, : 6.23e+7, : 6.97e+39, : 3.568e+7} ds = yt.load(fname,unit_base=unit_base,bounding_box=bbox,default_species_fields="ionized") ds.index ad= ds.all_data() yt : [INFO yt : [INFO yt : [INFO yt : [INFO yt : [INFO yt : [INFO yt : [WARNING yt : [INFO Loading particle index: 100%|██████████| 9/9 [00:00<00:00, 23.89it/s] yt : [WARNING ] 2022-04-16 15:46:20,442 Your Gadget-2 file may have extra columns or different precision! (16925336 diff => ['Bndry', 'Halo']?) yt : [INFO ] 2022-04-16 15:46:20,445 Allocating for 2.116e+06 particles Loading particle index: 100%|██████████| 9/9 [00:00<00:00, 24.75it/s] ] 2022-04-16 15:46:19,952 Omega Lambda is 0.0, so we are turning off Cosmology. ] 2022-04-16 15:46:20,002 Parameters: current_time ] 2022-04-16 15:46:20,003 Parameters: domain_dimensions ] 2022-04-16 15:46:20,004 Parameters: domain_left_edge ] 2022-04-16 15:46:20,005 Parameters: domain_right_edge ] 2022-04-16 15:46:20,006 Parameters: cosmological_simulation ] 2022-04-16 15:46:20,016 Your Gadget-2 file may have extra columns or different precision! (16925336 diff => ['Bndry', 'Halo']?) ] 2022-04-16 15:46:20,020 Allocating for 2.116e+06 particles = 17.999999932944775 = [1 1 1] = [-850. -850. -850.] = [850. 850. 850.] =0 ds.all_data() print ('left edge: ',ds.domain_left_edge) print ('right edge: ',ds.domain_right_edge) print ('center: ',ds.domain_center) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:873, in PWViewerMPL.__init__(self, *args, **kwargs) 871 872 self._plot_type = kwargs.pop("plot_type") self._splat_color = kwargs.pop("splat_color", None) PlotWindow.__init__(self, *args, **kwargs) --> 873 File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:257, in PlotWindow.__init__(self, data_source, bounds, buff_size, antialias, periodic, origin, 254 256 oblique, right_handed, window_size, fields, fontsize, aspect, setup) return if not self. data valid: self. () self._data_valid = True self._colorbar_valid = True self._units_config[field] self.setup callbacks() self. () _ _setup_plots --> 257 File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:965, in PWViewerMPL._setup_plots(self) 963 964 --> 965 966 967 __ _recreate_frb File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:317, in PlotWindow._recreate_frb(self) 315 316 --> 317 319 320 # At this point the frb has the valid bounds, size, aliasing, etc. if old_fields is None: self._frb._get_data_source_fields() # New frb, apply default units (if any) for field, field_unit in self._units_config.items(): File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/fixed_resolution.py:176, in FixedResolutionBuffer._get_data_source_fields(self) 174 for f in fields: 175 if f not in exclude and f[0] not in self.data_source.ds.particle_types: --> 176 self[f] File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/fixed_resolution.py:139, in FixedResolutionBuffer.__getitem__(self, item) 136 b = float(b.in_units("code_length")) 137 bounds.append(b) --> 139 .pixelize( .axis, self. item, bounds, self.buff_size, int(self.antialias), name, (args, kwargs) in self._filters: buff = filter_registry[name](*args[1:], **kwargs).apply(buff) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/coordinates/cartesian_coordinates.py:217, in CartesianCoordinateHandler.pixelize(self, dimension, data_sour ce, field, bounds, size, antialias, periodic) 214 216 --> 217 return np.squeeze(np.transpose(img, (yax, xax, ax))) elif self.axis_id.get(dimension, dimension) < 3: return self. ( , field, bounds, size, antialias, dimension, periodic ) 140 141 142 143 144 145 146 ) 148 for 149 , buff = self.ds. self. coordinates data_source data_source _ data_source ortho_pixelize 218 219 220 else: 221 return self._oblique_pixelize(data_source, field, bounds, size, antialias) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/coordinates/cartesian_coordinates.py:537, in CartesianCoordinateHandler._ortho_pixelize(self, data_source, field, bounds, size, antialias, dim, periodic) 533 buff = buff.transpose() 534 else: 535 pixelize_cartesian( 536 buff, 539 data_source["pdx"], 540 data_source["pdy"], 541 data_source[field], 542 bounds, 543 int(antialias), 544 period, 545 int(periodic), 546 ) 547 return buff File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/data_containers.py:255, in YTDataContainer.__getitem__(self, key) --> 537 data_source["px"], 538 data_source["py"], File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/selection_objects/slices.py:89, in YTSlice._generate_container_field(self, field) 87 88 ---> 89 90 91 253 254 257 tr = func(self) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/geometry_handler.py:292, in YTDataChunk._accumulate_values(self, method) 290 291 --> 292 293 294 for obj in self._fast_index or self.objs: f = getattr(obj, mname) arrs.append(f(self.dobj)) if method == "dtcoords": arrs = [arr[0] for arr in arrs] self.index._identify_base_chunk(self) if field == "px": return self._current_chunk.fcoords[:, xax] elif field == "py": return self._current_chunk.fcoords[:, yax] File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/geometry_handler.py:255, in cached_property.<locals>.cached_func(self) return getattr(self, n) if self.data_size is None: tr = self._accumulate_values(n[1:]) 256 else: --> 255 File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/index_subobjects/particle_container.py:17, in _non_indexed.<locals>._func_non_indexed(self, *args, **kw args) 16 def _func_non_indexed(self, *args, **kwargs): ---> 17 raise YTNonIndexedDataContainer(self) YTNonIndexedDataContainer: The data container (<class 'yt.data_objects.index_subobjects.particle_container.ParticleContainer'>) is an unindexed type. Operations such as ires, icoords, fcoords and fwidth will not work on it. In terms of the bbox_lim, you are right I forgot to change it, it should be 1.5 pc. In regards to the temperature the code to real value should be the UnitTemp_in_K * InternalEnergy, so I changed that. However the units are code_specific_energy not kelvin. These changes are below. import yt from yt.units import dimensions UnitTemp_in_K = 3.568e+7 def my_temperature(field, data): return (data["Gas","InternalEnergy"] * UnitTemp_in_K) yt.add_field( ("my_temperature"), function=my_temperature, sampling_type="local", units='code_specific_energy', force_override=True
I realize the traceback error was sent in a really weird format, sorry about that this should be easier to read. YTNonIndexedDataContainer Traceback (most recent call last) Input In [25], in <module> ----> 1 slc = yt.SlicePlot(ds, "z", [("my_temperature")], width=(1.5,"pc")) 3 slc.show() File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:2528, in SlicePlot(ds, normal, fields, axis, *args, **kwargs) 2522 mylog.warning( 2523 "Ignoring 'north_vector' keyword as it is ill-defined for " 2524 "an AxisAlignedSlicePlot object." 2525 ) 2526 del kwargs["north_vector"] -> 2528 return AxisAlignedSlicePlot(ds, normal, fields, *args, **kwargs) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:1658, in AxisAlignedSlicePlot.__init__(self, ds, axis, fields, center, width, axes_unit, origin, right_handed, fontsize, field_parameters, window_size, aspect, data_source, buff_size) 1656 slc.get_data(fields) 1657 validate_mesh_fields(slc, fields) -> 1658 PWViewerMPL.__init__( 1659 self, 1660 slc, 1661 bounds, 1662 origin=origin, 1663 fontsize=fontsize, 1664 fields=fields, 1665 window_size=window_size, 1666 aspect=aspect, 1667 right_handed=right_handed, 1668 buff_size=buff_size, 1669 ) 1670 if axes_unit is None: 1671 axes_unit = get_axes_unit(width, ds) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:873, in PWViewerMPL.__init__(self, *args, **kwargs) 871 self._plot_type = kwargs.pop("plot_type") 872 self._splat_color = kwargs.pop("splat_color", None) --> 873 PlotWindow.__init__(self, *args, **kwargs) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:257, in PlotWindow.__init__(self, data_source, bounds, buff_size, antialias, periodic, origin, oblique, right_handed, window_size, fields, fontsize, aspect, setup) 254 self._units_config[field] 256 self.setup_callbacks() --> 257 self._setup_plots() File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:965, in PWViewerMPL._setup_plots(self) 963 return 964 if not self._data_valid: --> 965 self._recreate_frb() 966 self._data_valid = True 967 self._colorbar_valid = True File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:317, in PlotWindow._recreate_frb(self) 315 # At this point the frb has the valid bounds, size, aliasing, etc. 316 if old_fields is None: --> 317 self._frb._get_data_source_fields() 319 # New frb, apply default units (if any) 320 for field, field_unit in self._units_config.items(): File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/fixed_resolution.py:176, in FixedResolutionBuffer._get_data_source_fields(self) 174 for f in fields: 175 if f not in exclude and f[0] not in self.data_source.ds.particle_types: --> 176 self[f] File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/fixed_resolution.py:139, in FixedResolutionBuffer.__getitem__(self, item) 136 b = float(b.in_units("code_length")) 137 bounds.append(b) --> 139 buff = self.ds.coordinates.pixelize( 140 self.data_source.axis, 141 self.data_source, 142 item, 143 bounds, 144 self.buff_size, 145 int(self.antialias), 146 ) 148 for name, (args, kwargs) in self._filters: 149 buff = filter_registry[name](*args[1:], **kwargs).apply(buff) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/coordinates/cartesian_coordinates.py:217, in CartesianCoordinateHandler.pixelize(self, dimension, data_source, field, bounds, size, antialias, periodic) 214 return np.squeeze(np.transpose(img, (yax, xax, ax))) 216 elif self.axis_id.get(dimension, dimension) < 3: --> 217 return self._ortho_pixelize( 218 data_source, field, bounds, size, antialias, dimension, periodic 219 ) 220 else: 221 return self._oblique_pixelize(data_source, field, bounds, size, antialias) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/coordinates/cartesian_coordinates.py:537, in CartesianCoordinateHandler._ortho_pixelize(self, data_source, field, bounds, size, antialias, dim, periodic) 533 buff = buff.transpose() 534 else: 535 pixelize_cartesian( 536 buff, --> 537 data_source["px"], 538 data_source["py"], 539 data_source["pdx"], 540 data_source["pdy"], 541 data_source[field], 542 bounds, 543 int(antialias), 544 period, 545 int(periodic), 546 ) 547 return buff File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/data_containers.py:255, in YTDataContainer.__getitem__(self, key) 253 if f not in self.field_data and key not in self.field_data: 254 if f in self._container_fields: --> 255 self.field_data[f] = self.ds.arr(self._generate_container_field(f)) 256 return self.field_data[f] 257 else: File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/selection_objects/slices.py:89, in YTSlice._generate_container_field(self, field) 87 self.index._identify_base_chunk(self) 88 if field == "px": ---> 89 return self._current_chunk.fcoords[:, xax] 90 elif field == "py": 91 return self._current_chunk.fcoords[:, yax] File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/geometry_handler.py:255, in cached_property.<locals>.cached_func(self) 253 return getattr(self, n) 254 if self.data_size is None: --> 255 tr = self._accumulate_values(n[1:]) 256 else: 257 tr = func(self) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/geometry_handler.py:292, in YTDataChunk._accumulate_values(self, method) 290 for obj in self._fast_index or self.objs: 291 f = getattr(obj, mname) --> 292 arrs.append(f(self.dobj)) 293 if method == "dtcoords": 294 arrs = [arr[0] for arr in arrs] File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/index_subobjects/particle_container.py:17, in _non_indexed.<locals>._func_non_indexed(self, *args, **kwargs) 16 def _func_non_indexed(self, *args, **kwargs): ---> 17 raise YTNonIndexedDataContainer(self) YTNonIndexedDataContainer: The data container (<class 'yt.data_objects.index_subobjects.particle_container.ParticleContainer'>) is an unindexed type. Operations such as ires, icoords, fcoords and fwidth will not work on it.
Hi! So I think in principle, what you'll want to do is compute a derived field: https://yt-project.org/doc/developing/creating_derived_fields.html I believe that it is safe to say that the input into these will be fields that are associated with units -- unless the fields are not known by yt (but the fields you list all should be) they will be assigned cgs units. So when you get them in, you can apply the appropriate formulation. If you want a bit more help, and can share the exact formulation you want to use to translate from the abundances and energies into temperature, I can go deeper into this. -Matt On Sun, Apr 10, 2022 at 3:26 PM <moonpigeona@gmail.com> wrote:
Hello! Thank you so much for your help! I am still having a slight issue though. I keep getting this error: YTNonIndexedDataContainer: The data container (<class 'yt.data_objects.index_subobjects.particle_container.ParticleContainer'>) is an unindexed type. Operations such as ires, icoords, fcoords and fwidth will not work on it. This is my code. (The way in which the temperature should be found is by the molecular weight and the internal energy) -------------------------------------------------------------------------------------------------------------------------------- Code Version 1: import yt from yt.units import dimensions UnitTemp_in_K = 3.568e+7 def my_temperature(field, data): return (data["Gas","Temperature"] * UnitTemp_in_K) yt.add_field( ("my_temperature"), function=_temperatures, sampling_type="local", units='code_mass*code_specific_energy'", force_override=True) fname = '/Users/Andrew/Downloads/output_00056/new_180' unit_base = {'UnitLength_in_cm' : 1.2e+17, 'UnitTime_in_yr' : 61.03633, 'UnitVelocity_in_cm_per_s' : 6.23e+7, 'UnitMass_in_g' : 6.97e+39} bbox_lim = 850 #kpc bbox = [[-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim]] ds = yt.load(fname,unit_base=unit_base,bounding_box=bbox,default_species_fields="ionized") ds.index ad= ds.all_data() slc = yt.SlicePlot(ds, "z", [("my_temperature")], width=(1.5,"pc")) slc.show() Then the error appears. -------------------------------------------------------------------------------------------------------------------------------- Code Version 2: import yt from yt.units import dimensions UnitTemp_in_K = 3.568e+7 def my_temperature(field, data): return (data["Gas","InternalEnergy"] * data["Gas","Mass"] * UnitTemp_in_K) yt.add_field( ("my_temperature"), function=_temperatures, sampling_type="local", units='code_mass*code_specific_energy'", force_override=True) fname = '/Users/Andrew/Downloads/output_00056/new_180' unit_base = {'UnitLength_in_cm' : 1.2e+17, 'UnitTime_in_yr' : 61.03633, 'UnitVelocity_in_cm_per_s' : 6.23e+7, 'UnitMass_in_g' : 6.97e+39} bbox_lim = 850 #kpc bbox = [[-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim]] ds = yt.load(fname,unit_base=unit_base,bounding_box=bbox,default_species_fields="ionized") ds.index ad= ds.all_data() slc = yt.SlicePlot(ds, "z", [("my_temperature")], width=(1.5,"pc")) slc.show() Then the error appears. -------------------------------------------------------------------------------------------------------------------------------- In addition, just in case this could help. The overall goal is to use the following code source_model = pyxsim.ThermalSourceModel("apec", 0.05, 11.0, 1000, Zmet=1.0, kT_min=0.025,temperature_field=None,emission_measure_field=("gas","emission_measure")) exp_time = (2000.0, "ks") # Exposure time area = (15.0, "kpc**2") # Collecting area redshift = 0.0 # because the source is local dist = (8.18, "kpc") # distance to the galactic center n_photons, n_cells = pyxsim.make_photons("gadg_photonstest",ad, redshift, area, exp_time, source_model,dist = dist) This should generate photons, however, none are being made. I know this part is pyxsim but if you have any information you think may help, that would be greatly appreciated. Thank you again!
Hi Andrew, I see a couple of issues here: 1. First, the derivation of the temperature field. In your first example, you have: UnitTemp_in_K = 3.568e+7 def my_temperature(field, data): return (data["Gas","Temperature"] * UnitTemp_in_K) yt.add_field( ("my_temperature"), function=_temperatures, sampling_type="local", units='code_mass*code_specific_energy'", force_override=True) But yt may recognize data["Gas","Temperature"] as a temperature field and return it in "K", or it may not (I’m not quite sure), and assume it is dimensionless. In either case, these are not the units of "code_mass*code_specific_energy" that you specify. The second field definition doesn’t have this problem, but I should note that pyxsim expects a temperature field to have units of K (or at least the same dimension). 2. The second potential issue, which likely causes your error, is: unit_base = {'UnitLength_in_cm' : 1.2e+17, 'UnitTime_in_yr' : 61.03633, 'UnitVelocity_in_cm_per_s' : 6.23e+7, 'UnitMass_in_g' : 6.97e+39} bbox_lim = 850 #kpc bbox = [[-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim]] Here in the unit_base, the code units of length are 1.2e17 cm, but is the bbox_lim really 850 kpc? Best, John
Here is the full traceback error: yt : [INFO yt : [INFO yt : [INFO yt : [INFO yt : [INFO ] 2022-04-16 17:25:18,009 xlim = -19.285485 19.285485 ] 2022-04-16 17:25:18,010 ylim = -19.285485 19.285485 ] 2022-04-16 17:25:18,011 xlim = -19.285485 19.285485 ] 2022-04-16 17:25:18,012 ylim = -19.285485 19.285485 ] 2022-04-16 17:25:18,014 Making a fixed resolution buffer of (('gas', 'my_temperature')) 800 by 800 3 slc.show() File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:2528, in SlicePlot(ds, normal, fields, axis, *args, **kwargs) 2522 2523 2524 2525 2526 -> 2528 mylog.warning( "Ignoring 'north_vector' keyword as it is ill-defined for " "an AxisAlignedSlicePlot object." ) del kwargs["north vector"] return (ds, normal, fields, *args, **kwargs) _ AxisAlignedSlicePlot File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:1658, in AxisAlignedSlicePlot.__init__(self, ds, axis, fields, center, width, axes_unit , origin, right_handed, fontsize, field_parameters, window_size, aspect, data_source, buff_size) 1656 1657 -> 1658 slc.get_data(fields) validate mesh_fields(slc, fields) slc, bounds, origin=origin, fontsize=fontsize, fields=fields, window_size=window_size, aspect=aspect, right_handed=right_handed, buff_size=buff_size, 253 254 --> 255 256 257 if f not in self.field_data and key not in self.field_data: if f in self._container_fields: self.field_data[f] = self.ds.arr(self._generate_container_field(f)) return self.field_data[f] else: _ PWViewerMPL ] 2022-04-16 15:46:16,784 Omega Lambda is 0.0, so we are turning off Cosmology. ] 2022-04-16 15:46:16,838 Parameters: current_time ] 2022-04-16 15:46:16,839 Parameters: domain_dimensions ] 2022-04-16 15:46:16,839 Parameters: domain_left_edge ] 2022-04-16 15:46:16,841 Parameters: domain_right_edge ] 2022-04-16 15:46:16,842 Parameters: cosmological_simulation ] 2022-04-16 15:46:16,851 Your Gadget-2 file may have extra columns or different precision! (16925336 diff => ['Bndry', 'Halo']?) ] 2022-04-16 15:46:16,857 Allocating for 2.116e+06 particles .__init__( self, 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 ) 1670 if axes_unit is None: 1671 axes_unit = get_axes_unit(width, ds) = 17.999999932944775 = [1 1 1] = [-850. -850. -850.] = [850. 850. 850.] =0 fname = '/Users/Andrew/Downloads/output_00056/new_180' unit_base = {'UnitLength_in_cm' 'UnitTime_in_yr' 'UnitVelocity_in_cm_per_s' 'UnitMass_in_g' 'UnitTemp_in_K' #3.568e+7 bbox_lim = 850 #kpc bbox = [[-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim], [-bbox_lim,bbox_lim]] : 1.2e+17, : 61.03633, : 6.23e+7, : 6.97e+39, : 3.568e+7} ds = yt.load(fname,unit_base=unit_base,bounding_box=bbox,default_species_fields="ionized") ds.index ad= ds.all_data() yt : [INFO yt : [INFO yt : [INFO yt : [INFO yt : [INFO yt : [INFO yt : [WARNING yt : [INFO Loading particle index: 100%|██████████| 9/9 [00:00<00:00, 23.89it/s] yt : [WARNING ] 2022-04-16 15:46:20,442 Your Gadget-2 file may have extra columns or different precision! (16925336 diff => ['Bndry', 'Halo']?) yt : [INFO ] 2022-04-16 15:46:20,445 Allocating for 2.116e+06 particles Loading particle index: 100%|██████████| 9/9 [00:00<00:00, 24.75it/s] ] 2022-04-16 15:46:19,952 Omega Lambda is 0.0, so we are turning off Cosmology. ] 2022-04-16 15:46:20,002 Parameters: current_time ] 2022-04-16 15:46:20,003 Parameters: domain_dimensions ] 2022-04-16 15:46:20,004 Parameters: domain_left_edge ] 2022-04-16 15:46:20,005 Parameters: domain_right_edge ] 2022-04-16 15:46:20,006 Parameters: cosmological_simulation ] 2022-04-16 15:46:20,016 Your Gadget-2 file may have extra columns or different precision! (16925336 diff => ['Bndry', 'Halo']?) ] 2022-04-16 15:46:20,020 Allocating for 2.116e+06 particles = 17.999999932944775 = [1 1 1] = [-850. -850. -850.] = [850. 850. 850.] =0 ds.all_data() print ('left edge: ',ds.domain_left_edge) print ('right edge: ',ds.domain_right_edge) print ('center: ',ds.domain_center) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:873, in PWViewerMPL.__init__(self, *args, **kwargs) 871 872 self._plot_type = kwargs.pop("plot_type") self._splat_color = kwargs.pop("splat_color", None) PlotWindow.__init__(self, *args, **kwargs) --> 873 File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:257, in PlotWindow.__init__(self, data_source, bounds, buff_size, antialias, periodic, origin, 254 256 oblique, right_handed, window_size, fields, fontsize, aspect, setup) return if not self. data valid: self. () self._data_valid = True self._colorbar_valid = True self._units_config[field] self.setup callbacks() self. () _ _setup_plots --> 257 File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:965, in PWViewerMPL._setup_plots(self) 963 964 --> 965 966 967 __ _recreate_frb File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:317, in PlotWindow._recreate_frb(self) 315 316 --> 317 319 320 # At this point the frb has the valid bounds, size, aliasing, etc. if old_fields is None: self._frb._get_data_source_fields() # New frb, apply default units (if any) for field, field_unit in self._units_config.items(): File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/fixed_resolution.py:176, in FixedResolutionBuffer._get_data_source_fields(self) 174 for f in fields: 175 if f not in exclude and f[0] not in self.data_source.ds.particle_types: --> 176 self[f] File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/fixed_resolution.py:139, in FixedResolutionBuffer.__getitem__(self, item) 136 b = float(b.in_units("code_length")) 137 bounds.append(b) --> 139 .pixelize( .axis, self. item, bounds, self.buff_size, int(self.antialias), name, (args, kwargs) in self._filters: buff = filter_registry[name](*args[1:], **kwargs).apply(buff) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/coordinates/cartesian_coordinates.py:217, in CartesianCoordinateHandler.pixelize(self, dimension, data_sour ce, field, bounds, size, antialias, periodic) 214 216 --> 217 return np.squeeze(np.transpose(img, (yax, xax, ax))) elif self.axis_id.get(dimension, dimension) < 3: return self. ( , field, bounds, size, antialias, dimension, periodic ) 140 141 142 143 144 145 146 ) 148 for 149 , buff = self.ds. self. coordinates data_source data_source _ data_source ortho_pixelize 218 219 220 else: 221 return self._oblique_pixelize(data_source, field, bounds, size, antialias) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/coordinates/cartesian_coordinates.py:537, in CartesianCoordinateHandler._ortho_pixelize(self, data_source, field, bounds, size, antialias, dim, periodic) 533 buff = buff.transpose() 534 else: 535 pixelize_cartesian( 536 buff, 539 data_source["pdx"], 540 data_source["pdy"], 541 data_source[field], 542 bounds, 543 int(antialias), 544 period, 545 int(periodic), 546 ) 547 return buff File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/data_containers.py:255, in YTDataContainer.__getitem__(self, key) --> 537 data_source["px"], 538 data_source["py"], File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/selection_objects/slices.py:89, in YTSlice._generate_container_field(self, field) 87 88 ---> 89 90 91 253 254 257 tr = func(self) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/geometry_handler.py:292, in YTDataChunk._accumulate_values(self, method) 290 291 --> 292 293 294 for obj in self._fast_index or self.objs: f = getattr(obj, mname) arrs.append(f(self.dobj)) if method == "dtcoords": arrs = [arr[0] for arr in arrs] self.index._identify_base_chunk(self) if field == "px": return self._current_chunk.fcoords[:, xax] elif field == "py": return self._current_chunk.fcoords[:, yax] File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/geometry_handler.py:255, in cached_property.<locals>.cached_func(self) return getattr(self, n) if self.data_size is None: tr = self._accumulate_values(n[1:]) 256 else: --> 255 File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/index_subobjects/particle_container.py:17, in _non_indexed.<locals>._func_non_indexed(self, *args, **kw args) 16 def _func_non_indexed(self, *args, **kwargs): ---> 17 raise YTNonIndexedDataContainer(self) YTNonIndexedDataContainer: The data container (<class 'yt.data_objects.index_subobjects.particle_container.ParticleContainer'>) is an unindexed type. Operations such as ires, icoords, fcoords and fwidth will not work on it. In terms of the bbox_lim, you are right I forgot to change it, it should be 1.5 pc. In regards to the temperature the code to real value should be the UnitTemp_in_K * InternalEnergy, so I changed that. However the units are code_specific_energy not kelvin. These changes are below. import yt from yt.units import dimensions UnitTemp_in_K = 3.568e+7 def my_temperature(field, data): return (data["Gas","InternalEnergy"] * UnitTemp_in_K) yt.add_field( ("my_temperature"), function=my_temperature, sampling_type="local", units='code_specific_energy', force_override=True
I realize the traceback error was sent in a really weird format, sorry about that this should be easier to read. YTNonIndexedDataContainer Traceback (most recent call last) Input In [25], in <module> ----> 1 slc = yt.SlicePlot(ds, "z", [("my_temperature")], width=(1.5,"pc")) 3 slc.show() File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:2528, in SlicePlot(ds, normal, fields, axis, *args, **kwargs) 2522 mylog.warning( 2523 "Ignoring 'north_vector' keyword as it is ill-defined for " 2524 "an AxisAlignedSlicePlot object." 2525 ) 2526 del kwargs["north_vector"] -> 2528 return AxisAlignedSlicePlot(ds, normal, fields, *args, **kwargs) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:1658, in AxisAlignedSlicePlot.__init__(self, ds, axis, fields, center, width, axes_unit, origin, right_handed, fontsize, field_parameters, window_size, aspect, data_source, buff_size) 1656 slc.get_data(fields) 1657 validate_mesh_fields(slc, fields) -> 1658 PWViewerMPL.__init__( 1659 self, 1660 slc, 1661 bounds, 1662 origin=origin, 1663 fontsize=fontsize, 1664 fields=fields, 1665 window_size=window_size, 1666 aspect=aspect, 1667 right_handed=right_handed, 1668 buff_size=buff_size, 1669 ) 1670 if axes_unit is None: 1671 axes_unit = get_axes_unit(width, ds) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:873, in PWViewerMPL.__init__(self, *args, **kwargs) 871 self._plot_type = kwargs.pop("plot_type") 872 self._splat_color = kwargs.pop("splat_color", None) --> 873 PlotWindow.__init__(self, *args, **kwargs) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:257, in PlotWindow.__init__(self, data_source, bounds, buff_size, antialias, periodic, origin, oblique, right_handed, window_size, fields, fontsize, aspect, setup) 254 self._units_config[field] 256 self.setup_callbacks() --> 257 self._setup_plots() File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:965, in PWViewerMPL._setup_plots(self) 963 return 964 if not self._data_valid: --> 965 self._recreate_frb() 966 self._data_valid = True 967 self._colorbar_valid = True File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/plot_window.py:317, in PlotWindow._recreate_frb(self) 315 # At this point the frb has the valid bounds, size, aliasing, etc. 316 if old_fields is None: --> 317 self._frb._get_data_source_fields() 319 # New frb, apply default units (if any) 320 for field, field_unit in self._units_config.items(): File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/fixed_resolution.py:176, in FixedResolutionBuffer._get_data_source_fields(self) 174 for f in fields: 175 if f not in exclude and f[0] not in self.data_source.ds.particle_types: --> 176 self[f] File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/visualization/fixed_resolution.py:139, in FixedResolutionBuffer.__getitem__(self, item) 136 b = float(b.in_units("code_length")) 137 bounds.append(b) --> 139 buff = self.ds.coordinates.pixelize( 140 self.data_source.axis, 141 self.data_source, 142 item, 143 bounds, 144 self.buff_size, 145 int(self.antialias), 146 ) 148 for name, (args, kwargs) in self._filters: 149 buff = filter_registry[name](*args[1:], **kwargs).apply(buff) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/coordinates/cartesian_coordinates.py:217, in CartesianCoordinateHandler.pixelize(self, dimension, data_source, field, bounds, size, antialias, periodic) 214 return np.squeeze(np.transpose(img, (yax, xax, ax))) 216 elif self.axis_id.get(dimension, dimension) < 3: --> 217 return self._ortho_pixelize( 218 data_source, field, bounds, size, antialias, dimension, periodic 219 ) 220 else: 221 return self._oblique_pixelize(data_source, field, bounds, size, antialias) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/coordinates/cartesian_coordinates.py:537, in CartesianCoordinateHandler._ortho_pixelize(self, data_source, field, bounds, size, antialias, dim, periodic) 533 buff = buff.transpose() 534 else: 535 pixelize_cartesian( 536 buff, --> 537 data_source["px"], 538 data_source["py"], 539 data_source["pdx"], 540 data_source["pdy"], 541 data_source[field], 542 bounds, 543 int(antialias), 544 period, 545 int(periodic), 546 ) 547 return buff File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/data_containers.py:255, in YTDataContainer.__getitem__(self, key) 253 if f not in self.field_data and key not in self.field_data: 254 if f in self._container_fields: --> 255 self.field_data[f] = self.ds.arr(self._generate_container_field(f)) 256 return self.field_data[f] 257 else: File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/selection_objects/slices.py:89, in YTSlice._generate_container_field(self, field) 87 self.index._identify_base_chunk(self) 88 if field == "px": ---> 89 return self._current_chunk.fcoords[:, xax] 90 elif field == "py": 91 return self._current_chunk.fcoords[:, yax] File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/geometry_handler.py:255, in cached_property.<locals>.cached_func(self) 253 return getattr(self, n) 254 if self.data_size is None: --> 255 tr = self._accumulate_values(n[1:]) 256 else: 257 tr = func(self) File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/geometry/geometry_handler.py:292, in YTDataChunk._accumulate_values(self, method) 290 for obj in self._fast_index or self.objs: 291 f = getattr(obj, mname) --> 292 arrs.append(f(self.dobj)) 293 if method == "dtcoords": 294 arrs = [arr[0] for arr in arrs] File ~/opt/miniconda3/envs/ciao-4.14/lib/python3.8/site-packages/yt/data_objects/index_subobjects/particle_container.py:17, in _non_indexed.<locals>._func_non_indexed(self, *args, **kwargs) 16 def _func_non_indexed(self, *args, **kwargs): ---> 17 raise YTNonIndexedDataContainer(self) YTNonIndexedDataContainer: The data container (<class 'yt.data_objects.index_subobjects.particle_container.ParticleContainer'>) is an unindexed type. Operations such as ires, icoords, fcoords and fwidth will not work on it.
participants (3)
-
John ZuHone -
Matthew Turk -
moonpigeona＠gmail.com