opencl poisson rotational

hysop/__init__.py

@@ -44,7 +44,7 @@ __TEST_ALL_OPENCL_PLATFORMS__ = get_env('TEST_ALL_OPENCL_PLATFORMS', False)
 __ENABLE_LONG_TESTS__ = get_env('ENABLE_LONG_TESTS', False)
 
 # OpenCL
-__DEFAULT_PLATFORM_ID__ = 2
+__DEFAULT_PLATFORM_ID__ = 1
 __DEFAULT_DEVICE_ID__   = 0
 
 if __MPI_ENABLED__:

hysop/backend/device/opencl/operator/poisson_rotational.py

+
+import primefac
+from hysop.tools.numpywrappers import npw
+from hysop.tools.decorators  import debug
+from hysop.tools.warning import HysopWarning
+from hysop.operator.base.poisson_rotational import PoissonRotationalOperatorBase
+from hysop.backend.device.opencl.opencl_symbolic import OpenClSymbolic
+from hysop.core.graph.graph import op_apply
+from hysop.core.memory.memory_request import MemoryRequest
+from hysop.backend.device.opencl.opencl_fft import OpenClFFT
+from hysop.backend.device.codegen.base.variables import dtype_to_ctype
+
+class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
+    '''
+    Solves the poisson-rotational equation using clFFT.
+    '''
+    
+    def generate_wave_numbers(self, dim, resolution, length, dtype, ctype, axes):
+        if (dim>3):
+            msg='clFFT only support 1D, 2D or 3D plans, got a {}D domain.'
+            msg=msg.format(dim)
+            raise ValueError(msg)
+
+        valid_factors = {2,3,5,7,11,13}
+        for Ni in resolution:
+            factors = tuple( primefac.primefac(int(Ni)) )
+            invalid_factors = set(factors) - valid_factors
+            if invalid_factors:
+                factorization = ' + '.join('{}^{}'.format(factor, factors.count(factor)) 
+                                                                for factor in set(factors))
+                candidates = ', '.join(str(vf) for vf in valid_factors)
+                msg ='\nInvalid transform shape {} for clFFT:'
+                msg+='\n   {} = {}'
+                msg+='\nOnly {} prime factors are available.'
+                msg+='\n'
+                msg=msg.format(resolution, Ni, factorization, candidates)
+                raise ValueError(msg)
+
+        K1 = ()
+        K2 = ()
+        for i in axes:
+            Ni = resolution[i]
+            Li = length[i]
+
+            if (i==dim-1):
+                k1 = 2*npw.pi*1j*npw.fft.rfftfreq(Ni,Li)*Ni
+            else:
+                k1 = 2*npw.pi*1j*npw.fft.fftfreq(Ni, Li)*Ni
+            k1 = k1.astype(dtype=ctype).copy()
+            k2 = (k1**2).real.astype(dtype=dtype).copy()
+            K1 += (k1,)
+            K2 += (k2,)
+        
+        shape = (sum(k.size for k in K1),)
+        K1d = self.backend.empty(shape=shape, dtype=ctype)
+        K2d = self.backend.empty(shape=shape, dtype=dtype)
+
+        Ksizes, Koffsets = (), ()
+        start, end = 0,0
+        for (k1,k2) in zip(K1,K2):
+            assert k1.size == k2.size
+            Koffsets += (start,)
+            Ksizes   += (k1.size,)
+            end      += k1.size
+            K1d[start:end] = k1
+            K2d[start:end] = k2
+            start = end
+        
+        self.K1d      = K1d
+        self.K2d      = K2d
+        self.Ksizes   = Ksizes
+        self.Koffsets = Koffsets
+    
+    @debug
+    def get_work_properties(self):
+        requests  = super(OpenClPoissonRotational,self).get_work_properties()
+        
+        axes  = self.axes
+        dW    = self.dW
+        ctype = self.ctype
+
+        shape = list(self.resolution)
+        shape[axes[0]] = shape[axes[0]] // 2 + 1
+        shape = tuple(shape)
+        assert npw.array_equal(shape, self.Ksizes[::-1])
+
+        request = MemoryRequest.empty_like(a=dW, dtype=ctype, 
+                shape = (dW.nb_components,)+shape,  # contiguous components
+                nb_components=1)
+        requests.push_mem_request('R2C_C2R', request)
+        
+        return requests
+
+    @debug
+    def setup(self, work):
+        super(OpenClPoissonRotational, self).setup(work)
+        self._build_fft_plans(work)
+
+    def _build_fft_plans(self, work):
+        axes    = self.axes
+        context = self.backend.cl_env.context
+        queue   = self.backend.cl_env.default_queue
+        
+        fft_buffer = work.get_buffer(self, 'R2C_C2R')[0]
+        forward_plans, backward_plans = [],[]
+        
+        for (i,Wi) in enumerate(self.W_buffers):
+            fp = OpenClFFT(context=context, queue=queue, 
+                        in_array=Wi, out_array=fft_buffer[0].handle,
+                        axes=axes, fast_math=False,
+                        real=True)
+            (forward_callbacks, forward_user_data) = self._build_forward_callbacks(i)
+            fp.plan.set_callback('post_callback', forward_callbacks['post'],
+                                    'post', user_data=forward_user_data)
+            fp.bake().allocate()
+            forward_plans.append(fp)
+
+        for (i,Ui) in enumerate(self.U_buffers):
+            (backward_callbacks, backward_user_data) = self._build_backward_callbacks(i)
+            bp = OpenClFFT(context=context, queue=queue, 
+                        in_array=fft_buffer[0].handle, out_array=Ui,
+                        axes=axes, fast_math=False,
+                        real=True)
+            bp.plan.set_callback('pre_callback', backward_callbacks['pre'],
+                                    'pre', user_data=backward_user_data)
+            bp.bake().allocate()
+            backward_plans.append(bp)
+
+        self.forward_plans  = forward_plans
+        self.backward_plans = backward_plans
+        self.fft_buffer     = fft_buffer
+
+    def _build_forward_callbacks(self, i):
+        fp = dtype_to_ctype(self.dtype)
+        (offsets, sizes) = (self.Koffsets, self.Ksizes)
+        N = npw.prod(sizes, dtype=npw.int64)
+        gencode = ''.join(
+                  '''
+                                i = (off % {Si});
+                                off /= {Si};
+                                C += K2[{Oi}+i];
+                  '''.format(Si=Si, Oi=Oi) for (Oi,Si) in zip(offsets, sizes))
+        forward_callbacks = {
+            'post': 
+                '''
+                void post_callback(__global void* output, 
+                                  const uint offset,                      
+                                  __global void* userdata,
+                                  {fp}2 res)
+                    {{                                                            
+                        __global {fp}2* out = (__global {fp}2*) output; 
+                        __global {fp}*  K2  = (__global {fp}*)  userdata;
+                        if (offset==0) {{
+                            res = ({fp}2)(0);
+                        }}
+                        else {{
+                            {fp} C = ({fp})(0);
+                            {{
+                                uint i;
+                                uint off = offset;
+                                {gencode}
+                            }}
+                            res /= C;
+                        }}
+                        out[{i}*{N}+offset] = res;
+                    }}
+                '''.format(fp=fp, i=i, N=N,
+                            gencode=gencode)
+        }
+        user_data = self.K2d.data
+        return (forward_callbacks, user_data)
+
+    def _build_backward_callbacks(self, i):
+        dim = self.dim
+        fp = dtype_to_ctype(self.dtype)
+        
+        if (dim == 2):
+            (Nx,Ny) = self.Ksizes
+            (Ox,Oy) = self.Koffsets
+            gencode = {
+                    0:
+                        '''
+                                uint iy = (offset / {Nx});
+                                res = cmul(-K1[{Oy}+iy], in[offset]);
+                        '''.format(Nx=Nx, Oy=Oy),
+                    1:
+                        '''
+                                uint ix = (offset % {Nx});
+                                res = cmul(+K1[{Ox}+ix], in[offset]);
+                        '''.format(Nx=Nx, Ox=Ox),
+            }
+        elif (dim==3):
+            (Nx,Ny,Nz) = self.Ksizes
+            (Ox,Oy,Oz) = self.Koffsets
+            N = Nx*Ny*Nz
+            gencode = {
+                    0:
+                        '''
+                                uint iy = ((offset / {Nx}) % {Ny});
+                                uint iz = (offset / ({Nx}*{Ny}));
+                                res += cmul(-K1[{Oy}+iy], in[2*{N}+offset]);
+                                res += cmul(+K1[{Oz}+iz], in[1*{N}+offset]);
+                        '''.format(Nx=Nx, Ny=Ny, N=N, Oy=Oy, Oz=Oz),
+                    1:
+                        '''
+                                uint ix = (offset % {Nx});
+                                uint iz = (offset / ({Nx}*{Ny}));
+                                res += cmul(-K1[{Oz}+iz], in[0*{N}+offset]);
+                                res += cmul(+K1[{Ox}+ix], in[2*{N}+offset]);
+                        '''.format(Nx=Nx, Ny=Ny, N=N, Ox=Ox, Oz=Oz),
+                    2:
+                        '''
+                                uint ix = (offset % {Nx});
+                                uint iy = ((offset / {Nx}) % {Ny});
+                                res += cmul(-K1[{Ox}+ix], in[1*{N}+offset]);
+                                res += cmul(+K1[{Oy}+iy], in[0*{N}+offset]);
+                        '''.format(Nx=Nx, Ny=Ny, N=N, Ox=Ox, Oy=Oy),
+            }
+        else:
+            msg='Unsupported dimension {}.'.format(dim)
+            raise RuntimeError(msg)
+
+        backward_callbacks = {
+            'pre': 
+                '''
+                {fp}2 cmul(const {fp}2 a, const {fp}2 b) {{
+                    return ({fp}2)(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);
+                }}
+
+                {fp}2 pre_callback(const __global void* input,
+                                   const uint offset,
+                                   const __global void* userdata)
+                    {{
+                        const __global {fp}2* in = (const __global {fp}2*) input;
+                        const __global {fp}2* K1 = (const __global {fp}2*) userdata;
+                        {fp}2 res = ({fp}2)(0);
+                        {{{gencode}}}
+                        return res;
+                    }}
+                '''.format(fp=fp, gencode=gencode[i])
+        }
+        print
+        print '{}D: axe {}'.format(dim,i)
+        print backward_callbacks['pre']
+        user_data = self.K1d.data
+        return (backward_callbacks, user_data)
+
+
+    @op_apply
+    def apply(self, **kwds):
+        '''Solve the PoissonRotational equation.'''
+        for fp in self.forward_plans:
+            evt, = fp.enqueue()
+        for bp in self.backward_plans:
+            evt, = bp.enqueue()
+        evt = self.dW.exchange_ghosts()
+        if evt:
+            evt.wait()

