diff --git a/ci/docker_images/ubuntu/bionic/Dockerfile b/ci/docker_images/ubuntu/bionic/Dockerfile
index dcd4f1cd110ce34a79048d31e69c72465265eacb..4c2cb4e8a5143cec09d74c8492d5ad82ac3ec774 100644
--- a/ci/docker_images/ubuntu/bionic/Dockerfile
+++ b/ci/docker_images/ubuntu/bionic/Dockerfile
@@ -65,7 +65,6 @@ RUN pip install --upgrade scipy
 RUN pip install --upgrade sympy
 RUN pip install --upgrade matplotlib
 RUN pip install --upgrade mpi4py
-RUN pip install --upgrade sphinx
 RUN pip install --upgrade h5py
 RUN pip install --upgrade gmpy2
 RUN pip install --upgrade psutil
@@ -82,6 +81,16 @@ RUN pip install --upgrade weave
 RUN pip install --upgrade argparse_color_formatter
 RUN pip install --upgrade numba
 
+# For documentation
+# RUN pip install --upgrade sphinx
+# RUN pip install --upgrade sphinxcontrib-bibtex
+# RUN pip install --upgrade sphinx_bootstrap_theme
+# RUN pip install --upgrade strip-hints
+# RUN cd /tmp && git clone https://github.com/sphinx-contrib/doxylink.git && cd doxylink/sphinxcontrib/doxylink \
+#  && mv doxylink.py doxylink.py3 && strip-hints doxylink.py3 > doxylink.py && rm doxylink.py3 \
+#  && mv parsing.py parsing.py3 && strip-hints parsing.py3 > parsing.py && rm parsing.py3 \
+#  && python setup.py install
+
 
 # scitools (python-scitools does not exist on ubuntu:bionic)
 RUN cd /tmp                                      \
diff --git a/cmake/hysop_tests.cmake b/cmake/hysop_tests.cmake
index 9978521f9a1799f500eae7b2636052a09edd1609..7c5e3c3e16e401e7201a8696015b295e1fdf4f68 100755
--- a/cmake/hysop_tests.cmake
+++ b/cmake/hysop_tests.cmake
@@ -22,7 +22,7 @@ file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/dataForTests)
 
 # Declaration of python test
 # Usage:
-#   Add_python_test(name file) with 'name' is the name of the test 
+#   Add_python_test(name file) with 'name' is the name of the test
 #   and file the source file for the test.
 macro(add_python_test test_name test_file)
   add_test(${test_name} py.test "${pytest_opt}" ${test_file})
@@ -31,7 +31,7 @@ macro(add_python_test test_name test_file)
   set_tests_properties(${test_name} PROPERTIES ENVIRONMENT "PYTHONPATH=$ENV{PYTHONPATH}:${HYSOP_BUILD_PYTHONPATH}")
   set_tests_properties(${test_name} PROPERTIES ENVIRONMENT "PYTHONPATH=$ENV{PYTHONPATH}:${HYSOP_BUILD_PYTHONPATH},LD_LIBRARY_PATH=$ENV{LD_LIBRARY_PATH}:${HYSOP_BUILD_PYTHONPATH}")
   set_tests_properties(${test_name} PROPERTIES TIMEOUT ${TEST_TIMEOUT})
-  
+
   if(WITH_MPI_TESTS)
     # Run the same test using mpi multi process run.
     # The number of processes used is set with NBPROCS_FOR_TESTS variable (user option for cmake, default=8)
@@ -44,7 +44,7 @@ macro(add_python_test test_name test_file)
     set_tests_properties(${test_name}_mpi PROPERTIES ENVIRONMENT "PYTHONPATH=$ENV{PYTHONPATH}:${HYSOP_BUILD_PYTHONPATH},LD_LIBRARY_PATH=$ENV{LD_LIBRARY_PATH}:${HYSOP_BUILD_PYTHONPATH}")
     set_tests_properties(${test_name}_mpi PROPERTIES TIMEOUT ${TEST_TIMEOUT})
  endif()
-  
+
 endmacro()
 
 # Choose python build dir as directory where tests will be run.
@@ -52,7 +52,7 @@ endmacro()
 set(testDir ${HYSOP_BUILD_PYTHONPATH})
 
 # === Set the list of all directories which may contain tests ===
-set(py_src_dirs core/arrays numerics fields operator)
+set(py_src_dirs core/arrays core/graph core/memory fields operator numerics numerics/splitting)
 
 # === Create the files list from all directories in py_src_dirs ===
 
@@ -68,7 +68,7 @@ endforeach()
 
 # Build doctests
 # Handling doctest in *.py files recursively for each directory of hysop/${py_src_dirs}
-# excluding  __init__ or test_ files. 
+# excluding  __init__ or test_ files.
 # Doctest are run for every line which contains '>>>'
 set(py_doctest_files)
 foreach(testdir ${py_src_dirs})
@@ -96,6 +96,11 @@ foreach(testfile ${py_doctest_files})
   add_test("doctest_${testName}" py.test -v --doctest-modules ${exename})
 endforeach()
 
+# Setup extra timeout for some specific test (used with ctest, not in CI)
+set_tests_properties("test_poisson" PROPERTIES TIMEOUT 600)
+set_tests_properties("test_poisson_curl" PROPERTIES TIMEOUT 600)
+set_tests_properties("test_solenoidal_projection" PROPERTIES TIMEOUT 600)
+set_tests_properties("test_spectral_derivative" PROPERTIES TIMEOUT 600)
+
 # Add custom target to launch continuous integration tests
 add_custom_target(ci COMMAND ${CMAKE_SOURCE_DIR}/ci/scripts/test.sh ${CMAKE_BINARY_DIR} ${CMAKE_SOURCE_DIR}/hysop)
-
diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt
index a571f89d7ee41266eb79179eb06b92688b35b5ae..4c58de609b76f0d85f428a0cd09da7d3b5293a6d 100644
--- a/docs/CMakeLists.txt
+++ b/docs/CMakeLists.txt
@@ -4,7 +4,7 @@
 # - automatic doxygen documentation, from C, Fortran, C++ sources
 # - automatic sphinx/apidoc documentaion from python sources
 # - sphinx documentation for manuals
-# 
+#
 # ===============================================================
 
 if(WITH_DOCUMENTATION)
@@ -12,7 +12,7 @@ if(WITH_DOCUMENTATION)
   if(NOT DOC_ROOT_DIR)
     set(DOC_ROOT_DIR ${CMAKE_BINARY_DIR}/docs/build/html CACHE INTERNAL "Path to html documentation.")
   endif()
-  
+
   # --------------------- Doxygen setup ---------------------
   find_package(Doxygen)
 
@@ -25,7 +25,7 @@ if(WITH_DOCUMENTATION)
   file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/build/html)
   # config file name (generated later from hysop.doxyfile.in)
   set(DOXY_CONFIG "${CMAKE_CURRENT_BINARY_DIR}/config/hysop.doxyfile")
-  
+
   # --------------------------------
   # Update doxy.config to take
   # into account sources/headers files
@@ -35,12 +35,12 @@ if(WITH_DOCUMENTATION)
         ${CMAKE_SOURCE_DIR}/src
         ${CMAKE_SOURCE_DIR}/hysop
     )
-  
+
   foreach(_dir ${DOXY_INPUTS})
     set(DOXYGEN_INPUTS "${DOXYGEN_INPUTS} ${_dir}")
   endforeach()
   list(REMOVE_DUPLICATES DOXYGEN_INPUTS)
-  
+
   set(GENERATE_HTML YES)
   set(GENERATE_XML YES)
   configure_file(${CMAKE_CURRENT_SOURCE_DIR}/config/hysop.doxyfile.in ${DOXY_CONFIG})
@@ -54,7 +54,7 @@ if(WITH_DOCUMENTATION)
   include(FindSphinxModule)
   find_python_module(sphinx REQUIRED)
   find_sphinx_module(sphinxcontrib bibtex REQUIRED)
-  find_sphinx_module(sphinxcontrib doxylink REQUIRED)
+#import ISSUE:  find_sphinx_module(sphinxcontrib doxylink REQUIRED)
 
   # generate sphinx conf file
   configure_file (
@@ -107,7 +107,7 @@ if(WITH_DOCUMENTATION)
       "${CMAKE_CURRENT_SOURCE_DIR}/sphinx/${filename}"
       "${CMAKE_CURRENT_BINARY_DIR}/sphinx/${filename}" @ONLY)
   endforeach()
-  
+
   # --> use make html or make latex ... to build documentation.
   set(SPHINX_EXECUTABLE sphinx-build )
   set(SPHINX_PARAMETERS
@@ -128,13 +128,13 @@ if(WITH_DOCUMENTATION)
     -d build/doctrees # path to doctree files
     ${CMAKE_CURRENT_BINARY_DIR}/sphinx  # path to rst source
     )
-  
+
   add_custom_target(latex
     COMMAND ${SPHINX_EXECUTABLE} ${SPHINX_LATEX_PARAMETERS}
     ${CMAKE_CURRENT_BINARY_DIR}/build/latex
     )
 
-  
+
   add_custom_target(doc DEPENDS html)
   # Doxygen css comes from sphinx setup
   # --> link to sphinx
@@ -148,5 +148,5 @@ if(WITH_DOCUMENTATION)
   # #      DEPENDS update_css)
   # add_dependencies(html doxygen)
   # add_dependencies(html apidoc)
-  
+
 endif()
diff --git a/hysop/backend/device/codegen/kernels/tests/test_directional_advection.py b/hysop/backend/device/codegen/kernels/tests/test_codegen_directional_advection.py
similarity index 100%
rename from hysop/backend/device/codegen/kernels/tests/test_directional_advection.py
rename to hysop/backend/device/codegen/kernels/tests/test_codegen_directional_advection.py
diff --git a/hysop/backend/device/codegen/kernels/tests/test_directional_remesh.py b/hysop/backend/device/codegen/kernels/tests/test_codegen_directional_remesh.py
similarity index 100%
rename from hysop/backend/device/codegen/kernels/tests/test_directional_remesh.py
rename to hysop/backend/device/codegen/kernels/tests/test_codegen_directional_remesh.py
diff --git a/hysop/backend/device/codegen/kernels/tests/test_directional_stretching.py b/hysop/backend/device/codegen/kernels/tests/test_codegen_directional_stretching.py
similarity index 100%
rename from hysop/backend/device/codegen/kernels/tests/test_directional_stretching.py
rename to hysop/backend/device/codegen/kernels/tests/test_codegen_directional_stretching.py
diff --git a/hysop/backend/device/codegen/kernels/tests/test_transpose.py b/hysop/backend/device/codegen/kernels/tests/test_codegen_transpose.py
similarity index 100%
rename from hysop/backend/device/codegen/kernels/tests/test_transpose.py
rename to hysop/backend/device/codegen/kernels/tests/test_codegen_transpose.py
diff --git a/hysop/core/graph/tests/test_graph.py b/hysop/core/graph/tests/test_graph.py
index 25746378e3a27f826ae22f947aaf2553c086f734..a359a69b1368d63ab28d1530216f5905673f8d10 100644
--- a/hysop/core/graph/tests/test_graph.py
+++ b/hysop/core/graph/tests/test_graph.py
@@ -1,4 +1,3 @@
-
 from hysop.domain.box        import Box
 from hysop.topology.cartesian_topology import CartesianTopology
 from hysop.tools.parameters  import CartesianDiscretization
@@ -9,7 +8,7 @@ class _ComputationalGraph(ComputationalGraph):
     @classmethod
     def supports_multiple_topologies(cls):
         return True
-    
+
 class _ComputationalGraphOperator(ComputationalGraphOperator):
     def apply(self):
         pass
@@ -18,71 +17,88 @@ class _ComputationalGraphOperator(ComputationalGraphOperator):
     def supports_multiple_topologies(cls):
         return True
 
-def test_graph_build(display=False):
-    box = Box(length=[1.0]*3, origin=[0.0]*3)
-    Vg    = Field(domain=box, name='Vg', is_vector=True)
-    Vp    = Field(domain=box, name='Vp', is_vector=True)
-    Wg    = Field(domain=box, name='Wg', is_vector=True)
-    Wp    = Field(domain=box, name='Wp', is_vector=True)
-    rho0g = Field(domain=box, name='rho0g')
-    rho0p = Field(domain=box, name='rho0p')
-    rho1g = Field(domain=box, name='rho1g')
-    rho1p = Field(domain=box, name='rho1p')
-    
-    d3d0 = CartesianDiscretization(resolution=(64,64,64), ghosts=None, default_boundaries=True)
-    d3d1 = CartesianDiscretization(resolution=(128,128,128), ghosts=None, default_boundaries=True)
-    t0  = CartesianTopology(domain=box, discretization=d3d0)
-    t1  = CartesianTopology(domain=box, discretization=d3d1)
-
-    ops = [
-                ('copyW',     [Wg],         [Wp]),
-                ('copyS0',    [rho0g],      [rho0p]),
-                ('copyS1',    [rho1g],      [rho1p]),
-                ('copyV',     [Vg],         [Vg]),
-            ]
-    g0 = _ComputationalGraph(name='g0')
-    for (opname,_in,_out) in ops:
-        invars  = {}
-        outvars = {}
-        for var in _in:
-            invars[var]  = t0 if (var.name.find('rho')) else t1
-        for var in _out:
-            outvars[var] = t0 if (var.name.find('rho')) else t1
-        g0.push_nodes(_ComputationalGraphOperator(name=opname,
-            input_fields=invars,output_fields=outvars))
-   
-    ops = [
-            ('advecWx',   [Vg, Wp],   [Wg]),
-            ('stretchWx', [Vg, Wg],   [Wg]),
-            ('diffWx',    [Wg],       [Wg]),
-            ('advecS0x',  [Vg,rho0p], [rho0g]),
-            ('advecS1x',  [Vg,rho1p], [rho1g]),
-            ('diffS0x',   [rho0g],    [rho0g]),
-            ('diffS1x',   [rho1g],    [rho1g]),
-            ('forceWx',   [rho0g,rho1g,Wg], [Wg]),
-            ]
-    g1 = _ComputationalGraph(name='g1')
-    for (opname,_in,_out) in ops:
-        invars  = {}
-        outvars = {}
-        for var in _in:
-            invars[var]  = t0 if (var.name.find('rho')) else t1
-        for var in _out:
-            outvars[var] = t0 if (var.name.find('rho')) else t1
-        g1.push_nodes(_ComputationalGraphOperator(name=opname,
-            input_fields=invars,output_fields=outvars))
-   
-    g = _ComputationalGraph(name='exec')
-    g.push_nodes(g0)
-    g.push_nodes(g1)
-    
-    g.initialize()
-    g.discretize()
-    g.setup(None)
-    g.apply()
-    
-    if display:
-        g.display()
+class TestGraph(object):
+
+    @classmethod
+    def setup_class(cls):
+        pass
+
+    @classmethod
+    def teardown_class(cls):
+        pass
+
+    def setup_method(self, method):
+        pass
+
+    def teardown_method(self, method):
+        pass
+
+    def test_graph_build(self, display=False):
+        box = Box(length=[1.0]*3, origin=[0.0]*3)
+        Vg    = Field(domain=box, name='Vg', is_vector=True)
+        Vp    = Field(domain=box, name='Vp', is_vector=True)
+        Wg    = Field(domain=box, name='Wg', is_vector=True)
+        Wp    = Field(domain=box, name='Wp', is_vector=True)
+        rho0g = Field(domain=box, name='rho0g')
+        rho0p = Field(domain=box, name='rho0p')
+        rho1g = Field(domain=box, name='rho1g')
+        rho1p = Field(domain=box, name='rho1p')
+
+        d3d0 = CartesianDiscretization(resolution=(64,64,64), ghosts=None, default_boundaries=True)
+        d3d1 = CartesianDiscretization(resolution=(128,128,128), ghosts=None, default_boundaries=True)
+        t0  = CartesianTopology(domain=box, discretization=d3d0)
+        t1  = CartesianTopology(domain=box, discretization=d3d1)
+
+        ops = [
+                    ('copyW',     [Wg],         [Wp]),
+                    ('copyS0',    [rho0g],      [rho0p]),
+                    ('copyS1',    [rho1g],      [rho1p]),
+                    ('copyV',     [Vg],         [Vg]),
+                ]
+        g0 = _ComputationalGraph(name='g0')
+        for (opname,_in,_out) in ops:
+            invars  = {}
+            outvars = {}
+            for var in _in:
+                invars[var]  = t0 if (var.name.find('rho')) else t1
+            for var in _out:
+                outvars[var] = t0 if (var.name.find('rho')) else t1
+            g0.push_nodes(_ComputationalGraphOperator(name=opname,
+                input_fields=invars,output_fields=outvars))
+
+        ops = [
+                ('advecWx',   [Vg, Wp],   [Wg]),
+                ('stretchWx', [Vg, Wg],   [Wg]),
+                ('diffWx',    [Wg],       [Wg]),
+                ('advecS0x',  [Vg,rho0p], [rho0g]),
+                ('advecS1x',  [Vg,rho1p], [rho1g]),
+                ('diffS0x',   [rho0g],    [rho0g]),
+                ('diffS1x',   [rho1g],    [rho1g]),
+                ('forceWx',   [rho0g,rho1g,Wg], [Wg]),
+                ]
+        g1 = _ComputationalGraph(name='g1')
+        for (opname,_in,_out) in ops:
+            invars  = {}
+            outvars = {}
+            for var in _in:
+                invars[var]  = t0 if (var.name.find('rho')) else t1
+            for var in _out:
+                outvars[var] = t0 if (var.name.find('rho')) else t1
+            g1.push_nodes(_ComputationalGraphOperator(name=opname,
+                input_fields=invars,output_fields=outvars))
+
+        g = _ComputationalGraph(name='exec')
+        g.push_nodes(g0)
+        g.push_nodes(g1)
+
+        g.initialize()
+        g.discretize()
+        g.setup(None)
+        g.apply()
+
+        if display:
+            g.display()
 
 if __name__ == '__main__':
-    test_graph_build(display=False)
+    test = TestGraph()
+    test.test_graph_build(display=False)
diff --git a/hysop/core/memory/memory_request.py b/hysop/core/memory/memory_request.py
index e6c65a33aa2b712a6c051e2f29f0b9f3eab8f18d..747f7b870597e37e1e78bfe46a380c972ee79bfb 100644
--- a/hysop/core/memory/memory_request.py
+++ b/hysop/core/memory/memory_request.py
@@ -494,12 +494,8 @@ class MultipleOperatorMemoryRequests(object):
 
             ss=''
             for (backend, backend_srequests) in all_requests.iteritems():
-<<<<<<< HEAD
-                kind = backend.kind
-=======
                 total = totals[backend]
                 kind = backend.kind
->>>>>>> master
                 if (kind == Backend.OPENCL):
                     precision = u' on device {}'.format(backend.device.name.strip())
                 else:
diff --git a/hysop/core/memory/tests/test_allocator.py b/hysop/core/memory/tests/test_allocator.py
index 919e65bc553e904bd082587053ae42408308ce63..eff44af2151080206bed02bf30522d351abe29af 100644
--- a/hysop/core/memory/tests/test_allocator.py
+++ b/hysop/core/memory/tests/test_allocator.py
@@ -1,13 +1,13 @@
-
 from hysop.backend import __HAS_OPENCL_BACKEND__, __HAS_CUDA_BACKEND__
 from hysop.backend.host.host_allocator import HostAllocator
+from hysop.core.mpi import default_mpi_params
 
 sizes = (1,89,7919,32*1024*1024)
 alignments = (1,2,4,8,16,32,64,128,256,512,1024)
 
 def test_host_allocator():
     allocator = HostAllocator()
-    
+
     for nbytes in sizes:
         for alignment in alignments:
             buf = allocator.allocate(nbytes)
@@ -27,16 +27,16 @@ if __HAS_OPENCL_BACKEND__:
                                                              OpenClDeferredAllocator
     from hysop.backend.device.opencl.opencl_tools import get_or_create_opencl_env
 
-    cl_env = get_or_create_opencl_env()
+    cl_env = get_or_create_opencl_env(mpi_params=default_mpi_params())
 
     def test_opencl_immediate_allocator():
-        allocator = OpenClImmediateAllocator(cl_env.default_queue)
+        allocator = OpenClImmediateAllocator(queue=cl_env.default_queue)
         for nbytes in sizes:
             buf = allocator.allocate(nbytes)
             assert buf.size == nbytes
-    
+
     def test_opencl_deferred_allocator():
-        allocator = OpenClDeferredAllocator(cl_env.default_queue)
+        allocator = OpenClDeferredAllocator(queue=cl_env.default_queue)
         for nbytes in sizes:
             buf = allocator.allocate(nbytes)
             cl.enqueue_barrier(cl_env.default_queue)
diff --git a/hysop/core/memory/tests/test_buffer.py b/hysop/core/memory/tests/test_buffer.py
index e7bcd97182a99643aaef5887f7a5668dbbbc27b1..abb00ff2d4b576a7d33ce724539a80fab01e5765 100644
--- a/hysop/core/memory/tests/test_buffer.py
+++ b/hysop/core/memory/tests/test_buffer.py
@@ -1,6 +1,6 @@
-
 from hysop.deps import np
 from hysop.backend.host.host_buffer import HostBuffer
+from hysop.core.mpi import default_mpi_params
 from hysop.backend import __HAS_OPENCL_BACKEND__, __HAS_CUDA_BACKEND__
 
 size = 32*1024*1024 # 32Mo
@@ -10,19 +10,19 @@ def test_host_buffer():
     assert buf1.get_int_ptr() == buf1.int_ptr
     assert buf1.get_size() == buf1.size
     assert buf1.get_buffer() is buf1.buf
-    
+
     buf2 = HostBuffer(size=size)
     assert buf1.int_ptr != buf2.int_ptr
     assert buf1.size == size
     assert buf1.size == buf2.size
-    
+
     arr1 = np.frombuffer(buf1, dtype=np.uint8)
     arr2 = np.frombuffer(buf2, dtype=np.uint8)
     arr1.fill(0)
     arr2.fill(1)
     assert (arr1==0).all()
     assert (arr2==1).all()
-    
+
     buf2 = HostBuffer.from_buffer(buf1)
     assert buf1.int_ptr == buf2.int_ptr
     assert buf1.size == buf2.size
@@ -30,7 +30,7 @@ def test_host_buffer():
     buf3 = HostBuffer.from_int_ptr(buf2.int_ptr, buf2.size)
     assert buf1.int_ptr == buf3.int_ptr
     assert buf1.size == buf3.size
