working poisson test

ci/scripts/test.sh

@@ -34,11 +34,12 @@ python -c 'import hysop; print hysop'
 
 python "$HYSOP_DIR/backend/device/codegen/kernels/tests/test_directional_advection.py"
 python "$HYSOP_DIR/backend/device/codegen/kernels/tests/test_directional_remesh.py"
+python "$HYSOP_DIR/operator/tests/test_custom_symbolic.py"
 python "$HYSOP_DIR/operator/tests/test_transpose.py"
 python "$HYSOP_DIR/operator/tests/test_directional_advection.py"
-python "$HYSOP_DIR/operator/tests/test_custom_symbolic.py"
 python "$HYSOP_DIR/operator/tests/test_directional_diffusion.py"
 python "$HYSOP_DIR/operator/tests/test_directional_stretching.py"
+python "$HYSOP_DIR/operator/tests/test_poisson_rotational.py"
 cp -r /root/.cache/* $CACHE_DIR/
 find $CACHE_DIR -name '*.lock' -delete
 

hysop/operator/tests/test_poisson_rotational.py

@@ -27,7 +27,7 @@ class TestPoissonRotationalOperator(object):
             cls.size_max = 9
         else:
             cls.size_min = 23
-            cls.size_max = 37
+            cls.size_max = 87
         
         cls.enable_extra_tests = enable_extra_tests
         cls.enable_debug_mode  = enable_debug_mode
@@ -44,7 +44,9 @@ class TestPoissonRotationalOperator(object):
         
         #at this moment we just test a periodic solver
         frame = SymbolicFrame(dim=3)
-        basis = tuple( (sm.cos(xi), sm.sin(xi)) for xi in 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) )
         def gen_psi():
             psi = sm.Integer(1)
             for i in xrange(dim):
@@ -65,10 +67,10 @@ class TestPoissonRotationalOperator(object):
                 Ws = Ws[-1:]
                 psis = psis[-1:]
                 coords = coords[:-1]
-            fWs = tuple(sm.lambdify(coords, W) for W in Ws)
-            fUs = tuple(sm.lambdify(coords, U) for U in Us)
+            fWs = tuple(sm.lambdify(coords[::-1], W) for W in Ws)
+            fUs = tuple(sm.lambdify(coords[::-1], U) for U in Us)
 
-            sols = {'Ws': Ws, Us: 'Us'}
+            sols = {'Ws': Ws, 'Us': Us}
             lambdified_sols = {'Ws': fWs, 'Us': fUs}
             analytic_solutions[dim] = sols
             analytic_functions[dim] = lambdified_sols
@@ -97,6 +99,7 @@ class TestPoissonRotationalOperator(object):
         
         self._test_one(shape=shape, dim=dim, dtype=dtype, 
                 domain=domain, U=U, W=W)
+        print
 
     @classmethod
     def __analytic_init(cls, data, coords, dtype, fns):
@@ -118,11 +121,18 @@ class TestPoissonRotationalOperator(object):
 
         print '\nTesting periodic {}D PoissonRotational: dtype={} shape={}'.format(
                 dim, dtype.__name__, shape)
+        print ' >Input analytic vorticity:'
+        for (i,wi) in enumerate(self.analytic_solutions[dim]['Ws']):
+            print '  *W{} = {}'.format(i, wi)
+        print ' >Expected output velocity is:'
+        for (i,ui) in enumerate(self.analytic_solutions[dim]['Us']):
+            print '  *U{} = {}'.format(i, ui)
+        print ' >Testing all implementations:'
 
         implementations = PoissonRotational.implementations()
        
         # Compute reference solution
-           
+        
         variables = { U:shape, W:shape }
 
         def iter_impl(impl):
@@ -130,7 +140,7 @@ class TestPoissonRotationalOperator(object):
                              implementation=impl, 
                              name='poisson_{}'.format(str(impl).lower()))
             if impl is Implementation.FORTRAN:
-                msg='  *Fortran FFTW: '
+                msg='   *Fortran FFTW: '
                 print msg,
                 yield PoissonRotational(**base_kwds)
             else:
@@ -185,10 +195,14 @@ class TestPoissonRotationalOperator(object):
                 iname = '{}{}'.format(name,i)
                 assert fout.dtype == fref.dtype, iname
                 assert fout.shape == fref.shape, iname
-
-                if npw.all(fout == fref):
-                    continue
                 
+                eps  = npw.finfo(fout.dtype).eps
+                dist = npw.abs(fout-fref)
+                dinf = npw.max(dist)
+                deps = int(npw.ceil(dinf/eps))
+                if (deps < 200):
+                    print '{}eps, '.format(deps),
+                    continue
                 has_nan = npw.any(npw.isnan(fout))
                 has_inf = npw.any(npw.isinf(fout))
                 
@@ -197,32 +211,34 @@ class TestPoissonRotationalOperator(object):
                 print 'Test output comparisson for {} failed for component {}:'.format(name, i)
                 print ' *has_nan: {}'.format(has_nan)
                 print ' *has_inf: {}'.format(has_inf)
+                print ' *dinf={} ({} eps)'.format(dinf, deps)
                 print
-                print 'REFERENCE INPUTS:'
-                for (i,w) in enumerate(Wref):
-                    print 'W{}'.format(i)
-                    print w
-                    print
-                if (name == 'U'):
-                    print 'REFERENCE OUTPUT:'
-                    for (i,u) in enumerate(Uref):
-                        print 'U{}'.format(i)
-                        print u
-                        print
-                    print
-                    print 'OPERATOR {} OUTPUT:'.format(op.name.upper())
-                    print
-                    for (i,u) in enumerate(Uout):
-                        print 'U{}'.format(i)
-                        print u
-                        print
-                else:
-                    print 'MODIFIED INPUTS:'
-                    for (i,w) in enumerate(Wout):
+                if cls.enable_debug_mode:
+                    print 'REFERENCE INPUTS:'
+                    for (i,w) in enumerate(Wref):
                         print 'W{}'.format(i)
                         print w
                         print
-                print
+                    if (name == 'U'):
+                        print 'REFERENCE OUTPUT:'
+                        for (i,u) in enumerate(Uref):
+                            print 'U{}'.format(i)
+                            print u
+                            print
+                        print
+                        print 'OPERATOR {} OUTPUT:'.format(op.name.upper())
+                        print
+                        for (i,u) in enumerate(Uout):
+                            print 'U{}'.format(i)
+                            print u
+                            print
+                    else:
+                        print 'MODIFIED INPUTS:'
+                        for (i,w) in enumerate(Wout):
+                            print 'W{}'.format(i)
+                            print w
+                            print
+                    print
                 
                 msg = 'Test failed for {} on component {} for implementation {}.'.format(name, i, impl)
                 raise RuntimeError(msg) 
@@ -250,7 +266,7 @@ class TestPoissonRotationalOperator(object):
     
 if __name__ == '__main__':
     TestPoissonRotationalOperator.setup_class(enable_extra_tests=False, 
-                                      enable_debug_mode=True)
+                                      enable_debug_mode=False)
     
     test = TestPoissonRotationalOperator()