hysop/backend/host/python/operator/poisson_rotational.py

@@ -4,160 +4,84 @@ from hysop.tools.decorators import debug
 from hysop.tools.numpywrappers import npw
 from hysop.backend.host.host_operator import HostOperator
 from hysop.core.graph.graph import op_apply
-from hysop.fields.continuous_field import Field
-from hysop.parameters.parameter import Parameter
-from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
-from hysop.constants import FieldProjection
+from hysop.operator.base.poisson_rotational import PoissonRotationalOperatorBase
 
-class PythonPoissonRotational(HostOperator):
+class PythonPoissonRotational(PoissonRotationalOperatorBase, HostOperator):
     """
     Solves the poisson rotational equation using numpy fftw.
     """
-
-    @debug
-    def __init__(self, vorticity, velocity, variables, 
-            projection=None, **kwds): 
-        """PoissonRotational operator to solve incompressible flows using FFTW in Fortran.
+   
+    def generate_wave_numbers(self, dim, resolution, length, dtype, ctype, axes):
+        Z  = (0,)*dim
+        K = npw.ix_(*tuple(2*npw.pi*1j*npw.fft.fftfreq(Ni,Li)*Ni for (Ni,Li) in zip(resolution,length)))[::-1]
         
-        Parameters
-        ----------
-        velocity : :class:`~hysop.fields.continuous_field.Field
-            solution field
-        vorticity:  :class:`~hysop.fields.continuous_field.Field`
-            right-hand side
-        variables: dict
-            dictionary of fields as keys and topologies as values.
-        projection: hysop.constants.FieldProjection or positive integer, optional
-            Project vorticity such that resolved velocity is divergence free (for 3D fields).
-            When active, projection is done prior to every solve, unless projection is 
-            an integer in which case it is done every n applies.
-            This parameter is ignored for 2D fields and defaults to no projection.
-        kwds : 
-            base class parameters.
-        """
+        k2 = (ki*ki for ki in K)
+        K2 = sum(k2).real.copy()
         