-    
+
     arr1 = np.frombuffer(buf1, dtype=np.uint8)
     arr2 = np.frombuffer(buf2, dtype=np.uint8)
     arr3 = np.frombuffer(buf3, dtype=np.uint8)
@@ -75,18 +75,18 @@ if __HAS_OPENCL_BACKEND__:
     from hysop.backend.device.opencl.opencl_buffer import OpenClBuffer
     from hysop.backend.device.opencl.opencl_tools import get_or_create_opencl_env
 
-    cl_env = get_or_create_opencl_env()
-    
+    cl_env = get_or_create_opencl_env(mpi_params=default_mpi_params())
+
     def test_opencl_buffer():
         mf = cl.mem_flags
-            
+
         buf  = cl.Buffer(cl_env.context, size=size, flags=mf.READ_ONLY)
 
         buf1 = OpenClBuffer(cl_env.context, size=size, mem_flags=mf.READ_WRITE)
         assert buf1.int_ptr == buf1.get_int_ptr()
-        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf1, 
+        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf1,
                 byte_count=size).wait()
-        
+
         assert buf1.ref_count() == 1
         buf1.release()
 
@@ -102,26 +102,26 @@ if __HAS_OPENCL_BACKEND__:
         assert buf3.ref_count() == 3
         assert buf1.int_ptr == buf2.int_ptr
         assert buf1.int_ptr == buf3.int_ptr
-        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf1, 
+        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf1,
                 byte_count=size).wait()
-        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf2, 
+        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf2,
                 byte_count=size).wait()
-        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf3, 
+        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf3,
                 byte_count=size).wait()
         buf1.release()
         assert buf2.ref_count() == 2
         assert buf3.ref_count() == 2
-        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf2, 
+        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf2,
                 byte_count=size).wait()
-        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf3, 
+        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf3,
                 byte_count=size).wait()
         buf3.release()
         assert buf2.ref_count() == 1
-        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf2, 
+        cl.enqueue_copy_buffer(queue=cl_env.default_queue, src=buf, dst=buf2,
                 byte_count=size).wait()
         buf2.release()
 
-       
+
     def test_opencl_host_device_buffer():
         mf = cl.mem_flags
 
@@ -131,11 +131,11 @@ if __HAS_OPENCL_BACKEND__:
         # try to use unified memory (or something like zero copy memory)
         hbuf = np.ndarray(shape=size, dtype=np.float32)
         buf = OpenClBuffer(cl_env.context, mem_flags=mf.USE_HOST_PTR, hostbuf=hbuf)
-        
+
         # copy host buf at creation
         hbuf = np.ndarray(shape=(1024,1024,), dtype=np.int64)
         buf = OpenClBuffer(cl_env.context, mem_flags=mf.COPY_HOST_PTR, hostbuf=hbuf)
-        
+
 
 if __name__ == '__main__':
     test_host_buffer()
diff --git a/hysop/core/memory/tests/test_mempool.py b/hysop/core/memory/tests/test_mempool.py
index 64e1f07fae9dbe11fa05420571ff164a54fabddd..ac43f5ee301aef79d4979afb4663b40d2c212704 100644
--- a/hysop/core/memory/tests/test_mempool.py
+++ b/hysop/core/memory/tests/test_mempool.py
@@ -1,9 +1,8 @@
-
 from hysop.deps import np
 from hysop.testsenv import opencl_failed, iter_clenv, \
                            __HAS_OPENCL_BACKEND__, __ENABLE_LONG_TESTS__
 from hysop.core.memory.mempool import MemoryPool
-    
+
 import random
 max_bytes_per_alloc = 1024*1024*128  # 128MB
 free   = lambda: bool(random.random()>0.1) # 80% probability of free
@@ -19,11 +18,11 @@ def test_mempool_opencl_immediate_allocator():
     from hysop.backend.device.opencl.opencl_allocator import OpenClImmediateAllocator
 
     for cl_env in iter_clenv():
-        allocator = OpenClImmediateAllocator(cl_env.default_queue)
+        allocator = OpenClImmediateAllocator(queue=cl_env.default_queue)
         _test_mempool_allocator(cl_env.platform.name, allocator)
 
 def _test_mempool_allocator(name, allocator):
-    pool = allocator.memory_pool(name=name, verbose=True)
+    pool = allocator.memory_pool(name=name)
     buffers = []
     try:
         while True:
diff --git a/hysop/fields/discrete_field.py b/hysop/fields/discrete_field.py
index e909ebbc733ff4ba59be47b8452ebf5651a55275..91be593970a969dd2899374a7763a4172ac11b43 100644
--- a/hysop/fields/discrete_field.py
+++ b/hysop/fields/discrete_field.py
@@ -25,7 +25,7 @@ from hysop.mesh.mesh import MeshView
 
 class DiscreteScalarFieldViewContainerI(object):
     """
-    Common abstract interface for scalar and tensor-like container of 
+    Common abstract interface for scalar and tensor-like container of
     discrete field views.
     """
 
@@ -38,7 +38,7 @@ class DiscreteScalarFieldViewContainerI(object):
     @abstractmethod
     def discrete_field_views(self):
         """
-        Return all unique discrete field views contained in this discrete 
+        Return all unique discrete field views contained in this discrete
         field view container.
         """
         pass
@@ -91,7 +91,7 @@ class DiscreteScalarFieldViewContainerI(object):
     @abstractmethod
     def clone(self, tstate=None):
         """
-        Create a new temporary DiscreteScalarField container and allocate it 
+        Create a new temporary DiscreteScalarField container and allocate it
         like the current object, with possibly a different topology state.
 
         This should only be used for debugging and testing purpose.
@@ -104,7 +104,7 @@ class DiscreteScalarFieldViewContainerI(object):
     def tmp_dfield_like(self, name, **kwds):
         """
         Create a new Field container and a new temporary CartesianDiscreteField.
-        like the current object, possibly on a different backend. 
+        like the current object, possibly on a different backend.
         /!\ The returned discrete field is not allocated.
         """
         pass
@@ -117,7 +117,7 @@ class DiscreteScalarFieldViewContainerI(object):
     @abstractmethod
     def build_ghost_exchanger(self, **kwds):
         """
-        Build a ghost exchanger, possibly on different data. 
+        Build a ghost exchanger, possibly on different data.
         Usefull for operator apply.
         """
         pass
@@ -149,7 +149,7 @@ class DiscreteScalarFieldViewContainerI(object):
         Return a view on this DiscreteScalarField using given topology state.
         """
         pass
-    
+
     @abstractmethod
     def as_any_dfield(self, memory_order, **kwds):
         """
@@ -201,7 +201,7 @@ class DiscreteScalarFieldViewContainerI(object):
     def _get_buffers():
         """Return all array data as a buffers as a tuple."""
         pass
-    
+
     @abstractmethod
     def _get_sdata(self):
         """Return contained array."""
@@ -216,7 +216,7 @@ class DiscreteScalarFieldViewContainerI(object):
     def _get_sbuffer():
         """Return contained buffer."""
         pass
-    
+
     def get_attributes(self, *attrs):
         """
         Return all matching attributes contained in self.discrete_field_views(),
@@ -301,9 +301,9 @@ class DiscreteScalarFieldViewContainerI(object):
     def get_unique_attribute(self, *attr):
         """
         Try to return the unique attribute common to all contained discrete fields.
-        Raise an AttributeError if a attribute is not unique accross contained 
+        Raise an AttributeError if a attribute is not unique accross contained
         discrete field views.
-        
+
         /!\ Can be slow to evaluate due to uniqueness check /!\
         """
         if self.has_unique_attribute(*attr):
