From fdc4d49e43e965507cdcb1752ff763d784222e91 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Sun, 25 Nov 2018 20:25:16 +0100
Subject: [PATCH] opencl fftw spectral operators

---
 hysop/__init__.py.in                          |  2 +-
 hysop/backend/host/host_operator.py           |  3 +++
 hysop/backend/host/python/operator/poisson.py |  4 ++--
 .../host/python/operator/poisson_curl.py      |  2 +-
 hysop/operator/base/derivative.py             |  2 +-
 hysop/operator/base/poisson_curl.py           |  2 +-
 hysop/operator/base/solenoidal_projection.py  |  2 +-
 hysop/operator/tests/test_poisson_curl.py     |  9 +++++++++
 .../tests/test_solenoidal_projection.py       |  9 +++++++++
 hysop/tools/numba_utils.py                    | 20 +++++++++++++++----
 10 files changed, 44 insertions(+), 11 deletions(-)

diff --git a/hysop/__init__.py.in b/hysop/__init__.py.in
index 0b1f41f2b..9e79793ef 100644
--- a/hysop/__init__.py.in
+++ b/hysop/__init__.py.in
@@ -55,7 +55,7 @@ __MAX_THREADS__ = int(get_env('MAX_THREADS', psutil.cpu_count(logical=True))
 set_env('OMP_NUM_THREADS',   __MAX_THREADS__)
 set_env('MKL_NUM_THREADS',   __MAX_THREADS__)
 set_env('NUMBA_NUM_THREADS', __MAX_THREADS__)
-set_env('NUMBA_THREADING_LAYER', 'omp')
+set_env('NUMBA_THREADING_LAYER', 'workqueue') # Use 'numba -s' to list support
 __DEFAULT_NUMBA_TARGET__ = ('parallel' if __ENABLE_THREADING__ else 'cpu')
 
 # OpenCL
diff --git a/hysop/backend/host/host_operator.py b/hysop/backend/host/host_operator.py
index 67ba01688..926e604b9 100644
--- a/hysop/backend/host/host_operator.py
+++ b/hysop/backend/host/host_operator.py
@@ -108,6 +108,9 @@ class OpenClMappable(object):
         self.__registered_getters = {}
         self.__mapped_objects     = {}
 
+    def __del__(self):
+        self.unmap_objects()
+
     @property
     def cl_env(self):
         return self.__cl_env
diff --git a/hysop/backend/host/python/operator/poisson.py b/hysop/backend/host/python/operator/poisson.py
index 51ee23512..2658e7a93 100644
--- a/hysop/backend/host/python/operator/poisson.py
+++ b/hysop/backend/host/python/operator/poisson.py
@@ -109,7 +109,7 @@ class PythonPoisson(PoissonOperatorBase, OpenClMappable, HostOperator):
                 Ft()
                 filter_poisson()
                 Bt()
-
+        
         for Fo in self.dFout.dfields:
-                Fo.exchange_ghosts()
+            Fo.exchange_ghosts()
 
diff --git a/hysop/backend/host/python/operator/poisson_curl.py b/hysop/backend/host/python/operator/poisson_curl.py
index 8685e74b4..a20e230bb 100644
--- a/hysop/backend/host/python/operator/poisson_curl.py
+++ b/hysop/backend/host/python/operator/poisson_curl.py
@@ -216,6 +216,6 @@ class PythonPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClMappable, HostOpe
                 Bt()
 
         self.dU.exchange_ghosts()
-        if (diffuse of project):
+        if (diffuse or project):
             self.dW.exchange_ghosts()
 
diff --git a/hysop/operator/base/derivative.py b/hysop/operator/base/derivative.py
index 6750d23c6..2f9541d16 100644
--- a/hysop/operator/base/derivative.py
+++ b/hysop/operator/base/derivative.py
@@ -179,7 +179,7 @@ class SpaceDerivativeBase(object):
         requests = super(SpaceDerivativeBase, self).get_work_properties()
         if self.require_tmp:
             request = MemoryRequest.empty_like(a=self.dFout, nb_components=1, 
-                    shape=self.dFout.compute_resolution)
+                    shape=self.dFout.compute_resolution, backend=self.backend)
             requests.push_mem_request('tmp', request)
         return requests
     