-        check_instance(velocity,  Field)
-        check_instance(vorticity, Field)
-        check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
-
-        assert velocity.domain is vorticity.domain, 'only one domain is supported'
-        assert variables[velocity] is variables[vorticity], 'only one topology is supported'
+        iK2 = npw.empty_like(K2)
+        iK2[K2!=0] = 1.0/K2[K2!=0]
+        iK2[Z]     = 0.0
+        
+        self.K = K
+        self.iK2 = iK2
 
-        vtopology = variables[velocity]
-        wtopology = variables[vorticity]
-        input_fields  = { vorticity: wtopology }
-        output_fields = { velocity:  vtopology }
+    @op_apply
+    def apply(self, simulation=None, **kwds):
+        """Apply analytic formula."""
+        super(PythonPoissonRotational, self).apply(**kwds)
         
-        if (projection == FieldProjection.NONE) or (projection is None) or (velocity.dim==2):
-            projection = FieldProjection.NONE
-            output_fields[vorticity] = wtopology
+        if (self.dim == 2):
+            self._solve2d(**kwds)
+        elif (self.dim == 3):
+            project=self._do_project(simu=simulation)
+            self._solve3d(project=project, **kwds)
         else:
-            assert (velocity.dim==3) and (vorticity.dim==3), 'Projection only available in 3D.'
-            output_fields[vorticity] = wtopology
-        
-        super(PythonPoissonRotational, self).__init__(input_fields=input_fields, 
-                output_fields=output_fields, **kwds)
-        
-        dim=velocity.dim
-        if (dim!=3) and (projection!=FieldProjection.NONE):
-            raise ValueError('Velocity reprojection only available in 3D.')
-        
-        if (projection == FieldProjection.NONE):
-            self._do_project = lambda simu: False
-        elif (projection == FieldProjection.EVERY_STEP):
-            self._do_project = lambda simu: True
-        else: # projection is an integer representing frenquency 
-            freq = projection
-            assert freq>=1
-            self._do_project = lambda simu: ((simu.current_iteration % freq)==0)
-        
-        self.W = vorticity
-        self.U = velocity
-        self.projection = projection
-        self.dim = vorticity.dim
-    
-        
-    @debug
-    def discretize(self):
-        if self.discretized:
-            return
-        super(PythonPoissonRotational, self).discretize()
-        dW = self.input_discrete_fields[self.W]
-        dU = self.output_discrete_fields[self.U]
-        resolution = dW.compute_resolution
-        length     = dW.domain.length
-        self.W_buffers = tuple(data[dW.compute_slices] for data in dW.buffers)
-        self.U_buffers = tuple(data[dU.compute_slices] for data in dU.buffers)
-        self.N = npw.prod(resolution)
-        k = npw.ix_(*tuple(2*npw.pi*1j*npw.fft.fftfreq(Ni,Li)*Ni for (Ni,Li) in zip(resolution,length)))[::-1]
-        k2 = (ki*ki for ki in k)
-        self.K2 = sum(k2).real
-        self.Z = (0,)*dW.dim
-        self.k = k
-        self.dW = dW
-        self.dU = dU
-
+            msg='Unsupported dimension {}.'.format(dim)
+            raise ValueError(msg)
 
     def _solve2d(self, **kwds):
         W_buffers, U_buffers = self.W_buffers, self.U_buffers