@@ -401,7 +401,7 @@ class DiscreteScalarFieldView(DiscreteScalarFieldViewContainerI, TaggedObjectVie
     def __new__(cls, dfield, topology_state, **kwds):
         check_instance(dfield, DiscreteScalarField, allow_none=issubclass(cls, DiscreteScalarField))
         check_instance(topology_state, TopologyState)
-        obj = super(DiscreteScalarFieldView, cls).__new__(cls, obj_view=dfield, 
+        obj = super(DiscreteScalarFieldView, cls).__new__(cls, obj_view=dfield,
                 variable_kind=Variable.DISCRETE_FIELD, **kwds)
         dfield = first_not_None(dfield, obj)
         obj._dfield = dfield
@@ -453,7 +453,7 @@ class DiscreteScalarFieldView(DiscreteScalarFieldViewContainerI, TaggedObjectVie
     def _get_is_read_only(self):
         """Return true if this view is a read-only."""
         return self._topology_state.is_read_only
-    
+
     def _get_name(self):
         """Get the name of the discrete field."""
         return self._dfield._name
@@ -466,7 +466,7 @@ class DiscreteScalarFieldView(DiscreteScalarFieldViewContainerI, TaggedObjectVie
     def _get_var_name(self):
         """Get the latex name of the discrete field."""
         return self._dfield._var_name
-    
+
     def _get_dtype(self):
         """Get the data type of the discrete field."""
         return self._dfield._field.dtype
@@ -506,10 +506,10 @@ class DiscreteScalarFieldView(DiscreteScalarFieldViewContainerI, TaggedObjectVie
             return not eq
         else:
             return eq
-    
+
     def honor_memory_request(self, *args, **kwds):
         self._dfield.honor_memory_request(*args, **kwds)
-    
+
     def _get_memory_request(self):
         """Get memory request that should be allocated for this TmpCartesianDiscreteField."""
         return getattr(self._dfield, '_memory_request', None)
@@ -518,7 +518,7 @@ class DiscreteScalarFieldView(DiscreteScalarFieldViewContainerI, TaggedObjectVie
         return getattr(self._dfield, '_memory_request_id', None)
 
     def __hash__(self):
-        h  =   id(self._dfield._topology) 
+        h  =   id(self._dfield._topology)
         h ^=   id(self._dfield._field)
         h ^= hash(self._dfield._name)
         h ^= hash(self._dfield._pretty_name)
@@ -536,12 +536,12 @@ class DiscreteScalarFieldView(DiscreteScalarFieldViewContainerI, TaggedObjectVie
     field = property(_get_field)
     topology_state = property(_get_topology_state)
     is_read_only = property(_get_is_read_only)
-    
+
     name  = property(_get_name)
     pretty_name  = property(_get_pretty_name)
     latex_name = property(_get_latex_name)
     var_name = property(_get_var_name)
-    
+
     dtype = property(_get_dtype)
     initial_values = property(_get_initial_values)
 
@@ -563,21 +563,21 @@ class DiscreteScalarField(NamedScalarContainerI, TaggedObject):
     wich are a collection of numpy like multidimensional arrays allocated using a
     specific backend (hysop.core.arrays.array.ArrayBackend).
 
-    A DiscreteScalarField is the result of discretizing a continuous Field 
-    (hysop.field.continuous_field.Field) defined on a specific domain 
-    (hysop.domain.domain.Domain) distributed accross processes through 
+    A DiscreteScalarField is the result of discretizing a continuous Field
+    (hysop.field.continuous_field.Field) defined on a specific domain
+    (hysop.domain.domain.Domain) distributed accross processes through
     a topology (hysop.topology.topology.Topology).
 
     Ghost exchangers are automatically built for all discrete fields.
-    Ghost exchangers may require additional memory buffers, 
+    Ghost exchangers may require additional memory buffers,
     depending on the discrete field topology backend and the ghost exchange strategy.
     """
 
     __metaclass__ = ABCMeta
 
     @debug
-    def __new__(cls, field, topology, register_discrete_field=True, 
-                name=None, pretty_name=None, 
+    def __new__(cls, field, topology, register_discrete_field=True,
+                name=None, pretty_name=None,
                 var_name=None, latex_name=None,
                 **kwds):
         """
@@ -609,12 +609,12 @@ class DiscreteScalarField(NamedScalarContainerI, TaggedObject):
         check_instance(pretty_name, (str,unicode), allow_none=True)
 
         _name, _pretty_name, _var_name, _latex_name = \
-                cls.format_discrete_names(field.name, 
-                                          field.pretty_name, 
+                cls.format_discrete_names(field.name,
+                                          field.pretty_name,
                                           field.var_name,
                                           field.latex_name,
                                           topology)
-        
+
         pretty_name = first_not_None(pretty_name, name, _pretty_name)
         var_name    = first_not_None(var_name,    name, _var_name)
         latex_name  = first_not_None(latex_name,  name, _latex_name)
@@ -668,8 +668,8 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
     A tensor discrete field is a collection of scalar discrete fields views.
 
     This object handles a numpy.ndarray of discrete scalar field views,
-    which may have different attributes (different data types for 
-    example) and different topology states. 
+    which may have different attributes (different data types for
+    example) and different topology states.
 
     A tensor field garanties that each different field objects have
     unique names and pretty names within the tensor field. A given
@@ -678,20 +678,20 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
 
     Some components may also be set to None.
 
-    Is also garanties that all fields shares the same domain, but contained 
+    Is also garanties that all fields shares the same domain, but contained
     discrete fields may be defined on different topologies.
     """
 
-    def __new__(cls, field, dfields, name=None, 
+    def __new__(cls, field, dfields, name=None,
                 pretty_name=None, latex_name=None, **kwds):
         check_instance(field, TensorField)
         check_instance(dfields, npw.ndarray, dtype=object, values=DiscreteScalarFieldView)
         assert npw.array_equal(field.shape, dfields.shape)
-        
+
         tensor_cls = cls.determine_tensor_cls(dfields)
         if (tensor_cls is not cls):
-            return tensor_cls.__new__(tensor_cls, 
-                                      field=field, dfields=dfields, 
+            return tensor_cls.__new__(tensor_cls,
+                                      field=field, dfields=dfields,
                                       name=name, pretty_name=pretty_name,
                                       **kwds)
 
@@ -701,9 +701,9 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
         pretty_name = first_not_None(pretty_name, _pretty_name)
         latex_name  = first_not_None(latex_name, _latex_name)
 
-        obj = super(DiscreteTensorField, cls).__new__(cls, name=name, 
+        obj = super(DiscreteTensorField, cls).__new__(cls, name=name,
                 pretty_name=pretty_name, latex_name=latex_name,
-                tag_prefix='tdf', tagged_cls=DiscreteTensorField, 
+                tag_prefix='tdf', tagged_cls=DiscreteTensorField,
                 contained_objects=dfields, **kwds)
         obj._field   = field
         obj._dfields = dfields
@@ -731,7 +731,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
     @classmethod
     def from_dfields(cls, name, dfields, shape, pretty_name=None, **kwds):
         """
-        Create a TensorField and a DiscreteTensorField 
+        Create a TensorField and a DiscreteTensorField
         from a list of discrete fields and a shape.
         """
         dfields = to_tuple(dfields)
@@ -747,7 +747,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
                                         fields=fields, shape=shape)
 
         dfields = npw.asarray(dfields, dtype=object).reshape(shape)
-        
+
         return cls(field=field, dfields=dfields, **kwds)
 
     @classmethod
@@ -773,7 +773,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
             msg='Tensor discrete field {} should at least contain a valid DiscreteScalarField.'
             msg=msg.format(name)
             raise ValueError(msg)
-    
+
     def _check_names(self):
         """Check that discrete fields names are unique."""
         names = {}
@@ -794,7 +794,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
 
     def discrete_field_views(self):
         """
-        Return all unique discrete field views contained in this discrete 
+        Return all unique discrete field views contained in this discrete
         field view container.
         """
         ordered_dfields = self._dfields.ravel().tolist()
@@ -836,7 +836,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
     def clone(self, tstate=None,
                 name=None, pretty_name=None, **kwds):
         """
-        Create a new temporary DiscreteScalarField container and allocate it 
+        Create a new temporary DiscreteScalarField container and allocate it
         like the current object, with possibly a different topology state.
 
         This should only be used for debugging and testing purpose.
@@ -844,19 +844,19 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
         field.
         """
         name = first_not_None(name, '{}__{}'.format(self.name, self._clone_id))
-        pretty_name = first_not_None(pretty_name, '{}__{}'.format(self.pretty_name, 
+        pretty_name = first_not_None(pretty_name, '{}__{}'.format(self.pretty_name,
                                                                   self._clone_id))
         dfields = npw.empty(shape=self.shape, dtype=object)
         for (idx, dfield) in self.nd_iter():
             dfields[idx] = dfield.clone(tstate=tstate, **kwds)
         self._clone_id += 1
-        return self.from_dfield_array(name=name, pretty_name=pretty_name, dfields=dfields, 
+        return self.from_dfield_array(name=name, pretty_name=pretty_name, dfields=dfields,
                                         **kwds)
 
     def tmp_dfield_like(self, name, pretty_name=None, **kwds):
         """
         Create a new Field container and a new temporary CartesianDiscreteField.
-        like the current object, possibly on a different backend. 
+        like the current object, possibly on a different backend.
         /!\ The returned discrete field is not allocated.
         """
         from hysop.core.memory.memory_request import OperatorMemoryRequests
@@ -867,9 +867,10 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
             _name = TensorField.default_name_formatter(basename=name, idx=idx)
             _pretty_name = TensorField.default_pretty_name_formatter(
                     basename=pretty_name, idx=idx)
-            _name, _pretty_name = DiscreteScalarField.format_discrete_names(_name, _pretty_name, 
-                    dfield.topology)
-            (dfield, request, request_id) = dfield.tmp_dfield_like(name=_name, 
+            _name, _pretty_name, _var_name, _latex_name = \
+                    DiscreteScalarField.format_discrete_names(_name, _pretty_name,
+                                                              _name, self.latex_name, dfield.topology)
+            (dfield, request, request_id) = dfield.tmp_dfield_like(name=_name,
                                                        pretty_name=_pretty_name,
                                                        **kwds)
             requests.push_mem_request(request_id, request)
@@ -877,7 +878,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
             dfields[idx] = dfield
         dfield = self.from_dfield_array(name=name, pretty_name=pretty_name, dfields=dfields)
         return (dfield, requests)
-   
+
     def honor_memory_request(self, work, op=None):
         """Honour memory requests for contained temporary discrete fields."""
         op = first_not_None(op, self.name)
@@ -894,7 +895,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
         """Check if two DiscreteScalarFieldViews container are equivalent."""
         check_instance(other, DiscreteTensorField)
         assert npw.array_equal(self.shape, other.shape)
-        eq = all(dfield.match(other[idx], invert=False) 
+        eq = all(dfield.match(other[idx], invert=False)
                             for (idx, dfield) in self.nd_iter())
         eq &= npw.array_equal(self.shape, other.shape)
         eq &= (self._name == other._name)
@@ -921,7 +922,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
             dfields[idx] = dfield.view(topology_state=topology_state)
         dfield = self.from_dfield_array(name=name, pretty_name=pretty_name, dfields=dfields)
         return dfield
-    
+
     def as_any_dfield(self, memory_order, name=None, pretty_name=None, **kwds):
         """
         Quickly take a view on this DiscreteScalarFieldView using self topology state
@@ -955,7 +956,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
     def short_description(self):
         """Short description of this discrete field container."""
         s  = '{}[name={}, pname={}, shape={}, nb_components={}]'
-        s  = s.format(self.full_tag, self.name, self.pretty_name, 
+        s  = s.format(self.full_tag, self.name, self.pretty_name,
                       self.shape, self.nb_components)
         return s
 
@@ -990,14 +991,14 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
             integrals[idx] = dfield.integrate(scale=scale, **kwds).tolist()[0]
         return integrals
 
-    def exchange_ghosts(self, build_exchanger=False, build_launcher=False, 
+    def exchange_ghosts(self, build_exchanger=False, build_launcher=False,
             evt=None, **kwds):
         """
         Exchange ghosts using cached ghost exchangers which are built at first use.
         ie. Exchange every ghosts components of self.data using current topology state.
         By default all ghosts are exchanged.
         """
-        if (build_launcher or build_exchanger): 
+        if (build_launcher or build_exchanger):
             assert (evt is None), 'Cannot spevify event while building a launcher.'
             from hysop.fields.ghost_exchangers import MultiGhostExchanger
             ghost_exchangers = MultiGhostExchanger(name='{}_ghost_exchange'.format(self.name))
@@ -1023,7 +1024,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
         ie. Exchange every ghosts components of self.data using current topology state.
         Specialization for ghost summation.
         """
-        if (build_launcher or build_exchanger): 
+        if (build_launcher or build_exchanger):
             assert (evt is None), 'Cannot spevify event while building a launcher.'
             from hysop.fields.ghost_exchangers import MultiGhostExchanger
             ghost_exchangers = MultiGhostExchanger(name='{}_ghost_exchange'.format(self.name))
@@ -1045,7 +1046,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
 
     def build_ghost_exchanger(self, **kwds):
         """
-        Build a ghost exchanger, possibly on different data. 
+        Build a ghost exchanger, possibly on different data.
         Usefull for operator apply.
         """
         from hysop.fields.ghost_exchangers import MultiGhostExchanger
@@ -1075,7 +1076,7 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
             dfield._set_data(data)
     def _get_buffers(self):
         return tuple(dfield.sbuffer for dfield in self.discrete_field_views())
-    
+
     def _get_sdata(self):
         self._raise_sdata()
     def _set_sdata(self, copy_data):
@@ -1088,10 +1089,10 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
         raise RuntimeError(msg)
 
     data     = property(_get_data, _set_data)
-    buffers  = property(_get_buffers) 
+    buffers  = property(_get_buffers)
     sdata    = property(_get_sdata, _set_sdata)
-    sbuffer  = property(_get_sbuffer) 
-    
+    sbuffer  = property(_get_sbuffer)
+
     def initialize(self, *args, **kwds):
         msg='This method is only defined for specific topologies '
         msg+='(see CartesianDiscreteTensorField).'
@@ -1109,4 +1110,3 @@ class DiscreteTensorField(NamedTensorContainerI, DiscreteScalarFieldViewContaine
 
 DiscreteField = (DiscreteScalarField, DiscreteTensorField)
 """A DiscreteField is either of DiscreteScalarField or a DiscreteTensorField"""
-
diff --git a/hysop/fields/tests/test_cartesian.py b/hysop/fields/tests/test_cartesian.py
index 3b8a5ed96122e4ebaca67f72b2d885430b955830..499b37484ce0bd0b85e7000a9a41528473d6e7c2 100644
--- a/hysop/fields/tests/test_cartesian.py
+++ b/hysop/fields/tests/test_cartesian.py
@@ -41,7 +41,7 @@ def test_serial_initialization_1d():
     dim = 1
     npts = (10,)
     nghosts = (2,)
-    
+
     for (lbd, rbd) in domain_boundary_iterator(dim):
         domain = Box(dim=dim, lboundaries=lbd,
                               rboundaries=rbd)
@@ -49,8 +49,8 @@ def test_serial_initialization_1d():
         F1 = Field(domain=domain, name='F1', nb_components=2)
         F2 = Field(domain=domain, name='F2', shape=(2,2))
         print '[{}]'.format(F0.format_boundaries())
-        
-        discretization = CartesianDiscretization(npts, nghosts, 
+
+        discretization = CartesianDiscretization(npts, nghosts,
                                             lboundaries=F0.lboundaries,
                                             rboundaries=F0.rboundaries)
 
@@ -124,8 +124,8 @@ def test_serial_initialization_2d():
         F0 = Field(domain=domain, name='F0', nb_components=1)
         F1 = Field(domain=domain, name='F1', nb_components=2)
         print '[{}]'.format(F0.format_boundaries())
-        
-        discretization = CartesianDiscretization(npts, nghosts, 
+
+        discretization = CartesianDiscretization(npts, nghosts,
                                             lboundaries=F0.lboundaries,
                                             rboundaries=F0.rboundaries)
 
@@ -147,7 +147,7 @@ def test_serial_initialization_2d():
             dfields = dF0.dfields + dF1.dfields
             assert all(np.all(d.local_lboundaries == F0.lboundaries) for d in dfields)
             assert all(np.all(d.local_rboundaries == F0.rboundaries) for d in dfields)
-            
+
             Ny,Nx = npts
             Gy,Gx = nghosts
             Ly,Lx = tuple(discretization.lboundaries)
@@ -212,7 +212,7 @@ def test_serial_initialization_2d():
                             assert b[Fy,Fx].shape == (Ny-2*Gy, Nx-2*Gx)
                         else:
                             raise NotImplementedError(ghost_mask)
-                        
+
                         if (Lx==BoundaryCondition.PERIODIC):
                             assert np.all(b[Fy,Xo[0]] == b[Fy,Xi[2]]), '\n'+str(d)
                         elif (Lx==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
@@ -222,7 +222,7 @@ def test_serial_initialization_2d():
                             assert (b[Fy,:Gx] == +b[Fy,Gx+1:2*Gx+1][Ix]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Lx))
-                        
+
                         if (Rx==BoundaryCondition.PERIODIC):
                             assert np.all(b[Fy,Xo[2]] == b[Fy,Xi[0]]), '\n'+str(d)
                         elif (Rx==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
@@ -232,7 +232,7 @@ def test_serial_initialization_2d():
                             assert (b[Fy,Nx-1:Nx+Gx-1] == +b[Fy,Nx+Gx:][Ix]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Rx))
-                        
+
                         if (Ly==BoundaryCondition.PERIODIC):
                             assert np.all(b[Yo[0],Fx] == b[Yi[2],Fx]), '\n'+str(d)
                         elif (Ly==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
@@ -242,7 +242,7 @@ def test_serial_initialization_2d():
                             assert (b[:Gy,Fx] == +b[Gy+1:2*Gy+1,Fx][Iy]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Ly))
-                        
+
                         if (Ry==BoundaryCondition.PERIODIC):
                             assert np.all(b[Yo[2],Fx] == b[Yi[0],Fx]), '\n'+str(d)
                         elif (Ry==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
@@ -285,8 +285,8 @@ def test_serial_initialization_3d():
         F0 = Field(domain=domain, name='F0', nb_components=1)
         F1 = Field(domain=domain, name='F1', nb_components=3)
         print '[{}]'.format(F0.format_boundaries())
-        
-        discretization = CartesianDiscretization(npts, nghosts, 
+
+        discretization = CartesianDiscretization(npts, nghosts,
                                             lboundaries=F0.lboundaries,
                                             rboundaries=F0.rboundaries)
 
@@ -308,7 +308,7 @@ def test_serial_initialization_3d():
             dfields = dF0.dfields + dF1.dfields
             assert all(np.all(d.local_lboundaries == F0.lboundaries) for d in dfields)
             assert all(np.all(d.local_rboundaries == F0.rboundaries) for d in dfields)
-            
+
             Nz,Ny,Nx = npts
             Gz,Gy,Gx = nghosts
             Lz,Ly,Lx = tuple(discretization.lboundaries)
@@ -389,7 +389,7 @@ def test_serial_initialization_3d():
                             assert b[Fz,Fy,Fx].shape == (Nz-2*Gz,Ny-2*Gy, Nx-2*Gx)
                         else:
                             raise NotImplementedError(ghost_mask)
-                        
+
                         if (Lx==BoundaryCondition.PERIODIC):
                             assert np.all(b[Fz,Fy,Xo[0]] == b[Fz,Fy,Xi[2]]), '\n'+str(d)
                         elif (Lx==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
@@ -399,7 +399,7 @@ def test_serial_initialization_3d():
                             assert (b[Fz,Fy,:Gx] == +b[Fz,Fy,Gx+1:2*Gx+1][Ix]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Lx))
-                        
+
                         if (Rx==BoundaryCondition.PERIODIC):
                             assert np.all(b[Fz,Fy,Xo[2]] == b[Fz,Fy,Xi[0]]), '\n'+str(d)
                         elif (Rx==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
@@ -409,7 +409,7 @@ def test_serial_initialization_3d():
                             assert (b[Fz,Fy,Nx-1:Nx+Gx-1] == +b[Fz,Fy,Nx+Gx:][Ix]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Rx))
-                        
+
                         if (Ly==BoundaryCondition.PERIODIC):
                             assert np.all(b[Fz,Yo[0],Fx] == b[Fz,Yi[2],Fx]), '\n'+str(d)
                         elif (Ly==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
@@ -419,7 +419,7 @@ def test_serial_initialization_3d():
                             assert (b[Fz,:Gy,Fx] == +b[Fz,Gy+1:2*Gy+1,Fx][Iy]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Ly))
-                        
+
                         if (Ry==BoundaryCondition.PERIODIC):
                             assert np.all(b[Fz,Yo[2],Fx] == b[Fz,Yi[0],Fx]), '\n'+str(d)
                         elif (Ry==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
@@ -429,7 +429,7 @@ def test_serial_initialization_3d():
                             assert (b[Fz,Ny-1:Ny+Gy-1,Fx] == +b[Fz,Ny+Gy:,Fx][Iy]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Ry))
-                        
+
                         if (Lz==BoundaryCondition.PERIODIC):
                             assert np.all(b[Zo[0],Fy,Fx] == b[Zi[2],Fy,Fx]), '\n'+str(d)
                         elif (Lz==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
@@ -439,7 +439,7 @@ def test_serial_initialization_3d():
                             assert (b[:Gz,Fy,Fx] == +b[Gz+1:2*Gz+1,Fy,Fx][Iz]).all(), '\n'+str(d)
                         else:
                             raise NotImplementedError('Unknown boundary condition {}.'.format(Lz))
-                        
+
                         if (Rz==BoundaryCondition.PERIODIC):
                             assert np.all(b[Zo[2],Fy,Fx] == b[Zi[0],Fy,Fx]), '\n'+str(d)
                         elif (Rz==BoundaryCondition.HOMOGENEOUS_DIRICHLET):
@@ -483,7 +483,10 @@ def iter_backends():
     for cl_env in iter_clenv():
         yield (Backend.OPENCL, cl_env)
 
-def test_mpi_ghost_exchange_periodic(comm):
+def test_mpi_ghost_exchange_periodic(comm=None):
+    if comm is None:
+        from mpi4py import MPI
+        comm = MPI.COMM_WORLD
     rank = comm.Get_rank()
     size = comm.Get_size()
     dtypes = (np.float32, np.float32, np.float64,
@@ -582,19 +585,22 @@ def test_mpi_ghost_exchange_periodic(comm):
                                         ldata = np.atleast_1d(ldata.get())
                                         target_vals = ghost_vals(ldata.shape, dtype, i, d,
                                                 left_rank,1)
-                                        assert np.allclose(ldata, target_vals), (rank, 
+                                        assert np.allclose(ldata, target_vals), (rank,
                                                                                     target_vals)
                                         rdata = np.atleast_1d(rdata.get())
-                                        target_vals = ghost_vals(rdata.shape, dtype, i, d, 
+                                        target_vals = ghost_vals(rdata.shape, dtype, i, d,
                                                                                 right_rank, 0)
-                                        assert np.allclose(rdata, target_vals), (rank, 
+                                        assert np.allclose(rdata, target_vals), (rank,
                                                                                         target_vals)
                             if rank==0:
                                 print
 
 
-def test_mpi_ghost_exchange_runtime(comm):
+def test_mpi_ghost_exchange_runtime(comm=None):
     """Just bruteforce all exchange possibilities to see if nothing crashes in 1D and 2D."""
+    if comm is None:
+        from mpi4py import MPI
+        comm = MPI.COMM_WORLD
     rank = comm.Get_rank()
     size = comm.Get_size()
     dtype = np.float32
@@ -614,13 +620,13 @@ def test_mpi_ghost_exchange_runtime(comm):
 
         npts = (17,16,19)[:dim]
         nghosts = (2,1,3)[:dim]
-        
+
         for shape in it.product(xrange(0,size+1), repeat=dim):
             if np.prod(shape, dtype=np.uint32)!=size:
                 continue
             if rank==0:
                 sys.stdout.write('  >CART SHAPE: {}\n'.format(shape))
-            
+
             for (backend, cl_env) in iter_backends():
                 if rank==0:
                     sys.stdout.write('    >BACKEND.{:<7} '.format(str(backend)+':'))
@@ -631,7 +637,7 @@ def test_mpi_ghost_exchange_runtime(comm):
                         sys.stdout.flush()
                         return True
                     return False
-                
+
                 i=0
                 brk = False
                 try:
@@ -639,10 +645,10 @@ def test_mpi_ghost_exchange_runtime(comm):
                         domain = Box(dim=dim, lboundaries=lbd,
                                               rboundaries=rbd)
 
-                        F = Field(domain=domain, name='F', nb_components=1, 
+                        F = Field(domain=domain, name='F', nb_components=1,
                                     dtype=dtype, _register=False)
 
-                        discretization = CartesianDiscretization(npts, nghosts, 
+                        discretization = CartesianDiscretization(npts, nghosts,
                                                 lboundaries=F.lboundaries,
                                                 rboundaries=F.rboundaries)
 
@@ -673,7 +679,10 @@ def test_mpi_ghost_exchange_runtime(comm):
                         sys.stdout.flush()
 
 
-def test_mpi_ghost_accumulate_periodic(comm):
+def test_mpi_ghost_accumulate_periodic(comm=None):
+    if comm is None:
+        from mpi4py import MPI
+        comm = MPI.COMM_WORLD
     rank = comm.Get_rank()
     size = comm.Get_size()
     if rank==0:
@@ -780,7 +789,7 @@ def test_mpi_ghost_accumulate_periodic(comm):
                                                 assert (ishape is not None)
                                                 data[oview] = ghost_vals(oshape, dtype, rank, directions,
                                                                                             displacements, i)
-                                    
+
                                     dF.accumulate_ghosts(directions=directions,
                                                          exchange_method=exchange_method)
 
@@ -836,9 +845,8 @@ if __name__ == '__main__':
             test_serial_initialization_2d()
             if __ENABLE_LONG_TESTS__:
                 test_serial_initialization_3d()
-        
+
         disable_hysop_warnings()
         test_mpi_ghost_exchange_runtime(comm=comm)
         test_mpi_ghost_exchange_periodic(comm=comm)
         test_mpi_ghost_accumulate_periodic(comm=comm)
-
diff --git a/hysop/fields/tests/test_fields.py b/hysop/fields/tests/test_fields.py
index 86314bb5b2fee90d41c6429238b8044d1566a7e9..e4bb0f94d412a516a24cf319c9a5519773e19b5f 100644
--- a/hysop/fields/tests/test_fields.py
+++ b/hysop/fields/tests/test_fields.py
@@ -34,10 +34,10 @@ def test_field():
     DF5 = F5.discretize(T0)
     DF6 = F6.discretize(T0)
     DF7 = F7.discretize(T0)
-    
-    requests  = DF5._dfield.memory_request(op=DF5.name, 
+
+    requests  = DF5._dfield.memory_request(op=DF5.name,
                     request_identifier=DF5.memory_request_id)
-    requests += DF7._dfield.memory_request(op=DF7.name, 
+    requests += DF7._dfield.memory_request(op=DF7.name,
                     request_identifier=DF7.memory_request_id)
     work = requests.allocate(True)
     DF5.honor_memory_request(work)
@@ -89,19 +89,19 @@ def test_field():
     assert F4.dtype == np.uint16
     assert F4.name == 'F4'
     assert F4.pretty_name == 'F4'
-    
+
     def func(data, coords):
         for (d, coord) in zip(data, coords):
             x,y,z = coord
             d[...] = x*y*z
-    
+
     DF7.initialize(func)
     DF7.initialize(DF5.data)
     DF7.fill(4)
     DF7.randomize()
     DF7.copy(DF0)
     DF7.integrate()
-    
+
     DF7.has_ghosts()
     DF7.exchange_ghosts()
     DF7.accumulate_ghosts()
@@ -123,9 +123,10 @@ def test_tensor_field():
     T0 = TensorField(domain, 'T0', (3,3))
     T0_bis = Field(domain, 'T0b', shape=(3,3))
     T1 = T0.field_like('T1', dtype=np.float16)
-    T2 = T0[1:,1:].rename(name='T2')
+    T2 = T0[1:,1:]
+    T2.rename(name='T2')
     T3 = T0[1,1]
-    T4 = TensorField.from_fields(name='T4', 
+    T4 = TensorField.from_fields(name='T4',
             fields=(T0[0,0], T0[1,1], T1[0,1], T1[1,0]), shape=(2,2))
     T5 = TensorField.from_field_array('T5', T0._fields)
     T6 = T5.field_like(name='T6')
@@ -178,10 +179,10 @@ def test_tensor_field():
     DT6 = T6.discretize(topo)
     DT7 = T7.discretize(topo)
 
-    DT8 = DiscreteTensorField.from_dfields(name='DT8', 
-            dfields=(DT0.dfields[0], 
-                     DT1.dfields[0], 
-                     DT2.dfields[0], 
+    DT8 = DiscreteTensorField.from_dfields(name='DT8',
+            dfields=(DT0.dfields[0],
+                     DT1.dfields[0],
+                     DT2.dfields[0],
                      DT3.dfields[0]),
             shape=(2,2))
     DT9 = DT0[:2,:2]
@@ -207,7 +208,7 @@ def test_tensor_field():
     DT9.fill(4)
     DT9.randomize()
     DT9.copy(DT0[1:,1:])
-    
+
     DT9.has_ghosts()
     DT9.exchange_ghosts()
     DT9.accumulate_ghosts()
@@ -329,8 +330,8 @@ def test_boundaries():
     """This test checks that all boundaries are compatible for velocity and vorticity."""
     for dim in (1,2,):
         i=0
-        for (lbd, rbd) in domain_boundary_iterator(dim): 
-            domain = Box(dim=dim, lboundaries=lbd, 
+        for (lbd, rbd) in domain_boundary_iterator(dim):
+            domain = Box(dim=dim, lboundaries=lbd,
                                   rboundaries=rbd)
             V = VelocityField(domain)
             S = ScalarField(name='S0', domain=domain)
@@ -358,7 +359,7 @@ def test_boundaries():
                 print '{} (VORTICITY) BOUNDARIES:'.format(rotV.pretty_name)
                 for Wi in rotV.fields:
                     print ' *{} boundaries=[{}]'.format(Wi.pretty_name, Wi.format_boundaries())
-    
+
 
 
 if __name__ == '__main__':
diff --git a/hysop/numerics/splitting/test/test_strang.py b/hysop/numerics/splitting/test/test_strang.py
index fb86d81c959aaedfb75a9050fbbeb90ba42c93fd..ac6048a007526bf67185159abb49987678dbf953 100644
--- a/hysop/numerics/splitting/test/test_strang.py
+++ b/hysop/numerics/splitting/test/test_strang.py
@@ -1,6 +1,5 @@
-
 from hysop.problem          import Problem
-from hysop.operators        import DirectionalAdvection, Poisson
+from hysop.operators        import DirectionalAdvection, PoissonCurl
 from hysop.numerics.splitting.strang import StrangSplitting
 
 from hysop.domain.box         import Box
@@ -10,8 +9,6 @@ from hysop.tools.parameters   import CartesianDiscretization
 from hysop.simulation import Simulation
 
 from hysop.methods import *
-from hysop.backend.device.kernel_autotuner import AutotunerConfig, AutotunerFlags
-# from hysop.constants import callback_profiler 
 
 class TestStrang(object):
 
@@ -22,7 +19,7 @@ class TestStrang(object):
     @classmethod
     def teardown_class(cls):
         pass
-    
+
     def setup_method(self, method):
         pass
 
@@ -49,43 +46,45 @@ class TestStrang(object):
 
     def _test_strang_nd(self, dim, n):
         ## Simulation
-        simu = Simulation(start=0.0, end=1.0, max_iter=100, time_step=0.1)
+        dt0=0.1
+        simu = Simulation(start=0.0, end=1.0, max_iter=100, dt0=dt0)
         simu.initialize()
-        
+
         ## Domain and continuous fields
         box, velo, vorti, scalars = self.make_fields(dim,2)
 
         ## CartesianDiscretizations and topologies
         resolution     = (n,)*dim
-        dnd            = CartesianDiscretization(resolution=resolution, ghosts=None)
-        topo           = CartesianTopology(box,dnd,dim)
+        dnd            = CartesianDiscretization(resolution=resolution, ghosts=None, default_boundaries=True)
+        topo           = CartesianTopology(domain=box,discretization=dnd)
 
         ## Operators
         advec = DirectionalAdvection(
-                name='advec',
-                velocity = velo,       
-                advected_fields = (vorti,)+scalars,
-                variables = { velo: topo, vorti: topo,
-                              scalars[0]: topo, scalars[1]: topo},
-                split_fields=(vorti, scalars)
+            name='advec',
+            dt = simu.dt,
+            velocity = velo,velocity_cfl=1.*dt0*n,
+            advected_fields = (vorti,)+scalars,
+            variables = { velo: topo, vorti: topo,
+                          scalars[0]: topo, scalars[1]: topo},
             )
 
-        poisson = Poisson(name='poisson', 
-                velocity=velo, vorticity=vorti, 
-                topology=topo)
+        poisson = PoissonCurl(name='poisson',
+                              velocity=velo, vorticity=vorti,
+                              variables={velo: topo, vorti: topo})
 
         ## Graph
-        splitting = StrangSplitting(dim=dim, order=StrangOrder.STRANG_FIRST_ORDER)
-        splitting.push_operator(advec)
+        splitting = StrangSplitting(splitting_dim=dim,
+                                    order=StrangOrder.STRANG_FIRST_ORDER)
+        splitting.push_operators(advec)
 
         problem = Problem()
         problem.insert(splitting)
         problem.insert(poisson)
         problem.build()
-        
+
         problem.display()
         problem.finalize()
-    
+
     def test_strang_2d(self, n=33):
         self._test_strang_nd(dim=2, n=n)
     def test_strang_3d(self, n=33):
diff --git a/hysop/numerics/tests/test_fft.py b/hysop/numerics/tests/test_fft.py
index 4b42e4e001d74ed1e7f9c4d6631e9ea3ea768f9a..58f4112b21ea92ac1bf0900a2e5fcaa086cc7b6c 100644
--- a/hysop/numerics/tests/test_fft.py
+++ b/hysop/numerics/tests/test_fft.py
@@ -26,14 +26,14 @@ class TestFFT(object):
 
     implementations = {
         Implementation.PYTHON: {
-            'numpy': NumpyFFT(warn_on_allocation=False), 
-            'scipy': ScipyFFT(warn_on_allocation=False), 
+            'numpy': NumpyFFT(warn_on_allocation=False),
+            'scipy': ScipyFFT(warn_on_allocation=False),
             'fftw':  FftwFFT(warn_on_allocation=False,
                              warn_on_misalignment=False)
         },
         Implementation.OPENCL: {}
     }
-    
+
     print
     print ':: STARTING FFT BACKEND TESTS ::'
     for (i,cl_env) in enumerate(iter_clenv()):
@@ -59,7 +59,7 @@ class TestFFT(object):
         print '::Testing 1D transform, precision {}::'.format(dtype.__name__)
         eps = np.finfo(dtype).eps
         ctype = float_to_complex_dtype(dtype)
-    
+
         def check_distances(results, eps, report_eps, tag, failures):
             if len(results.keys())==0:
                 print 'no support'
@@ -111,9 +111,9 @@ class TestFFT(object):
                         Bin = mk_buffer(backend=impl.backend, shape=shape, dtype=ctype)
                         plan = impl.fft(a=Bin).setup()
                         Bout = plan.output_array
-                        assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape, 
+                        assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape,
                                 name)
-                        assert Bout.dtype == ctype, self.msg_dtype.format(ctype, Bout.dtype, 
+                        assert Bout.dtype == ctype, self.msg_dtype.format(ctype, Bout.dtype,
                                 name)
                         Bin[...] = Href
                         plan.execute()
@@ -130,9 +130,9 @@ class TestFFT(object):
                     except HysopFFTDataLayoutError as e:
                         results[name] = e
             check_distances(results, eps, self.report_eps, 'forward C2C', failures)
-        
+
         print '\n BACKWARD C2C: complex to complex backward transform'
-        for (shape, cshape, rshape, N, Nc, Nr, 
+        for (shape, cshape, rshape, N, Nc, Nr,
                 ghosts, mk_buffer) in self.iter_shapes():
             print '   IFFT shape={:12s} ghosts={:12s} '.format(shape, str(ghosts)+':'),
             Href = np.random.rand(2*N).astype(dtype).view(dtype=ctype).reshape(shape)
@@ -145,9 +145,9 @@ class TestFFT(object):
                         Bin = mk_buffer(backend=impl.backend, shape=shape, dtype=ctype)
                         plan = impl.ifft(a=Bin).setup()
                         Bout = plan.output_array
-                        assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape, 
+                        assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape,
                                 name)
-                        assert Bout.dtype == ctype, self.msg_dtype.format(ctype, Bout.dtype, 
+                        assert Bout.dtype == ctype, self.msg_dtype.format(ctype, Bout.dtype,
                                 name)
                         Bin[...] = Href
                         plan.execute()
@@ -164,9 +164,9 @@ class TestFFT(object):
                     except HysopFFTDataLayoutError as e:
                         results[name] = e
             check_distances(results, eps, self.report_eps, 'backward C2C', failures)
-            
+
         print '\n FORWARD R2C: real to hermitian complex transform'
-        for (shape, cshape, rshape, N, Nc, Nr, 
+        for (shape, cshape, rshape, N, Nc, Nr,
                 ghosts, mk_buffer) in self.iter_shapes():
             print '   RFFT shape={:12s} ghosts={:12s} '.format(shape, str(ghosts)+':'),
             Href = np.random.rand(*shape).astype(dtype).reshape(shape)
@@ -179,9 +179,9 @@ class TestFFT(object):
                         Bin = mk_buffer(backend=impl.backend, shape=shape, dtype=dtype)
                         plan = impl.rfft(a=Bin).setup()
                         Bout = plan.output_array
-                        assert Bout.shape == cshape, self.msg_shape.format(cshape, Bout.shape, 
+                        assert Bout.shape == cshape, self.msg_shape.format(cshape, Bout.shape,
                                 name)
-                        assert Bout.dtype == ctype, self.msg_dtype.format(ctype, Bout.dtype, 
+                        assert Bout.dtype == ctype, self.msg_dtype.format(ctype, Bout.dtype,
                                 name)
                         Bin[...] = Href
                         plan.execute()
@@ -198,9 +198,9 @@ class TestFFT(object):
                     except HysopFFTDataLayoutError as e:
                         results[name] = e
             check_distances(results, eps, self.report_eps, 'R2C', failures)
-        
+
         print '\n BACKWARD C2R: real to hermitian complex transform'
-        for (shape, cshape, rshape, N, Nc, Nr, 
+        for (shape, cshape, rshape, N, Nc, Nr,
                 ghosts, mk_buffer) in self.iter_shapes():
             print '   IRFFT shape={:12s} ghosts={:12s} '.format(shape, str(ghosts)+':'),
             Href = np.random.rand(2*Nc).astype(dtype).view(dtype=ctype).reshape(cshape)
@@ -213,9 +213,9 @@ class TestFFT(object):
                         Bin = mk_buffer(backend=impl.backend, shape=cshape, dtype=ctype)
                         plan = impl.irfft(a=Bin).setup()
                         Bout = plan.output_array
-                        assert Bout.shape == rshape, self.msg_shape.format(rshape, Bout.shape, 
+                        assert Bout.shape == rshape, self.msg_shape.format(rshape, Bout.shape,
                                 name)
-                        assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype, 
+                        assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype,
                                 name)
                         Bin[...] = Href
                         plan.execute()
@@ -235,7 +235,7 @@ class TestFFT(object):
 
         print ('\n BACKWARD FORCED C2R: real to hermitian complex transform with specified '
                 +'shape')
-        for (shape, cshape, rshape, N, Nc, Nr, 
+        for (shape, cshape, rshape, N, Nc, Nr,
                 ghosts, mk_buffer) in self.iter_shapes():
             print '   IRFFT shape={:12s} ghosts={:12s} '.format(shape, str(ghosts)+':'),
             Href = np.random.rand(2*Nc).astype(dtype).view(dtype=ctype).reshape(cshape)
@@ -248,9 +248,9 @@ class TestFFT(object):
                         Bin = mk_buffer(backend=impl.backend, shape=cshape, dtype=ctype)
                         plan = impl.irfft(a=Bin, n=shape[-1]).setup()
                         Bout = plan.output_array
-                        assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape, 
+                        assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape,
                                 name)
-                        assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype, 
+                        assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype,
                                 name)
                         Bin[...] = Href
                         plan.execute()
@@ -267,12 +267,12 @@ class TestFFT(object):
                     except HysopFFTDataLayoutError as e:
                         results[name] = e
             check_distances(results, eps, self.report_eps, 'forced C2R', failures)
-    
+
         types = ['I  ','II ','III','IV ']
         for (itype,stype) in enumerate(types, 1):
             print '\n DCT-{}: real to real discrete cosine transform {}'.format(
                     stype.strip(), itype)
-            for (shape, cshape, rshape, N, Nc, Nr, 
+            for (shape, cshape, rshape, N, Nc, Nr,
                     ghosts, mk_buffer) in self.iter_shapes():
                 print '   DCT-{} shape={:12s} ghosts={:12s} '.format(stype, shape, str(ghosts)+':'),
                 if (itype==1): # real size is 2*(N-1)
@@ -289,9 +289,9 @@ class TestFFT(object):
                             Bin = mk_buffer(backend=impl.backend, shape=shape, dtype=dtype)
                             plan = impl.dct(a=Bin, type=itype).setup()
                             Bout = plan.output_array
-                            assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape, 
+                            assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape,
                                     name)
-                            assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype, 
+                            assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype,
                                     name)
                             Bin[...] = Href
                             plan.execute()
@@ -307,14 +307,14 @@ class TestFFT(object):
                             results[name] = H1 / shape[-1]
                         except HysopFFTDataLayoutError as e:
                             results[name] = e
-                check_distances(results, eps, self.report_eps, 
+                check_distances(results, eps, self.report_eps,
                         'DCT-{}'.format(stype), failures)
-        
+
         for (itype,stype) in enumerate(types, 1):
             iitype = [1,3,2,4][itype-1]
             print '\n IDCT-{}: real to real inverse discrete cosine transform {}'.format(
                     stype.strip(), itype)
-            for (shape, cshape, rshape, N, Nc, Nr, 
+            for (shape, cshape, rshape, N, Nc, Nr,
                     ghosts, mk_buffer) in self.iter_shapes():
                 print '   IDCT-{} shape={:12s} ghosts={:12s} '.format(stype, shape, str(ghosts)+':'),
                 if (iitype==1): # real size is 2*(N-1)
@@ -331,9 +331,9 @@ class TestFFT(object):
                             Bin = mk_buffer(backend=impl.backend, shape=shape, dtype=dtype)
                             plan = impl.idct(a=Bin, type=itype).setup()
                             Bout = plan.output_array
-                            assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape, 
+                            assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape,
                                     name)
-                            assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype, 
+                            assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype,
                                     name)
                             Bin[...] = Href
                             plan.execute()
@@ -349,14 +349,14 @@ class TestFFT(object):
                             results[name] = H1
                         except HysopFFTDataLayoutError as e:
                             results[name] = e
-                check_distances(results, eps, self.report_eps, 
+                check_distances(results, eps, self.report_eps,
                         'IDCT-{}'.format(stype), failures)
 
         types = ['I  ','II ','III','IV ']
         for (itype,stype) in enumerate(types, 1):
             print '\n DST-{}: real to real discrete sine transform {}'.format(
                     stype.strip(), itype)
-            for (shape, cshape, rshape, N, Nc, Nr, 
+            for (shape, cshape, rshape, N, Nc, Nr,
                     ghosts, mk_buffer) in self.iter_shapes():
                 print '   DST-{} shape={:12s} ghosts={:12s} '.format(stype, shape, str(ghosts)+':'),
                 if (itype==1): # real size will be 2*(N+1)
@@ -373,9 +373,9 @@ class TestFFT(object):
                             Bin = mk_buffer(backend=impl.backend, shape=shape, dtype=dtype)
                             plan = impl.dst(a=Bin, type=itype).setup()
                             Bout = plan.output_array
-                            assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape, 
+                            assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape,
                                     name)
-                            assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype, 
+                            assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype,
                                     name)
                             Bin[...] = Href
                             plan.execute()
@@ -391,14 +391,14 @@ class TestFFT(object):
                             results[name] = H1 / shape[-1]
                         except HysopFFTDataLayoutError as e:
                             results[name] = e
-                check_distances(results, eps, self.report_eps, 
+                check_distances(results, eps, self.report_eps,
                         'DST-{}'.format(stype), failures)
-        
+
         for (itype,stype) in enumerate(types, 1):
             iitype = [1,3,2,4][itype-1]
             print '\n IDST-{}: real to real inverse discrete sine transform {}'.format(
                     stype.strip(), itype)
-            for (shape, cshape, rshape, N, Nc, Nr, 
+            for (shape, cshape, rshape, N, Nc, Nr,
                     ghosts, mk_buffer) in self.iter_shapes():
                 print '   IDST-{} shape={:12s} ghosts={:12s} '.format(stype, shape, str(ghosts)+':'),
                 if (iitype==1): # real size will be 2*(N+1)
@@ -415,9 +415,9 @@ class TestFFT(object):
                             Bin = mk_buffer(backend=impl.backend, shape=shape, dtype=dtype)
                             plan = impl.idst(a=Bin, type=itype).setup()
                             Bout = plan.output_array
-                            assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape, 
+                            assert Bout.shape == shape, self.msg_shape.format(shape, Bout.shape,
                                     name)
-                            assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype, 
+                            assert Bout.dtype == dtype, self.msg_dtype.format(dtype, Bout.dtype,
                                     name)
                             Bin[...] = Href
                             plan.execute()
@@ -433,9 +433,9 @@ class TestFFT(object):
                             results[name] = H1
                         except HysopFFTDataLayoutError as e:
                             results[name] = e
-                check_distances(results, eps, self.report_eps, 
+                check_distances(results, eps, self.report_eps,
                         'IDST-{}'.format(stype), failures)
-   
+
 
 
 
@@ -467,7 +467,7 @@ class TestFFT(object):
                 raise RuntimeError(msg)
 
         print '\n C2C-C2C transform'
-        for (shape, cshape, rshape, N, Nc, Nr, 
+        for (shape, cshape, rshape, N, Nc, Nr,
                 ghosts, mk_buffer) in self.iter_shapes():
             print '   X - IFFT(FFT(X)) shape={:12s} ghosts={:12s}'.format(shape, str(ghosts)+':'),
             Href = np.random.rand(2*N).astype(dtype).view(dtype=ctype).reshape(shape)
@@ -525,7 +525,7 @@ class TestFFT(object):
                     except HysopFFTDataLayoutError as e:
                         results[name] = e
             check_distances(results)
-        
+
         print '\n R2R-R2R transforms'
 
         types = ['I  ','II ','III','IV ']
@@ -568,7 +568,7 @@ class TestFFT(object):
                         except HysopFFTDataLayoutError as e:
                             results[name] = e
                 check_distances(results)
-        
+
         for (itype,stype) in enumerate(types, 1):
             print '\n DST-{}: real to real discrete sinine transform {}'.format(
                     stype.strip(), itype)
@@ -608,7 +608,7 @@ class TestFFT(object):
                         except HysopFFTDataLayoutError as e:
                             results[name] = e
                 check_distances(results)
-    
+
     def iter_shapes(self):
         minj=(3,3)
         maxj=(6,6)
@@ -643,7 +643,7 @@ class TestFFT(object):
                         assert buf.dtype == dtype
                         return buf
                     yield (shape, cshape, rshape, N, Nc, Nr, ghosts, mk_buffer)
-    
+
     def report_failures(self, failures):
         print
         print '== TEST FAILURES REPORT =='
@@ -682,7 +682,7 @@ class TestFFT(object):
             msg+='\n***************************************************'
             msg+='\n'
             print msg
-    
+
 
     def perform_tests(self):
         # not testing np.longdouble because only one implementation supports it
@@ -696,11 +696,18 @@ class TestFFT(object):
             self._test_forward_backward_1d(dtype=dtype)
             self._test_1d(dtype=dtype, failures=failures.setdefault(dtype.__name__, {}))
         self.report_failures(failures)
-    
+
+    def test_1d(self):
+        failures = {}
+        self._test_forward_backward_1d(dtype=HYSOP_REAL)
+        self._test_1d(dtype=HYSOP_REAL, failures=failures.setdefault(HYSOP_REAL.__name__, {}))
+        self.report_failures(failures)
+
+
 if __name__ == '__main__':
     test = TestFFT()
 
-    with printoptions(threshold=10000, linewidth=240, 
-                      nanstr='nan', infstr='inf', 
+    with printoptions(threshold=10000, linewidth=240,
+                      nanstr='nan', infstr='inf',
                       formatter={'float': lambda x: '{:>0.2f}'.format(x)}):
         test.perform_tests()
diff --git a/hysop/operator/base/custom_symbolic_operator.py b/hysop/operator/base/custom_symbolic_operator.py
index 913530cd267ad84b2cf5a6d984a3cd31abdcc830..2c5562461b370a4a26de37399204e529ee68b2e0 100644
--- a/hysop/operator/base/custom_symbolic_operator.py
+++ b/hysop/operator/base/custom_symbolic_operator.py
@@ -344,12 +344,8 @@ class SymbolicExpressionInfo(object):
         else:
             msg='No discrete fields found in custom symbolic operator.'
             raise RuntimeError(msg)
-<<<<<<< HEAD
-
-=======
         self.axes = axes
 
->>>>>>> master
         SymbolicExpressionParser.discretize_expressions(self)
 
         self.check_dfield_sizes()
diff --git a/hysop/operator/tests/test_analytic.py b/hysop/operator/tests/test_analytic.py
index 7822a30cab1283a89b42c740950dfd21effe457f..47b66cd7e059252f98f42af617768e832b8cc876 100644
--- a/hysop/operator/tests/test_analytic.py
+++ b/hysop/operator/tests/test_analytic.py
@@ -19,26 +19,26 @@ from hysop import Field, Box
 class TestAnalyticField(object):
 
     @classmethod
-    def setup_class(cls, 
+    def setup_class(cls,
             enable_extra_tests=__ENABLE_LONG_TESTS__,
             enable_debug_mode=False):
 
         IO.set_default_path('/tmp/hysop_tests/test_analytic')
-        
+
         if enable_debug_mode:
             cls.size_min = 15
             cls.size_max = 16
         else:
             cls.size_min = 23
             cls.size_max = 87
-        
+
         cls.enable_extra_tests = enable_extra_tests
         cls.enable_debug_mode  = enable_debug_mode
-        
+
         cls.t = ScalarParameter(name='t', dtype=HYSOP_REAL)
-        
+
         cls.build_analytic_expressions()
-    
+
     @classmethod
     def build_analytic_expressions(cls):
         from hysop.tools.sympy_utils import enable_pretty_printing
@@ -46,7 +46,7 @@ class TestAnalyticField(object):
         from hysop.symbolic.frame import SymbolicFrame
         from hysop.symbolic.field import curl, laplacian
         enable_pretty_printing()
-        
+
         #at this moment we just test a periodic solver
         analytic_expressions = {}
         analytic_functions = {}
@@ -57,13 +57,13 @@ class TestAnalyticField(object):
             def gen_F0():
                 kis = tuple(random.randint(1,5) for _ in xrange(dim))
                 qis = tuple(npw.random.rand(dim).round(decimals=3).tolist())
-                basis = tuple( (sm.cos(ki*xi+qi), sm.sin(ki*xi+qi)) 
+                basis = tuple( (sm.cos(ki*xi+qi), sm.sin(ki*xi+qi))
                         for (ki,qi,xi) in zip(kis, qis, coords) )
                 F0 = sm.Integer(1) / (sm.Integer(1) + cls.t.s)
                 for i in xrange(dim):
                     F0 *= random.choice(basis[i])
                 return F0
-            
+
             def gen_F1():
                 base = tuple(tuple( coords[i]**k for k in xrange(5) ) for i in xrange(dim))
                 base = tuple(set(npw.prod(x) for x in it.product(*base)))
@@ -71,7 +71,7 @@ class TestAnalyticField(object):
                 F1 = sm.Integer(1) / (sm.Integer(1) + cls.t.s)
                 F1 *= (base * (npw.random.rand(len(base)).round(2))).sum()
                 return F1
-            
+
             params = coords + (cls.t.s,)
             Fs = npw.asarray([gen_F0(), gen_F1()], dtype=object).view(TensorBase)
             fFs = tuple(sm.lambdify(params, F) for F in Fs)
@@ -110,16 +110,16 @@ class TestAnalyticField(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())
-        
+
         domain = Box(length=(1,)*dim)
         F = Field(domain=domain, name='F', dtype=dtype,
                 nb_components=2, register_object=False)
-        
-        self._test_one(shape=shape, dim=dim, dtype=dtype, 
+
+        self._test_one(shape=shape, dim=dim, dtype=dtype,
                 domain=domain, F=F)
 
-    
-    def _test_one(self, shape, dim, dtype, 
+
+    def _test_one(self, shape, dim, dtype,
             domain, F):
 
         print '\nTesting {}D AnalyticField: dtype={} shape={}'.format(
@@ -132,19 +132,19 @@ class TestAnalyticField(object):
         print ' >Testing all implementations:'
 
         implementations = AnalyticField.implementations()
-       
+
         variables = { F:shape }
         fns   = self.analytic_functions[dim]['F']
         exprs = self.analytic_expressions[dim]['F']
 
         def iter_impl(impl):
             base_kwds = dict(field=F, variables=variables,
-                             implementation=impl, 
+                             implementation=impl,
                              name='analytic_{}'.format(str(impl).lower()))
             if impl is Implementation.PYTHON:
                 msg='   *Python: '
                 print msg,
-                yield AnalyticField(formula=self.__analytic_init, 
+                yield AnalyticField(formula=self.__analytic_init,
                                     extra_input_kwds={'fns':fns, 't':self.t},
                                     **base_kwds)
                 print
@@ -153,7 +153,7 @@ class TestAnalyticField(object):
                 msg='   *Opencl: '
                 print msg
                 for cl_env in iter_clenv():
-                    msg='      >platform {}, device {}:'.format(cl_env.platform.name.strip(), 
+                    msg='      >platform {}, device {}:'.format(cl_env.platform.name.strip(),
                                                           cl_env.device.name.strip())
                     print msg,
                     yield AnalyticField(cl_env=cl_env, formula=exprs, **base_kwds)
@@ -162,24 +162,24 @@ class TestAnalyticField(object):
             else:
                 msg='Unknown implementation to test {}.'.format(impl)
                 raise NotImplementedError(msg)
-        
+
         # Compare to analytic solution
         Fref = None
         for impl in implementations:
             for op in iter_impl(impl):
                 op = op.build()
                 dF = op.get_output_discrete_field(F)
-                
+
                 if (Fref is None):
                     dF.initialize(self.__analytic_init, fns=fns, t=self.t)
                     Fref = tuple( data.get().handle.copy() for data in dF.data )
-                
+
                 dF.initialize(self.__random_init)
                 op.apply()
-                
+
                 Fout   = tuple( data.get().handle.copy() for data in dF.data )
                 self._check_output(impl, op, Fref, Fout)
-    
+
     @classmethod
     def _check_output(cls, impl, op, Fref, Fout):
         check_instance(Fref, tuple,   values=npw.ndarray)
@@ -197,12 +197,12 @@ class TestAnalyticField(object):
                 print
                 msg = msg0.format(iname)
                 raise ValueError(msg)
-        
+
         for i, (fout,fref) in enumerate(zip(Fout, Fref)):
             iname = '{}{}'.format('F',i)
             assert fout.dtype == fref.dtype, iname
             assert fout.shape == fref.shape, iname
-            
+
             eps  = npw.finfo(fout.dtype).eps
             dist = npw.abs(fout-fref)
             dinf = npw.max(dist)
@@ -212,7 +212,7 @@ class TestAnalyticField(object):
                 continue
             has_nan = npw.any(npw.isnan(fout))
             has_inf = npw.any(npw.isinf(fout))
-            
+
             print
             print
             print 'Test output comparisson for {} failed for component {}:'.format(name, i)
@@ -246,9 +246,9 @@ class TestAnalyticField(object):
                         print w
                         print
                 print
-            
+
             msg = 'Test failed for {} on component {} for implementation {}.'.format(name, i, impl)
-            raise RuntimeError(msg) 
+            raise RuntimeError(msg)
 
 
     def test_1d_float32(self):
@@ -274,15 +274,15 @@ class TestAnalyticField(object):
             self.test_1d_float64()
             self.test_2d_float64()
             self.test_3d_float64()
-    
+
 if __name__ == '__main__':
-    TestAnalyticField.setup_class(enable_extra_tests=False, 
+    TestAnalyticField.setup_class(enable_extra_tests=False,
                                       enable_debug_mode=False)
-    
+
     test = TestAnalyticField()
 
-    with printoptions(threshold=10000, linewidth=240, 
-                      nanstr='nan', infstr='inf', 
+    with printoptions(threshold=10000, linewidth=240,
+                      nanstr='nan', infstr='inf',
                       formatter={'float': lambda x: '{:>6.2f}'.format(x)}):
         test.perform_tests()
 
diff --git a/hysop/operator/tests/test_custom_symbolic.py b/hysop/operator/tests/test_custom_symbolic.py
index 62dd3e00bc4f95851103d5d99e362545b41a3c12..8f3ca661b1b9769411fac469b4d352aa9fd19ab5 100644
--- a/hysop/operator/tests/test_custom_symbolic.py
+++ b/hysop/operator/tests/test_custom_symbolic.py
@@ -1,4 +1,3 @@
-
 from hysop import Field, Box
 from hysop.deps import np, it, sm
 from hysop.constants import Implementation, ComputeGranularity, SpaceDiscretization, HYSOP_REAL
@@ -23,12 +22,12 @@ from hysop.numerics.odesolvers.runge_kutta import TimeIntegrator, Euler, RK2, RK
 class TestCustomSymbolic(object):
 
     @classmethod
-    def setup_class(cls, 
+    def setup_class(cls,
             enable_extra_tests=__ENABLE_LONG_TESTS__,
             enable_debug_mode=False):
-        
+
         IO.set_default_path('/tmp/hysop_tests/test_custom_symbolic')
-        
+
         if enable_debug_mode:
             cls.size_min0 = 20
             cls.size_max0 = 20
@@ -39,10 +38,10 @@ class TestCustomSymbolic(object):
             cls.size_max0 = 4096
             cls.size_min = 3
             cls.size_max = 16
-        
+
         cls.enable_extra_tests = enable_extra_tests
         cls.enable_debug_mode  = enable_debug_mode
-            
+
         cls.dtypes = [np.int8, np.int16, np.int32, np.int64,
                      np.uint8, np.uint16, np.uint32, np.uint64,
                      np.float32, np.float64]
@@ -56,17 +55,17 @@ class TestCustomSymbolic(object):
         shape = data[0].shape
         if is_integer(dtype):
             for d in data:
-                d[...] = np.random.random_integers(low=0, high=255, size=shape) 
+                d[...] = np.random.random_integers(low=0, high=255, size=shape)
         elif is_fp(dtype):
             for d in data:
                 if pollute:
                     d[...] = np.nan
                 else:
-                    d[...] = np.random.random(size=d.shape) 
+                    d[...] = np.random.random(size=d.shape)
         else:
             msg='Unknown dtype {}.'.format(dtype)
             raise NotImplementedError(msg)
-   
+
     @staticmethod
     def iter_implementations(op_cls, base_kwds):
         for impl in op_cls.implementations():
@@ -74,7 +73,7 @@ class TestCustomSymbolic(object):
             base_kwds['implementation'] = impl
             if impl is Implementation.OPENCL:
                 for cl_env in iter_clenv():
-                    print '      *platform {}, device {}: '.format(cl_env.platform.name.strip(), 
+                    print '      *platform {}, device {}: '.format(cl_env.platform.name.strip(),
                                                           cl_env.device.name.strip()),
                     yield impl, op_cls(cl_env=cl_env, **base_kwds)
             else:
@@ -82,7 +81,7 @@ class TestCustomSymbolic(object):
                 yield impl, op_cls(**base_kwds)
 
     @classmethod
-    def _check_output(cls, impl, op, in_names, refin_buffers, out_names, refout_buffers, 
+    def _check_output(cls, impl, op, in_names, refin_buffers, out_names, refout_buffers,
             out_buffers):
         check_instance(out_buffers, tuple, values=np.ndarray)
         check_instance(refout_buffers, tuple, values=np.ndarray)
@@ -91,7 +90,7 @@ class TestCustomSymbolic(object):
         for i, (oname, out, refout) in enumerate(zip(out_names, out_buffers, refout_buffers)):
             assert refout.dtype == out.dtype, '{} vs {}'.format(refout.dtype, out.dtype)
             assert refout.shape == out.shape, '{} vs {}'.format(refout.shape, out.shape)
-            
+
             has_nan = np.any(np.isnan(refout))
             has_inf = np.any(np.isinf(refout))
             if (out.dtype != refout.dtype):
@@ -103,7 +102,7 @@ class TestCustomSymbolic(object):
 
             has_nan = np.any(np.isnan(out))
             has_inf = np.any(np.isinf(out))
-            
+
             distances = np.abs(out - refout) / np.max(refout)
             if is_integer(out.dtype):
                 mask = (out == refout)
@@ -123,7 +122,7 @@ class TestCustomSymbolic(object):
             if (not has_nan) and (not has_inf) and mask.all():
                 print msg,
                 continue
-            print 
+            print
             print 'Failed to match output of {}:'.format(oname)
             print
             print 'Test output comparisson failed for component {}:'.format(i)
@@ -131,7 +130,7 @@ class TestCustomSymbolic(object):
             print ' *has_inf: {}'.format(has_inf)
             print
             if cls.enable_debug_mode:
-                mask[...] = False 
+                mask[...] = False
             print 'REFERENCE INPUTS:'
             for name, _in in zip(in_names, refin_buffers):
                 print name
@@ -154,16 +153,16 @@ class TestCustomSymbolic(object):
                 print '{} DISTANCES (EPS):'.format(oname)
                 print eps_distances[~mask]
                 print
-            
+
             msg = 'Test failed on component {} for implementation {}.'.format(i, impl)
-            raise RuntimeError(msg) 
+            raise RuntimeError(msg)
         print
-    
+
     def _test_simple(self, dim, _size_min=None, _size_max=None):
         enable_extra_tests = self.enable_extra_tests
         assert dim >= 1
-            
-        print 'TEST SIMPLE' 
+
+        print 'TEST SIMPLE'
         print ' DIM {}'.format(dim)
 
         size_min = first_not_None(_size_min, self.size_min)
@@ -173,34 +172,34 @@ class TestCustomSymbolic(object):
 
         domain = Box(length=(1.0,)*dim)
 
-        discretization  = tuple(np.random.randint(low=size_min, high=size_max+1, 
+        discretization  = tuple(np.random.randint(low=size_min, high=size_max+1,
                                 size=dim-1).tolist())
-        discretization0 = tuple(np.random.randint(low=size_min0, high=size_max0+1, 
+        discretization0 = tuple(np.random.randint(low=size_min0, high=size_max0+1,
                                 size=1).tolist())
         discretization = discretization + discretization0
 
         print ' DISCRETIZATION {}'.format(discretization)
-        
+
         if __ENABLE_LONG_TESTS__:
             dtypes = (np.float32, np.float64)
         else:
             dtypes = (HYSOP_REAL,)
-    
+
         for dtype in dtypes:
             print '  DTYPE {}'.format(dtype.__name__)
-            A = Field(domain=domain, name='A', dtype=dtype,  
+            A = Field(domain=domain, name='A', dtype=dtype,
                     nb_components=1, register_object=False)
-            B = Field(domain=domain, name='B', dtype=dtype, 
+            B = Field(domain=domain, name='B', dtype=dtype,
                     nb_components=1, register_object=False)
-            C = Field(domain=domain, name='C', dtype=dtype, 
+            C = Field(domain=domain, name='C', dtype=dtype,
                     nb_components=3, register_object=False)
-            
+
             P0 = ScalarParameter('P0', dtype=np.float32, initial_value=1.0)
             P1 = ScalarParameter('P1', dtype=np.float64, initial_value=2.0, const=True)
-            T0 = TensorParameter('T0', shape=(13,), dtype=np.int32, 
+            T0 = TensorParameter('T0', shape=(13,), dtype=np.int32,
                     initial_value=np.arange(13, dtype=np.int32))
             T1 = TensorParameter('T1', shape=(3,3), dtype=np.int32)
-            
+
             As = A.s()
             Bs = B.s()
             Cs = C.s()
@@ -209,60 +208,60 @@ class TestCustomSymbolic(object):
 
             for granularity in xrange(dim):
                 print '   GRANULARITY {}'.format(granularity)
-                
+
                 self._test_affect((A,), (42,),      discretization, granularity)
                 self._test_affect((A,B,C), (1,2,3), discretization, granularity)
-                
+
                 if __ENABLE_LONG_TESTS__:
                     expr = Assignment(As, 1+P0s)
                     compute_outputs = lambda ifields, iparams, dfields, ovar, i: \
                             np.full(fill_value=1+iparams['P0'], shape=ovar.resolution, dtype=ovar.dtype)
-                    self._test_expr(expr, compute_outputs, 
+                    self._test_expr(expr, compute_outputs,
                             variables={A:discretization}, method={ComputeGranularity:granularity})
                     P0.value = 4.0
                     compute_outputs = lambda ifields, iparams, dfields, ovar, i: np.full(fill_value=5, shape=ovar.resolution, dtype=ovar.dtype)
-                    self._test_expr(expr, compute_outputs, 
+                    self._test_expr(expr, compute_outputs,
                             variables={A:discretization}, method={ComputeGranularity:granularity})
-                    
+
                     expr = Assignment(Bs, np.sum(T0s[1::2]))
                     compute_outputs = lambda ifields, iparams, dfields, ovar, i: np.full(fill_value=1+3+5+7+9+11, shape=ovar.resolution, dtype=ovar.dtype)
-                    self._test_expr(expr, compute_outputs, 
+                    self._test_expr(expr, compute_outputs,
                             variables={B:discretization}, method={ComputeGranularity:granularity})
 
                     compute_outputs = lambda ifields, iparams, dfields, ovar, i: ifields[C[0]]
-                    self._test_expr(Assignment(As, C.s[0]()), compute_outputs, 
-                            variables={A:discretization, C:discretization}, method={ComputeGranularity:granularity}) 
-                
+                    self._test_expr(Assignment(As, C.s[0]()), compute_outputs,
+                            variables={A:discretization, C:discretization}, method={ComputeGranularity:granularity})
+
                 compute_outputs = lambda ifields, iparams, dfields, ovar, i: 2*ifields[C[0]] + 8
-                self._test_expr(Assignment(As, 2*C.s[0]()+8), compute_outputs, 
-                        variables={A:discretization, C:discretization}, method={ComputeGranularity:granularity}) 
-                
+                self._test_expr(Assignment(As, 2*C.s[0]()+8), compute_outputs,
+                        variables={A:discretization, C:discretization}, method={ComputeGranularity:granularity})
+
                 if (dim==3):
-                    compute_outputs = lambda ifields, iparams, dfields, ovar, i: ifields[C][0]*ifields[C][1]*ifields[C][2]
-                    self._test_expr(Assignment(As, C.s[0]()*C.s[1]()*C.s[2]()), compute_outputs, 
+                    compute_outputs = lambda ifields, iparams, dfields, ovar, i: ifields[C[0]]*ifields[C[1]]*ifields[C[2]]
+                    self._test_expr(Assignment(As, C.s[0]()*C.s[1]()*C.s[2]()), compute_outputs,
                             variables={A:discretization, C:discretization}, method={ComputeGranularity:granularity})
-                    
+
                     x0  = domain.frame.coords[0]
-                    compute_outputs = lambda ifields, iparams, dfields, ovar, i: 2*ifields[A] - 3*ifields[C][0]*ifields[C][1]*ifields[C][2]
-                    self._test_expr(Assignment(As, 2*As-3*C.s[0]()*C.s[1]()*C.s[2]()), compute_outputs, 
+                    compute_outputs = lambda ifields, iparams, dfields, ovar, i: 2*ifields[A] - 3*ifields[C[0]]*ifields[C[1]]*ifields[C[2]]
+                    self._test_expr(Assignment(As, 2*As-3*C.s[0]()*C.s[1]()*C.s[2]()), compute_outputs,
                             variables={A:discretization, C:discretization}, method={ComputeGranularity:granularity})
-                    
-                    compute_outputs = lambda ifields, iparams, dfields, ovar, i: (2*np.cos(ifields[C][0])*np.sin(ifields[C][1])*np.tan(ifields[C][2])).astype(ovar.dtype)
-                    self._test_expr(Assignment(As, 2*sm.cos(C.s[0]())*sm.sin(C.s[1]())*sm.tan(C.s[2]())), compute_outputs, 
+
+                    compute_outputs = lambda ifields, iparams, dfields, ovar, i: (2*np.cos(ifields[C[0]])*np.sin(ifields[C[1]])*np.tan(ifields[C[2]])).astype(ovar.dtype)
+                    self._test_expr(Assignment(As, 2*sm.cos(C.s[0]())*sm.sin(C.s[1]())*sm.tan(C.s[2]())), compute_outputs,
                             variables={A:discretization, C:discretization}, method={ComputeGranularity:granularity})
-    
-    def _test_expr(self, exprs, compute_outputs, variables, method, 
+
+    def _test_expr(self, exprs, compute_outputs, variables, method,
             apply_kwds=None, no_ref_view=False, dt=None):
         exprs = to_tuple(exprs)
         print '    CustomExpr: {}'.format(' || '.join(str(e) for e in exprs))
-        
+
         assert ComputeGranularity in method
         base_kwds = dict(name='array_affect', exprs=exprs, variables=variables, method=method, dt=dt)
         apply_kwds = first_not_None(apply_kwds, {})
 
         for impl, op in self.iter_implementations(CustomSymbolicOperator, base_kwds):
             problem = op.build()
-    
+
             for field, dfield in problem.iter_output_discrete_fields():
                 dfield.initialize(self.__field_init, dtype=dfield.dtype, pollute=True, only_finite=False)
             for field, dfield in problem.iter_input_discrete_fields():
@@ -295,7 +294,7 @@ class TestCustomSymbolic(object):
                 out_names.append(pname)
 
             problem.apply(**apply_kwds)
-            
+
             out = {'f': {}, 'p': {}, 'dfields': {}}
             for field, ofield in problem.iter_output_discrete_fields(as_scalars=True):
                 out['f'][field] = ()
@@ -322,7 +321,7 @@ class TestCustomSymbolic(object):
                 for _ in x['p'].values():
                     out += (np.asarray(_),)
                 return out
-            
+
             self._check_output(impl, op, in_names, flatten(refin), out_names, flatten(refout), flatten(out))
 
     def _test_affect(self, fields, rhs, discretization, granularity):
@@ -333,9 +332,9 @@ class TestCustomSymbolic(object):
             for i in xrange(f.nb_components):
                 e = Assignment(f.s[i](),c)
                 exprs += (e,)
-        
+
         print '    CustomExpr: {}'.format(' || '.join(str(e) for e in exprs))
-        
+
         method={ComputeGranularity: granularity}
         base_kwds = dict(name='custom_affect', exprs=exprs, variables=variables, method=method)
         for impl, op in self.iter_implementations(CustomSymbolicOperator, base_kwds):
@@ -350,9 +349,9 @@ class TestCustomSymbolic(object):
                 for i in xrange(dfield.nb_components):
                     refin += (dfield.data[i].get().handle[view],)
                     in_names.append(ofield.name+'::{}'.format(i))
-            
+
             problem.apply()
-            
+
             out = ()
             refout = ()
             out_names = []
@@ -374,7 +373,7 @@ class TestCustomSymbolic(object):
     def _test_stencil(self, dim, _size_min=None, _size_max=None):
         enable_extra_tests = self.enable_extra_tests
         assert dim >= 1
-        print 'TEST STENCIL' 
+        print 'TEST STENCIL'
         print ' DIM {}'.format(dim)
 
         size_min = first_not_None(_size_min, self.size_min)
@@ -385,23 +384,23 @@ class TestCustomSymbolic(object):
         domain = Box(length=(1.0,)*dim)
         frame = domain.frame
 
-        discretization  = tuple(np.random.randint(low=size_min, high=size_max+1, 
+        discretization  = tuple(np.random.randint(low=size_min, high=size_max+1,
                                 size=dim-1).tolist())
-        discretization0 = tuple(np.random.randint(low=size_min0, high=size_max0+1, 
+        discretization0 = tuple(np.random.randint(low=size_min0, high=size_max0+1,
                                 size=1).tolist())
         discretization = discretization + discretization0
 
         print ' DISCRETIZATION {}'.format(discretization)
-        A = Field(domain=domain, name='A', dtype=np.float32,  
+        A = Field(domain=domain, name='A', dtype=np.float32,
                 nb_components=1, register_object=False)
-        B = Field(domain=domain, name='B', dtype=np.float32, 
+        B = Field(domain=domain, name='B', dtype=np.float32,
                 nb_components=2, register_object=False)
-        C = Field(domain=domain, name='C', dtype=np.float64, 
+        C = Field(domain=domain, name='C', dtype=np.float64,
                 nb_components=3, register_object=False)
 
         P0 = ScalarParameter('P0', dtype=np.float32, initial_value=3.14)
         T0 = TensorParameter('T',  shape=(3,), dtype=np.int32, initial_value=[4,3,2])
-        
+
         As = A.s()
         Bs = B.s()
         Cs = C.s()
@@ -409,7 +408,7 @@ class TestCustomSymbolic(object):
         T0s = T0.s
 
         x0  = frame.coords[0]
-        
+
         csg = CenteredStencilGenerator()
         csg.configure(dtype=MPQ, dim=1)
 
@@ -417,31 +416,31 @@ class TestCustomSymbolic(object):
             orders = (2,4)
         else:
             orders = (4,)
-        
+
         for order in orders:
             print '  ORDER {}'.format(order)
             for granularity in xrange(dim):
                 print '   GRANULARITY {}'.format(granularity)
-                
+
                 if __ENABLE_LONG_TESTS__:
                     expr    = Assignment(As, B.s[0]().diff(x0,x0))
                     stencil = csg.generate_exact_stencil(derivative=2, order=order)
                     def compute_outputs(fields, iparams, ifields, ovar, i):
                         return stencil.apply(fields[B[0]], symbols={stencil.dx:ovar.space_step[-1]}, axis=-1)[ifields[B[0]].compute_slices]
-                    self._test_expr(expr, compute_outputs, 
-                            variables={A:discretization, B:discretization}, 
+                    self._test_expr(expr, compute_outputs,
+                            variables={A:discretization, B:discretization},
                             method={ComputeGranularity:granularity, SpaceDiscretization:order})
-                    
-                    
+
+
                     expr    = Assignment(As, 1+2*B.s[0]().diff(x0))
                     stencil = csg.generate_exact_stencil(derivative=1, order=order)
                     def compute_outputs(fields, iparams, ifields, ovar, i):
                         res = stencil.apply(fields[B[0]], symbols={stencil.dx:ovar.space_step[-1]}, axis=-1)[ifields[B[0]].compute_slices]
                         return 1+2*res
-                    self._test_expr(expr, compute_outputs, 
-                            variables={A:discretization, B:discretization}, 
+                    self._test_expr(expr, compute_outputs,
+                            variables={A:discretization, B:discretization},
                             method={ComputeGranularity:granularity, SpaceDiscretization:order})
-                    
+
                     expr    = Assignment(Bs[0](), (1+np.dot([3,-2], Bs.diff(x0)))*As.diff(x0))
 
                     stencil = csg.generate_exact_stencil(derivative=1, order=order)
@@ -454,10 +453,10 @@ class TestCustomSymbolic(object):
                         else:
                             return ovar.data[i].get().handle
 
-                    self._test_expr(expr, compute_outputs, 
-                            variables={A:discretization, B:discretization}, 
+                    self._test_expr(expr, compute_outputs,
+                            variables={A:discretization, B:discretization},
                             method={ComputeGranularity:granularity, SpaceDiscretization:order})
-                
+
                 expr0    = Assignment(B.s[0](), (1+np.dot([3,-2], Bs.diff(x0)))*As.diff(x0))
                 expr1    = Assignment(B.s[1](), (1-np.dot([3,-2], Bs.diff(x0)))*As.diff(x0))
                 expr2    = Assignment(A.s[0](), 0.27 + 1.57*sm.cos(Bs[0].diff(x0))*sm.sin(Bs[1].diff(x0)))
@@ -473,10 +472,10 @@ class TestCustomSymbolic(object):
                             return (1-(3*res0)+(2*res1))*res2
                     else:
                         return 0.27 + 1.57*np.cos(res0)*np.sin(res1)
-                self._test_expr(expr, compute_outputs, 
-                        variables={A:discretization, B:discretization}, 
+                self._test_expr(expr, compute_outputs,
+                        variables={A:discretization, B:discretization},
                         method={ComputeGranularity:granularity, SpaceDiscretization:order})
-                
+
                 expr = Assignment(B.s[0](), T0s[2]*((4+P0s)*As*B.s[0]().diff(x0)+T0s[1]*Bs[1]).diff(x0) + P0s*T0s[0])
                 stencil0 = csg.generate_exact_stencil(derivative=1, order=order)
                 stencil1 = csg.generate_exact_stencil(derivative=2, order=order)
@@ -494,15 +493,15 @@ class TestCustomSymbolic(object):
                         return P0[0]*T0[0] + T0[2]*((4+P0[0])*(A0*d2B0+dA0*dB0) + T0[1]*dB1)
                     else:
                         return ovar.dfield.data[i].get().handle
-                self._test_expr(expr, compute_outputs, 
-                        variables={A:discretization, B:discretization}, 
+                self._test_expr(expr, compute_outputs,
+                        variables={A:discretization, B:discretization},
                         method={ComputeGranularity:granularity, SpaceDiscretization:order})
-    
-                
+
+
     def _test_time_integrator(self, dim, _size_min=None, _size_max=None):
         enable_extra_tests = self.enable_extra_tests
         assert dim >= 1
-        print 'TEST INTEGRATOR' 
+        print 'TEST INTEGRATOR'
         print ' DIM {}'.format(dim)
 
         size_min = first_not_None(_size_min, self.size_min)
@@ -513,38 +512,38 @@ class TestCustomSymbolic(object):
         domain = Box(length=(1.0,)*dim)
         frame = domain.frame
 
-        discretization  = tuple(np.random.randint(low=size_min, high=size_max+1, 
+        discretization  = tuple(np.random.randint(low=size_min, high=size_max+1,
                                 size=dim-1).tolist())
-        discretization0 = tuple(np.random.randint(low=size_min0, high=size_max0+1, 
+        discretization0 = tuple(np.random.randint(low=size_min0, high=size_max0+1,
                                 size=1).tolist())
         discretization = discretization + discretization0
 
         print ' DISCRETIZATION {}'.format(discretization)
-        A = Field(domain=domain, name='A', dtype=np.float32,  
+        A = Field(domain=domain, name='A', dtype=np.float32,
                 nb_components=1, register_object=False)
-        B = Field(domain=domain, name='B', dtype=np.float32, 
+        B = Field(domain=domain, name='B', dtype=np.float32,
                 nb_components=2, register_object=False)
-        C = Field(domain=domain, name='C', dtype=np.float32, 
+        C = Field(domain=domain, name='C', dtype=np.float32,
                 nb_components=3, register_object=False)
-            
+
         P0 = ScalarParameter('P0', dtype=np.float32, initial_value=3.14)
         T0 = TensorParameter('T',  shape=(3,), dtype=np.int32, initial_value=[4,3,2])
-        
+
         As = A.s()
         Bs = B.s()
         Cs = C.s()
         P0s = P0.s
         T0s = T0.s
-        
+
         T  = ScalarParameter(name=frame.time,   dtype=np.float32, initial_value=0.0)
         DT = ScalarParameter(name=dtime_symbol, dtype=np.float32, initial_value=0.1)
 
         t = T.s
         x0  = frame.coords[0]
-        
+
         csg = CenteredStencilGenerator()
         csg.configure(dtype=MPQ, dim=1)
-        
+
         order = 4
         granularity = 0
         print '  GRANULARITY {}'.format(granularity)
@@ -553,7 +552,7 @@ class TestCustomSymbolic(object):
         D1 = csg.generate_exact_stencil(derivative=1, order=order)
         D2 = csg.generate_exact_stencil(derivative=2, order=order)
         D3 = csg.generate_exact_stencil(derivative=3, order=order)
-        
+
         if __ENABLE_LONG_TESTS__:
             integrators = [Euler, RK2, RK4, RK4_38]
         else:
@@ -561,7 +560,7 @@ class TestCustomSymbolic(object):
 
         for integrator in integrators:
             print '   INTEGRATOR {}'.format(integrator)
-            
+
             if __ENABLE_LONG_TESTS__:
                 expr = Assignment(As.diff(t), 0)
                 def compute_outputs(fields, iparams, ifields, ovar, i):
@@ -570,9 +569,9 @@ class TestCustomSymbolic(object):
                         return fields[A].copy()
                     else:
                         return ovar.dfield.data[i].get().handle
-                self._test_expr(expr, compute_outputs, 
-                        variables={A:discretization, B:discretization, C:discretization}, 
-                        method={ComputeGranularity:granularity, SpaceDiscretization:order, 
+                self._test_expr(expr, compute_outputs,
+                        variables={A:discretization, B:discretization, C:discretization},
+                        method={ComputeGranularity:granularity, SpaceDiscretization:order,
                                 TimeIntegrator:integrator}, dt=DT)
 
                 expr = Assignment(As.diff(t), 1)
@@ -587,12 +586,12 @@ class TestCustomSymbolic(object):
                         return Xout['a']
                     else:
                         return ovar.dfield.data[i].get().handle
-                self._test_expr(expr, compute_outputs, 
-                        variables={A:discretization, B:discretization, C:discretization}, 
-                        method={ComputeGranularity:granularity, SpaceDiscretization:order, 
+                self._test_expr(expr, compute_outputs,
+                        variables={A:discretization, B:discretization, C:discretization},
+                        method={ComputeGranularity:granularity, SpaceDiscretization:order,
                                 TimeIntegrator:integrator}, dt=DT)
-                
-                
+
+
                 expr = Assignment(As.diff(t), np.dot([2,-3], Bs))
                 def compute_outputs(fields, iparams, ifields, ovar, i):
                     A0 = fields[A]
@@ -607,11 +606,11 @@ class TestCustomSymbolic(object):
                         return Xout['a']
                     else:
                         return ovar.dfield.data[i].get().handle
-                self._test_expr(expr, compute_outputs, 
-                        variables={A:discretization, B:discretization, C:discretization}, 
-                        method={ComputeGranularity:granularity, SpaceDiscretization:order, 
+                self._test_expr(expr, compute_outputs,
+                        variables={A:discretization, B:discretization, C:discretization},
+                        method={ComputeGranularity:granularity, SpaceDiscretization:order,
                                 TimeIntegrator:integrator}, dt=DT)
-                
+
                 expr = Assignment(As.diff(t), As*Bs[0] + Bs[1])
                 def compute_outputs(fields, iparams, ifields, ovar, i):
                     A0 = fields[A]
@@ -626,11 +625,11 @@ class TestCustomSymbolic(object):
                         return Xout['a']
                     else:
                         return ovar.dfield.data[i].get().handle
-                self._test_expr(expr, compute_outputs, 
-                        variables={A:discretization, B:discretization, C:discretization}, 
-                        method={ComputeGranularity:granularity, SpaceDiscretization:order, 
+                self._test_expr(expr, compute_outputs,
+                        variables={A:discretization, B:discretization, C:discretization},
+                        method={ComputeGranularity:granularity, SpaceDiscretization:order,
                                 TimeIntegrator:integrator}, dt=DT)
-                
+
                 expr0 = Assignment(Bs[0].diff(t), 1+Bs[1])
                 expr1 = Assignment(Bs[1].diff(t), 1-Bs[0])
                 expr = (expr0, expr1)
@@ -649,12 +648,12 @@ class TestCustomSymbolic(object):
                         return Xout[varname]
                     else:
                         return ovar.dfield.data[i].get().handle
-                self._test_expr(expr, compute_outputs, 
-                        variables={A:discretization, B:discretization, C:discretization}, 
-                        method={ComputeGranularity:granularity, SpaceDiscretization:order, 
+                self._test_expr(expr, compute_outputs,
+                        variables={A:discretization, B:discretization, C:discretization},
+                        method={ComputeGranularity:granularity, SpaceDiscretization:order,
                                 TimeIntegrator:integrator}, dt=DT)
-                
-                
+
+
                 expr = Assignment(As.diff(t), Bs[0].diff(x0))
                 def compute_outputs(fields, iparams, ifields, ovar, i):
                     A0 = fields[A]
@@ -669,12 +668,12 @@ class TestCustomSymbolic(object):
                         return Xout['a']
                     else:
                         return ovar.dfield.data[i].get().handle
-                self._test_expr(expr, compute_outputs, 
-                        variables={A:discretization, B:discretization}, 
-                        method={ComputeGranularity:granularity, SpaceDiscretization:order, 
+                self._test_expr(expr, compute_outputs,
+                        variables={A:discretization, B:discretization},
+                        method={ComputeGranularity:granularity, SpaceDiscretization:order,
                                 TimeIntegrator:integrator}, dt=DT)
-                
-                
+
+
                 expr = Assignment(As.diff(t), As.diff(x0))
                 def compute_outputs(fields, iparams, ifields, ovar, i):
                     A0 = fields[A]
@@ -689,11 +688,11 @@ class TestCustomSymbolic(object):
                         return Xout['a']
                     else:
                         return ovar.dfield.data[i].get().handle
-                self._test_expr(expr, compute_outputs, 
-                        variables={A:discretization}, 
-                        method={ComputeGranularity:granularity, SpaceDiscretization:order, 
+                self._test_expr(expr, compute_outputs,
+                        variables={A:discretization},
+                        method={ComputeGranularity:granularity, SpaceDiscretization:order,
                                 TimeIntegrator:integrator}, dt=DT, no_ref_view=True)
-                
+
                 expr = Assignment(As.diff(t), P0s*Bs[0]*As.diff(x0) + Bs[1].diff(x0,x0)*As)
                 def compute_outputs(fields, iparams, ifields, ovar, i):
                     A0 = fields[A]
@@ -715,11 +714,11 @@ class TestCustomSymbolic(object):
                         return Xout['a']
                     else:
                         return ovar.dfield.data[i].get().handle
-                self._test_expr(expr, compute_outputs, 
-                        variables={A:discretization, B:discretization}, 
-                        method={ComputeGranularity:granularity, SpaceDiscretization:order, 
+                self._test_expr(expr, compute_outputs,
+                        variables={A:discretization, B:discretization},
+                        method={ComputeGranularity:granularity, SpaceDiscretization:order,
                                 TimeIntegrator:integrator}, dt=DT, no_ref_view=True)
-            
+
             expr0 = Assignment(Bs[0].diff(t), 1 - As*Bs[1].diff(x0))
             expr1 = Assignment(Bs[1].diff(t), 1 + As*Bs[0].diff(x0))
             expr = (expr0, expr1)
@@ -727,13 +726,13 @@ class TestCustomSymbolic(object):
                 A0 = fields[A]
                 B0 = fields[B[0]]
                 B1 = fields[B[1]]
-                b0_view = ifields[B[0]].compute_slices 
-                b1_view = ifields[B[1]].compute_slices 
+                b0_view = ifields[B[0]].compute_slices
+                b1_view = ifields[B[1]].compute_slices
                 assert order%2==0
                 bg = order//2
                 v = (Ellipsis, slice(bg, -bg))
                 if (ovar.dfield._field in B.fields):
-                    Xin   = { 'b0': B0, 
+                    Xin   = { 'b0': B0,
                               'b1': B1 }
                     Xout  = { 'b0': np.empty_like(Xin['b0'][ifields[B[0]].compute_slices]),
                               'b1': np.empty_like(Xin['b1'][ifields[B[1]].compute_slices]) }
@@ -753,9 +752,9 @@ class TestCustomSymbolic(object):
                     return Xout[varname]
                 else:
                     return ovar.dfield.data[i].get().handle
-            self._test_expr(expr, compute_outputs, 
-                    variables={A:discretization, B:discretization}, 
-                    method={ComputeGranularity:granularity, SpaceDiscretization:order, 
+            self._test_expr(expr, compute_outputs,
+                    variables={A:discretization, B:discretization},
+                    method={ComputeGranularity:granularity, SpaceDiscretization:order,
                             TimeIntegrator:integrator}, dt=DT, no_ref_view=True)
 
     def test_simple_1d(self):
@@ -782,7 +781,7 @@ class TestCustomSymbolic(object):
         self.test_simple_2d()
         if __ENABLE_LONG_TESTS__:
             self.test_simple_3d()
-       
+
         self.test_stencil_1d()
         self.test_stencil_2d()
         if __ENABLE_LONG_TESTS__:
@@ -792,13 +791,13 @@ class TestCustomSymbolic(object):
         self.test_time_integrator_2d()
         if __ENABLE_LONG_TESTS__:
             self.test_time_integrator_3d()
-    
+
 if __name__ == '__main__':
-    TestCustomSymbolic.setup_class(enable_extra_tests=False, 
+    TestCustomSymbolic.setup_class(enable_extra_tests=False,
                                       enable_debug_mode=False)
-    
+
     enable_pretty_printing()
-    
+
     test = TestCustomSymbolic()
     test.perform_tests()
 
diff --git a/hysop/operator/tests/test_diffusion.py b/hysop/operator/tests/test_diffusion.py
index ddf2cbf9a4eaacb57ef553548c9a1d703454b449..4ac8706b20d697db6816486d8787ce969c060398 100644
--- a/hysop/operator/tests/test_diffusion.py
+++ b/hysop/operator/tests/test_diffusion.py
@@ -43,7 +43,7 @@ class TestDiffusionOperator(object):
         size_max = first_not_None(size_max, self.size_max) + 1
 
         shape = tuple(npw.random.randint(low=size_min, high=size_max, size=dim).tolist())
-    
+
         if __ENABLE_LONG_TESTS__:
             flt_types = (npw.float32, npw.float64)
         else:
@@ -105,10 +105,10 @@ class TestDiffusionOperator(object):
 
         implementations = Diffusion.implementations()
         ref_impl = Implementation.PYTHON # ref impl is always the first
-        
+
         def iter_impl(impl):
             base_kwds = dict(Fin=fin, nu=nu, dt=dt,
-                    variables=variables, implementation=impl, 
+                    variables=variables, implementation=impl,
                     name='test_diffusion_{}'.format(str(impl).lower()))
             if not is_inplace:
                 base_kwds['Fout'] = fout
@@ -224,7 +224,7 @@ class TestDiffusionOperator(object):
                             if not mask.all():
                                 msg='\nFATAL ERROR: Output is not finite on axis {}.\n'.format(i)
                                 print msg
-                                npw.fancy_print(output[i], 
+                                npw.fancy_print(output[i],
                                         replace_values={(lambda a: npw.isfinite(a)): '.'})
                                 raise RuntimeError(msg)
                             di = npw.abs(reference[i] - output[i])
@@ -269,7 +269,8 @@ class TestDiffusionOperator(object):
     def test_diffusion_2D_inplace(self):
         self._test(dim=2, is_inplace=True)
     def test_diffusion_3D_inplace(self):
-        self._test(dim=3, is_inplace=True)
+        if __ENABLE_LONG_TESTS__:
+            self._test(dim=3, is_inplace=True)
 
     def perform_tests(self):
         self.test_diffusion_1D_inplace()
diff --git a/hysop/operator/tests/test_directional_advection.py b/hysop/operator/tests/test_directional_advection.py
index 9463d3745899f5668bdcf916e8a5a3e58d0d7449..05e81ccff195130237328d68888f2b34d5cd4271 100644
--- a/hysop/operator/tests/test_directional_advection.py
+++ b/hysop/operator/tests/test_directional_advection.py
@@ -124,14 +124,14 @@ class TestDirectionalAdvectionOperator(object):
         implementations = [ref_impl] + implementations
 
         method = {TimeIntegrator: time_integrator, Remesh: remesh_kernel}
-       
+
         def iter_impl(impl):
             base_kwds = dict(velocity=vin, advected_fields=sin, advected_fields_out=sout, dt=dt,
-                    velocity_cfl=velocity_cfl, variables=variables, implementation=impl, 
+                    velocity_cfl=velocity_cfl, variables=variables, implementation=impl,
                     method=method, name='advection_{}'.format(str(impl).lower()))
             da = DirectionalAdvection(**base_kwds)
             if impl is Implementation.PYTHON:
-                split = StrangSplitting(splitting_dim=dim, 
+                split = StrangSplitting(splitting_dim=dim,
                                         order=StrangOrder.STRANG_SECOND_ORDER)
                 split.push_operators(da)
                 force_tstate = ForceTopologyState(fields=variables.keys(),
@@ -142,9 +142,9 @@ class TestDirectionalAdvectionOperator(object):
                 yield 'default', graph
             elif impl is Implementation.OPENCL:
                 for cl_env in iter_clenv():
-                    msg='platform {}, device {}'.format(cl_env.platform.name.strip(), 
+                    msg='platform {}, device {}'.format(cl_env.platform.name.strip(),
                                                         cl_env.device.name.strip())
-                    split = StrangSplitting(splitting_dim=dim, 
+                    split = StrangSplitting(splitting_dim=dim,
                                            extra_kwds=dict(cl_env=cl_env),
                                            order=StrangOrder.STRANG_SECOND_ORDER)
                     split.push_operators(da)
@@ -178,7 +178,7 @@ class TestDirectionalAdvectionOperator(object):
                     #print 'SWITCHING TO AXES {}'.format(axes)
                     ref_outputs = reference_fields.setdefault(axes, {})
                     napplies = 10
-                    Vi = npw.asarray([+1.0 if (i in axes) else +0.0 
+                    Vi = npw.asarray([+1.0 if (i in axes) else +0.0
                                         for i in xrange(dim)], dtype=dtype)
 
                     dvin  = graph.get_input_discrete_field(vin)
@@ -225,7 +225,7 @@ class TestDirectionalAdvectionOperator(object):
                                     print 'DSREF - DSOUT'
                                     for (output, ref) in zip(dsout, dsref):
                                         print (output.sdata[output.compute_slices].get() - ref.sdata[ref.compute_slices].get())
-                                    msg='Test failed with V={}, k={}, dxk={}, inter-field L2 distances are {}.' 
+                                    msg='Test failed with V={}, k={}, dxk={}, inter-field L2 distances are {}.'
                                     msg=msg.format(Vi, k, to_tuple(dxk, cast=float), to_tuple(d, cast=float))
                                     raise RuntimeError(msg)
                                 ref_outputs[k] = output
@@ -272,25 +272,25 @@ class TestDirectionalAdvectionOperator(object):
                         sys.stdout.flush()
                         raise
 
-    def test_advec_1D_out_of_place(self):
-        self._test(dim=1, is_inplace=False)
-    def test_advec_2D_out_of_place(self):
-        self._test(dim=2, is_inplace=False)
-    def test_advec_3D_out_of_place(self):
-        self._test(dim=3, is_inplace=False)
-    
+    # def test_advec_1D_out_of_place(self):
+    #     self._test(dim=1, is_inplace=False)
+    # def test_advec_2D_out_of_place(self):
+    #     self._test(dim=2, is_inplace=False)
+    # def test_advec_3D_out_of_place(self):
+    #     self._test(dim=3, is_inplace=False)
+
     def test_advec_1D_inplace(self):
         self._test(dim=1, is_inplace=True)
     def test_advec_2D_inplace(self):
         self._test(dim=2, is_inplace=True)
     def test_advec_3D_inplace(self):
         self._test(dim=3, is_inplace=True)
-    
+
     def perform_tests(self):
         self.test_advec_1D_inplace()
         self.test_advec_2D_inplace()
         self.test_advec_3D_inplace()
-        
+
         #self.test_advec_1D_out_of_place()
         #self.test_advec_2D_out_of_place()
         #self.test_advec_3D_out_of_place()
diff --git a/hysop/operator/tests/test_directional_diffusion.py b/hysop/operator/tests/test_directional_diffusion.py
index 9d70d201c3c636e90371e5530616237684794585..4949daad99774bbfb3261ceb694c8b2828041bcb 100644
--- a/hysop/operator/tests/test_directional_diffusion.py
+++ b/hysop/operator/tests/test_directional_diffusion.py
@@ -152,10 +152,10 @@ class TestDirectionalDiffusionOperator(object):
                 if (not is_inplace):
                     split.push_copy(dst=fin, src=fout)
                 split = split.build()
-                
+
                 dfin  = split.get_input_discrete_field(fin)
                 dfout = split.get_output_discrete_field(fout)
-                
+
                 view = dfin.compute_slices
 
                 reference = None
@@ -170,15 +170,15 @@ class TestDirectionalDiffusionOperator(object):
                     if not npw.all(npw.isfinite(S0)):
                         msg='Integral is not finite. Got {}.'.format(S0)
                         raise RuntimeError(msg)
-                    
-                    _input = tuple(dfin.data[i].get().handle.copy() 
+
+                    _input = tuple(dfin.data[i].get().handle.copy()
                                     for i in xrange(dfin.nb_components))
 
                     if (reference is None):
                         reference  = self._compute_reference(fin=dfin, dfin=_input, dt=dt(), dim=dim,
                             time_integrator=time_integrator, order=order, coeffs=coeffs, view=view)
                     split.apply()
-                    
+
                     output = tuple(dfout.data[i].get().handle.copy()[view]
                                 for i in xrange(dfout.nb_components))
 
@@ -272,12 +272,12 @@ class TestDirectionalDiffusionOperator(object):
         return out
 
 
-    def test_diffusion_1D_out_of_place(self):
-        self._test(dim=1, is_inplace=False)
-    def test_diffusion_2D_out_of_place(self):
-        self._test(dim=2, is_inplace=False)
-    def test_diffusion_3D_out_of_place(self):
-        self._test(dim=3, is_inplace=False)
+    # def test_diffusion_1D_out_of_place(self):
+    #     self._test(dim=1, is_inplace=False)
+    # def test_diffusion_2D_out_of_place(self):
+    #     self._test(dim=2, is_inplace=False)
+    # def test_diffusion_3D_out_of_place(self):
+    #     self._test(dim=3, is_inplace=False)
 
     def test_diffusion_1D_inplace(self):
         self._test(dim=1, is_inplace=True)
@@ -298,9 +298,9 @@ class TestDirectionalDiffusionOperator(object):
 
 
 if __name__ == '__main__':
-    TestDirectionalDiffusionOperator.setup_class(enable_extra_tests=False, 
+    TestDirectionalDiffusionOperator.setup_class(enable_extra_tests=False,
                                       enable_debug_mode=False)
-    
+
     test = TestDirectionalDiffusionOperator()
 
     enable_pretty_printing()
diff --git a/hysop/operator/tests/test_directional_stretching.py b/hysop/operator/tests/test_directional_stretching.py
index 99f8006e251d9677312a3ee403957621932efb78..326327d3acf0c7ce3945347ea431ee5b211f064e 100644
--- a/hysop/operator/tests/test_directional_stretching.py
+++ b/hysop/operator/tests/test_directional_stretching.py
@@ -119,8 +119,8 @@ class TestDirectionalStretchingOperator(object):
             msg='Out of place stretching has not been implemented yet.'
             raise NotImplementedError(msg)
 
-        implementations = (Implementation.PYTHON, Implementation.OPENCL) 
-        # currently we have to different implementations: 
+        implementations = (Implementation.PYTHON, Implementation.OPENCL)
+        # currently we have to different implementations:
         #  DirectionalStretching and StaticDirectionalStretching
 
         method = {TimeIntegrator: time_integrator, SpaceDiscretization: order}
@@ -158,11 +158,11 @@ class TestDirectionalStretchingOperator(object):
                 if (not is_inplace):
                     split.push_copy(dst=win, src=wout)
                 split = split.build()
-                
+
                 dvin  = split.get_input_discrete_field(vin)
                 dwin  = split.get_input_discrete_field(win)
                 dwout = split.get_output_discrete_field(wout)
-                
+
                 sys.stdout.write('.')
                 sys.stdout.flush()
                 try:
@@ -170,9 +170,9 @@ class TestDirectionalStretchingOperator(object):
                         dwin.initialize(self.__random_init)
                         dvin.initialize(self.__random_init)
 
-                        V_input = tuple(df.sdata.get()[df.compute_slices].handle.copy() 
+                        V_input = tuple(df.sdata.get()[df.compute_slices].handle.copy()
                                         for df in dvin)
-                        W_input = tuple(df.sdata.get()[df.compute_slices].handle.copy() 
+                        W_input = tuple(df.sdata.get()[df.compute_slices].handle.copy()
                                         for df in dwin)
 
                         V0 = dvin.integrate()
@@ -194,8 +194,8 @@ class TestDirectionalStretchingOperator(object):
                     dvin.initialize(V_input, only_finite=False)
 
                     split.apply(dt=dt)
-                    
-                    W_output = tuple(df.sdata.get()[df.compute_slices].handle.copy() 
+
+                    W_output = tuple(df.sdata.get()[df.compute_slices].handle.copy()
                                      for df in dwout)
 
                     V1 = dvin.integrate()
@@ -273,7 +273,7 @@ class TestDirectionalStretchingOperator(object):
         wis   = tuple('W{}'.format(i) for i in xrange(3))
         Xin   = dict(zip(wis, Win))
         views = { 'W{}'.format(i): Wi.compute_slices for (i,Wi) in enumerate(dwin) }
-        
+
         ghost_exchanger = dwin.build_ghost_exchanger(data=Win)
 
         for (j, dt_coeff) in zip(directions, dt_coeffs):
@@ -291,7 +291,7 @@ class TestDirectionalStretchingOperator(object):
                     D1_VW = ()
                     for (i,Vi) in enumerate(V):
                         ghosts = [0,]*3
-                        ghosts[2-j] = dvin.ghosts[-1] 
+                        ghosts[2-j] = dvin.ghosts[-1]
                         a=Vi*W[j]
                         d1_vw = D1(a=a, axis=ndir, symbols={D1.dx:dwin.space_step[ndir]})
                         D1_VW += (d1_vw,)
@@ -340,23 +340,23 @@ class TestDirectionalStretchingOperator(object):
 
     def test_stretching_3D_inplace_conservative(self):
         self._test(formulation=StretchingFormulation.CONSERVATIVE, is_inplace=True)
-    def test_stretching_3D_inplace_gradU_W(self):
-        self._test(formulation=StretchingFormulation.GRAD_UW, is_inplace=True)
-    def test_stretching_3D_inplace_gradU_T_W(self):
-        self._test(formulation=StretchingFormulation.GRAD_UW_T, is_inplace=True)
-    def test_stretching_3D_inplace_mixed_gradient(self):
-        As = (None, sm.Rational(1,3))
-        self._test(formulation=StretchingFormulation.MIXED_GRAD_UW, is_inplace=True, As=As)
-
-    def test_stretching_3D_out_of_place_conservative(self):
-        self._test(formulation=StretchingFormulation.CONSERVATIVE, is_inplace=False)
-    def test_stretching_3D_out_of_place_gradU_W(self):
-        self._test(formulation=StretchingFormulation.GRAD_UW, is_inplace=False)
-    def test_stretching_3D_out_of_place_gradU_T_W(self):
-        self._test(formulation=StretchingFormulation.GRAD_UW_T, is_inplace=False)
-    def test_stretching_3D_out_of_place_mixed_gradient(self):
-        As = (None, sm.Rational(1,3))
-        self._test(formulation=StretchingFormulation.MIXED_GRAD_UW, is_inplace=False, As=As)
+    # def test_stretching_3D_inplace_gradU_W(self):
+    #     self._test(formulation=StretchingFormulation.GRAD_UW, is_inplace=True)
+    # def test_stretching_3D_inplace_gradU_T_W(self):
+    #     self._test(formulation=StretchingFormulation.GRAD_UW_T, is_inplace=True)
+    # def test_stretching_3D_inplace_mixed_gradient(self):
+    #     As = (None, sm.Rational(1,3))
+    #     self._test(formulation=StretchingFormulation.MIXED_GRAD_UW, is_inplace=True, As=As)
+
+    # def test_stretching_3D_out_of_place_conservative(self):
+    #     self._test(formulation=StretchingFormulation.CONSERVATIVE, is_inplace=False)
+    # def test_stretching_3D_out_of_place_gradU_W(self):
+    #     self._test(formulation=StretchingFormulation.GRAD_UW, is_inplace=False)
+    # def test_stretching_3D_out_of_place_gradU_T_W(self):
+    #     self._test(formulation=StretchingFormulation.GRAD_UW_T, is_inplace=False)
+    # def test_stretching_3D_out_of_place_mixed_gradient(self):
+    #     As = (None, sm.Rational(1,3))
+    #     self._test(formulation=StretchingFormulation.MIXED_GRAD_UW, is_inplace=False, As=As)
 
     def perform_tests(self):
         # self.test_stretching_3D_inplace_gradU_W()
diff --git a/hysop/operator/tests/test_fd_derivative.py b/hysop/operator/tests/test_fd_derivative.py
index 334f8bed9b79800878d7b2c700bdd2d5f317d1bf..6bca7fbe37e20fb9af0b529a08e95b6cfccf175a 100644
--- a/hysop/operator/tests/test_fd_derivative.py
+++ b/hysop/operator/tests/test_fd_derivative.py
@@ -25,22 +25,22 @@ from hysop import Field, Box
 class TestFiniteDifferencesDerivative(object):
 
     @classmethod
-    def setup_class(cls, 
+    def setup_class(cls,
             enable_extra_tests=__ENABLE_LONG_TESTS__,
             enable_debug_mode=False):
 
         IO.set_default_path('/tmp/hysop_tests/test_fd_derivative')
-        
+
         cls.size_min = 23
         cls.size_max = 35
-        
+
         cls.enable_extra_tests = enable_extra_tests
         cls.enable_debug_mode  = enable_debug_mode
-        
+
         cls.t = ScalarParameter(name='t', dtype=HYSOP_REAL)
-        
+
         cls.build_analytic_expressions()
-    
+
     @classmethod
     def build_analytic_expressions(cls):
         from hysop.tools.sympy_utils import enable_pretty_printing
@@ -48,7 +48,7 @@ class TestFiniteDifferencesDerivative(object):
         from hysop.symbolic.frame import SymbolicFrame
         from hysop.symbolic.field import curl, laplacian
         enable_pretty_printing()
-        
+
         #at this moment we just test a periodic solver
         analytic_expressions = {}
         analytic_functions = {}
@@ -59,7 +59,7 @@ class TestFiniteDifferencesDerivative(object):
             def gen_Fi():
                 kis = tuple(random.randint(1,5)*2*sm.pi for _ in xrange(dim))
                 qis = tuple(npw.random.rand(dim).round(decimals=3).tolist())
-                basis = tuple( (sm.cos(ki*xi+qi), sm.sin(ki*xi+qi)) 
+                basis = tuple( (sm.cos(ki*xi+qi), sm.sin(ki*xi+qi))
                         for (ki,qi,xi) in zip(kis, qis, coords) )
                 F0 = sm.Integer(1) / (sm.Integer(1) + npw.random.randint(1,5)*cls.t.s)
                 for i in xrange(dim):
@@ -69,18 +69,18 @@ class TestFiniteDifferencesDerivative(object):
             def gen_F():
                 F = npw.round(random.random(),2)*gen_Fi() - npw.round(random.random(),2)*gen_Fi()
                 return F
-            
+
             params = coords + (cls.t.s,)
             Fs = npw.asarray([gen_F(), gen_F(), gen_F()], dtype=object).view(TensorBase)
             fFs = tuple(sm.lambdify(params, F) for F in Fs)
-            
+
             dFs = npw.empty(shape=(3, dim), dtype=object)
             fdFs = npw.empty_like(dFs)
             for (i, Fi) in enumerate(Fs):
                 for (j, xj) in enumerate(coords):
                     dFij = Fi.diff(xj)
                     dFs[i,j] = dFij
-                    fdFs[i,j] = sm.lambdify(params, dFij) 
+                    fdFs[i,j] = sm.lambdify(params, dFij)
 
             analytic_expressions[dim] = {'F':Fs, 'gradF':dFs}
             analytic_functions[dim]   = {'F':fFs, 'gradF':fdFs}
@@ -116,17 +116,17 @@ class TestFiniteDifferencesDerivative(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())
-        
+
         domain = Box(length=(1,)*dim)
         F = Field(domain=domain, name='F', dtype=dtype,
                 nb_components=3, register_object=False)
         gradF = F.gradient(register_object=False)
-        
-        self._test_one(shape=shape, dim=dim, dtype=dtype, 
+
+        self._test_one(shape=shape, dim=dim, dtype=dtype,
                 domain=domain, F=F, gradF=gradF)
 
-    
-    def _test_one(self, shape, dim, dtype, 
+
+    def _test_one(self, shape, dim, dtype,
             domain, F, gradF):
 
         print '\nTesting {}D Gradient: dtype={} shape={}'.format(
@@ -139,15 +139,15 @@ class TestFiniteDifferencesDerivative(object):
         for (i,j) in npw.ndindex(*self.analytic_expressions[dim]['gradF'].shape):
             dFij = self.analytic_expressions[dim]['gradF'][i,j]
             print '   *{p}F{}/{p}x{}(x,t) = {}'.format(subscript(i).encode('utf-8'),
-                                                       subscript(j).encode('utf-8'), 
-                                                       dFij, 
+                                                       subscript(j).encode('utf-8'),
+                                                       dFij,
                                                        p=partial.encode('utf-8'))
         self.t.value = random.random()
         print ' >Parameter t has been set to {}.'.format(self.t())
         print ' >Testing all implementations:'
 
         implementations = FiniteDifferencesSpaceDerivative.implementations()
-       
+
         variables = { F:shape, gradF: shape }
         fns    = self.analytic_functions[dim]['F']
         dfns   = tuple(self.analytic_functions[dim]['gradF'].ravel().tolist())
@@ -156,7 +156,7 @@ class TestFiniteDifferencesDerivative(object):
 
         def iter_impl(impl):
             base_kwds = dict(F=F, gradF=gradF, variables=variables,
-                             implementation=impl, 
+                             implementation=impl,
                              method={SpaceDiscretization: 8})
             if impl is Implementation.PYTHON:
                 msg='   *Python: '
@@ -168,12 +168,12 @@ class TestFiniteDifferencesDerivative(object):
                 msg='   *Opencl: '
                 print msg
                 for cl_env in iter_clenv():
-                    msg='      >platform {}, device {}:'.format(cl_env.platform.name.strip(), 
+                    msg='      >platform {}, device {}:'.format(cl_env.platform.name.strip(),
                                                           cl_env.device.name.strip())
                     print msg,
                     op0 = Gradient(cl_env=cl_env, **base_kwds)
                     op1 = ForceTopologyState(fields=gradF, variables={gradF: shape},
-                                                  backend=Backend.OPENCL, 
+                                                  backend=Backend.OPENCL,
                                                   extra_kwds={'cl_env':cl_env})
                     op = op0.to_graph()
                     op.push_nodes(op1)
@@ -183,8 +183,8 @@ class TestFiniteDifferencesDerivative(object):
             else:
                 msg='Unknown implementation to test {}.'.format(impl)
                 raise NotImplementedError(msg)
-        
-        
+
+
         # Compare to analytic solution
         Fref = None
         for impl in implementations:
@@ -204,12 +204,12 @@ class TestFiniteDifferencesDerivative(object):
 
                 dgradF.initialize(self.__random_init)
                 op.apply()
-                
+
                 Fout     = tuple( data.get().handle.copy() for data in dF.data )
                 gradFout = tuple( data.get().handle.copy() for data in dgradF.data )
-                
+
                 self._check_output(impl, op, Fref, gradFref, Fout, gradFout)
-    
+
     @classmethod
     def _check_output(cls, impl, op, Fref, gradFref, Fout, gradFout):
         check_instance(Fref, tuple, values=npw.ndarray)
@@ -222,7 +222,7 @@ class TestFiniteDifferencesDerivative(object):
                 iname = '{}{}'.format(name,i)
                 assert fout.dtype == fref.dtype, iname
                 assert fout.shape == fref.shape, iname
-                
+
                 eps  = npw.finfo(fout.dtype).eps
                 dist = npw.abs(fout-fref)
                 dL1 = npw.nansum(dist) / dist.size
@@ -231,7 +231,7 @@ class TestFiniteDifferencesDerivative(object):
                     continue
                 has_nan = npw.any(npw.isnan(fout))
                 has_inf = npw.any(npw.isinf(fout))
-                
+
                 print
                 print
                 print 'Test output comparisson for {} failed for component {}:'.format(name, i)
@@ -265,9 +265,9 @@ class TestFiniteDifferencesDerivative(object):
                             print w
                             print
                     print
-                
+
                 msg = 'Test failed for {} on component {} for implementation {}.'.format(name, i, impl)
-                raise RuntimeError(msg) 
+                raise RuntimeError(msg)
 
 
 
@@ -297,15 +297,15 @@ class TestFiniteDifferencesDerivative(object):
             self.test_2d_float64()
             if __ENABLE_LONG_TESTS__:
                 self.test_3d_float64()
-    
+
 if __name__ == '__main__':
-    TestFiniteDifferencesDerivative.setup_class(enable_extra_tests=False, 
+    TestFiniteDifferencesDerivative.setup_class(enable_extra_tests=False,
                              enable_debug_mode=False)
-    
+
     test = TestFiniteDifferencesDerivative()
 
-    with printoptions(threshold=10000, linewidth=1000, 
-                      nanstr='nan', infstr='inf', 
+    with printoptions(threshold=10000, linewidth=1000,
+                      nanstr='nan', infstr='inf',
                       formatter={'float': lambda x: '{:>6.2f}'.format(x)}):
         test.perform_tests()
 
diff --git a/hysop/operator/tests/test_poisson.py b/hysop/operator/tests/test_poisson.py
index d4cf58ea095dc926c07c905c3fc2d392fb5ebc41..c874a5e5b4ebf6f9fbe9fa18ceef4876aafec785 100644
--- a/hysop/operator/tests/test_poisson.py
+++ b/hysop/operator/tests/test_poisson.py
@@ -30,7 +30,7 @@ class TestPoissonOperator(object):
 
         cls.enable_extra_tests = enable_extra_tests
         cls.enable_debug_mode  = enable_debug_mode
-        
+
         from hysop.tools.sympy_utils import enable_pretty_printing
         enable_pretty_printing()
 
@@ -39,23 +39,23 @@ class TestPoissonOperator(object):
         pass
 
     @classmethod
-    def build_analytic_solutions(cls, polynomial, 
+    def build_analytic_solutions(cls, polynomial,
                                       dim, nb_components,
                                       lboundaries, rboundaries,
                                       origin, end):
         from hysop.symbolic.base  import TensorBase
         from hysop.symbolic.frame import SymbolicFrame
         from hysop.symbolic.field import laplacian
-        
+
         frame = SymbolicFrame(dim=dim)
         coords = frame.coords
-                
+
         def gen_psi():
             psis = ()
             for i in xrange(nb_components):
                 if polynomial:
                     psi, y = make_multivariate_polynomial(origin, end,
-                                                        lboundaries, rboundaries, 
+                                                        lboundaries, rboundaries,
                                                         10, 4)
                 else:
                     psi, y = make_multivariate_trigonometric_polynomial(origin, end,
@@ -109,7 +109,7 @@ class TestPoissonOperator(object):
             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))) )
-       
+
         domain_boundaries = list(domain_boundary_iterator(dim=dim))
         periodic = domain_boundaries[0]
         domain_boundaries = domain_boundaries[1:]
@@ -127,7 +127,7 @@ class TestPoissonOperator(object):
             W = Field(domain=domain, name='W', dtype=dtype,
                     nb_components=nb_components, register_object=False)
 
-            self._test_one(shape=shape, dim=dim, dtype=dtype, 
+            self._test_one(shape=shape, dim=dim, dtype=dtype,
                     domain=domain, Psi=Psi, W=W,
                     polynomial=polynomial, nb_components=nb_components)
             if (max_runs is not None) and (i==max_runs):
@@ -147,15 +147,15 @@ class TestPoissonOperator(object):
 
     def _test_one(self, shape, dim, dtype,
             domain, Psi, W, polynomial, nb_components):
-        
+
         (analytic_expressions, analytic_functions) = \
             self.build_analytic_solutions(
-                dim=dim, nb_components=nb_components, polynomial=polynomial, 
+                dim=dim, nb_components=nb_components, polynomial=polynomial,
                 lboundaries=W.lboundaries[::-1], # => boundaries in variable order x0,...,xn
                 rboundaries=W.rboundaries[::-1],
                 origin=domain.origin[::-1],
                 end=domain.end[::-1])
-        
+
         def format_expr(e):
             return truncate_expr(round_expr(e, 3), 80)
 
@@ -341,13 +341,13 @@ class TestPoissonOperator(object):
         self._test(dim=3, dtype=npw.float64, **kwds)
     def test_4d_float64(self, **kwds):
         self._test(dim=4, dtype=npw.float64, **kwds)
-    
-    def test_polynomial_1d_float32(self, **kwds):
-        self._test(dim=1, dtype=npw.float32, polynomial=True, **kwds)
-    def test_polynomial_2d_float32(self, **kwds):
-        self._test(dim=2, dtype=npw.float32, polynomial=True, **kwds)
-    def test_polynomial_3d_float32(self, **kwds):
-        self._test(dim=3, dtype=npw.float32, polynomial=True, **kwds)
+
+    # def test_polynomial_1d_float32(self, **kwds):
+    #     self._test(dim=1, dtype=npw.float32, polynomial=True, **kwds)
+    # def test_polynomial_2d_float32(self, **kwds):
+    #     self._test(dim=2, dtype=npw.float32, polynomial=True, **kwds)
+    # def test_polynomial_3d_float32(self, **kwds):
+    #     self._test(dim=3, dtype=npw.float32, polynomial=True, **kwds)
 
 
     def perform_tests(self):
diff --git a/hysop/operator/tests/test_poisson_curl.py b/hysop/operator/tests/test_poisson_curl.py
index 33e779441918ce96aa44add9969a752c5ef8bb09..82d7bdd3de2bdee6dc3949c98cb527056dbcbb44 100644
--- a/hysop/operator/tests/test_poisson_curl.py
+++ b/hysop/operator/tests/test_poisson_curl.py
@@ -33,17 +33,17 @@ class TestPoissonCurlOperator(object):
 
         cls.enable_extra_tests = enable_extra_tests
         cls.enable_debug_mode  = enable_debug_mode
-        
+
         from hysop.tools.sympy_utils import enable_pretty_printing
         enable_pretty_printing()
-    
+
     @classmethod
     def teardown_class(cls):
         pass
 
 
     @classmethod
-    def build_analytic_solutions(cls, polynomial, 
+    def build_analytic_solutions(cls, polynomial,
                                       dim, nb_components,
                                       lboundaries, rboundaries,
                                       origin, end):
@@ -53,7 +53,7 @@ class TestPoissonCurlOperator(object):
 
         assert len(lboundaries)==nb_components
         assert len(rboundaries)==nb_components
-        
+
         frame = SymbolicFrame(dim=dim)
         coords = frame.coords
         def gen_psi():
@@ -61,7 +61,7 @@ class TestPoissonCurlOperator(object):
             for i in xrange(nb_components):
                 if polynomial:
                     psi, y = make_multivariate_polynomial(origin, end,
-                                                        lboundaries[i], rboundaries[i], 
+                                                        lboundaries[i], rboundaries[i],
                                                         10, 4)
                 else:
                     psi, y = make_multivariate_trigonometric_polynomial(origin, end,
@@ -69,7 +69,7 @@ class TestPoissonCurlOperator(object):
                 psi = psi.xreplace({yi: xi for (yi,xi) in zip(y, coords)})
                 psis += (psi,)
             return npw.asarray(psis).view(TensorBase)
-        
+
         Psis  = gen_psi()
         Ws    = npw.atleast_1d(-laplacian(Psis, frame))
         Us    = curl(Psis, frame)
@@ -82,7 +82,7 @@ class TestPoissonCurlOperator(object):
         analytic_functions   = {'Psi':fPsis, 'W':fWs, 'U':fUs}
         return (analytic_expressions, analytic_functions)
 
-    
+
     @staticmethod
     def __random_init(data, coords, dtype):
         for d in data:
@@ -91,7 +91,7 @@ class TestPoissonCurlOperator(object):
             else:
                 msg = 'Unknown dtype {}.'.format(d.dtype)
                 raise NotImplementedError(msg)
-    
+
     @staticmethod
     def __analytic_init(data, coords, dtype, fns):
         assert len(fns) == len(data)
@@ -115,7 +115,7 @@ class TestPoissonCurlOperator(object):
             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))) )
-        
+
         domain_boundaries = list(domain_boundary_iterator(dim=dim))
         periodic = domain_boundaries[0]
         domain_boundaries = domain_boundaries[1:]
@@ -130,7 +130,7 @@ class TestPoissonCurlOperator(object):
             U = VelocityField(domain=domain, dtype=dtype)
             W = VorticityField(velocity=U, dtype=dtype)
 
-            self._test_one(shape=shape, dim=dim, dtype=dtype, 
+            self._test_one(shape=shape, dim=dim, dtype=dtype,
                     domain=domain, W=W, U=U, polynomial=polynomial)
             if (max_runs is not None) and (i==max_runs):
                 missing = ((4**(dim+1) - 1) / 3) - i
@@ -148,15 +148,15 @@ class TestPoissonCurlOperator(object):
 
     def _test_one(self, shape, dim, dtype,
             domain, U, W, polynomial):
-        
+
         (analytic_expressions, analytic_functions) = \
             self.build_analytic_solutions(
-                dim=dim, nb_components=W.nb_components, polynomial=polynomial, 
+                dim=dim, nb_components=W.nb_components, polynomial=polynomial,
                 lboundaries=[Wi.lboundaries[::-1] for Wi in W.fields], # => boundaries in variable order x0,...,xn
                 rboundaries=[Wi.rboundaries[::-1] for Wi in W.fields],
                 origin=domain.origin[::-1],
                 end=domain.end[::-1])
-        
+
         def format_expr(e):
             return truncate_expr(round_expr(e, 3), 80)
 
@@ -221,10 +221,10 @@ class TestPoissonCurlOperator(object):
                 name='{}_{}'.format(impl, i)
 
                 op = op.build()
-                
+
                 dw = op.get_input_discrete_field(W).as_contiguous_dfield()
                 du = op.get_output_discrete_field(U).as_contiguous_dfield()
-                
+
                 dw.initialize(self.__analytic_init, dtype=dtype,
                                     fns=analytic_functions['W'])
 
diff --git a/hysop/operator/tests/test_scales_advection.py b/hysop/operator/tests/test_scales_advection.py
index 77672ac4354ddddbf653c06a3a981618c4a43d48..3b69adde2a16a9de2a0f3087bd300ac9f18a7b19 100644
--- a/hysop/operator/tests/test_scales_advection.py
+++ b/hysop/operator/tests/test_scales_advection.py
@@ -169,7 +169,7 @@ class TestScalesAdvectionOperator(object):
                     #print 'SWITCHING TO AXES {}'.format(axes)
                     ref_outputs = reference_fields.setdefault(axes, {})
                     napplies = 10
-                    Vi = npw.asarray([+1.0 if (i in axes) else +0.0 
+                    Vi = npw.asarray([+1.0 if (i in axes) else +0.0
                                         for i in xrange(dim)], dtype=dtype)
 
                     dvin  = graph.get_input_discrete_field(vin).as_contiguous_dfield()
@@ -213,7 +213,7 @@ class TestScalesAdvectionOperator(object):
                                     print dsref.sdata[dsref.compute_slices]
                                     print 'DSREF - DSOUT'
                                     print (dsout.sdata[dsout.compute_slices].get() - dsref.sdata[dsref.compute_slices].get())
-                                    msg='Test failed with V={}, k={}, dxk={}, inter-field L2 distances are {}.' 
+                                    msg='Test failed with V={}, k={}, dxk={}, inter-field L2 distances are {}.'
                                     msg=msg.format(Vi, k, to_tuple(dxk, cast=float), to_tuple(d, cast=float))
                                     raise RuntimeError(msg)
                                 ref_outputs[k] = output
@@ -265,6 +265,9 @@ class TestScalesAdvectionOperator(object):
         self._test(dim=3, is_inplace=True)
         print
 
+    def test_3D(self):
+        self._test(dim=3, is_inplace=True)
+
 
 if __name__ == '__main__':
     import hysop
diff --git a/hysop/operator/tests/test_solenoidal_projection.py b/hysop/operator/tests/test_solenoidal_projection.py
index 49da261b081919cedf36884e4f156c74a8ec1664..234a1e8459a7459b8b8fa9b201627ab092d3930c 100644
--- a/hysop/operator/tests/test_solenoidal_projection.py
+++ b/hysop/operator/tests/test_solenoidal_projection.py
@@ -1,4 +1,3 @@
-
 import random, primefac, scipy
 from hysop.deps import it, sm, random, np
 from hysop.constants import HYSOP_REAL, Implementation, BoxBoundaryCondition
@@ -20,18 +19,18 @@ from hysop import Field, Box
 class TestSolenoidalProjectionOperator(object):
 
     @classmethod
-    def setup_class(cls, 
+    def setup_class(cls,
             enable_extra_tests=__ENABLE_LONG_TESTS__,
             enable_debug_mode=False):
 
         IO.set_default_path('/tmp/hysop_tests/test_solenoidal_projection')
-        
+
         cls.size_min = 8
         cls.size_max = 16
-        
+
         cls.enable_extra_tests = enable_extra_tests
         cls.enable_debug_mode  = enable_debug_mode
-        
+
         from hysop.tools.sympy_utils import enable_pretty_printing
         enable_pretty_printing()
 
@@ -43,17 +42,17 @@ class TestSolenoidalProjectionOperator(object):
 
         assert len(lboundaries)==4
         assert len(rboundaries)==4
-        
+
         frame = SymbolicFrame(dim=3)
         coords = frame.coords
-        
+
         def gen_fn(nb_components, left_boundaries, right_boundaries):
             assert len(left_boundaries)==len(right_boundaries)==nb_components
             fns = ()
             for i in xrange(nb_components):
                 if polynomial:
                     fn, y = make_multivariate_polynomial(origin, end,
-                                                        left_boundaries[i], right_boundaries[i], 
+                                                        left_boundaries[i], right_boundaries[i],
                                                         10, 2)
                 else:
                     fn, y = make_multivariate_trigonometric_polynomial(origin, end,
@@ -64,11 +63,11 @@ class TestSolenoidalProjectionOperator(object):
 
         analytic_solutions = {}
         analytic_functions = {}
-        
-        psis = gen_fn(nb_components=3, 
+
+        psis = gen_fn(nb_components=3,
                       left_boundaries=lboundaries[:3],
                       right_boundaries=rboundaries[:3])
-        phis = gen_fn(nb_components=1, 
+        phis = gen_fn(nb_components=1,
                       left_boundaries=lboundaries[3:],
                       right_boundaries=rboundaries[3:])
 
@@ -85,7 +84,7 @@ class TestSolenoidalProjectionOperator(object):
         fUs    = tuple(sm.lambdify(coords, Ui) for Ui in Us)
         fdivUs = tuple(sm.lambdify(coords, Ui) for Ui in divUs)
 
-        sols = {'U0s': U0s, 'U1s': U1s, 'Us':Us, 
+        sols = {'U0s': U0s, 'U1s': U1s, 'Us':Us,
                 'divUs':divUs, 'divU0s':divU0s, 'divU1s':divU1s}
         lambdified_sols = {'U0s': fU0s, 'U1s': fU1s, 'Us':fUs, 'divUs':fdivUs}
 
@@ -98,14 +97,14 @@ class TestSolenoidalProjectionOperator(object):
     def teardown_class(cls):
         pass
 
-    
+
     def _test(self, dtype, max_runs=5,
             polynomial=False, size_min=None, size_max=None):
         enable_extra_tests = self.enable_extra_tests
 
         size_min = first_not_None(size_min, self.size_min)
         size_max = first_not_None(size_max, self.size_max)
-        
+
         dim = 3
         valid_factors = {2,3,5,7,11,13}
         factors = {1}
@@ -114,19 +113,19 @@ class TestSolenoidalProjectionOperator(object):
             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))) )
-        
+
         domain_boundaries = list(domain_boundary_iterator(dim=dim))
         periodic = domain_boundaries[0]
         domain_boundaries = domain_boundaries[1:]
         random.shuffle(domain_boundaries)
         domain_boundaries.insert(0, periodic)
-        
+
         for i, (lboundaries, rboundaries) in enumerate(domain_boundaries, 1):
             domain = Box(origin=(npw.random.rand(dim)-0.5),
                          length=(0.5+npw.random.rand(dim))*2*npw.pi,
                          lboundaries=lboundaries,
                          rboundaries=rboundaries)
-        
+
             U     = VelocityField(domain=domain, name='U', dtype=dtype)
             U0    = VelocityField(domain=domain, name='U0', dtype=dtype)
             U1    = VelocityField(domain=domain, name='U1', dtype=dtype)
@@ -156,12 +155,12 @@ class TestSolenoidalProjectionOperator(object):
         assert len(fns) == len(data)
         for (d,fn,coord) in zip(data,fns,coords):
             d[...] = npw.asarray(fn(*coord)).astype(dtype)
-    
+
     @classmethod
     def __zero_init(cls, data, coords, dtype):
         for d in data:
             d[...] = 0
-    
+
     @staticmethod
     def __random_init(data, coords, dtype):
         shape = data[0].shape
@@ -174,36 +173,36 @@ class TestSolenoidalProjectionOperator(object):
 
     def _test_one(self, shape, dtype, polynomial,
             domain, U, U0, U1, divU, divU0, divU1):
-        
+
         print '\nTesting {}D SolenoidalProjection: dtype={} shape={}, polynomial={}, bc=[{}]'.format(
                 3, dtype.__name__, shape, polynomial, domain.format_boundaries())
         print ' >Building U = U0 + U1 = curl(Psi) + grad(Phi)...'
-            
+
         # U = curl(Psi) + grad(Phi)
         Psi = VorticityField(velocity=U, name='Psi')
         assert all(all(psic.lboundaries == u.lboundaries) for (psic,u) in zip(Psi.curl().fields, U.fields))
         assert all(all(psic.rboundaries == u.rboundaries) for (psic,u) in zip(Psi.curl().fields, U.fields))
-        
+
         Phi = U.div(name='Phi')
         assert all(all(phig.lboundaries == u.lboundaries) for (phig,u) in zip(Phi.gradient().fields, U.fields))
         assert all(all(phig.rboundaries == u.rboundaries) for (phig,u) in zip(Phi.gradient().fields, U.fields))
-            
+
         lboundaries = tuple(psi.lboundaries[::-1] for psi in Psi.fields) + (Phi.lboundaries[::-1],)
         rboundaries = tuple(psi.rboundaries[::-1] for psi in Psi.fields) + (Phi.rboundaries[::-1],)
 
         (analytic_solutions, analytic_functions) = \
             self.build_analytic_solutions(
-                polynomial=polynomial, 
+                polynomial=polynomial,
                 lboundaries=lboundaries, # => boundaries in variable order x0,...,xn
                 rboundaries=rboundaries,
                 origin=domain.origin[::-1],
                 end=domain.end[::-1])
-        
+
         del Psi
         del Phi
         del lboundaries
         del rboundaries
-        
+
         def format_expr(e):
             return truncate_expr(round_expr(e, 3), 80)
 
@@ -222,14 +221,14 @@ class TestSolenoidalProjectionOperator(object):
         print ' >Testing all available implementations:'
 
         implementations = SolenoidalProjection.implementations()
-       
+
         # Compute reference solution
         variables = { U:shape, divU: shape, U0:shape, divU0:shape}
 
         def iter_impl(impl):
             base_kwds = dict(input_field=U, input_field_div=divU,
                              output_field=U0, output_field_div=divU0,
-                             variables=variables, implementation=impl, 
+                             variables=variables, implementation=impl,
                              name='projection_{}'.format(str(impl).lower()))
             if impl is Implementation.PYTHON:
                 msg='   *Python FFTW: '
@@ -239,7 +238,7 @@ class TestSolenoidalProjectionOperator(object):
                 msg='   *OpenCl CLFFT: '
                 print msg
                 for cl_env in iter_clenv():
-                    msg='     |platform {}, device {}'.format(cl_env.platform.name.strip(), 
+                    msg='     |platform {}, device {}'.format(cl_env.platform.name.strip(),
                                                           cl_env.device.name.strip())
                     print msg,
                     yield SolenoidalProjection(cl_env=cl_env, **base_kwds)
@@ -255,7 +254,7 @@ class TestSolenoidalProjectionOperator(object):
             else:
                 msg='Unknown implementation to test {}.'.format(impl)
                 raise NotImplementedError(msg)
-        
+
         # Compare to analytic solution
         Uref, U0ref = None, None
         for impl in implementations:
@@ -266,14 +265,14 @@ class TestSolenoidalProjectionOperator(object):
                 du0     = op.get_output_discrete_field(U0)
                 du_div  = op.get_output_discrete_field(divU)
                 du0_div = op.get_output_discrete_field(divU0)
-                
-                du.initialize(self.__analytic_init, dtype=dtype, 
+
+                du.initialize(self.__analytic_init, dtype=dtype,
                                     fns=analytic_functions['Us'])
 
                 if (Uref is None):
-                    du0.initialize(self.__analytic_init, dtype=dtype, 
+                    du0.initialize(self.__analytic_init, dtype=dtype,
                                     fns=analytic_functions['U0s'])
-                    du_div.initialize(self.__analytic_init, dtype=dtype, 
+                    du_div.initialize(self.__analytic_init, dtype=dtype,
                                     fns=analytic_functions['divUs'])
                     du0_div.initialize(self.__zero_init, dtype=dtype)
                     Uref     = tuple( data.get().handle.copy() for data in du.data )
@@ -284,9 +283,9 @@ class TestSolenoidalProjectionOperator(object):
                 du0.initialize(self.__random_init, dtype=dtype)
                 du_div.initialize(self.__random_init, dtype=dtype)
                 du0_div.initialize(self.__random_init, dtype=dtype)
-                
+
                 op.apply(simulation=None)
-                
+
                 Uout     = tuple( data.get().handle.copy() for data in du.data )
                 U0out    = tuple( data.get().handle.copy() for data in du0.data )
                 divUout  = tuple( data.get().handle.copy() for data in du_div.data )
@@ -297,13 +296,13 @@ class TestSolenoidalProjectionOperator(object):
                         npw.sqrt(npw.sum(divUout[0]**2)*s),
                         npw.sqrt(npw.sum(divU0out[0]**2)*s)),
 
-                self._check_output(impl, op, Uref, divUref, U0ref, divU0ref, 
+                self._check_output(impl, op, Uref, divUref, U0ref, divU0ref,
                                              Uout, divUout, U0out, divU0out)
                 print
-    
+
     @classmethod
-    def _check_output(cls, impl, op, 
-                      Uref, divUref, U0ref, divU0ref, 
+    def _check_output(cls, impl, op,
+                      Uref, divUref, U0ref, divU0ref,
                       Uout, divUout, U0out, divU0out):
         check_instance(Uref,     tuple, values=npw.ndarray, size=3)
         check_instance(U0ref,    tuple, values=npw.ndarray, size=3)
@@ -328,16 +327,16 @@ class TestSolenoidalProjectionOperator(object):
                     print
                     msg = msg0.format(iname)
                     raise ValueError(msg)
-        
-        for (out_buffers, ref_buffers, name) in zip((Uout, U0out, divUout, divU0out), 
-                                                    (Uref, U0ref, divUref, divU0ref), 
+
+        for (out_buffers, ref_buffers, name) in zip((Uout, U0out, divUout, divU0out),
+                                                    (Uref, U0ref, divUref, divU0ref),
                                                     ('U', 'U0', 'divU', 'divU0')):
             print '| {}=('.format(name),
             for i, (fout,fref) in enumerate(zip(out_buffers, ref_buffers)):
                 iname = '{}{}'.format(name,i)
                 assert fout.dtype == fref.dtype, iname
                 assert fout.shape == fref.shape, iname
-                
+
                 eps  = npw.finfo(fout.dtype).eps
                 dist = npw.abs(fout-fref)
                 dinf = npw.max(dist)
@@ -385,10 +384,10 @@ class TestSolenoidalProjectionOperator(object):
                             print
                         print
                     print
-                
+
                 msg = 'Test failed for {} on component {} for implementation {}.'
                 msg = msg.format(name, i, impl)
-                raise RuntimeError(msg) 
+                raise RuntimeError(msg)
             print ')',
 
 
@@ -403,15 +402,15 @@ class TestSolenoidalProjectionOperator(object):
             self.test_3d_float32(max_runs=max_3d_runs)
         if __ENABLE_LONG_TESTS__ or (HYSOP_REAL==npw.float64):
             self.test_3d_float64(max_runs=max_3d_runs)
-    
+
 if __name__ == '__main__':
-    TestSolenoidalProjectionOperator.setup_class(enable_extra_tests=False, 
+    TestSolenoidalProjectionOperator.setup_class(enable_extra_tests=False,
                                       enable_debug_mode=False)
-    
+
     test = TestSolenoidalProjectionOperator()
 
-    with printoptions(threshold=10000, linewidth=240, 
-                      nanstr='nan', infstr='inf', 
+    with printoptions(threshold=10000, linewidth=240,
+                      nanstr='nan', infstr='inf',
                       formatter={'float': lambda x: '{:>6.2f}'.format(x)}):
         test.perform_tests()
 
diff --git a/hysop/operator/tests/test_spectral_derivative.py b/hysop/operator/tests/test_spectral_derivative.py
index 5702c7f73ec2256b14f911062cc23b97f4a072aa..d5004381ad6c3c985b3c38bb230217fab62d280d 100644
--- a/hysop/operator/tests/test_spectral_derivative.py
+++ b/hysop/operator/tests/test_spectral_derivative.py
@@ -24,20 +24,20 @@ from hysop import Field, Box
 class TestSpectralDerivative(object):
 
     @classmethod
-    def setup_class(cls, 
+    def setup_class(cls,
             enable_extra_tests=__ENABLE_LONG_TESTS__,
             enable_debug_mode=False):
 
         IO.set_default_path('/tmp/hysop_tests/test_spectral_derivative')
-        
+
         cls.size_min = 8
         cls.size_max = 16
-        
+
         cls.enable_extra_tests = enable_extra_tests
         cls.enable_debug_mode  = enable_debug_mode
-        
+
         cls.t = ScalarParameter(name='t', dtype=HYSOP_REAL)
-    
+
     @classmethod
     def build_analytic_expressions(cls, polynomial, dim, max_derivative,
                                         lboundaries, rboundaries, origin, end):
@@ -46,15 +46,15 @@ class TestSpectralDerivative(object):
         from hysop.symbolic.frame import SymbolicFrame
         from hysop.symbolic.field import curl, laplacian
         enable_pretty_printing()
-        
+
         frame = SymbolicFrame(dim=dim)
         coords = frame.coords
         params = coords + (cls.t.s,)
-        
+
         def gen_F():
             if polynomial:
                 f, y = make_multivariate_polynomial(origin, end,
-                                                    lboundaries, rboundaries, 
+                                                    lboundaries, rboundaries,
                                                     10, 4)
             else:
                 f, y = make_multivariate_trigonometric_polynomial(origin, end,
@@ -65,7 +65,7 @@ class TestSpectralDerivative(object):
 
         F  = gen_F()
         fF = sm.lambdify(params, F)
-    
+
         dFs  = {}
         fdFs = {}
         symbolic_dvars = {}
@@ -80,7 +80,7 @@ class TestSpectralDerivative(object):
                     continue
                 dF = dF.diff(ci,i)
             dFs[idx]  = dF
-            fdFs[idx] = sm.lambdify(params, dF) 
+            fdFs[idx] = sm.lambdify(params, dF)
 
         analytic_expressions = {'F':F,  'dF':dFs}
         analytic_functions = {'F':fF, 'dF':fdFs}
@@ -110,18 +110,18 @@ class TestSpectralDerivative(object):
     def _test(self, dim, dtype, polynomial, max_derivative=2,
             size_min=None, size_max=None, max_runs=None):
         enable_extra_tests = self.enable_extra_tests
-        
+
         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())
-        
+
         domain_boundaries = list(domain_boundary_iterator(dim=dim))
         periodic = domain_boundaries[0]
         domain_boundaries = domain_boundaries[1:]
         random.shuffle(domain_boundaries)
         domain_boundaries.insert(0, periodic)
-        
+
         for i, (lboundaries, rboundaries) in enumerate(domain_boundaries, 1):
             domain = Box(origin=(npw.random.rand(dim)-0.5),
                          length=(npw.random.rand(dim)+0.5)*2*npw.pi,
@@ -130,11 +130,11 @@ class TestSpectralDerivative(object):
 
             F = Field(domain=domain, name='F', dtype=dtype)
 
-            self._test_one(shape=shape, dim=dim, dtype=dtype, 
-                    domain=domain, F=F, 
-                    polynomial=polynomial, 
+            self._test_one(shape=shape, dim=dim, dtype=dtype,
+                    domain=domain, F=F,
+                    polynomial=polynomial,
                     max_derivative=max_derivative)
-            
+
             if (max_runs is not None) and (i==max_runs):
                 missing = ((4**(dim+1) - 1) / 3) - i
                 print
@@ -149,24 +149,24 @@ class TestSpectralDerivative(object):
             print
             print
 
-    
-    def _test_one(self, shape, dim, dtype, 
+
+    def _test_one(self, shape, dim, dtype,
             domain, F, polynomial, max_derivative):
-        
+
         implementations = SpectralSpaceDerivative.implementations()
 
         (symbolic_dvars, analytic_expressions, analytic_functions) = \
             self.build_analytic_expressions(
-                dim=dim, polynomial=polynomial, 
+                dim=dim, polynomial=polynomial,
                 max_derivative=max_derivative,
                 lboundaries=F.lboundaries[::-1], # => boundaries in variable order x0,...,xn
                 rboundaries=F.rboundaries[::-1],
                 origin=domain.origin[::-1],
                 end=domain.end[::-1])
-        
+
         Fs  = analytic_expressions['F']
         fFs = analytic_functions['F']
-        
+
         def format_expr(e):
             return truncate_expr(round_expr(e, 3), 80)
 
@@ -190,11 +190,11 @@ class TestSpectralDerivative(object):
             dFs   = analytic_expressions['dF'][idx]
             fdFs  = analytic_functions['dF'][idx]
             print '   *{}'.format(dF.pretty_name)
-            
+
             variables = { F:shape, dF: shape }
-        
+
             def iter_impl(impl):
-                base_kwds = dict(F=F, dF=dF, derivative=idx, 
+                base_kwds = dict(F=F, dF=dF, derivative=idx,
                                  variables=variables,
                                  implementation=impl,
                                  testing=True)
@@ -209,7 +209,7 @@ class TestSpectralDerivative(object):
                     print msg
                     for cl_env in iter_clenv():
                         msg='        >platform {}, device {}:'.format(
-                                                              cl_env.platform.name.strip(), 
+                                                              cl_env.platform.name.strip(),
                                                               cl_env.device.name.strip())
                         print msg,
                         op = SpectralSpaceDerivative(cl_env=cl_env, **base_kwds)
@@ -219,7 +219,7 @@ class TestSpectralDerivative(object):
                 else:
                     msg='Unknown implementation to test {}.'.format(impl)
                     raise NotImplementedError(msg)
-        
+
             # Compare to analytic solution
             Fref = None
             for impl in implementations:
@@ -239,12 +239,12 @@ class TestSpectralDerivative(object):
 
                     dFd.initialize(self.__random_init)
                     op.apply()
-                    
+
                     Fout  = tuple( data.get().handle.copy() for data in Fd.data )
                     dFout = tuple( data.get().handle.copy() for data in dFd.data )
-                    
+
                     self._check_output(impl, op, Fref, dFref, Fout, dFout, idx)
-    
+
     @classmethod
     def _check_output(cls, impl, op, Fref, dFref, Fout, dFout, idx):
         nidx = sum(idx)
@@ -253,19 +253,19 @@ class TestSpectralDerivative(object):
         check_instance(Fout,  tuple, values=npw.ndarray, size=len(Fref))
         check_instance(dFout, tuple, values=npw.ndarray, size=len(dFref))
 
-        for j,(out_buffers, ref_buffers, name) in enumerate(zip((Fout, dFout), 
-                                                                (Fref, dFref), 
+        for j,(out_buffers, ref_buffers, name) in enumerate(zip((Fout, dFout),
+                                                                (Fref, dFref),
                                                                 ('F', 'dF'))):
             for i, (fout,fref) in enumerate(zip(out_buffers, ref_buffers)):
                 iname = '{}{}'.format(name,i)
                 assert fout.dtype == fref.dtype, iname
                 assert fout.shape == fref.shape, iname
-                
+
                 assert not npw.any(npw.isnan(fref))
                 assert not npw.any(npw.isinf(fref))
                 has_nan = npw.any(npw.isnan(fout))
                 has_inf = npw.any(npw.isinf(fout))
-                
+
                 if (has_nan or has_inf):
                     pass
                 else:
@@ -279,7 +279,7 @@ class TestSpectralDerivative(object):
                         else:
                             print '{}eps, '.format(deps),
                         continue
-                
+
                 print
                 print
                 print 'Test output comparisson for {} failed for component {}:'.format(name, i)
@@ -314,31 +314,31 @@ class TestSpectralDerivative(object):
                             print w
                             print
                     print
-                
-                msg = 'Test failed for {} on component {} for implementation {}.'.format(name, 
+
+                msg = 'Test failed for {} on component {} for implementation {}.'.format(name,
                         i, impl)
-                raise RuntimeError(msg) 
+                raise RuntimeError(msg)
 
 
 
-    def test_1d_trigonometric_float32(self, **kwds):
-        self._test(dim=1, dtype=npw.float32, polynomial=False, **kwds)
-    def test_2d_trigonometric_float32(self, **kwds):
-        self._test(dim=2, dtype=npw.float32, polynomial=False, **kwds)
+    # def test_1d_trigonometric_float32(self, **kwds):
+    #     self._test(dim=1, dtype=npw.float32, polynomial=False, **kwds)
+    # def test_2d_trigonometric_float32(self, **kwds):
+    #     self._test(dim=2, dtype=npw.float32, polynomial=False, **kwds)
     def test_3d_trigonometric_float32(self, **kwds):
         self._test(dim=3, dtype=npw.float32, polynomial=False, **kwds)
 
-    def test_1d_trigonometric_float64(self, **kwds):
-        self._test(dim=1, dtype=npw.float64, polynomial=False, **kwds)
-    def test_2d_trigonometric_float64(self, **kwds):
-        self._test(dim=2, dtype=npw.float64, polynomial=False, **kwds)
+    # def test_1d_trigonometric_float64(self, **kwds):
+    #     self._test(dim=1, dtype=npw.float64, polynomial=False, **kwds)
+    # def test_2d_trigonometric_float64(self, **kwds):
+    #     self._test(dim=2, dtype=npw.float64, polynomial=False, **kwds)
     def test_3d_trigonometric_float64(self, **kwds):
         self._test(dim=3, dtype=npw.float64, polynomial=False, **kwds)
-    
-    def test_1d_polynomial_float32(self, **kwds):
-        self._test(dim=1, dtype=npw.float32, polynomial=True, **kwds)
-    def test_2d_polynomial_float32(self, **kwds):
-        self._test(dim=2, dtype=npw.float32, polynomial=True, **kwds)
+
+    # def test_1d_polynomial_float32(self, **kwds):
+    #     self._test(dim=1, dtype=npw.float32, polynomial=True, **kwds)
+    # def test_2d_polynomial_float32(self, **kwds):
+    #     self._test(dim=2, dtype=npw.float32, polynomial=True, **kwds)
     def test_3d_polynomial_float32(self, **kwds):
         self._test(dim=3, dtype=npw.float32, polynomial=True, **kwds)
 
@@ -358,14 +358,14 @@ class TestSpectralDerivative(object):
             self.test_1d_polynomial_float32(max_derivative=3)
             self.test_2d_polynomial_float32(max_derivative=2)
             self.test_3d_polynomial_float32(max_derivative=1)
-    
+
 if __name__ == '__main__':
-    TestSpectralDerivative.setup_class(enable_extra_tests=False, 
+    TestSpectralDerivative.setup_class(enable_extra_tests=False,
                                      enable_debug_mode=False)
-    
+
     test = TestSpectralDerivative()
 
-    with printoptions(threshold=10000, linewidth=1000, 
+    with printoptions(threshold=10000, linewidth=1000,
                       nanstr='nan', infstr='inf',
                       formatter={'float': lambda x: '{:>6.2f}'.format(x)}):
         test.perform_tests()