Fix NumPy copy functions when using GPU backend

pytest/test_numpy_transform.py

@@ -0,0 +1,104 @@
+"""
+Direct unit tests for numpy_to_omega_h_read and omega_h_write_to_numpy.
+
+Verifies that data is properly deep-copied through host memory
+(not aliased from numpy buffer or device pointer).
+
+Run with: pytest test_numpy_transform.py -v
+"""
+
+import numpy as np
+import PyOmega_h as omega_h
+
+
+def test_numpy_to_read_copies_data():
+    """numpy_to_omega_h_read must copy data through host, not alias."""
+    size = 10
+    original = np.arange(size, dtype=np.float64)
+    read_view = omega_h.numpy_to_read_float64(original)
+
+    # Access Read data via HostRead (buffer protocol)
+    host_read = omega_h.HostRead_float64(read_view)
+    result = np.array(host_read, copy=False)
+    np.testing.assert_array_equal(result, original)
+
+    # Mutate original — result must be unchanged (deep copy)
+    original[:] = -999.0
+    result_after = np.array(host_read, copy=False)
+    expected = np.arange(size, dtype=np.float64)
+    np.testing.assert_array_equal(result_after, expected)
+
+
+def test_write_to_numpy_copies_data():
+    """omega_h_write_to_numpy must copy device data to host numpy array."""
+    size = 10
+    expected = np.arange(size, dtype=np.float64)
+
+    host_write = omega_h.HostWrite_float64(size)
+    # Fill via buffer protocol
+    view = np.array(host_write, copy=False)
+    view[:] = expected
+
+    # HostWrite.write() returns Write; convert to numpy
+    write_view = host_write.write()
+    result = omega_h.write_to_numpy_float64(write_view)
+    np.testing.assert_array_equal(result, expected)
+
+    # Mutate host — result must be unchanged (deep copy)
+    view[:] = -999.0
+    np.testing.assert_array_equal(result, expected)
+
+
+def test_numpy_to_read_sliced_step():
+    """numpy_to_omega_h_read must handle non-contiguous sliced arrays."""
+    base = np.arange(20, dtype=np.float64)
+    sliced = base[::3]  # [0, 3, 6, 9, 12, 15, 18] — non-contiguous
+    assert not sliced.flags["C_CONTIGUOUS"]
+
+    read_view = omega_h.numpy_to_read_float64(sliced)
+    host_read = omega_h.HostRead_float64(read_view)
+    result = np.array(host_read, copy=False)
+    np.testing.assert_array_equal(result, sliced)
+
+    # Verify deep copy — mutating base doesn't affect result
+    base[:] = -999.0
+    result_after = np.array(host_read, copy=False)
+    expected = np.array([0, 3, 6, 9, 12, 15, 18], dtype=np.float64)
+    np.testing.assert_array_equal(result_after, expected)
+
+
+def test_numpy_to_read_sliced_reversed():
+    """numpy_to_omega_h_read must handle reversed (negative stride) arrays."""
+    base = np.array([1, 2, 3, 4, 5], dtype=np.int32)
+    sliced = base[::-1]  # [5, 4, 3, 2, 1] — negative stride
+    assert not sliced.flags["C_CONTIGUOUS"]
+
+    read_view = omega_h.numpy_to_read_int32(sliced)
+    host_read = omega_h.HostRead_int32(read_view)
+    result = np.array(host_read, copy=False)
+    np.testing.assert_array_equal(result, sliced)
+
+    # Verify deep copy
+    base[:] = 0
+    result_after = np.array(host_read, copy=False)
+    expected = np.array([5, 4, 3, 2, 1], dtype=np.int32)
+    np.testing.assert_array_equal(result_after, expected)
+
+
+def test_numpy_to_read_sliced_offset():
+    """numpy_to_omega_h_read must handle sliced arrays with offset + step."""
+    base = np.array([10, 11, 12, 13, 14, 15, 16], dtype=np.int64)
+    sliced = base[2:6:2]  # [12, 14] — offset by 2, step 2
+    assert not sliced.flags["C_CONTIGUOUS"]
+
+    read_view = omega_h.numpy_to_read_int64(sliced)
+    host_read = omega_h.HostRead_int64(read_view)
+    result = np.array(host_read, copy=False)
+    np.testing.assert_array_equal(result, sliced)
+
+    # Verify deep copy
+    base[2] = 999
+    result_after = np.array(host_read, copy=False)
+    expected = np.array([12, 14], dtype=np.int64)
+    np.testing.assert_array_equal(result_after, expected)
+

src/PyOmega_h_array.cpp