-        N, K2, Z, k = self.N, self.K2, self.Z, self.k
+        iK2, K = self.iK2, self.K
 
         W0 = W_buffers[0]
         U0 = U_buffers[0]
         U1 = U_buffers[1]
         
         W0_hat = npw.fft.fftn(W0)
-        W0_hat[...] /= K2
-        W0_hat[Z] = 0
+        W0_hat[...] *= iK2
         Ui_hat = npw.empty_like(W0_hat)
         Ui_hat[...] = W0_hat
-        Ui_hat[...] *= -k[1]
+        Ui_hat[...] *= -K[1]
         U0[...] = npw.fft.ifftn(Ui_hat).real
         Ui_hat[...] = W0_hat
-        Ui_hat[...] *= k[0]
+        Ui_hat[...] *= K[0]
         U1[...] = npw.fft.ifftn(Ui_hat).real
         self.dU.exchange_ghosts()
 
     def _solve3d(self, project=False, **kwds):
         W_buffers, U_buffers = self.W_buffers, self.U_buffers
-        N, K2, Z, k = self.N, self.K2, self.Z, self.k
+        iK2, K = self.iK2, self.K
 
         W0,W1,W2 = W_buffers
         U0,U1,U2 = U_buffers
         
-        W0_hat = npw.fft.fftn(W0)
-        W1_hat = npw.fft.fftn(W1)
-        W2_hat = npw.fft.fftn(W2)        
+        W0_hat = npw.fft.fftn(W0) * iK2
+        W1_hat = npw.fft.fftn(W1) * iK2
+        W2_hat = npw.fft.fftn(W2) * iK2   
 
         if project:
-            divW = k[0]*W0_hat + k[1]*W1_hat + k[2]*W2_hat
-            W0_hat -= k[0]*divW / K2
-            W1_hat -= k[1]*divW / K2
-            W2_hat -= k[2]*divW / K2
+            divW = K[0]*W0_hat + K[1]*W1_hat + K[2]*W2_hat
+            W0_hat -= K[0]*divW
+            W1_hat -= K[1]*divW
+            W2_hat -= K[2]*divW
         
         Ui_hat = npw.empty_like(W0_hat)
-        Ui_hat[...] = (-k[1]*W2_hat + k[2]*W1_hat) / K2
-        Ui_hat[Z] = 0
+        Ui_hat[...] = (-K[1]*W2_hat + K[2]*W1_hat)
         U0[...] = npw.fft.ifftn(Ui_hat).real
-        Ui_hat[...] = (-k[2]*W0_hat + k[0]*W2_hat) / K2
-        Ui_hat[Z] = 0
+        Ui_hat[...] = (-K[2]*W0_hat + K[0]*W2_hat)
         U1[...] = npw.fft.ifftn(Ui_hat).real
-        Ui_hat[...] = (-k[0]*W1_hat + k[1]*W0_hat) / K2
-        Ui_hat[Z] = 0
+        Ui_hat[...] = (-K[0]*W1_hat + K[1]*W0_hat)
         U2[...] = npw.fft.ifftn(Ui_hat).real
+
         self.dU.exchange_ghosts()
-         
 
-    @op_apply
-    def apply(self, simulation=None, **kwds):
-        """Apply analytic formula."""
-        super(PythonPoissonRotational, self).apply(**kwds)
-        
-        if (self.dim == 2):
-            self._solve2d(**kwds)
-        elif (self.dim == 3):
-            project=self._do_project(simu=simulation)
-            self._solve3d(project=project, **kwds)
-        else:
-            msg='Unsupported dimension {}.'.format(dim)
-            raise ValueError(msg)

hysop/operator/base/poisson_rotational.py

