working opencl poisson rotational

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

@@ -151,7 +151,7 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
                         __global {fp}2* out = (__global {fp}2*) output; 
                         __global {fp}*  K2  = (__global {fp}*)  userdata;
                         if (offset==0u) {{
-                            res = ({fp}2)(0);
+                             res = ({fp}2)(0);
                         }}
                         else {{
                             {fp} C = ({fp})(0);
@@ -162,14 +162,12 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
                             }}
                             res /= C;
                         }}
-                        printf("\\nOFFSET {i}: %u", offset);
-                        out[{i}*{N}+offset] = ({fp}2)(0,0);
+                        out[{i}*{N}+offset] = res;
                     }}
                 '''.format(fp=fp, i=i, N=N,
                             gencode=gencode)
         }
         user_data = self.K2d.data
-        print forward_callbacks['post']
         return (forward_callbacks, user_data)
 
     def _build_backward_callbacks(self, i):
@@ -241,9 +239,6 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
                     }}
                 '''.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)
 
@@ -251,10 +246,10 @@ class OpenClPoissonRotational(PoissonRotationalOperatorBase, OpenClSymbolic):
     @op_apply
     def apply(self, **kwds):
         '''Solve the PoissonRotational equation.'''
-        for fp in self.forward_plans:
+        # /!\ clFFT use the destination buffer as a scratch
+        # so we reverse the order of forward transforms.
+        for fp in self.forward_plans[::-1]:
             evt, = fp.enqueue()
-            evt.wait()
-        print self.fft_buffer
         for bp in self.backward_plans:
             evt, = bp.enqueue()
         evt = self.dW.exchange_ghosts()

hysop/operator/tests/test_poisson_rotational.py

 
-import random
+import random, primefac
 from hysop.deps import it, sm, random
 from hysop.constants import HYSOP_REAL
 from hysop.testsenv import __ENABLE_LONG_TESTS__, __HAS_OPENCL_BACKEND__
@@ -45,9 +45,10 @@ class TestPoissonRotationalOperator(object):
         #at this moment we just test a periodic solver
         frame = SymbolicFrame(dim=3)
         def gen_psi():
-            kis = (1,)*3 #tuple(random.randint(1,5) for _ in xrange(3))
-            qis = (0,)*3 #tuple(npw.random.rand(3).round(decimals=3).tolist())
-            basis = tuple( (sm.cos(ki*xi+qi), sm.sin(ki*xi+qi)) for (ki,qi,xi) in zip(kis, qis, frame.coords) )
+            kis = tuple(random.randint(1,5) for _ in xrange(3))
+            qis = tuple(npw.random.rand(3).round(decimals=3).tolist())
+            basis = tuple( (sm.cos(ki*xi+qi), sm.sin(ki*xi+qi)) \
+                            for (ki,qi,xi) in zip(kis, qis, frame.coords) )
             psi = sm.Integer(1)
             for i in xrange(dim):
                 psi *= random.choice(basis[i])
@@ -56,7 +57,8 @@ class TestPoissonRotationalOperator(object):
         analytic_solutions = {}
         analytic_functions = {}
         for dim in (2,3):
-            psis = npw.asarray([gen_psi() for _ in xrange(3)], dtype=object).view(TensorBase)
+            psis = npw.asarray([gen_psi() for _ in xrange(3)], 
+                                        dtype=object).view(TensorBase)
             if (dim==2):
                 psis[:2] = 0
             Ws = -laplacian(psis, frame)
@@ -90,7 +92,13 @@ class TestPoissonRotationalOperator(object):
         size_min = first_not_None(size_min, self.size_min)
         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())
+        valid_factors = {2,3,5,7,11,13}
+        factors = {1}
+        while (factors-valid_factors):
+            factors.clear()
+            shape = tuple(npw.random.randint(low=size_min, high=size_max+1, size=dim).tolist())
+            for Si in shape:
+                factors.update( set(primefac.primefac(int(Si-1))) )
         
         domain = Box(length=(2*npw.pi,)*dim)
         U = Field(domain=domain, name='U', dtype=dtype, 
@@ -106,8 +114,7 @@ class TestPoissonRotationalOperator(object):
     def __analytic_init(cls, data, coords, dtype, fns):
         assert len(fns) == len(data)
         for (d,fn) in zip(data,fns):
-            print fn(*coords)
-            d[...] = fn(*coords).astype(dtype)
+            d[...] = npw.asarray(fn(*coords)).astype(dtype)
 
     @staticmethod
     def __random_init(data, coords, dtype):
@@ -170,9 +177,11 @@ class TestPoissonRotationalOperator(object):
                 dw = op.input_discrete_fields[W]
                 du = op.output_discrete_fields[U]
                 
-                dw.initialize(self.__analytic_init, dtype=dtype, fns=self.analytic_functions[dim]['Ws'])
+                dw.initialize(self.__analytic_init, dtype=dtype, 
+                                    fns=self.analytic_functions[dim]['Ws'])
                 if (Uref is None):
-                    du.initialize(self.__analytic_init, dtype=dtype, fns=self.analytic_functions[dim]['Us'])
+                    du.initialize(self.__analytic_init, dtype=dtype, 
+                                    fns=self.analytic_functions[dim]['Us'])
                     Wref = tuple( data.get().handle.copy() for data in dw.data )
                     Uref = tuple( data.get().handle.copy() for data in du.data )
                 du.initialize(self.__random_init, dtype=dtype)
@@ -215,7 +224,7 @@ class TestPoissonRotationalOperator(object):
                 dist = npw.abs(fout-fref)
                 dinf = npw.max(dist)
                 deps = int(npw.ceil(dinf/eps))
-                if (deps < 200):
+                if (deps < 500):
                     print '{}eps, '.format(deps),
                     continue
                 has_nan = npw.any(npw.isnan(fout))
@@ -255,7 +264,8 @@ class TestPoissonRotationalOperator(object):
                             print
                     print
                 
-                msg = 'Test failed for {} on component {} for implementation {}.'.format(name, i, impl)
+                msg = 'Test failed for {} on component {} for implementation {}.'
+                msg = msg.format(name, i, impl)
                 raise RuntimeError(msg) 
 
 
@@ -270,15 +280,15 @@ class TestPoissonRotationalOperator(object):
         self._test(dim=3, dtype=npw.float64)
 
     def perform_tests(self):
-        # self.test_2d_float32()
+        self.test_2d_float32()
         self.test_3d_float32()
 
-        # self.test_2d_float64()
-        # self.test_3d_float64()
+        self.test_2d_float64()
+        self.test_3d_float64()
     
 if __name__ == '__main__':
     TestPoissonRotationalOperator.setup_class(enable_extra_tests=False, 
-                                      enable_debug_mode=True)
+                                      enable_debug_mode=False)
     
     test = TestPoissonRotationalOperator()