Backport PR #4547 on branch yt-4.2.x (BUG: Fix initialization of max level in load_octree)

yt/frontends/stream/data_structures.py

@@ -866,6 +866,8 @@ def _initialize_oct_handler(self):
             "partial_coverage": self.ds.partial_coverage,
         }
         self.oct_handler = OctreeContainer.load_octree(header)
+        # We do now need to get the maximum level set, as well.
+        self.ds.max_level = self.oct_handler.max_level
 
     def _identify_base_chunk(self, dobj):
         if getattr(dobj, "_chunk_info", None) is None:

yt/frontends/stream/tests/test_stream_octree.py

@@ -1,4 +1,5 @@
 import numpy as np
+from numpy.testing import assert_equal
 
 import yt
 
@@ -36,7 +37,7 @@ def test_octree():
     octree_mask = np.array(OCT_MASK_LIST, dtype=np.uint8)
 
     quantities = {}
-    quantities[("gas", "density")] = np.ones((22, 1), dtype="float64")
+    quantities[("gas", "density")] = np.random.random((22, 1))
 
     bbox = np.array([[-10.0, 10.0], [-10.0, 10.0], [-10.0, 10.0]])
 
@@ -50,3 +51,10 @@ def test_octree():
 
     proj = ds.proj(("gas", "density"), "x")
     proj[("gas", "density")]
+
+    assert_equal(ds.r[:]["ones"].size, 22)
+    rho1 = quantities["gas", "density"].ravel()
+    rho2 = ds.r[:]["density"].copy()
+    rho1.sort()
+    rho2.sort()
+    assert_equal(rho1, rho2)

yt/geometry/oct_container.pxd

@@ -80,6 +80,7 @@ cdef class OctreeContainer:
     # The fill_style is the ordering, C or F, of the octs in the file.  "o"
     # corresponds to C, and "r" is for Fortran.
     cdef public object fill_style
+    cdef public int max_level
 
 cdef class SparseOctreeContainer(OctreeContainer):
     cdef OctKey *root_nodes

yt/geometry/oct_container.pyx

@@ -92,6 +92,7 @@ cdef class OctreeContainer:
         visitor.global_index = -1
         visitor.level = 0
         visitor.nz = visitor.nzones = 1
+        visitor.max_level = 0
         assert(ref_mask.shape[0] / float(visitor.nzones) ==
             <int>(ref_mask.shape[0]/float(visitor.nzones)))
         obj.allocate_domains([ref_mask.shape[0] / visitor.nzones])
@@ -135,6 +136,7 @@ cdef class OctreeContainer:
         if obj.nocts * visitor.nz != ref_mask.size:
             raise KeyError(ref_mask.size, obj.nocts, obj.nz,
                 obj.partial_coverage, visitor.nzones)
+        obj.max_level = visitor.max_level
         return obj
 
     def __dealloc__(self):

yt/geometry/oct_visitors.pxd

@@ -134,6 +134,7 @@ cdef class LoadOctree(OctVisitor):
     cdef Oct* octs
     cdef np.uint64_t *nocts
     cdef np.uint64_t *nfinest
+    cdef np.uint64_t max_level
 
 cdef class MortonIndexOcts(OctVisitor):
     cdef np.uint8_t[:] level_arr

yt/geometry/oct_visitors.pyx

@@ -294,6 +294,8 @@ cdef class LoadOctree(OctVisitor):
     cdef void visit(self, Oct* o, np.uint8_t selected):
         cdef int i, ii
         ii = cind(self.ind[0], self.ind[1], self.ind[2])
+        if self.level > self.max_level:
+            self.max_level = self.level
         if self.ref_mask[self.index] == 0:
             # We only want to do this once.  Otherwise we end up with way too many
             # nfinest for our tastes.

yt/geometry/particle_oct_container.pyx

@@ -75,7 +75,6 @@ cdef class ParticleOctreeContainer(OctreeContainer):
     cdef Oct** oct_list
     #The starting oct index of each domain
     cdef np.int64_t *dom_offsets
-    cdef public int max_level
     #How many particles do we keep before refining
     cdef public int n_ref
 
@@ -1945,7 +1944,6 @@ cdef class ParticleBitmapSelector:
 
 cdef class ParticleBitmapOctreeContainer(SparseOctreeContainer):
     cdef Oct** oct_list
-    cdef public int max_level
     cdef public int n_ref
     cdef int loaded # Loaded with load_octree?
     cdef np.uint8_t* _ptr_index_base_roots