+
+from abc import abstractmethod
+from hysop.tools.numpywrappers import npw
+from hysop.tools.types import check_instance, first_not_None
+from hysop.tools.decorators import debug
+from hysop.tools.numerics import float_to_complex_dtype
+from hysop.fields.continuous_field import Field
+from hysop.operator.base.fft_operator import FftOperatorBase
+from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
+from hysop.constants import FieldProjection
+from hysop.fields.continuous_field import Field
+
+class PoissonRotationalOperatorBase(FftOperatorBase):
+    """
+    Solves the poisson-rotational equation using a specific implementation.
+    """
+    
+    @debug
+    def __init__(self, vorticity, velocity, variables, 
+            projection=None, **kwds): 
+        """PoissonRotational operator to solve incompressible flows using FFTW in Fortran.
+        
+        Parameters
+        ----------
+        velocity : :class:`~hysop.fields.continuous_field.Field
+            solution field
+        vorticity:  :class:`~hysop.fields.continuous_field.Field`
+            right-hand side
+        variables: dict
+            dictionary of fields as keys and topologies as values.
+        projection: hysop.constants.FieldProjection or positive integer, optional
+            Project vorticity such that resolved velocity is divergence free (for 3D fields).
+            When active, projection is done prior to every solve, unless projection is 
+            an integer in which case it is done every n applies.
+            This parameter is ignored for 2D fields and defaults to no projection.
+        kwds : 
+            base class parameters.
+        """
+        
+        check_instance(velocity,  Field)
+        check_instance(vorticity, Field)
+        check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
+
+        assert velocity.domain is vorticity.domain, 'only one domain is supported'
+        assert variables[velocity] is variables[vorticity], 'only one topology is supported'
+
+        vtopology = variables[velocity]
+        wtopology = variables[vorticity]
+        input_fields  = { vorticity: wtopology }
+        output_fields = { velocity:  vtopology }
+        
+        if (projection == FieldProjection.NONE) or (projection is None) or (velocity.dim==2):
+            projection = FieldProjection.NONE
+            output_fields[vorticity] = wtopology
+        else:
+            assert (velocity.dim==3) and (vorticity.dim==3), 'Projection only available in 3D.'
+            output_fields[vorticity] = wtopology
+        
+        super(PoissonRotationalOperatorBase, self).__init__(input_fields=input_fields, 
+                output_fields=output_fields, **kwds)
+        
+        dim=velocity.dim
+        wcomp = vorticity.nb_components
+        if (dim==2) and (wcomp!=1):
+            msg='Vorticity component mistmach, got {} components but expected 1.'.format(wcomp)
+            raise RuntimeError(msg)
+        if (dim==3) and (wcomp!=3):
+            msg='Vorticity component mistmach, got {} components but expected 3.'.format(wcomp)
+            raise RuntimeError(msg)
+        if (dim!=3) and (projection!=FieldProjection.NONE):
+            raise ValueError('Velocity reprojection only available in 3D.')
+        
+        if (projection == FieldProjection.NONE):
+            self._do_project = lambda simu: False
+        elif (projection == FieldProjection.EVERY_STEP):
+            self._do_project = lambda simu: True
+        else: # projection is an integer representing frenquency 
+            freq = projection
+            assert freq>=1
+            self._do_project = lambda simu: ((simu.current_iteration % freq)==0)
+        
+        self.W = vorticity
+        self.U = velocity
+        self.projection = projection
+    
+    @debug
+    def discretize(self):
+        if self.discretized:
+            return
+        super(PoissonRotationalOperatorBase, self).discretize()
+        dW = self.input_discrete_fields[self.W]
+        dU = self.output_discrete_fields[self.U]
+        
+        W_buffers = tuple(data[dW.compute_slices] for data in dW.buffers)
+        U_buffers = tuple(data[dU.compute_slices] for data in dU.buffers)
+        
+        dim        = dU.dim
+        resolution = dU.compute_resolution
+        length     = dU.domain.length
+        dtype      = dU.dtype
+        ctype      = float_to_complex_dtype(dtype)
+        axes       = npw.arange(dim)[::-1]
+        for b in (W_buffers + U_buffers):
+            assert b.dtype == dtype, b.dtype
+            assert npw.array_equal(npw.argsort(b.strides), axes)
+        
+        self.dW = dW
+        self.dU = dU
+        self.W_buffers = W_buffers
+        self.U_buffers = U_buffers
+        
+        self.dim = dim
+        self.dtype = dtype
+        self.ctype = ctype
+        self.axes  = axes
+        self.resolution = resolution
+        self.dtype = dtype
+        self.backend = dW.backend
+
+        self.generate_wave_numbers(dim, resolution, length, dtype, ctype, axes)
+
+    @abstractmethod
+    def generate_wave_numbers(self, dim, resolution, length, dtype, ctype, axes):
+        pass