@@ -1,5 +1,6 @@
 #include <Omega_h_array.hpp>
 #include <PyOmega_h.hpp>
+#include <PyOmega_h_numpy_transform.hpp>
 
 namespace Omega_h {
 
@@ -48,6 +49,17 @@ static void pybind11_array_type(py::module& m,
       &deep_copy;
   m.def(deepcopy_name.c_str(), deep_copy_type, py::arg("a"),
       py::arg("name") = "");
+  // Expose numpy transform functions for direct testing
+  m.def(("numpy_to_read_" + py_scalar).c_str(),
+      [](py::array_t<Scalar> arr) { return numpy_to_omega_h_read<Scalar>(arr); },
+      py::arg("arr"),
+      "Convert 1D numpy array to Read array (copies through host)");
+  m.def(("write_to_numpy_" + py_scalar).c_str(),
+      [](Write<Scalar> w) {
+        return omega_h_write_to_numpy<Scalar>(w);
+      },
+      py::arg("write_view"),
+      "Convert Write array to numpy array (copies through host)");
 }
 
 void pybind11_array(py::module& m) {

src/PyOmega_h_numpy_transform.hpp

@@ -17,10 +17,17 @@ Omega_h::Read<T> numpy_to_omega_h_read(py::array_t<T> arr)
   if (buf.ndim != 1) {
     throw std::runtime_error("Number of dimensions must be 1");
   }
-  Kokkos::View<T*, Kokkos::DefaultExecutionSpace::memory_space,
-               Kokkos::MemoryTraits<Kokkos::Unmanaged>>
-    view(reinterpret_cast<T*>(buf.ptr), buf.shape[0]);
-  Omega_h::Write<T> write_view(view);
+  // Copy element-by-element through HostWrite instead of wrapping
+  // the numpy pointer directly in a Kokkos view. This ensures:
+  // 1. Proper memory ownership (Write owns its allocation even with kokkos CPU backend)
+  // 2. Correct handling of non-contiguous (sliced) numpy arrays via stride
+  auto write_view_host = Omega_h::HostWrite<T>(buf.shape[0]);
+  T* ptr = static_cast<T*>(buf.ptr);
+  ssize_t stride = buf.strides[0] / sizeof(T);
+  for (Omega_h::LO i = 0; i < buf.shape[0]; ++i) {
+    write_view_host[i] = ptr[i * stride];
+  }
+  auto write_view = Omega_h::Write<T>(write_view_host);
   Omega_h::Read<T> read_view(write_view);
   return read_view;
 }
@@ -47,26 +54,34 @@ Omega_h::Write<T> numpy_to_omega_h_write(py::array_t<T> arr)
   if (buf.ndim != 1) {
     throw std::runtime_error("Number of dimensions must be 1");
   }
-  // Get host mirror and copy data
+  // Copy element-by-element through HostWrite instead of wrapping
+  // the numpy pointer directly in a Kokkos view. This ensures:
+  // 1. Proper memory ownership (Write owns its allocation even with kokkos CPU backend)
+  // 2. Correct handling of non-contiguous (sliced) numpy arrays via stride
   auto write_view_host = Omega_h::HostWrite<T>(buf.shape[0]);
   T* ptr = static_cast<T*>(buf.ptr);
+  ssize_t stride = buf.strides[0] / sizeof(T);
   for (Omega_h::LO i = 0; i < buf.shape[0]; ++i) {
-    write_view_host[i] = ptr[i];
+    write_view_host[i] = ptr[i * stride];
   }
   auto write_view = Omega_h::Write<T>(write_view_host);
   return write_view;
 }
 
-// Helper to convert Omega_h::Write to numpy array (creates a reference)
+// Helper to convert Omega_h::Write to numpy array (creates a copy)
 template <typename T>
 py::array_t<T> omega_h_write_to_numpy(Omega_h::Write<T> write_view)
 {
-  return py::array_t<T>({static_cast<py::ssize_t>(write_view.size())}, // shape
-                        {sizeof(T)},         // strides
-                        write_view.data(),   // data pointer
-                        py::cast(write_view) // base object to manage lifetime
-  );
+  auto write_view_host = Omega_h::HostWrite<T>(write_view);
+  py::array_t<T> result(write_view_host.size());
+  py::buffer_info buf = result.request();
+  T* ptr = static_cast<T*>(buf.ptr);
+  for (Omega_h::LO i = 0; i < write_view_host.size(); ++i) {
+    ptr[i] = write_view_host[i];
+  }
+  return result;
 }
 } // namespace Omega_h
 
-#endif
+#endif
+