[BUGFIX] Fix particle selection for corner cases of morton index values

.travis.yml

@@ -1,5 +1,5 @@
 language: python
-dist: xenial
+dist: bionic
 cache:
   pip: true
   directories:

setup.py

@@ -106,7 +106,8 @@ def _compile(
               include_dirs=["yt/utilities/lib/",
                             "yt/utilities/lib/ewahboolarray"],
               language="c++",
-              libraries=std_libs),
+              libraries=std_libs,
+              extra_compile_args=["-std=c++11"]),
     Extension("yt.geometry.selection_routines",
               ["yt/geometry/selection_routines.pyx"],
               include_dirs=["yt/utilities/lib/"],

yt/frontends/gadget/tests/test_outputs.py

@@ -5,7 +5,8 @@
 import tempfile
 
 import yt
-from yt.testing import requires_file
+from yt.testing import requires_file, \
+    ParticleSelectionComparison
 from yt.utilities.answer_testing.framework import \
     data_dir_load, \
     requires_ds, \
@@ -107,6 +108,73 @@ def test_multifile_read():
     assert isinstance(data_dir_load(snap_33), GadgetDataset)
     assert isinstance(data_dir_load(snap_33_dir), GadgetDataset)
 
+@requires_file(snap_33)
+def test_particle_subselection():
+    #This checks that we correctly subselect from a dataset, first by making
+    #sure we get all the particles, then by comparing manual selections against
+    #them.
+    ds = data_dir_load(snap_33)
+    psc = ParticleSelectionComparison(ds)
+
+    sp1 = ds.sphere("c", (0.1, "unitary"))
+    psc.compare_dobj_selection(sp1)
+
+    sp2 = ds.sphere("c", (0.2, "unitary"))
+    psc.compare_dobj_selection(sp2)
+
+    # Test wrapping around each axis individually: x
+    sp3_x = ds.sphere((1.0, 12.5, 12.5), (2.0, "code_length"))
+    psc.compare_dobj_selection(sp3_x)
+
+    # Test wrapping around each axis individually: y
+    sp3_y = ds.sphere((12.5, 1.0, 12.5), (2.0, "code_length"))
+    psc.compare_dobj_selection(sp3_y)
+
+    # Test wrapping around each axis individually: z
+    sp3_z = ds.sphere((12.5, 12.5, 1.0), (2.0, "code_length"))
+    psc.compare_dobj_selection(sp3_z)
+
+    # Test wrapping around all three axes simultaneously on left
+    sp3_all = ds.sphere((1.0, 1.0, 1.0), (2.0, "code_length"))
+    psc.compare_dobj_selection(sp3_all)
+
+    # Test wrapping around each axis individually on right: x
+    sp4_x = ds.sphere((24.0, 12.5, 12.5), (2.0, "code_length"))
+    psc.compare_dobj_selection(sp4_x)
+
+    # Test wrapping around each axis individually on right: y
+    sp4_y = ds.sphere((12.5, 24.0, 12.5), (2.0, "code_length"))
+    psc.compare_dobj_selection(sp4_y)
+
+    # Test wrapping around each axis individually on right: z
+    sp4_z = ds.sphere((12.5, 12.5, 24.0), (2.0, "code_length"))
+    psc.compare_dobj_selection(sp4_z)
+
+    # Test wrapping around all three axes simultaneously on right
+    sp4_all = ds.sphere((24.0, 24.0, 24.0), (2.0, "code_length"))
+    psc.compare_dobj_selection(sp4_all)
+
+    sp5 = ds.sphere("c", (0.5, "unitary"))
+    psc.compare_dobj_selection(sp5)
+
+    dd = ds.all_data()
+    psc.compare_dobj_selection(dd)
+
+    reg1 = ds.r[ (0.1, 'unitary'):(0.9, 'unitary'),
+                 (0.1, 'unitary'):(0.9, 'unitary'),
+                 (0.1, 'unitary'):(0.9, 'unitary')]
+    psc.compare_dobj_selection(reg1)
+
+    reg2 = ds.r[ (0.8, 'unitary'):(0.85, 'unitary'),
+                 (0.8, 'unitary'):(0.85, 'unitary'),
+                 (0.8, 'unitary'):(0.85, 'unitary')]
+    psc.compare_dobj_selection(reg2)
+
+    reg3 = ds.r[ (0.3, 'unitary'):(0.6, 'unitary'),
+                 (0.2, 'unitary'):(0.8, 'unitary'),
+                 (0.0, 'unitary'):(0.1, 'unitary')]
+    psc.compare_dobj_selection(reg3)
+
 @requires_ds(BE_Gadget)
 def test_bigendian_field_access():
     ds = data_dir_load(BE_Gadget)

yt/geometry/particle_geometry_handler.py

@@ -167,21 +167,30 @@ def _initialize_refined_index(self):
         sub_mi1 = np.zeros(max_npart, "uint64")
         sub_mi2 = np.zeros(max_npart, "uint64")
         pb = get_pbar("Initializing refined index", len(self.data_files))
+        mask_threshold = getattr(self, '_index_mask_threshold', 2)
+        count_threshold = getattr(self, '_index_count_threshold', 256)
+        mylog.debug("Using estimated thresholds of %s and %s for refinement", mask_threshold, count_threshold)
+        total_refined = 0
+        total_coarse_refined = ((mask >= 2) & (self.regions.particle_counts > count_threshold)).sum()
+        mylog.debug("This should produce roughly %s zones, for %s of the domain",
+                    total_coarse_refined, 100 * total_coarse_refined / mask.size)
         for i, data_file in enumerate(self.data_files):
+            coll = None
             pb.update(i)
             nsub_mi = 0
             for ptype, pos in self.io._yield_coordinates(data_file):
+                if pos.size == 0: continue
                 if hasattr(self.ds, '_sph_ptypes') and ptype == self.ds._sph_ptypes[0]:
                     hsml = self.io._get_smoothing_length(
                         data_file, pos.dtype, pos.shape)
                 else:
                     hsml = None
-                nsub_mi = self.regions._refined_index_data_file(
-                    pos, hsml, mask, sub_mi1, sub_mi2,
-                    data_file.file_id, nsub_mi)
-            self.regions._set_refined_index_data_file(
-                sub_mi1, sub_mi2,
-                data_file.file_id, nsub_mi)
+                nsub_mi, coll = self.regions._refined_index_data_file(
+                    coll, pos, hsml, mask, sub_mi1, sub_mi2,
+                    data_file.file_id, nsub_mi, count_threshold = count_threshold,
+                    mask_threshold = mask_threshold)
+                total_refined += nsub_mi
+            self.regions.bitmasks.append(data_file.file_id, coll)
         pb.finish()
         self.regions.find_collisions_refined()