yt
yt copied to clipboard
yt-project
error with demeshing yt in accessing deposited fields in arbitrary grids
Bug report
Bug summary
When creating an arbitrary grid with demeshed yt, one cannot access the deposited gas mass from an SPH simulation.
Code for reproduction
from __future__ import print_function
import yt
fname = '/users/desika.narayanan/Dropbox/yt_datasets/GadgetDiskGalaxy/snapshot_200.hdf5'
ds = yt.load(fname)
#from the notebook
center = [ 31995.63476562,31473.6640625,28969.88671875]
arb_center = ds.arr(center,'code_length')
left_edge = arb_center-ds.quan(100,'kpc')
right_edge = arb_center+ds.quan(100,'kpc')
arb = ds.arbitrary_grid(left_edge, right_edge,dims=[128, 128, 128],field_parameters={'center': center})
mass = arb[('deposit', 'PartType0_mass')].in_units('Msun')
Actual outcome
In [8]: run arb
yt : [INFO ] 2018-09-04 06:30:37,956 Calculating time from 3.448e-01 to be 1.108e+17 seconds
yt : [INFO ] 2018-09-04 06:30:37,956 Assuming length units are in kpc/h (comoving)
yt : [INFO ] 2018-09-04 06:30:37,968 Parameters: current_time = 1.1075810732534829e+17 s
yt : [INFO ] 2018-09-04 06:30:37,968 Parameters: domain_dimensions = [1 1 1]
yt : [INFO ] 2018-09-04 06:30:37,968 Parameters: domain_left_edge = [0. 0. 0.]
yt : [INFO ] 2018-09-04 06:30:37,969 Parameters: domain_right_edge = [64000. 64000. 64000.]
yt : [INFO ] 2018-09-04 06:30:37,969 Parameters: cosmological_simulation = 1
yt : [INFO ] 2018-09-04 06:30:37,969 Parameters: current_redshift = 1.8999965286929705
yt : [INFO ] 2018-09-04 06:30:37,969 Parameters: omega_lambda = 0.721
yt : [INFO ] 2018-09-04 06:30:37,969 Parameters: omega_matter = 0.279
yt : [INFO ] 2018-09-04 06:30:37,969 Parameters: hubble_constant = 0.7
yt : [INFO ] 2018-09-04 06:30:38,150 Allocating for 1.191e+07 particles
Loading particle index: 100%|██████████████████| 19/19 [00:00<00:00, 216.84it/s]
Interpolating SPH field ('deposit', 'PartType0_mass'): 0it [00:00, ?it/s]---------------------------------------------------------------------------
NeedsGridType Traceback (most recent call last)
~/yt/yt/data_objects/data_containers.py in _generate_fluid_field(self, field)
309 try:
--> 310 finfo.check_available(gen_obj)
311 except NeedsGridType as ngt_exception:
~/yt/yt/fields/derived_field.py in check_available(self, data)
205 for validator in self.validators:
--> 206 validator(data)
207 # If we don't get an exception, we're good to go
~/yt/yt/fields/derived_field.py in __call__(self, data)
451 if not getattr(data, '_spatial', False):
--> 452 raise NeedsGridType(self.ghost_zones,self.fields)
453 if self.ghost_zones <= data._num_ghost_zones:
NeedsGridType: (0, None)
During handling of the above exception, another exception occurred:
YTNonIndexedDataContainer Traceback (most recent call last)
~/Desktop/scratch/arb.py in <module>()
14
15 arb = ds.arbitrary_grid(left_edge, right_edge,dims=[128, 128, 128],field_parameters={'center': center})
---> 16 mass = arb[('deposit', 'PartType0_mass')].in_units('Msun')
~/yt/yt/data_objects/data_containers.py in __getitem__(self, key)
255 return self.field_data[f]
256 else:
--> 257 self.get_data(f)
258 # fi.units is the unit expression string. We depend on the registry
259 # hanging off the dataset to define this unit object.
~/yt/yt/data_objects/construction_data_containers.py in get_data(self, fields)
684 if len(part) > 0 and len(alias) == 0:
685 if(is_sph_field):
--> 686 self._fill_sph_particles(fields)
687 for field in fields:
688 if field in gen:
~/yt/yt/data_objects/construction_data_containers.py in _fill_sph_particles(self, fields)
756
757 pbar = tqdm(desc="Interpolating SPH field {}".format(field))
--> 758 for chunk in self._data_source.chunks([field],"io"):
759 px = chunk[(ptype,'particle_position_x')].in_base("code").d
760 py = chunk[(ptype,'particle_position_y')].in_base("code").d
~/yt/yt/data_objects/data_containers.py in chunks(self, fields, chunking_style, **kwargs)
1297 continue
1298 with self._chunked_read(chunk):
-> 1299 self.get_data(fields)
1300 # NOTE: we yield before releasing the context
1301 yield self
~/yt/yt/data_objects/data_containers.py in get_data(self, fields)
1398
1399 fields_to_generate += gen_fluids + gen_particles
-> 1400 self._generate_fields(fields_to_generate)
1401 for field in list(self.field_data.keys()):
1402 if field not in ofields:
~/yt/yt/data_objects/data_containers.py in _generate_fields(self, fields_to_generate)
1418 fi = self.ds._get_field_info(*field)
1419 try:
-> 1420 fd = self._generate_field(field)
1421 if fd is None:
1422 raise RuntimeError
~/yt/yt/data_objects/data_containers.py in _generate_field(self, field)
292 tr = self._generate_particle_field(field)
293 else:
--> 294 tr = self._generate_fluid_field(field)
295 if tr is None:
296 raise YTCouldNotGenerateField(field, self.ds)
~/yt/yt/data_objects/data_containers.py in _generate_fluid_field(self, field)
310 finfo.check_available(gen_obj)
311 except NeedsGridType as ngt_exception:
--> 312 rv = self._generate_spatial_fluid(field, ngt_exception.ghost_zones)
313 else:
314 rv = finfo(gen_obj)
~/yt/yt/data_objects/data_containers.py in _generate_spatial_fluid(self, field, ngz)
337 o = self._current_chunk.objs[0]
338 if accumulate:
--> 339 rv = self.ds.arr(np.empty(o.ires.size, dtype="float64"),
340 units)
341 outputs.append(rv)
~/yt/yt/data_objects/data_containers.py in ires(self)
1575 if self._current_chunk is None:
1576 self.index._identify_base_chunk(self)
-> 1577 return self._current_chunk.ires
1578
1579 @property
~/yt/yt/geometry/geometry_handler.py in cached_func(self)
267 return getattr(self, n)
268 if self.data_size is None:
--> 269 tr = self._accumulate_values(n[1:])
270 else:
271 tr = func(self)
~/yt/yt/geometry/geometry_handler.py in _accumulate_values(self, method)
295 for obj in self._fast_index or self.objs:
296 f = getattr(obj, mname)
--> 297 arrs.append(f(self.dobj))
298 if method == "dtcoords":
299 arrs = [arr[0] for arr in arrs]
~/yt/yt/data_objects/particle_container.py in _func_non_indexed(self, *args, **kwargs)
28 def _non_indexed(name):
29 def _func_non_indexed(self, *args, **kwargs):
---> 30 raise YTNonIndexedDataContainer(self)
31 return _func_non_indexed
32
YTNonIndexedDataContainer: The data container (<class 'yt.data_objects.particle_container.ParticleContainer'>) is an unindexed type. Operations such as ires, icoords, fcoords and fwidth will not work on it.
Expected outcome
Version Information
- Operating System: OSX and/or Linux
- Python Version:
In [2]: sys.version Out[2]: '3.6.5 |Anaconda, Inc.| (default, Apr 26 2018, 08:42:37) \n[GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]' - yt version: 4.0
- Other Libraries (if applicable):
installed from source
I'm interested in looking into this. I know what the issue is, but I don't know the solution. I'm happy to chat with someone about this
I just pulled this and it gives me a different error now:
File "/home/mturk/yt/yt/yt/data_objects/construction_data_containers.py", line 890, in _fill_sph_particles
raise KeyError(f"{ptype} is not a SPH particle type!")
KeyError: 'deposit is not a SPH particle type!'
but, if I ask for:
arb["PartType0", "mass"]
it works. So I think this is ready to be closed, if you're OK with that, @dnarayanan .
I think arb["PartType0", "mass"] returns the SPH smoothed mass. But, sometimes you might want to use some of the existing deposit machinery, but the SPH changes prevent that from happening. I assume arb[('deposit', 'PartType1_mass')] still works. This is directly related to #2004