diff --git a/hysop/operator/base/poisson_curl.py b/hysop/operator/base/poisson_curl.py
index 896e4fb87..9555cce2d 100644
--- a/hysop/operator/base/poisson_curl.py
+++ b/hysop/operator/base/poisson_curl.py
@@ -268,7 +268,7 @@ class SpectralPoissonCurlOperatorBase(PoissonCurlOperatorBase, SpectralOperatorB
                 assert (Ft.output_shape == Bt.input_shape), (Ft.output_shape, Bt.input_shape)
             shape = Ft.output_shape
             dtype = Ft.output_dtype
-            request = MemoryRequest(backend=Ft.backend, dtype=dtype, 
+            request = MemoryRequest(backend=self.backend, dtype=dtype, 
                                     shape=shape, nb_components=1,
                                     alignment=self.min_fft_alignment)
             requests.push_mem_request('fft_buffer_{}'.format(i), request)
diff --git a/hysop/operator/base/solenoidal_projection.py b/hysop/operator/base/solenoidal_projection.py
index 93f66ec73..8e45b23f8 100644
--- a/hysop/operator/base/solenoidal_projection.py
+++ b/hysop/operator/base/solenoidal_projection.py
@@ -189,7 +189,7 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
             assert (Ft.output_shape == Bt.input_shape), (Ft.output_shape, Bt.input_shape)
             shape = Ft.output_shape
             dtype = Ft.output_dtype
-            request = MemoryRequest(backend=Ft.backend, dtype=dtype, 
+            request = MemoryRequest(backend=self.backend, dtype=dtype, 
                                     shape=shape, nb_components=1,
                                     alignment=self.min_fft_alignment)
             requests.push_mem_request('fft_buffer_{}'.format(i), request)
diff --git a/hysop/operator/tests/test_poisson_curl.py b/hysop/operator/tests/test_poisson_curl.py
index 0efce527f..50906379e 100644
--- a/hysop/operator/tests/test_poisson_curl.py
+++ b/hysop/operator/tests/test_poisson_curl.py
@@ -194,6 +194,15 @@ class TestPoissonCurlOperator(object):
                                                           cl_env.device.name.strip())
                     print msg,
                     yield PoissonCurl(cl_env=cl_env, projection=0, **base_kwds)
+                cpu_envs = tuple(iter_clenv(device_type='cpu'))
+                if cpu_envs:
+                    msg='   *OpenCl FFTW: '
+                    print msg
+                    for cl_env in cpu_envs:
+                        msg='     |platform {}, device {}'.format(cl_env.platform.name.strip(),
+                                                                  cl_env.device.name.strip())
+                        print msg,
+                        yield PoissonCurl(cl_env=cl_env, enforce_implementation=False, **base_kwds)
             else:
                 msg='Unknown implementation to test {}.'.format(impl)
                 raise NotImplementedError(msg)
diff --git a/hysop/operator/tests/test_solenoidal_projection.py b/hysop/operator/tests/test_solenoidal_projection.py
index 6b1990535..9e7f4fd24 100644
--- a/hysop/operator/tests/test_solenoidal_projection.py
+++ b/hysop/operator/tests/test_solenoidal_projection.py
@@ -243,6 +243,15 @@ class TestSolenoidalProjectionOperator(object):
                                                           cl_env.device.name.strip())
                     print msg,
                     yield SolenoidalProjection(cl_env=cl_env, **base_kwds)
+                cpu_envs = tuple(iter_clenv(device_type='cpu'))
+                if cpu_envs:
+                    msg='   *OpenCl FFTW: '
+                    print msg
+                    for cl_env in cpu_envs:
+                        msg='     |platform {}, device {}'.format(cl_env.platform.name.strip(),
+                                                                  cl_env.device.name.strip())
+                        print msg,
+                        yield SolenoidalProjection(cl_env=cl_env, enforce_implementation=False, **base_kwds)
             else:
                 msg='Unknown implementation to test {}.'.format(impl)
                 raise NotImplementedError(msg)
diff --git a/hysop/tools/numba_utils.py b/hysop/tools/numba_utils.py
index ff2e524d2..527f7270d 100644
--- a/hysop/tools/numba_utils.py
+++ b/hysop/tools/numba_utils.py
@@ -3,7 +3,11 @@ import numba as nb
 import numpy as np
 from hysop.core.arrays.array import Array
 
-def make_numba_signature(*args):
+def make_numba_signature(*args, **kwds):
+    raise_on_cl_array = kwds.pop('raise_on_cl_array', True)
+    if kwds:
+        msg='Unknown kwds {}.'.forma(kwds.keys())
+        raise RuntimeError(kwds)
     dtype_to_ntype = {
             int:        np.int32,
             long:       np.int64,
@@ -44,13 +48,21 @@ def make_numba_signature(*args):
     
     numba_args = ()
     numba_layout = ()
-    for a in args:
+    for i,a in enumerate(args):
+        from hysop.backend.device.opencl import clArray
         if isinstance(a, Array):
+            a = a.handle
+        if isinstance(a, clArray.Array):
+            # some opencl arrays can be mapped to host
+            if raise_on_cl_array:
+                msg='Numba signature: Got a cl.Array or hysop.OpenClArray for argument {} (shape={}, dtype={}).'
+                msg=msg.format(i, a.shape, a.dtype)
+                raise ValueError(msg)
             assert a.dtype.type in dtype_to_ntype, a.dtype.type
             dtype = dtype_to_ntype[a.dtype.type]
-            if a.is_c_contiguous:
+            if a.flags.c_contiguous:
                 na = nb.types.Array(dtype=dtype, ndim=a.ndim, layout='C')
-            elif a.is_fortran_contiguous:
+            elif a.flags.f_contiguous:
                 na = nb.types.Array(dtype=dtype, ndim=a.ndim, layout='F')
             else:
                 na = nb.types.Array(dtype=dtype, ndim=a.ndim, layout='A')
-- 
GitLab