diff --git a/pystencils/datahandling/pycuda.py b/pystencils/datahandling/pycuda.py
new file mode 100644
index 0000000000000000000000000000000000000000..7309c54aee2a49e7c96dd4e1d044113edfbeb28a
--- /dev/null
+++ b/pystencils/datahandling/pycuda.py
@@ -0,0 +1,42 @@
+try:
+ import pycuda.gpuarray as gpuarray
+except ImportError:
+ gpuarray = None
+import numpy as np
+
+import pystencils
+
+
+class PyCudaArrayHandler:
+
+ def __init__(self):
+ import pycuda.autoinit # NOQA
+
+ def zeros(self, shape, dtype=np.float64, order='C'):
+ return gpuarray.zeros(shape, dtype, order)
+
+ def ones(self, shape, dtype, order='C'):
+ return gpuarray.ones(shape, dtype, order)
+
+ def empty(self, shape, dtype=np.float64, layout=None):
+ if layout:
+ cpu_array = pystencils.field.create_numpy_array_with_layout(shape, dtype, layout)
+ return self.from_numpy(cpu_array)
+ else:
+ return gpuarray.empty(shape, dtype)
+
+ def from_numpy(self, array, dtype=np.float64, layout=None):
+ return gpuarray.to_gpu(array)
+
+ def to_gpu(self, shape, dtype=np.float64, order='C'):
+ gpuarray.ones(shape, dtype, order)
+
+ def upload(gpuarray, numpy_array):
+ gpuarray.set(numpy_array)
+
+ def download(gpuarray, numpy_array):
+ gpuarray.get(numpy_array)
+
+ def randn(self, shape, dtype=np.float64):
+ cpu_array = np.random.randn(*shape).astype(dtype)
+ return self.from_numpy(cpu_array)
diff --git a/pystencils/datahandling/pyopencl.py b/pystencils/datahandling/pyopencl.py
new file mode 100644
index 0000000000000000000000000000000000000000..acc99ef298a949f4998fece481cee2450f923a74
--- /dev/null
+++ b/pystencils/datahandling/pyopencl.py
@@ -0,0 +1,40 @@
+try:
+ import pyopencl.array as gpuarray
+except ImportError:
+ gpuarray = None
+
+import numpy as np
+
+import pystencils
+
+
+class PyOpenClArrayHandler:
+
+ def __init__(self, queue):
+ self.queue = queue
+
+ def zeros(self, shape, dtype=np.float64, order='C'):
+ return gpuarray.zeros(shape, dtype, order)
+
+ def ones(self, shape, dtype, order='C'):
+ return gpuarray.ones(self.queue, shape, dtype, order)
+
+ def empty(self, shape, dtype=np.float64, layout=None):
+ if layout:
+ cpu_array = pystencils.field.create_numpy_array_with_layout(shape, dtype, layout)
+ return self.from_numpy(cpu_array)
+ else:
+ return gpuarray.empty(shape, dtype)
+
+ def from_numpy(self, array):
+ return gpuarray.to_device(self.queue, array)
+
+ def upload(self, gpuarray, numpy_array):
+ gpuarray.set(self.queue, numpy_array)
+
+ def download(self, gpuarray, numpy_array):
+ gpuarray.get(numpy_array)
+
+ def randn(self, shape, dtype=np.float64):
+ cpu_array = np.random.randn(*shape).astype(dtype)
+ return self.from_numpy(cpu_array)
diff --git a/pystencils/datahandling/serial_datahandling.py b/pystencils/datahandling/serial_datahandling.py
index 8c79a0931f06fad7c7ab58143cf7058a3f52dbac..8287a36f210c7d313f1f0def5f4ca9b76c9998dc 100644
--- a/pystencils/datahandling/serial_datahandling.py
+++ b/pystencils/datahandling/serial_datahandling.py
@@ -6,22 +6,24 @@ import numpy as np
from pystencils.datahandling.blockiteration import SerialBlock
from pystencils.datahandling.datahandling_interface import DataHandling
+from pystencils.datahandling.pycuda import PyCudaArrayHandler
+from pystencils.datahandling.pyopencl import PyOpenClArrayHandler
from pystencils.field import (
Field, create_numpy_array_with_layout, layout_string_to_tuple, spatial_layout_string_to_tuple)
from pystencils.slicing import normalize_slice, remove_ghost_layers
from pystencils.utils import DotDict
-try:
- import pycuda.gpuarray as gpuarray
- import pycuda.autoinit # NOQA
-except ImportError:
- gpuarray = None
-
class SerialDataHandling(DataHandling):
- def __init__(self, domain_size: Sequence[int], default_ghost_layers: int = 1, default_layout: str = 'SoA',
- periodicity: Union[bool, Sequence[bool]] = False, default_target: str = 'cpu') -> None:
+ def __init__(self,
+ domain_size: Sequence[int],
+ default_ghost_layers: int = 1,
+ default_layout: str = 'SoA',
+ periodicity: Union[bool, Sequence[bool]] = False,
+ default_target: str = 'cpu',
+ opencl_queue=None,
+ array_handler=None) -> None:
"""
Creates a data handling for single node simulations.
@@ -43,6 +45,19 @@ class SerialDataHandling(DataHandling):
self.custom_data_gpu = DotDict()
self._custom_data_transfer_functions = {}
+ if array_handler:
+ self.array_handler = array_handler
+ else:
+ try:
+ self.array_handler = PyCudaArrayHandler()
+ except Exception:
+ self.array_handler = None
+
+ if opencl_queue:
+ default_target = 'gpu'
+ assert opencl_queue
+ self.array_handler = PyOpenClArrayHandler(opencl_queue)
+
if periodicity is None or periodicity is False:
periodicity = [False] * self.dim
if periodicity is True:
@@ -125,7 +140,7 @@ class SerialDataHandling(DataHandling):
if gpu:
if name in self.gpu_arrays:
raise ValueError("GPU Field with this name already exists")
- self.gpu_arrays[name] = gpuarray.to_gpu(cpu_arr)
+ self.gpu_arrays[name] = self.array_handler.to_gpu(cpu_arr)
assert all(f.name != name for f in self.fields.values()), "Symbolic field with this name already exists"
self.fields[name] = Field.create_from_numpy_array(name, cpu_arr, index_dimensions=index_dimensions)
@@ -234,14 +249,14 @@ class SerialDataHandling(DataHandling):
transfer_func = self._custom_data_transfer_functions[name][1]
transfer_func(self.custom_data_gpu[name], self.custom_data_cpu[name])
else:
- self.gpu_arrays[name].get(self.cpu_arrays[name])
+ self.array_handler.download(self.gpu_arrays[name], self.cpu_arrays[name])
def to_gpu(self, name):
if name in self._custom_data_transfer_functions:
transfer_func = self._custom_data_transfer_functions[name][0]
transfer_func(self.custom_data_gpu[name], self.custom_data_cpu[name])
else:
- self.gpu_arrays[name].set(self.cpu_arrays[name])
+ self.array_handler.upload(self.gpu_arrays[name], self.cpu_arrays[name])
def is_on_gpu(self, name):
return name in self.gpu_arrays
diff --git a/pystencils/gpucuda/kernelcreation.py b/pystencils/gpucuda/kernelcreation.py
index 4430abba6e8aaacfd4d2b00286af46d9eb879741..eecfe278e4378cdbed8738275d2cfe813d350d7b 100644
--- a/pystencils/gpucuda/kernelcreation.py
+++ b/pystencils/gpucuda/kernelcreation.py
@@ -20,7 +20,7 @@ def create_cuda_kernel(assignments,
all_fields = fields_read.union(fields_written)
read_only_fields = set([f.name for f in fields_read - fields_written])
- buffers = set([f for f in all_fields if FieldType.is_buffer(f)])
+ buffers = set([f for f in all_fields if FieldType.is_buffer(f) or FieldType.is_custom(f)])
fields_without_buffers = all_fields - buffers
field_accesses = set()