Issue regarding the mesh_sampling_particle_field with Ramses

Hi all,

I tried to use the mesh_sampling_particle_field with ramses following https://yt-project.org/doc/analyzing/fields.html#mesh-sampling-particle-fields coupled with the bbox from the same page. I faced some issues, which can be reproduced with output_00080 from https://yt-project.org/data/.

here is the code I tried:

import numpy as np
import yt
import pyrats
#mesh particle sampling OK
ds = yt.load('output_00080/info_00080.txt')
ds.add_mesh_sampling_particle_field(('gas', 'temperature'), ptype='all')
ds.r['all', 'cell_index']
ds.r['all', 'cell_gas_temperature']

=> YTArray([ 932633.36800546, 11671881.80641624, 2462189.62059189, ..., 55344.84089513, 23934.30031332, 67203.92686486]) K

#bug with bbox
center = [0.5,0.5,0.5]
ds = yt.load('output_00080/info_00080.txt')
radius = 5
w = ds.arr(radius, 'Mpc').to('code_length').v
bbox = [list(center-w), list(center+w)]
ds._bbox = bbox
ds.add_mesh_sampling_particle_field(('gas', 'temperature'), ptype='all')
ds.r['all', 'cell_index']
ds.r['all', 'cell_gas_temperature']

yt : [INFO     ] 2020-05-03 17:39:09,185 Parameters: current_time              = 1.119216564055017
yt : [INFO     ] 2020-05-03 17:39:09,189 Parameters: domain_dimensions         = [64 64 64]
yt : [INFO     ] 2020-05-03 17:39:09,192 Parameters: domain_left_edge          = [0. 0. 0.]
yt : [INFO     ] 2020-05-03 17:39:09,195 Parameters: domain_right_edge         = [1. 1. 1.]
yt : [INFO     ] 2020-05-03 17:39:09,199 Parameters: cosmological_simulation   = 1
yt : [INFO     ] 2020-05-03 17:39:09,202 Parameters: current_redshift          = 0.14255728632206321
yt : [INFO     ] 2020-05-03 17:39:09,205 Parameters: omega_lambda              = 0.723999977111816
yt : [INFO     ] 2020-05-03 17:39:09,207 Parameters: omega_matter              = 0.276000022888184
yt : [INFO     ] 2020-05-03 17:39:09,210 Parameters: omega_radiation           = 0.0
yt : [INFO     ] 2020-05-03 17:39:09,211 Parameters: hubble_constant           = 0.703000030517578
yt : [WARNING  ] 2020-05-03 17:39:09,227 Detected 2 extra particle fields assuming kind `double`. Consider using the `extra_particle_fields` keyword argument if you have unexpected behavior.
yt : [WARNING  ] 2020-05-03 17:39:09,281 This cooling file format is no longer supported. Cooling field loading skipped.
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-18-9fd7755fd3ae> in <module>()
      8 ds.add_mesh_sampling_particle_field(('gas', 'temperature'), ptype='all')
      9 ds.r['all', 'cell_index']
---> 10 ds.r['all', 'cell_gas_temperature']

~/Codes/yt/yt/data_objects/region_expression.py in __getitem__(self, item)
     40             return self.all_data[item]
     41         if isinstance(item, tuple) and isinstance(item[1], string_types):
---> 42             return self.all_data[item]
     43         if isinstance(item, slice):
     44             if obj_length(item.start) == 3 and obj_length(item.stop) == 3:

~/Codes/yt/yt/data_objects/data_containers.py in __getitem__(self, key)
    254                 return self.field_data[f]
    255             else:
--> 256                 self.get_data(f)
    257         # fi.units is the unit expression string. We depend on the registry
    258         # hanging off the dataset to define this unit object.

~/Codes/yt/yt/data_objects/data_containers.py in get_data(self, fields)
   1582 
   1583         fields_to_generate += gen_fluids + gen_particles
-> 1584         self._generate_fields(fields_to_generate)
   1585         for field in list(self.field_data.keys()):
   1586             if field not in ofields:

~/Codes/yt/yt/data_objects/data_containers.py in _generate_fields(self, fields_to_generate)
   1602                 fi = self.ds._get_field_info(*field)
   1603                 try:
-> 1604                     fd = self._generate_field(field)
   1605                     if fd is None:
   1606                         raise RuntimeError

~/Codes/yt/yt/data_objects/data_containers.py in _generate_field(self, field)
    289                 tr = self._generate_container_field(field)
    290             if finfo.particle_type: # This is a property now
--> 291                 tr = self._generate_particle_field(field)
    292             else:
    293                 tr = self._generate_fluid_field(field)

~/Codes/yt/yt/data_objects/data_containers.py in _generate_particle_field(self, field)
    374         else:
    375             with self._field_type_state(ftype, finfo, gen_obj):
--> 376                 rv = self.ds._get_field_info(*field)(gen_obj)
    377         return rv
    378 

~/Codes/yt/yt/fields/derived_field.py in __call__(self, data)
    254                 "for %s" % (self.name,))
    255         with self.unit_registry(data):
--> 256             dd = self._function(self, data)
    257         for field_name in data.keys():
    258             if field_name not in original_fields:

~/Codes/yt/yt/geometry/oct_geometry_handler.py in _mesh_sampling_particle_field(field, data)
    100                 cell_data = subset[ftype, deposit_field].T.reshape(-1)
    101 
--> 102                 ret[mask] = cell_data[icell[mask]]
    103 
    104             return data.ds.arr(ret, input_units=cell_data.units)

~/Codes/yt/yt/units/yt_array.py in __getitem__(self, item)
   1063 
   1064     def __getitem__(self, item):
-> 1065         ret = super(YTArray, self).__getitem__(item)
   1066         if ret.shape == ():
   1067             return YTQuantity(ret, self.units, bypass_validation=True)

IndexError: index -4611686018427387904 is out of bounds for axis 1 with size 194632

It's not a big issue, as it works when the whole box is loaded, but it would be much faster if we could use both options together!

Thanks for your work! Hugo

yt-project / yt

Issue regarding the mesh_sampling_particle_field with Ramses #2573