hysop/operator/poisson_rotational.py

@@ -11,8 +11,9 @@ from hysop.fields.continuous_field import Field
 from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
 from hysop.core.graph.computational_node_frontend import ComputationalGraphNodeFrontend
 
-from hysop.backend.host.fortran.operator.poisson_rotational import PoissonRotationalFFTW
-from hysop.backend.host.python.operator.poisson_rotational import PythonPoissonRotational
+from hysop.backend.host.fortran.operator.poisson_rotational  import PoissonRotationalFFTW
+from hysop.backend.host.python.operator.poisson_rotational   import PythonPoissonRotational
+from hysop.backend.device.opencl.operator.poisson_rotational import OpenClPoissonRotational
 
 class PoissonRotational(ComputationalGraphNodeFrontend):
     """
@@ -23,6 +24,7 @@ class PoissonRotational(ComputationalGraphNodeFrontend):
     
     __implementations = {
             Implementation.PYTHON: PythonPoissonRotational,
+            Implementation.OPENCL_CODEGEN: OpenClPoissonRotational,
             #Implementation.FORTRAN: PoissonRotationalFFTW
     }
 

hysop/operator/tests/test_poisson_rotational.py

@@ -23,8 +23,8 @@ class TestPoissonRotationalOperator(object):
         IO.set_default_path('/tmp/hysop_tests/test_poisson_rotational')
         
         if enable_debug_mode:
-            cls.size_min = 8
-            cls.size_max = 9
+            cls.size_min = 4
+            cls.size_max = 4
         else:
             cls.size_min = 23
             cls.size_max = 87
@@ -91,7 +91,6 @@ class TestPoissonRotationalOperator(object):
         size_max = first_not_None(size_max, self.size_max)
 
         shape = tuple(npw.random.randint(low=size_min, high=size_max+1, size=dim).tolist())
-        shape = (11, 21, 31)[:dim]
         
         domain = Box(length=(2*npw.pi,)*dim)
         U = Field(domain=domain, name='U', dtype=dtype, 
@@ -149,6 +148,14 @@ class TestPoissonRotationalOperator(object):
                 msg='   *Python FFTW: '
                 print msg,
                 yield PoissonRotational(**base_kwds)
+            elif impl is Implementation.OPENCL_CODEGEN:
+                msg='   *OpenCl CLFFT: '
+                print msg
+                for cl_env in iter_clenv():
+                    msg='     |platform {}, device {}'.format(cl_env.platform.name.strip(), 
+                                                          cl_env.device.name.strip())
+                    print msg,
+                    yield PoissonRotational(cl_env=cl_env, **base_kwds)
             else:
                 msg='Unknown implementation to test {}.'.format(impl)
                 raise NotImplementedError(msg)
@@ -174,6 +181,7 @@ class TestPoissonRotationalOperator(object):
                 Wout = tuple( data.get().handle.copy() for data in dw.data )
                 Uout = tuple( data.get().handle.copy() for data in du.data )
                 self._check_output(impl, op, Wref, Uref, Wout, Uout)
+                print
     
     @classmethod
     def _check_output(cls, impl, op, Wref, Uref, Wout, Uout):
@@ -269,7 +277,7 @@ class TestPoissonRotationalOperator(object):
     
 if __name__ == '__main__':
     TestPoissonRotationalOperator.setup_class(enable_extra_tests=False, 
-                                      enable_debug_mode=False)
+                                      enable_debug_mode=True)
     
     test = TestPoissonRotationalOperator()