Gadget snapshot is determined to be 'cell-based' not 'particle-based'
Closed this issue · 10 comments
I'm having an issue where trident determines that my gadget is a 'cell based' simulation where it is actually a 'particle-based' simulation. We thing this has a significant affect on the density calculations of our rays.
I have tried running this code to diagnose the problems and force the code to use a sampling type of 'particle' but I get an error of 'ValueError: ('gas', 'smoothing_length') not in derived_field_list'.
How do I get trident to recognize that this is a SPH simulation and not a grid simulation.
Python Version: '3.6.4'
YT version: '3.5.dev0'
Trident Version: '1.1'
import yt
import trident
fn = "/Users/abatten/PhD/data/AURORA/L012N0128/Aurora_L012N0128_FSN1.0_FESC0.5/data/snapshot_060/snap_060.0.hdf5"
ds = yt.load(fn)
trident.add_ion_fields(ds, ions=['H I', 'He I', 'He II'], sampling_type='particle')
yt : [INFO ] 2018-09-06 11:51:46,130 Parameters: current_time = 3.035652263793517e+16 s
yt : [INFO ] 2018-09-06 11:51:46,131 Parameters: domain_dimensions = [2 2 2]
yt : [INFO ] 2018-09-06 11:51:46,132 Parameters: domain_left_edge = [0. 0. 0.]
yt : [INFO ] 2018-09-06 11:51:46,132 Parameters: domain_right_edge = [12.5 12.5 12.5]
yt : [INFO ] 2018-09-06 11:51:46,133 Parameters: cosmological_simulation = 1
yt : [INFO ] 2018-09-06 11:51:46,133 Parameters: current_redshift = 5.9999999999999964
yt : [INFO ] 2018-09-06 11:51:46,134 Parameters: omega_lambda = 0.735
yt : [INFO ] 2018-09-06 11:51:46,135 Parameters: omega_matter = 0.265
yt : [INFO ] 2018-09-06 11:51:46,135 Parameters: hubble_constant = 0.71
yt : [INFO ] 2018-09-06 11:51:46,152 Allocating for 4.194e+06 particles (index particle type 'all')
yt : [INFO ] 2018-09-06 11:51:46,594 Identified 2.952e+05 octs
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-3-44efd89e07c2> in <module>()
5 fn = "/Users/abatten/PhD/data/AURORA/L012N0128/Aurora_L012N0128_FSN1.0_FESC0.5/data/snapshot_060/snap_060.0.hdf5"
6 ds = yt.load(fn)
----> 7 trident.add_ion_fields(ds, ions=['H I', 'He I', 'He II'], sampling_type='particle')
~/anaconda3/lib/python3.6/site-packages/trident/ion_balance.py in add_ion_fields(ds, ions, ftype, ionization_table, field_suffix, line_database, force_override, sampling_type, particle_type)
328 add_ion_mass_field(atom, ion, ds, ftype, ionization_table,
329 field_suffix=field_suffix, force_override=force_override,
--> 330 sampling_type=sampling_type)
331
332 def add_ion_fraction_field(atom, ion, ds, ftype="gas",
~/anaconda3/lib/python3.6/site-packages/trident/ion_balance.py in add_ion_mass_field(atom, ion, ds, ftype, ionization_table, field_suffix, force_override, sampling_type, particle_type)
875 field_suffix=field_suffix,
876 force_override=force_override,
--> 877 sampling_type=sampling_type)
878 ds.add_field((ftype, field), function=_ion_mass, units=r"g",
879 sampling_type=sampling_type,
~/anaconda3/lib/python3.6/site-packages/trident/ion_balance.py in add_ion_density_field(atom, ion, ds, ftype, ionization_table, field_suffix, force_override, sampling_type, particle_type)
736 field_suffix=field_suffix,
737 force_override=force_override,
--> 738 sampling_type=sampling_type)
739 ds.add_field((ftype, field), function=_ion_density,
740 units="g/cm**3", sampling_type=sampling_type,
~/anaconda3/lib/python3.6/site-packages/trident/ion_balance.py in add_ion_number_density_field(atom, ion, ds, ftype, ionization_table, field_suffix, force_override, sampling_type, particle_type)
598 field_suffix=field_suffix,
599 force_override=force_override,
--> 600 sampling_type=sampling_type)
601 ds.add_field((ftype, field),function=_ion_number_density,
602 units="cm**-3", sampling_type=sampling_type,
~/anaconda3/lib/python3.6/site-packages/trident/ion_balance.py in add_ion_fraction_field(atom, ion, ds, ftype, ionization_table, field_suffix, force_override, sampling_type, particle_type)
471 # if ion particle field, add a smoothed deposited version to gas fields
472 if sampling_type == 'particle':
--> 473 new_field = ds.add_smoothed_particle_field((ftype, field))
474 if ftype != "gas":
475 ds.field_info.alias(("gas", field), new_field)
~/anaconda3/lib/python3.6/site-packages/yt/frontends/sph/data_structures.py in add_smoothed_particle_field(self, smooth_field, method, nneighbors, kernel_name)
71 return super(SPHDataset, self).add_smoothed_particle_field(
72 smooth_field=smooth_field, method=method, nneighbors=nneighbors,
---> 73 kernel_name=kernel_name
74 )
~/anaconda3/lib/python3.6/site-packages/yt/data_objects/static_output.py in add_smoothed_particle_field(self, smooth_field, method, nneighbors, kernel_name)
1347 if (ptype, smoothing_length_name) not in self.derived_field_list:
1348 raise ValueError("%s not in derived_field_list" %
-> 1349 ((ptype, smoothing_length_name),))
1350 density_name = "density"
1351 registry = self.field_info
ValueError: ('gas', 'smoothing_length') not in derived_field_list
As supplementary, I tried to add the smoothing length of the particles as a field by creating a ray (in the hdf5 file it is SmoothingLength) however I encounter the same error. It also here tells me that there is a mismatch between the SAMPLING_TYPE and FTYPE.
ray = trident.make_simple_ray(fn, start_position=[0.2,0.5, 0],
end_position=[0.2, 0.5, 12.5],
data_filename="ray.h5",
lines=["H", "He"],
fields=[('PartType0', 'SmoothingLength'),
('PartType0', 'Density')])
yt : [INFO ] 2018-09-06 12:01:36,594 Parameters: current_time = 3.035652263793517e+16 s
yt : [INFO ] 2018-09-06 12:01:36,597 Parameters: domain_dimensions = [2 2 2]
yt : [INFO ] 2018-09-06 12:01:36,599 Parameters: domain_left_edge = [0. 0. 0.]
yt : [INFO ] 2018-09-06 12:01:36,600 Parameters: domain_right_edge = [12.5 12.5 12.5]
yt : [INFO ] 2018-09-06 12:01:36,603 Parameters: cosmological_simulation = 1
yt : [INFO ] 2018-09-06 12:01:36,604 Parameters: current_redshift = 5.9999999999999964
yt : [INFO ] 2018-09-06 12:01:36,604 Parameters: omega_lambda = 0.735
yt : [INFO ] 2018-09-06 12:01:36,605 Parameters: omega_matter = 0.265
yt : [INFO ] 2018-09-06 12:01:36,605 Parameters: hubble_constant = 0.71
yt : [INFO ] 2018-09-06 12:01:36,612 Allocating for 4.194e+06 particles (index particle type 'all')
yt : [INFO ] 2018-09-06 12:01:37,061 Identified 2.952e+05 octs
yt : [WARNING ] 2018-09-06 12:01:39,780 ===================================================
yt : [WARNING ] 2018-09-06 12:01:39,781 MISMATCH BETWEEN SAMPLING_TYPE OF FTYPE AND DATASET
yt : [WARNING ] 2018-09-06 12:01:39,781 sampling_type of (gas, 'density') = cell
yt : [WARNING ] 2018-09-06 12:01:39,782 sampling_type of dataset = particle
yt : [WARNING ] 2018-09-06 12:01:39,782 THIS IS PROBABLY UNDESIRED BEHAVIOR. PLEASE CHOOSE A DIFFERENT FTYPE.
yt : [WARNING ] 2018-09-06 12:01:39,782 ===================================================
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-82-26d709e94ff5> in <module>()
4 lines=["H", "He"],
5 fields=[('PartType0', 'SmoothingLength'),
----> 6 ('PartType0', 'Density')])
~/anaconda3/lib/python3.6/site-packages/trident/ray_generator.py in make_simple_ray(dataset_file, start_position, end_position, lines, ftype, fields, solution_filename, data_filename, trajectory, redshift, setup_function, load_kwargs, line_database, ionization_table)
248 ftype=ftype,
249 ionization_table=ionization_table,
--> 250 sampling_type=sampling_type)
251
252 # To assure there are no fields that are double specified or that collide
~/anaconda3/lib/python3.6/site-packages/trident/ion_balance.py in add_ion_number_density_field(atom, ion, ds, ftype, ionization_table, field_suffix, force_override, sampling_type, particle_type)
598 field_suffix=field_suffix,
599 force_override=force_override,
--> 600 sampling_type=sampling_type)
601 ds.add_field((ftype, field),function=_ion_number_density,
602 units="cm**-3", sampling_type=sampling_type,
~/anaconda3/lib/python3.6/site-packages/trident/ion_balance.py in add_ion_fraction_field(atom, ion, ds, ftype, ionization_table, field_suffix, force_override, sampling_type, particle_type)
471 # if ion particle field, add a smoothed deposited version to gas fields
472 if sampling_type == 'particle':
--> 473 new_field = ds.add_smoothed_particle_field((ftype, field))
474 if ftype != "gas":
475 ds.field_info.alias(("gas", field), new_field)
~/anaconda3/lib/python3.6/site-packages/yt/frontends/sph/data_structures.py in add_smoothed_particle_field(self, smooth_field, method, nneighbors, kernel_name)
71 return super(SPHDataset, self).add_smoothed_particle_field(
72 smooth_field=smooth_field, method=method, nneighbors=nneighbors,
---> 73 kernel_name=kernel_name
74 )
~/anaconda3/lib/python3.6/site-packages/yt/data_objects/static_output.py in add_smoothed_particle_field(self, smooth_field, method, nneighbors, kernel_name)
1347 if (ptype, smoothing_length_name) not in self.derived_field_list:
1348 raise ValueError("%s not in derived_field_list" %
-> 1349 ((ptype, smoothing_length_name),))
1350 density_name = "density"
1351 registry = self.field_info
ValueError: ('gas', 'smoothing_length') not in derived_field_list
Hi @abatten . Thanks for the message and for joining our community!
You probably want to specify the ftype keyword in add_ion_fields() as the field type in your dataset that pertains to the actual gas particles. Probably for Gadget data, this is PartType0. This is mentioned a few times in the docs (in the API docs for add_ion_fields(): https://trident.readthedocs.io/en/latest/generated/trident.add_ion_fields.html#trident.add_ion_fields, and in the docs on adding ion fields: https://trident.readthedocs.io/en/latest/ion_balance.html), but I should probably add an entry to the FAQ about this. But I believe if you set ftype correctly, what you're trying to do should work.
It's also noting that in the current stable release of Trident, particle-based datasets like Gadget are deposited to an octree grid structure and then treated as grid-based datasets for sampling. This is described at length in the method paper, and it can sometimes introduce artifacts in regions of poor resolution. However, there has been substantial work from both the yt and trident communities to treat particle-based datasets natively without depositing to a grid structure called the "demeshening." There are full instructions on installing and using the demeshened versions of yt and trident here (it's linked on our installation page): https://nbviewer.jupyter.org/url/trident-project.org/notebooks/trident_demesh_install.ipynb
I hope this helps with your issue, but please let us know if you continue to have problems.
Hi @chummels .
Thanks for the tip about the demeshed versions of YT and Trident. I was able to download and install them. It seems to be almost working as I would like.
I had a few questions still:
-
What is the difference between 'H_p1_number_density' and 'H_nuclei_density'? When I try to generate a ray with "H" in the line list it only creates 'H_nuclei_number_density' and 'H_number_density' in the output ray.h5, however when I use the example gizmo script on the FIRE data it creates both 'H_nuclei_number_density' and 'H_p1_number_density' with different values.
-
Secondly I'm having problems getting the ray to output 'He_p2_number_density' in the ray. It seems to work on the FIRE data but, for my simulations it doesn't create the the field in my simulations.
import yt import trident fn = "/Users/abatten/PhD/data/AURORA/L012N0128/Aurora_L012N0128_FSN1.0_FESC0.5/data/snapshot_060/snap_060.0.hdf5" ds = yt.load(fn) xi, yi, zi = 10, 1.00, 0.0 xf, yf, zf = 10, 1.00, 12.5 ray_start = [xi, yi, zi] ray_end = [xf, yf, zf] ray = trident.make_simple_ray(fn, start_position=ray_start, end_position=ray_end, lines=['H', 'He', 'H I', 'H II', 'He I', 'He II', 'He III'], ftype='PartType0', data_filename='test_ray.h5')
I have also tried using the following, but this leads to a "RecursionError: maximum recursion depth exceeded while calling a Python object"
trident.add_ion_number_density_field("He", 3, ds, ftype="PartType0", force_override=True)
yt : [WARNING ] 2018-09-10 16:31:28,313 Field ('gas', 'He_p2_ion_fraction') already exists. To override use force_override=True.
---------------------------------------------------------------------------
RecursionError Traceback (most recent call last)
<ipython-input-50-98226b082339> in <module>()
----> 1 trident.add_ion_number_density_field("He", 3, ds, ftype="PartType0", force_override=True)
2 ray = trident.make_simple_ray(fn, start_position=[0.2,0.5, 0],
3 end_position=[0.2, 0.5, 12.5],
4 data_filename="ray.h5",
5 lines=["H", "He", "C"],
~/Packages/trident/trident/ion_balance.py in add_ion_number_density_field(atom, ion, ds, ftype, ionization_table, field_suffix, force_override, sampling_type, particle_type)
508
509 add_ion_fraction_field(atom, ion, ds, ionization_table,
--> 510 field_suffix=field_suffix)
511 ds.add_field(("gas", field),function=_ion_number_density,
512 units="cm**-3", sampling_type='local',
~/Packages/trident/trident/ion_balance.py in add_ion_fraction_field(atom, ion, ds, ftype, ionization_table, field_suffix, force_override, sampling_type, particle_type)
405
406 ds.add_field(("gas", field), function=_ion_fraction_field, units="",
--> 407 sampling_type="local", force_override=force_override)
408 if ion == 1: # add aliased field too
409 ds.field_info.alias(("gas", alias_field), ("gas", field))
~/Packages/yt/yt/data_objects/static_output.py in add_field(self, name, function, sampling_type, **kwargs)
1252 self.field_info.add_field(name, sampling_type, function=function, **kwargs)
1253 self.field_info._show_field_errors.append(name)
-> 1254 deps, _ = self.field_info.check_derived_fields([name])
1255 self.field_dependencies.update(deps)
1256
~/Packages/yt/yt/fields/field_info_container.py in check_derived_fields(self, fields_to_check)
377 fi = self[field]
378 try:
--> 379 fd = fi.get_dependencies(ds = self.ds)
380 except Exception as e:
381 if field in self._show_field_errors:
~/Packages/yt/yt/fields/derived_field.py in get_dependencies(self, *args, **kwargs)
216 e.requested.append(self.name)
217 else:
--> 218 e[self.name]
219 return e
220
~/Packages/yt/yt/fields/field_detector.py in __missing__(self, item)
106 vv = finfo(self)
107 if not permute_params:
--> 108 vv = finfo(self)
109 except NeedsGridType as exc:
110 ngz = exc.ghost_zones
~/Packages/yt/yt/fields/derived_field.py in __call__(self, data)
256 "for %s" % (self.name,))
257 with self.unit_registry(data):
--> 258 dd = self._function(self, data)
259 for field_name in data.keys():
260 if field_name not in original_fields:
~/Packages/trident/trident/ion_balance.py in _ion_fraction_field(field, data)
851 raise RuntimeError("This data file format is not supported.")
852
--> 853 fraction = np.power(10, interp(data))
854 fraction[fraction <= fraction_zero_point] = 0.0
855 if not isinstance(data, FieldDetector) and (fraction > 1.0).any():
~/Packages/yt/yt/utilities/linear_interpolators.py in __call__(self, data_object)
207
208 def __call__(self, data_object):
--> 209 orig_shape = data_object[self.x_name].shape
210 x_vals = data_object[self.x_name].ravel().astype('float64')
211 y_vals = data_object[self.y_name].ravel().astype('float64')
~/Packages/yt/yt/fields/field_detector.py in __missing__(self, item)
106 vv = finfo(self)
107 if not permute_params:
--> 108 vv = finfo(self)
109 except NeedsGridType as exc:
110 ngz = exc.ghost_zones
~/Packages/yt/yt/fields/derived_field.py in __call__(self, data)
256 "for %s" % (self.name,))
257 with self.unit_registry(data):
--> 258 dd = self._function(self, data)
259 for field_name in data.keys():
260 if field_name not in original_fields:
~/Packages/trident/trident/ion_balance.py in _log_nH(field, data)
101 """
102 if ("gas", "H_nuclei_density") in data.ds.derived_field_list:
--> 103 log_nH_field = np.log10(data["gas", "H_nuclei_density"])
104 else:
105 log_nH_field = np.log10(data["gas", "density"] * to_nH)
~/Packages/yt/yt/fields/field_detector.py in __missing__(self, item)
106 vv = finfo(self)
107 if not permute_params:
--> 108 vv = finfo(self)
109 except NeedsGridType as exc:
110 ngz = exc.ghost_zones
~/Packages/yt/yt/fields/derived_field.py in __call__(self, data)
256 "for %s" % (self.name,))
257 with self.unit_registry(data):
--> 258 dd = self._function(self, data)
259 for field_name in data.keys():
260 if field_name not in original_fields:
~/Packages/yt/yt/fields/species_fields.py in _nuclei_density(field, data)
190
191 field_data = np.zeros_like(data[ftype, "%s_number_density" %
--> 192 data.ds.field_info.species_names[0]])
193 for species in data.ds.field_info.species_names:
194 nucleus = species
~/Packages/yt/yt/fields/field_detector.py in __missing__(self, item)
106 vv = finfo(self)
107 if not permute_params:
--> 108 vv = finfo(self)
109 except NeedsGridType as exc:
110 ngz = exc.ghost_zones
~/Packages/yt/yt/fields/derived_field.py in __call__(self, data)
256 "for %s" % (self.name,))
257 with self.unit_registry(data):
--> 258 dd = self._function(self, data)
259 for field_name in data.keys():
260 if field_name not in original_fields:
~/Packages/trident/trident/ion_balance.py in _ion_number_density(field, data)
801
802 if atom == 'H' or atom == 'He':
--> 803 number_density = solar_abundance[atom] * data[fraction_field_name]
804 else:
805 number_density = data.ds.quan(solar_abundance[atom], "1.0/Zsun") * \
... last 13 frames repeated, from the frame below ...
~/Packages/yt/yt/fields/field_detector.py in __missing__(self, item)
106 vv = finfo(self)
107 if not permute_params:
--> 108 vv = finfo(self)
109 except NeedsGridType as exc:
110 ngz = exc.ghost_zones
RecursionError: maximum recursion depth exceeded while calling a Python object
Also thirdly, I have a question about the dl,
I created a frequency distribution of the various lengths of dl from 163 rays. We were confused as to why the dl is sharply peaked near zero, as we were expecting something more gaussian like. How is the dl calculated for the demeshed Trident?
We were guessing that the dl is the length piercing the through individual particles (i.e so dl = 2*r_s if the ray goes through the centre of the particle, with r_s being the smoothing length and dl ~ 0, when the ray clips the edge of a particle.) However it would be unlikely for a ray to clip the edge of many particles as opposed to be somewhere in the middle. Which is why we didn't expect a peak near zero for dl.
Finally is the sph kernel integrated for weighting densities when creating rays?
Hi @abatten . Glad to hear things are almost working for your use case!
To answer your questions, H_p1_number_density refers to the number density of H II. (p1 stands for plus one). H_nuclei_density is the number density of all hydrogen (both neutral and ionized). For reference, H_p0_number_density and H_number_density are both the neutral hydrogen (H I) number density. It's unfortunate, but in the initial implementation in yt, H_number_density was just set up to be neutral hydrogen, not total hydrogen, so in order to retain backwards compatibility it has thus remained as such, and we just alias H_p0_number_density (H plus 0 = H I) to it. Confusing, I know. Maybe in the next full release we'll break backwards compatibility to fix this problem...
Hmm, that's odd about H II fields not being generated for your dataset. Have you just tried explicitly adding them with:
import trident
trident.add_ion_fields(ds, 'H II')
As for your problem with He III, again, this seems odd. Can you again try running the code above but replacing H II with He III and see if it works?
To address your final question about the dl values peaking close to zero for a sightline passing through your SPH volume, I think this makes sense. In a grid-based dataset, the sightline just passes through different cells, and so dl is simply the geometric pathlength over which it intersects those cells. However, in the demeshening version of trident/yt, those path lengths are not simple geometric lengths representing the amount of the SPH particle that has been intersected by the path. It also folds in information about the SPH kernel itself, since that gives the weighting to whatever particle-field you're sampling. So to get the integrated column density of H I along a line of sight, it is still sum(n_HI * dl), but now the dl includes the line integral through the kernel at whatever impact parameter the sightline passes relative to the particle center.
So it's like what you suggest, but weighted with by the smoothing kernel function, which falls off with radius. So I think this will bias sightlines hitting a particle off-center to have smaller and smaller dl values, effectively making it look like a lot of stuff is grazing.
But perhaps @ngoldbaum can speak more to this result, since he was one of the main architects of the demeshening in yt.
Thanks @chummels , that explaination of dl was what we wanted. However I'm still having the issue with HII and HeIII. I tried adding them explicitly with what you suggested, however it gives me erros that the ion fields already exist but doesn't add them to the ray. I also tried adding fields for C I , C II instead of H II and He III which seems to work, just not for H II and HeIII
import yt
import trident
xi, yi, zi = 10, 1.00, 0.0
xf, yf, zf = 10, 1.00, 12.5
ray_start = [xi, yi, zi]
ray_end = [xf, yf, zf]
trident.add_ion_fields(ds, ions=['H II'])
trident.add_ion_fields(ds, ions=['He III'])
sray = trident.make_simple_ray(fn,
start_position=ray_start,
end_position=ray_end,
lines=['H', 'He', 'H I', 'H II', 'He I', 'He II', 'He III'],
ftype='PartType0',
data_filename='test_ray.h5')
These are the warnings it keeps giving me about the ion fields already existing. And using the force_override=True gives that maximum recursion depth exceeded while calling a Python object error from before.
yt : [WARNING ] 2018-09-12 15:47:28,145 Field ('gas', 'H_p1_ion_fraction') already exists. To override use force_override=True.
yt : [WARNING ] 2018-09-12 15:47:28,156 Field ('gas', 'H_p1_number_density') already exists. To override use force_override=True.
yt : [WARNING ] 2018-09-12 15:47:28,164 Field ('gas', 'H_p1_density') already exists. To override use force_override=True.
yt : [WARNING ] 2018-09-12 15:47:28,171 Field ('gas', 'H_p1_mass') already exists. To override use force_override=True.
yt : [WARNING ] 2018-09-12 15:47:28,178 Field ('gas', 'He_p2_ion_fraction') already exists. To override use force_override=True.
yt : [WARNING ] 2018-09-12 15:47:28,184 Field ('gas', 'He_p2_number_density') already exists. To override use force_override=True.
yt : [WARNING ] 2018-09-12 15:47:28,191 Field ('gas', 'He_p2_density') already exists. To override use force_override=True.
yt : [WARNING ] 2018-09-12 15:47:28,201 Field ('gas', 'He_p2_mass') already exists. To override use force_override=True.
It looks like from your example, you're not actually using the dataset, ds, to generate your sray object, since you pass it the filename of the dataset, not the actual dataset that you've processed. Try changing this line from:
sray = trident.make_simple_ray(fn,
to
sray = trident.make_simple_ray(ds,
so trident uses the dataset that you already added the H II and He III lines to. If this doesn't address the problem for you, then I'm not sure how to fix it. I may have to play with your dataset. Alternatively, you could try to reproduce your bug with one of the public datasets, like the FIRE snapshot referenced in the docs. But first let me know if the solution above fixes your issue.
@abatten : Did my suggestions address your problem? If so, can you close this issue? If not, let me know and we can further address things.
Yes I think this worked! Thank you.