WIP: Try hotfixing

Author: cphyc

yt/data_objects/index_subobjects/grid_patch.py

@@ -373,16 +373,25 @@ def deposit(
         # octree deposit routines, so we reverse it here to match the
         # convention there
 
-        # append a dummy dimension because we are only depositing onto
-        # one grid
-        nvals = tuple(self.ActiveDimensions[::-1]) + (1,)
-        if vector_field:
+        nvals = tuple(self.ActiveDimensions[::-1])
+        scalar_as_vector = (
+            not vector_field
+            and method != "count"
+            and fields is not None
+            and len(fields) > 0
+            and getattr(fields[0], "ndim", 0) == 1
+        )
+        if scalar_as_vector:
+            fields = [np.ascontiguousarray(f.reshape(f.shape[0], 1)) for f in fields]
+
+        use_vector_path = vector_field or scalar_as_vector
+        # Append a dummy grid dimension because we only deposit onto one grid.
+        # For vector fields, append an additional vector-component dimension.
+        if use_vector_path:
             vec_size = fields[0].shape[-1] if fields else 1
-            nvals = nvals + (vec_size,)
+            op = cls(nvals + (1, vec_size), kernel_name)
         else:
-            vec_size = 1
-
-        op = cls(nvals, kernel_name)
+            op = cls(nvals + (1,), kernel_name)
 
         op.initialize()
         if positions.size > 0:
@@ -392,12 +401,16 @@ def deposit(
             return
         # Fortran-ordered, so transpose.
         vals = vals.transpose()
-        # squeeze dummy dimension(s) we appended above
-        if vector_field:
-            new_shape = (vec_size, *self.ActiveDimensions)
+        if use_vector_path:
+            # (vec, 1, Nx, Ny, Nz) -> (Nx, Ny, Nz, vec)
+            vals = np.squeeze(vals, axis=1)
+            vals = np.moveaxis(vals, 0, -1)
+            if scalar_as_vector:
+                vals = np.squeeze(vals, axis=-1)
         else:
-            new_shape = self.ActiveDimensions
-        return vals.reshape(new_shape)
+            # (1, Nx, Ny, Nz) -> (Nx, Ny, Nz)
+            vals = np.squeeze(vals, axis=0)
+        return vals
 
     def select_blocks(self, selector):
         mask = self._get_selector_mask(selector)

yt/data_objects/static_output.py

@@ -1884,7 +1884,6 @@ def add_deposited_particle_field(
             ptype, deposit_field = deposit_field[0], deposit_field[1]
         else:
             raise RuntimeError
-
         if weight_field is None:
             weight_field = (ptype, "particle_mass")
         units = self.field_info[ptype, deposit_field].output_units
@@ -1923,13 +1922,16 @@ def _deposit_field(data):
             fields = [data[ptype, deposit_field]]
             if method == "weighted_mean":
                 fields.append(data[ptype, weight_field])
-            fields = [np.ascontiguousarray(f) for f in fields]
+            # Cython deposition kernels operate on float64 buffers.
+            fields = [np.ascontiguousarray(f, dtype="float64") for f in fields]
             d = data.deposit(
                 pos,
                 fields,
                 method=method,
                 kernel_name=kernel_name,
-                vector_field=vector_field,
+                # Count deposition tracks number of particles per cell and
+                # should remain scalar even if the source field is vector-valued.
+                vector_field=vector_field and method != "count",
             )
 
             d = data.ds.arr(d, units=units)

yt/geometry/particle_deposit.pyx

@@ -26,8 +26,7 @@ from yt.utilities.lib.misc_utilities import OnceIndirect
 cdef append_axes(np.ndarray arr, int naxes):
     if arr.ndim == naxes:
         return arr
-    for _ in range((naxes - arr.ndim)):
-        arr = np.expand_dims(arr)
+    arr = np.expand_dims(arr, axis=tuple(range(arr.ndim, naxes)))
     return arr
 
 cdef class ParticleDepositOperation:
@@ -58,7 +57,7 @@ cdef class ParticleDepositOperation:
         nf = len(fields)
 
         # Convert scalar fields into vector fields of length 1
-        fields = [append_axes(field, 2) for field in fields]
+        fields = [append_axes(field, 1) for field in fields]
 
         # Make sure all of them have the same length
         nvec = 1
@@ -140,7 +139,7 @@ cdef class ParticleDepositOperation:
         nf = len(fields)
 
         # Convert scalar fields into vector fields of length 1
-        fields = [append_axes(field, 2) for field in fields]
+        fields = [append_axes(field, 1) for field in fields]
 
         # Make sure all of them have the same length
         nvec = 1
@@ -382,11 +381,11 @@ cdef class CICDeposit(ParticleDepositOperation):
     cdef np.float64_t[:,:,:,:,:] field
     cdef public object ofield
     def initialize(self):
-        if not all(_ > 1 for _ in self.nvals[:-1]):
+        if not all(_ > 1 for _ in self.nvals[:3]):
             from yt.utilities.exceptions import YTBoundsDefinitionError
             raise YTBoundsDefinitionError(
                 "CIC requires minimum of 2 zones in all spatial dimensions.",
-                self.nvals[:-1])
+                self.nvals[:3])
         self.field = append_axes(
             np.zeros(self.nvals, dtype="float64", order='F'), 5)