From 53f50dd0327b403afa4ecc2b48219332633f8a27 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Sun, 1 Nov 2020 21:11:44 +0100
Subject: [PATCH] many fixes

---
 README.md                                     |   3 +-
 ci/scripts/test.sh                            |  24 ++---
 hysop/__init__.py.in                          |   7 +-
 hysop/backend/device/autotunable_kernel.py    |   6 +-
 hysop/backend/device/codegen/base/codegen.py  |   2 +-
 .../device/codegen/base/struct_codegen.py     |   4 +-
 .../backend/device/codegen/base/variables.py  |  14 +--
 .../device/codegen/functions/apply_stencil.py |   4 +-
 .../device/codegen/functions/cache_load.py    |   2 +-
 .../device/codegen/functions/compute_index.py |   2 +-
 .../device/codegen/kernels/custom_symbolic.py |  27 ++---
 .../codegen/kernels/directional_advection.py  |   4 +-
 .../test_codegen_directional_stretching.py    |   2 +-
 .../kernels/tests/test_codegen_transpose.py   |   2 +-
 .../device/codegen/kernels/transpose.py       |   2 +-
 hysop/backend/device/codegen/symbolic/expr.py |   6 +-
 .../kernels/custom_symbolic_time_integrate.py |   4 +-
 hysop/backend/device/kernel_autotuner.py      |   4 +-
 .../opencl/autotunable_kernels/transpose.py   |   2 +-
 hysop/backend/device/opencl/clpeak.py         |   2 +-
 .../opencl/opencl_autotunable_kernel.py       |   2 +-
 .../backend/device/opencl/opencl_operator.py  |   2 +-
 .../backend/device/opencl/opencl_symbolic.py  |   3 +
 hysop/backend/device/opencl/opencl_types.py   |   6 +-
 .../device/opencl/operator/diffusion.py       |  12 ++-
 .../operator/directional/advection_dir.py     |  14 +--
 .../device/opencl/operator/external_force.py  |   2 +-
 hysop/backend/hardware/pci.py                 |   4 +-
 .../host/fortran/operator/diffusion.py        |   9 +-
 .../host/fortran/operator/fortran_fftw.py     |  15 ++-
 hysop/backend/host/host_array_backend.py      |   4 +-
 hysop/backend/host/host_operator.py           |  39 +++----
 hysop/backend/host/python/operator/curl.py    |   2 +-
 .../backend/host/python/operator/diffusion.py |  12 +--
 .../operator/directional/stretching_dir.py    |   5 +
 .../host/python/operator/poisson_curl.py      |  26 ++---
 hysop/core/checkpoints.py                     |   2 +-
 hysop/core/graph/computational_graph.py       |  26 ++---
 hysop/core/graph/computational_node.py        |   4 +-
 hysop/core/graph/graph_builder.py             |   8 +-
 hysop/core/graph/node_requirements.py         |  37 ++++---
 hysop/core/memory/memory_request.py           |   2 +-
 hysop/core/mpi/topo_tools.py                  |   1 -
 hysop/fields/continuous_field.py              |   2 +-
 hysop/fields/discrete_field.py                |   2 +-
 hysop/fields/ghost_exchangers.py              |   2 +-
 hysop/fields/tests/test_cartesian.py          |  27 +++--
 hysop/mesh/cartesian_mesh.py                  |   2 +-
 hysop/numerics/fft/gpyfft_fft.py              |   4 +-
 hysop/numerics/fft/host_fft.py                |  34 +++---
 hysop/numerics/interpolation/polynomial.py    |  14 +--
 hysop/numerics/odesolvers/runge_kutta.py      |   8 +-
 hysop/numerics/remesh/kernel_generator.py     |  14 +--
 hysop/numerics/stencil/stencil.py             |   2 +-
 hysop/numerics/stencil/stencil_generator.py   |   2 +-
 hysop/numerics/tests/bench_fft.py             |   2 +-
 hysop/numerics/tests/test_fft.py              |   2 +-
 hysop/operator/analytic.py                    |  24 ++---
 hysop/operator/base/convergence.py            |   2 +-
 hysop/operator/base/curl.py                   |  40 ++++---
 .../operator/base/custom_symbolic_operator.py |   8 +-
 hysop/operator/base/diffusion.py              |  14 ++-
 hysop/operator/base/memory_reordering.py      |   2 +-
 hysop/operator/base/poisson.py                |  10 ++
 hysop/operator/base/poisson_curl.py           |  25 ++++-
 hysop/operator/base/solenoidal_projection.py  |  58 +++++-----
 hysop/operator/base/spectral_operator.py      |  66 ++++++++----
 hysop/operator/base/stretching_dir.py         |  13 ++-
 hysop/operator/directional/diffusion_dir.py   |  13 +++
 hysop/operator/directional/stretching_dir.py  |  26 +++++
 hysop/operator/hdf_io.py                      |  22 +++-
 hysop/operator/memory_reordering.py           |  70 ++++++------
 .../operator/tests/test_bilevel_advection.py  |   5 +-
 hysop/operator/tests/test_custom_symbolic.py  |  18 ++--
 hysop/operator/tests/test_diffusion.py        |   4 +-
 .../tests/test_directional_diffusion.py       |   6 +-
 .../tests/test_directional_stretching.py      |   8 +-
 hysop/operator/tests/test_poisson.py          |   4 +-
 .../tests/test_solenoidal_projection.py       |   6 +-
 hysop/operator/transpose.py                   |   4 +-
 hysop/problem.py                              |   2 +-
 hysop/symbolic/spectral.py                    |   6 ++
 hysop/tools/callback.py                       |   8 +-
 hysop/tools/debug_dumper.py                   |   4 +-
 hysop/tools/decorators.py                     |   2 +-
 hysop/tools/enum.py                           |   2 +-
 hysop/tools/field_utils.py                    |   2 +-
 hysop/tools/hptt_utils.py                     |   6 +-
 hysop/tools/misc.py                           |   2 +-
 hysop/tools/mpi_utils.py                      |   8 +-
 hysop/tools/numba_utils.py                    |   4 +-
 hysop/tools/profiler.py                       |   4 +-
 hysop/tools/spectral_utils.py                 |   2 +-
 hysop/tools/transposition_states.py           |   2 +-
 hysop/topology/cartesian_descriptor.py        |   3 +-
 hysop_examples/example_utils.py               |   3 +-
 hysop_examples/examples/analytic/analytic.py  |   2 -
 .../particles_above_salt_bc.py                | 102 +++++++++---------
 .../particles_above_salt_periodic.py          |   2 +-
 .../particles_above_salt_symmetrized.py       |   2 +-
 setup.py.in                                   |   2 +-
 101 files changed, 631 insertions(+), 464 deletions(-)

diff --git a/README.md b/README.md
index 6a37b8f96..14ed56a2c 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,6 @@
 [![Platform](https://img.shields.io/badge/platform-linux--64%20%7C%C2%A0%20osx--64-lightgrey.svg)]()
-[![Python 2.7](https://img.shields.io/badge/python-2.7-blue.svg)](https://www.python.org/downloads/release/python-270/)
+[![Python 3.8](https://img.shields.io/badge/python-3.8-blue.svg)](https://www.python.org/downloads/release/python-380/)
+[![Python 3.9](https://img.shields.io/badge/python-3.9-blue.svg)](https://www.python.org/downloads/release/python-390/)
 [![Licence](https://img.shields.io/badge/licence-APLv2-blue.svg)](https://www.apache.org/licenses/LICENSE-2.0)
 [![Pipeline Status](https://gricad-gitlab.univ-grenoble-alpes.fr/particle_methods/hysop/badges/master/pipeline.svg)](https://gricad-gitlab.univ-grenoble-alpes.fr/particle_methods/hysop/commits/master)
 [![Docker Pulls](https://img.shields.io/docker/pulls/keckj/hysop.svg)](https://hub.docker.com/r/keckj/hysop/tags)
diff --git a/ci/scripts/test.sh b/ci/scripts/test.sh
index 1d67f1da2..df3cb5b42 100755
--- a/ci/scripts/test.sh
+++ b/ci/scripts/test.sh
@@ -68,7 +68,7 @@ fi
 
 export PYTHONPATH="${INSTALL_DIR}/lib/python3.8/site-packages:${INSTALL_DIR}"
 export MPLBACKEND='cairo'
-export HYSOP_VERBOSE=0
+export HYSOP_VERBOSE=1
 export HYSOP_DEBUG=0
 export HYSOP_PROFILE=0
 export HYSOP_KERNEL_DEBUG=0
@@ -126,16 +126,16 @@ if [ "$RUN_TESTS" = true ]; then
     #hysop_test "operator/tests/test_penalization.py"
     #hysop_test "operator/tests/test_velocity_correction.py"
     #hysop_test "operator/tests/test_restriction_filter.py"
-    hysop_test "operator/tests/test_directional_advection.py"
-    hysop_test "operator/tests/test_directional_diffusion.py"
-    hysop_test "operator/tests/test_directional_stretching.py"
-    hysop_test "operator/tests/test_custom_symbolic.py"
-    hysop_test "operator/tests/test_spectral_derivative.py"
-    hysop_test "operator/tests/test_spectral_curl.py"
-    hysop_test "operator/tests/test_diffusion.py"
-    hysop_test "operator/tests/test_poisson.py"
-    hysop_test "operator/tests/test_solenoidal_projection.py"
-    hysop_test "operator/tests/test_poisson_curl.py"
+    #hysop_test "operator/tests/test_directional_advection.py"
+    #hysop_test "operator/tests/test_directional_diffusion.py"
+    #hysop_test "operator/tests/test_directional_stretching.py"
+    #hysop_test "operator/tests/test_custom_symbolic.py"
+    #hysop_test "operator/tests/test_spectral_derivative.py"
+    #hysop_test "operator/tests/test_spectral_curl.py"
+    #hysop_test "operator/tests/test_diffusion.py"
+    #hysop_test "operator/tests/test_poisson.py"
+    #hysop_test "operator/tests/test_solenoidal_projection.py"
+    #hysop_test "operator/tests/test_poisson_curl.py"
     ${HYSOP_DIR}/fields/tests/test_cartesian.sh
     ${HYSOP_DIR}/core/tests/test_checkpoint.sh
 fi
@@ -163,7 +163,7 @@ if [ "${RUN_EXAMPLES}" = true ]; then
 fi
 
 if [ "${HAS_CACHE_DIR}" = true ]; then
-    rsync -rtvu "${HYSOP_CACHE_DIR}/" "${CACHE_DIR}/"
+    rsync -rtu "${HYSOP_CACHE_DIR}/" "${CACHE_DIR}/"
     find "${CACHE_DIR}" -name '*.lock' -delete
 fi
 
diff --git a/hysop/__init__.py.in b/hysop/__init__.py.in
index 0da4c6dfb..997c1b758 100644
--- a/hysop/__init__.py.in
+++ b/hysop/__init__.py.in
@@ -133,14 +133,17 @@ if '__formatwarning' not in locals():
 def __new_formatwarning(warning, category, filename, lineno, line=None):
     # filter out useless intel opencl build "warnings" stacktraces
     message = str(warning)
-    if ('Compilation started' in message) and ('Done.' in message):
+    warning_msg = __formatwarning(warning, category, filename, lineno, line=line)
+    if not warning_msg:
+        pass
+    elif ('Compilation started' in message):
         pass
     elif ('Non-empty compiler output encountered.' in message):
         pass
     elif __TRACE_WARNINGS__:
         print('::WARNING STACK TRACE::')
         traceback.print_stack()
-    return __formatwarning(warning, category, filename, lineno, line=line)
+    return warning_msg
 warnings.formatwarning = __new_formatwarning
 
 if __TRACE_CALLS__:
diff --git a/hysop/backend/device/autotunable_kernel.py b/hysop/backend/device/autotunable_kernel.py
index 9fdf35af4..155126ca5 100644
--- a/hysop/backend/device/autotunable_kernel.py
+++ b/hysop/backend/device/autotunable_kernel.py
@@ -145,7 +145,7 @@ class AutotunableKernel(object, metaclass=ABCMeta):
                         s += ss
                         s += '\nHASHED ARGUMENT {}: {}'.format(i, h)
             if kwds:
-                items = sorted(kwds.items(), key=lambda x: x[0])
+                items = tuple(sorted(kwds.items(), key=lambda x: x[0]))
                 if (h is None):
                     h, s  = _hash_karg(*items[0])
                 else:
@@ -336,7 +336,7 @@ class AutotunableKernel(object, metaclass=ABCMeta):
                 msg+='and not its instance id.'
                 msg=msg.format(t, type(v), str(type(v)))
                 raise RuntimeError(msg)
-        items = sorted(extra_parameters.items(), key=lambda x: x[0])
+        items = tuple(sorted(extra_parameters.items(), key=lambda x: x[0]))
         return hash(frozenset(items))
 
     @abstractmethod
@@ -707,7 +707,7 @@ class AutotunerWorkConfiguration(object):
 
         for fname in self.filter_names:
             fn = self._filters[fname]
-            candidates = filter(fn, candidates)
+            candidates = tuple(filter(fn, candidates))
             if self.__debug_filters:
                 candidates, _ = it.tee(candidates)
                 msg=' *Filter {}:\n {}\n'.format(fname,
diff --git a/hysop/backend/device/codegen/base/codegen.py b/hysop/backend/device/codegen/base/codegen.py
index e2c844801..cc92c7d72 100644
--- a/hysop/backend/device/codegen/base/codegen.py
+++ b/hysop/backend/device/codegen/base/codegen.py
@@ -667,7 +667,7 @@ class CodeGenerator(object):
         self.generate(blocks,genset)
         if blocks:
             entire_code = ''
-            blocks = sorted(blocks.items(), key=operator.itemgetter(1))
+            blocks = tuple(sorted(blocks.items(), key=operator.itemgetter(1)))
             for bn,(priority,code,comment) in blocks:
                 if not code:
                     continue
diff --git a/hysop/backend/device/codegen/base/struct_codegen.py b/hysop/backend/device/codegen/base/struct_codegen.py
index e6b0efee2..3a218a9de 100644
--- a/hysop/backend/device/codegen/base/struct_codegen.py
+++ b/hysop/backend/device/codegen/base/struct_codegen.py
@@ -35,8 +35,8 @@ class StructCodeGenerator(OpenClCodeGenerator):
 
     def c_decl(self):
         assert (self.context is not None)
-        dtype,cdecl = cl.tools.match_dtype_to_c_struct( \
-                self.device,self.ctype.replace('struct',''),self.dtype,context=self.context)
+        (dtype, cdecl) = cl.tools.match_dtype_to_c_struct( \
+                self.device, self.ctype.replace('struct',''), self.dtype, context=self.context)
         return cdecl
 
     def gencode(self, comments, ctype_overrides):
diff --git a/hysop/backend/device/codegen/base/variables.py b/hysop/backend/device/codegen/base/variables.py
index 910fc0cd3..a46430f0d 100644
--- a/hysop/backend/device/codegen/base/variables.py
+++ b/hysop/backend/device/codegen/base/variables.py
@@ -430,9 +430,9 @@ class CodegenVariable(object):
 
         if (init is not None):
             if compact:
-                code = '{}={}'.format(code,init)
+                code = '{}={}'.format(code, init)
             else:
-                code = '{} $= {}'.format(code,init)
+                code = '{} $= {}'.format(code, init)
         elif self.known():
             self.force_symbolic(False)
             sval = self.sval()
@@ -769,13 +769,13 @@ class CodegenVectorClBuiltin(CodegenVector):
         return super(CodegenVectorClBuiltin, self).newvar(name=name, btype=btype, dim=dim,
                 **kwds)
 
-    def view(self,name,components,const=False):
-        if isinstance(components,slice):
+    def view(self, name, components, const=False):
+        if isinstance(components, slice):
             start, stop, step = components.indices(self.dim)
-            it = range(start,stop,step)
+            it = tuple(range(start, stop, step))
             dim = len(it)
         elif components is Ellipsis:
-            it = range(self.dim)
+            it = tuple(range(self.dim))
             dim = self.dim
         else:
             raise ValueError('Unknown components type {}.'.format(components.__class__.__name__))
@@ -831,7 +831,7 @@ class CodegenVectorClBuiltin(CodegenVector):
         if isinstance(key, range):
             key = tuple(key)
         if isinstance(key,slice) :
-            ids = range(*key.indices(dim))
+            ids = tuple(range(*key.indices(dim)))
             if self.declared and key.indices(dim)==(0,dim,1):
                 return self.name
             else:
diff --git a/hysop/backend/device/codegen/functions/apply_stencil.py b/hysop/backend/device/codegen/functions/apply_stencil.py
index d9d83534e..54e8f1e5f 100644
--- a/hysop/backend/device/codegen/functions/apply_stencil.py
+++ b/hysop/backend/device/codegen/functions/apply_stencil.py
@@ -31,7 +31,7 @@ class ApplyStencilFunction(OpenClFunctionCodeGenerator):
         symbol2vars = first_not_None(symbol2vars, {})
         assert set(symbol2vars.keys())==stencil.variables()
 
-        extra_inputs  = set(extra_inputs + symbol2vars.values())
+        extra_inputs  = set(extra_inputs).union(symbol2vars.values())
         scalar_inputs = set(scalar_inputs)
         vector_inputs = set(vector_inputs)
 
@@ -40,7 +40,7 @@ class ApplyStencilFunction(OpenClFunctionCodeGenerator):
         vtype = typegen.vtype(ftype,components)
 
         if (vector_suffixes is None):
-            vector_suffixes = range(components)
+            vector_suffixes = tuple(range(components))
 
         has_custom_id = (custom_id is not None)
 
diff --git a/hysop/backend/device/codegen/functions/cache_load.py b/hysop/backend/device/codegen/functions/cache_load.py
index 2cdd8c664..f54332125 100644
--- a/hysop/backend/device/codegen/functions/cache_load.py
+++ b/hysop/backend/device/codegen/functions/cache_load.py
@@ -179,7 +179,7 @@ class CacheLoadFunction(OpenClFunctionCodeGenerator):
 
             s.barrier(_local=True)
             with s._block_():
-                nested_loops = [_cache_iterate_(i) for i in range(work_dim-1,-1,-1)]
+                nested_loops = [_cache_iterate_(i) for i in range(work_dim-1, -1, -1)]
                 with nested(*nested_loops):
                     with s._block_():
                         s.append(LID.declare(const=True,init=compute_index(idx=local_pos,
diff --git a/hysop/backend/device/codegen/functions/compute_index.py b/hysop/backend/device/codegen/functions/compute_index.py
index a6188ae6b..f349df121 100644
--- a/hysop/backend/device/codegen/functions/compute_index.py
+++ b/hysop/backend/device/codegen/functions/compute_index.py
@@ -36,7 +36,7 @@ class ComputeIndexFunction(OpenClFunctionCodeGenerator):
             with s._if_(wrap()):
                 self.append('{idx} = ({idx}+{size}) % {size};'.format(idx=idx(),size=size))
             ss = '{}'.format(idx[dim-1])
-            for i in range(dim-2,-1,-1):
+            for i in range(dim-2, -1, -1):
                 ss = '({}*{}+{})'.format(ss,size[i],idx[i])
             ret = 'return {};'.format(ss)
             self.append(ret)
diff --git a/hysop/backend/device/codegen/kernels/custom_symbolic.py b/hysop/backend/device/codegen/kernels/custom_symbolic.py
index 30af9a64f..47eb23757 100644
--- a/hysop/backend/device/codegen/kernels/custom_symbolic.py
+++ b/hysop/backend/device/codegen/kernels/custom_symbolic.py
@@ -300,8 +300,7 @@ class SymbolicCodegenContext(object):
         array_contiguous_ghosts = SortedDict()
         array_sizes = SortedDict()
 
-        arrays = set(f for f in (expr_info.input_arrays.values() +
-                                            expr_info.output_arrays.values()))
+        arrays = set(expr_info.input_arrays.values()).union(expr_info.output_arrays.values())
 
         for a in arrays:
             ctype = a.ctype
@@ -710,7 +709,7 @@ class CustomSymbolicKernelGenerator(KernelCodeGenerator, metaclass=ABCMeta):
             xmin, dx, inv_dx = None, None, None
         else:
             declare_mesh_properties = True
-            mesh_info_0 = field_mesh_infos.values()[0]
+            mesh_info_0 = next(iter(field_mesh_infos.values()))
             dx     = mesh_info_0['dx'].alias('dx', const=True)
             inv_dx = mesh_info_0['inv_dx'].alias('inv_dx', const=True)
             xmin   = mesh_info_0['local_mesh']['xmin'].alias('xmin', const=True)
@@ -1048,7 +1047,7 @@ class CustomSymbolicKernelGenerator(KernelCodeGenerator, metaclass=ABCMeta):
 
         if (granularity>0):
             self.jumpline()
-        self.decl_vars(*tuple([loop_id]+array_vids.values()))
+        self.decl_vars(*tuple([loop_id]+list(array_vids.values())))
 
         if self.field_mesh_infos:
             x0 = csc.space_symbols[symbolic_space_symbols[0]]
@@ -1079,11 +1078,11 @@ class CustomSymbolicKernelGenerator(KernelCodeGenerator, metaclass=ABCMeta):
                 code = idx_i.affect(self, init=loop_id[i])
 
             if self.field_mesh_infos:
-                fmi = self.field_mesh_infos[array_vids.keys()[0]]
+                fmi = self.field_mesh_infos[next(iter(array_vids))]
                 for i in range(kdim, kdim+granularity):
                     xi = csc.space_symbols[symbolic_space_symbols[i]]
                     code = '{xi} = {x0} + {vid}*{dx};'.format(
-                            xi=xi, vid=array_vids.values()[0][i],
+                            xi=xi, vid=next(iter(array_vids.values()))[i],
                             voffset=self.voffset,
                             x0=fmi['local_mesh']['xmin'][i], dx=fmi['dx'][i])
                     self.append(code)
@@ -1123,10 +1122,10 @@ class CustomSymbolicKernelGenerator(KernelCodeGenerator, metaclass=ABCMeta):
                         code = idx_i.affect(self, init=loop_id[i])
 
                         if self.field_mesh_infos:
-                            fmi = self.field_mesh_infos[array_vids.keys()[0]]
+                            fmi = self.field_mesh_infos[next(iter(array_vids))]
                             xi = csc.space_symbols[symbolic_space_symbols[i]]
                             code = '{xi} = {x0} + {vid}*{dx};'.format(
-                                    xi=xi, vid=array_vids.values()[0][i],
+                                    xi=xi, vid=next(iter(array_vids.values()))[i],
                                     voffset=self.voffset,
                                     x0=fmi['local_mesh']['xmin'][i], dx=fmi['dx'][i])
                             self.append(code)
@@ -1156,11 +1155,11 @@ class CustomSymbolicKernelGenerator(KernelCodeGenerator, metaclass=ABCMeta):
                                             align=True)
 
                         if self.field_mesh_infos:
-                            fmi = self.field_mesh_infos[array_vids.keys()[0]]
+                            fmi = self.field_mesh_infos[next(iter(array_vids))]
                             xi = csc.space_symbols[symbolic_space_symbols[i]]
                             code = '{xi} = {x0} + convert_{vftype}({vid}+{voffset}+{lo})*{dx};'
                             code=code.format(
-                                    xi=xi, vid=array_vids.values()[0][i], lo=local_offset,
+                                    xi=xi, vid=next(iter(array_vids.values()))[i], lo=local_offset,
                                     voffset=self.voffset, vftype=csc.vftype,
                                     x0=fmi['local_mesh']['xmin'][i], dx=fmi['dx'][i])
                             self.append(code)
@@ -1209,7 +1208,7 @@ class CustomSymbolicKernelGenerator(KernelCodeGenerator, metaclass=ABCMeta):
             except:
                 raise
 
-        nested_loops = [work_iterate(i) for i in range(kdim-1,-1,-1)]
+        nested_loops = [work_iterate(i) for i in range(kdim-1, -1, -1)]
         return nested_loops
 
     def gencode(self):
@@ -1256,8 +1255,10 @@ class CustomSymbolicKernelGenerator(KernelCodeGenerator, metaclass=ABCMeta):
             s.decl_aligned_vars(*(aij for ai in s.array_args.values() for aij in ai.values()))
             s.decl_aligned_vars(*csc.buffer_args.values())
             s.comment('Local and private memory arrays')
-            s.decl_aligned_vars(*(aij for ai in array_local_data.values() +
-                                    array_local_rdata.values()+array_private_data.values()
+            s.decl_aligned_vars(*(aij for ai in
+                                    tuple(array_local_data.values())  +
+                                    tuple(array_local_rdata.values()) +
+                                    tuple(array_private_data.values())
                                     for aij in filter(lambda x:x, ai.values())))
             s.comment('Iterating over array lines')
             nested_loops = self._generate_loop_context()
diff --git a/hysop/backend/device/codegen/kernels/directional_advection.py b/hysop/backend/device/codegen/kernels/directional_advection.py
index f59792772..85a941c1f 100644
--- a/hysop/backend/device/codegen/kernels/directional_advection.py
+++ b/hysop/backend/device/codegen/kernels/directional_advection.py
@@ -545,7 +545,7 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
                     yield ctx
             except:
                 raise
-        nested_loops = [_work_iterate_(i) for i in range(work_dim-1,-1,-1)]
+        nested_loops = [_work_iterate_(i) for i in range(work_dim-1, -1, -1)]
 
         with s._kernel_():
             field_infos.declare(s)
@@ -575,7 +575,7 @@ class DirectionalAdvectionKernelGenerator(KernelCodeGenerator):
 
             with s._align_() as al:
                 npart.declare(al, const=True, align=True)
-                poffset.declare(al, init=range(nparticles), const=True, align=True)
+                poffset.declare(al, init=tuple(range(nparticles)), const=True, align=True)
             s.jumpline()
 
             if is_cached:
diff --git a/hysop/backend/device/codegen/kernels/tests/test_codegen_directional_stretching.py b/hysop/backend/device/codegen/kernels/tests/test_codegen_directional_stretching.py
index 1780cb3de..adbb6a169 100644
--- a/hysop/backend/device/codegen/kernels/tests/test_codegen_directional_stretching.py
+++ b/hysop/backend/device/codegen/kernels/tests/test_codegen_directional_stretching.py
@@ -172,7 +172,7 @@ class TestDirectionalStretching(object):
 
         def deriv(field):
             res = np.zeros_like(field)
-            for i in range(-order//2,order//2+1,1):
+            for i in range(-order//2, order//2+1, 1):
                 res += stencil[i+order//2]*np.roll(field,-i,axis=2)
             return res*self.inv_dx[0]
 
diff --git a/hysop/backend/device/codegen/kernels/tests/test_codegen_transpose.py b/hysop/backend/device/codegen/kernels/tests/test_codegen_transpose.py
index 5d6444b4b..f9de9412f 100644
--- a/hysop/backend/device/codegen/kernels/tests/test_codegen_transpose.py
+++ b/hysop/backend/device/codegen/kernels/tests/test_codegen_transpose.py
@@ -249,7 +249,7 @@ class TestTranspose(object):
                 local_work_size[i] = 1
             max_work = device.max_work_group_size
             while (prod(local_work_size) > max_work):
-                for i in range(work_dim-1,-1,-1):
+                for i in range(work_dim-1, -1, -1):
                     if local_work_size[i] > 1:
                         break
                 local_work_size[i] //= 2
diff --git a/hysop/backend/device/codegen/kernels/transpose.py b/hysop/backend/device/codegen/kernels/transpose.py
index cdc6d0abc..7f399b71b 100644
--- a/hysop/backend/device/codegen/kernels/transpose.py
+++ b/hysop/backend/device/codegen/kernels/transpose.py
@@ -177,7 +177,7 @@ class TransposeKernelGenerator(KernelCodeGenerator):
         if len(_axes) != pdim:
             raise ValueError(msg)
 
-        _permutation = (axes != range(pdim))
+        _permutation = (axes != set(range(pdim)))
         _naxes = sum(_permutation)
         if (_naxes == 0):
             msg='There is nothing to transpose with given axes {}.'
diff --git a/hysop/backend/device/codegen/symbolic/expr.py b/hysop/backend/device/codegen/symbolic/expr.py
index 507ac6a40..1ee0e5a7e 100644
--- a/hysop/backend/device/codegen/symbolic/expr.py
+++ b/hysop/backend/device/codegen/symbolic/expr.py
@@ -368,9 +368,9 @@ class UpdateVars(Expr):
                 for (src, dst) in self.private_stores:
                     dst.affect(al, init=src, align=True)
         if self.local_stores:
-            srcs    = map(lambda x: x[0], self.local_stores)
-            ptrs    = map(lambda x: x[1], self.local_stores)
-            offsets = map(lambda x: x[2], self.local_stores)
+            srcs    = tuple(map(lambda x: x[0], self.local_stores))
+            ptrs    = tuple(map(lambda x: x[1], self.local_stores))
+            offsets = tuple(map(lambda x: x[2], self.local_stores))
             codegen.multi_vstore_if(csc.is_last_active,
                     lambda i: '{}+{} < {}'.format(csc.full_offset, i, csc.compute_grid_size[0]),
                     csc.vectorization, csc.local_offset,
diff --git a/hysop/backend/device/codegen/symbolic/kernels/custom_symbolic_time_integrate.py b/hysop/backend/device/codegen/symbolic/kernels/custom_symbolic_time_integrate.py
index f07fd2571..7f17339ac 100644
--- a/hysop/backend/device/codegen/symbolic/kernels/custom_symbolic_time_integrate.py
+++ b/hysop/backend/device/codegen/symbolic/kernels/custom_symbolic_time_integrate.py
@@ -65,8 +65,8 @@ class CustomSymbolicTimeIntegrateKernelGenerator(CustomSymbolicKernelGenerator):
 
                     dval = CustomSymbolicFunction.default_out_of_bounds_value(ctype_to_dtype(ctype))
                     if is_local:
-                        size = csc.array_size(lhs.field, lhs.index)
-                        assert (size > 2*ghosts)
+                        size = int(csc.array_size(lhs.field, lhs.index))
+                        assert (size > 2*ghosts), (size, 2*ghosts)
                         vload = Vload(tg, ctype, vectorization, default_val=dval,
                                 itype=itype, restrict=True, storage=var.storage,
                                 known_args=dict(size=size))
diff --git a/hysop/backend/device/kernel_autotuner.py b/hysop/backend/device/kernel_autotuner.py
index 41a1c12ea..12dd091a1 100644
--- a/hysop/backend/device/kernel_autotuner.py
+++ b/hysop/backend/device/kernel_autotuner.py
@@ -462,7 +462,7 @@ class KernelAutotuner(object, metaclass=ABCMeta):
             keep_only = max(previous_pow2(kept_count),1)
             self._print_first_step_results(total_count, kept_count, pruned_count,
                     failed_count, keep_only)
-            candidates = sorted(bench_results.items(), key=lambda x: x[1][self.kernel_statistics_idx])
+            candidates = tuple(sorted(bench_results.items(), key=lambda x: x[1][self.kernel_statistics_idx]))
             candidates = candidates[:keep_only]
             while(len(candidates)>1):
                 step_count += 1
@@ -478,7 +478,7 @@ class KernelAutotuner(object, metaclass=ABCMeta):
                                              best_stats=best_stats,
                                              global_work_size=global_work_size,
                                              local_work_size=local_work_size)
-                candidates = sorted(candidates, key=lambda x: x[1][self.kernel_statistics_idx])
+                candidates = tuple(sorted(candidates, key=lambda x: x[1][self.kernel_statistics_idx]))
                 self._print_step_results(candidates, self.kernel_statistics_idx)
                 candidates = candidates[:max(previous_pow2(len(candidates)),1)]
                 ks.push_step(step_count, candidates)
diff --git a/hysop/backend/device/opencl/autotunable_kernels/transpose.py b/hysop/backend/device/opencl/autotunable_kernels/transpose.py
index b39d0742f..c64306d84 100644
--- a/hysop/backend/device/opencl/autotunable_kernels/transpose.py
+++ b/hysop/backend/device/opencl/autotunable_kernels/transpose.py
@@ -154,7 +154,7 @@ class OpenClAutotunableTransposeKernel(OpenClAutotunableKernel):
         imax = int(math.log(max_tile_size, 2))
         jmax = int(math.log(max_tile_size, 3)) if flag in (AutotunerFlags.EXHAUSTIVE,) else 0
         tile_sizes = tuple( int((2**i)*(3**j))
-                for (i,j) in it.product(range(0,imax+1), range(0,jmax+1)))
+                for (i,j) in it.product(range(0, imax+1), range(0, jmax+1)))
         tile_sizes = (max_tile_size,) + tuple(sorted(tile_sizes, reverse=True))
         tile_sizes = tuple(filter(lambda x: (x>=max_tile_size//8) and (x<=max_tile_size), tile_sizes))
 
diff --git a/hysop/backend/device/opencl/clpeak.py b/hysop/backend/device/opencl/clpeak.py
index 46d322f6a..c1504b8c2 100644
--- a/hysop/backend/device/opencl/clpeak.py
+++ b/hysop/backend/device/opencl/clpeak.py
@@ -251,7 +251,7 @@ class ClPeakInfo(object):
             def criteria(key):
                 values = transfer_bandwidth[key]
                 return -(values['enqueuewritebuffer'] + values['enqueuereadbuffer'])
-            params = sorted(transfer_bandwidth.keys(), key=criteria)
+            params = tuple(sorted(transfer_bandwidth.keys(), key=criteria))
             self._set_exec_params(*params[0])
             self._preferred_exec_params = params[0]
             self._transfer_bandwidth_values = transfer_bandwidth
diff --git a/hysop/backend/device/opencl/opencl_autotunable_kernel.py b/hysop/backend/device/opencl/opencl_autotunable_kernel.py
index d25c49df9..6f5e2e6f6 100644
--- a/hysop/backend/device/opencl/opencl_autotunable_kernel.py
+++ b/hysop/backend/device/opencl/opencl_autotunable_kernel.py
@@ -246,7 +246,7 @@ class OpenClAutotunableKernel(AutotunableKernel):
             global_work_size = tuple(x for x in global_work_size) + (1,)*(3-dim)
             local_work_size  = tuple(x for x in local_work_size)  + (1,)*(3-dim)
 
-        sorted_args = sorted(args_mapping.items(), key=lambda x:x[1][0])
+        sorted_args = tuple(sorted(args_mapping.items(), key=lambda x:x[1][0]))
         assert len(sorted_args) == len(args_list)
 
         dump_folder = self.autotuner_config.dump_folder
diff --git a/hysop/backend/device/opencl/opencl_operator.py b/hysop/backend/device/opencl/opencl_operator.py
index 4c942f37f..3dcff415c 100644
--- a/hysop/backend/device/opencl/opencl_operator.py
+++ b/hysop/backend/device/opencl/opencl_operator.py
@@ -90,7 +90,7 @@ class OpenClOperator(ComputationalGraphOperator, metaclass=ABCMeta):
                 if not _vars:
                     msg='Cannot deduce domain without input or output fields.'
                     raise RuntimeError(msg)
-                domain = _vars.keys()[0].domain
+                domain = next(iter(_vars)).domain
                 mpi_params = MPIParams(comm=domain.task_comm,
                                        task_id=domain.current_task())
                 cl_env = get_or_create_opencl_env(mpi_params)
diff --git a/hysop/backend/device/opencl/opencl_symbolic.py b/hysop/backend/device/opencl/opencl_symbolic.py
index a57ec6a0c..84cda9eeb 100644
--- a/hysop/backend/device/opencl/opencl_symbolic.py
+++ b/hysop/backend/device/opencl/opencl_symbolic.py
@@ -50,6 +50,9 @@ class OpenClSymbolic(OpenClOperator):
         am.update(cls.__available_methods)
         return am
 
+    @debug
+    def __new__(cls, **kwds):
+        return super(OpenClSymbolic, cls).__new__(cls, **kwds)
 
     @debug
     def __init__(self, **kwds):
diff --git a/hysop/backend/device/opencl/opencl_types.py b/hysop/backend/device/opencl/opencl_types.py
index 000dd1515..f62930e1c 100644
--- a/hysop/backend/device/opencl/opencl_types.py
+++ b/hysop/backend/device/opencl/opencl_types.py
@@ -447,9 +447,9 @@ class OpenClTypeGen(TypeGen):
             msg += 'and regular expression \'{}\'.'
             msg=msg.format(sversion,_regexp)
             raise RuntimeError(msg)
-        major = match.group(1)
-        minor = match.group(2)
-        return (major,minor)
+        major = int(match.group(1))
+        minor = int(match.group(2))
+        return (major, minor)
 
     def dtype_from_str(self,stype):
         stype = stype.replace('ftype', self.fbtype).replace('fbtype',self.fbtype)
diff --git a/hysop/backend/device/opencl/operator/diffusion.py b/hysop/backend/device/opencl/operator/diffusion.py
index 9d7bbdc2c..48f2b697e 100644
--- a/hysop/backend/device/opencl/operator/diffusion.py
+++ b/hysop/backend/device/opencl/operator/diffusion.py
@@ -17,7 +17,11 @@ class OpenClDiffusion(DiffusionOperatorBase, OpenClSymbolic):
     """
     Solves the diffusion equation using an OpenCL FFT backend.
     """
-    
+
+    @debug
+    def __new__(cls, **kwds):
+        return super(OpenClDiffusion, cls).__new__(cls, **kwds)
+
     @debug
     def __init__(self, **kwds):
         super(OpenClDiffusion, self).__init__(**kwds)
@@ -29,7 +33,7 @@ class OpenClDiffusion(DiffusionOperatorBase, OpenClSymbolic):
         for (i,(Ft,Wn)) in enumerate(zip(self.forward_transforms, self.wave_numbers)):
             Fhs = Ft.output_symbolic_array('F{}_hat'.format(i))
             indices = local_indices_symbols[:Fhs.dim]
-        
+
             kname = 'filter_diffusion_{}d_{}'.format(Fhs.dim, i)
             kernel_names += (kname,)
 
@@ -38,11 +42,11 @@ class OpenClDiffusion(DiffusionOperatorBase, OpenClSymbolic):
                 indexed_Wi = self.tg._indexed_wave_numbers[Wi]
                 F += indexed_Wi
             expr = Assignment(Fhs, Fhs / (1 - nu*dt*F))
-            
+
             self.require_symbolic_kernel(kname, expr)
 
         self._kernel_names = kernel_names
-    
+
     @debug
     def setup(self, work):
         super(OpenClDiffusion, self).setup(work)
diff --git a/hysop/backend/device/opencl/operator/directional/advection_dir.py b/hysop/backend/device/opencl/operator/directional/advection_dir.py
index 45ec39661..346e57752 100644
--- a/hysop/backend/device/opencl/operator/directional/advection_dir.py
+++ b/hysop/backend/device/opencl/operator/directional/advection_dir.py
@@ -175,26 +175,28 @@ class OpenClDirectionalAdvection(DirectionalAdvectionBase, OpenClDirectionalOper
         if self.DEBUG:
             queue.flush()
             queue.finish()
+            dfin  = tuple(self.dadvected_fields_in.values())
+            dfout = tuple(self.dadvected_fields_out.values())
             print('OPENCL_DT= {}'.format(dt))
             self.advection_kernel_launcher(queue=queue, dt=dt).wait()
             print('OPENCL_P')
             self.dposition.print_with_ghosts()
             print('OPENCL_Sin (before remesh)')
-            print(self.dadvected_fields_in.values()[0].data[0].get(queue=queue))
+            print(dfin[0].data[0].get(queue=queue))
             print('OPENCL_Sout (before remesh)')
-            print(self.dadvected_fields_out.values()[0].data[0].get(queue=queue))
+            print(dfout.data[0].get(queue=queue))
             print()
             self.remesh_kernel_launcher(queue=queue).wait()
             print('OPENCL_Sout (before accumulation)')
-            data = self.dadvected_fields_out.values()[0].data[0]
+            data = dfout[0].data[0]
             print(data.get(queue=queue))
             print('OPENCL_Sout (before accumulation, no ghosts)  ID={}'.format(id(data)))
-            self.dadvected_fields_out.values()[0].print_with_ghosts()
+            dfout[0].print_with_ghosts()
             self.accumulate_and_exchange(queue=queue).wait()
             print('OPENCL_Sout (after accumulation)')
-            print(self.dadvected_fields_out.values()[0].data[0])
+            print(dfout[0].data[0])
             print('OPENCL_Sout (after accumulation, no ghosts)')
-            self.dadvected_fields_out.values()[0].print_with_ghosts()
+            dfout[0].print_with_ghosts()
         else:
             self.all_kernels(queue=queue, dt=dt)
 
diff --git a/hysop/backend/device/opencl/operator/external_force.py b/hysop/backend/device/opencl/operator/external_force.py
index f5e66c6df..dcc2690d3 100644
--- a/hysop/backend/device/opencl/operator/external_force.py
+++ b/hysop/backend/device/opencl/operator/external_force.py
@@ -108,7 +108,7 @@ class SymbolicExternalForce(ExternalForce):
                     if field in self._diffusion:
                         assert field in backward_transforms, field.name
                 replace = {sf:forward_transforms[sf.field].s for sf in efields}
-                frame = replace.values()[0].frame
+                frame = next(iter(replace.values)).frame
                 e = e.xreplace(replace)
             fft_expressions += (e,)
 
diff --git a/hysop/backend/hardware/pci.py b/hysop/backend/hardware/pci.py
index 483c7cc87..fef7bd6a2 100644
--- a/hysop/backend/hardware/pci.py
+++ b/hysop/backend/hardware/pci.py
@@ -299,7 +299,7 @@ class PciBridge(TopologyObject):
         super(PciBridge, self).__init__(parent, bridge)
 
     def pci_devices(self, split=False):
-        devs = sorted(self._pci_devices, key=lambda x: x.os_index())
+        devs = tuple(sorted(self._pci_devices, key=lambda x: x.os_index()))
         if split:
             devices = [dev for dev in devs if isinstance(dev,PciBridge)]
             devices +=[dev for dev in devs if isinstance(dev,PciDevice)]
@@ -346,7 +346,7 @@ class PciBridge(TopologyObject):
             pci_device = pci_device.to_string(expand_pci_tree)
             pci_device = pci_device.split('\n')
             pci_device[0] = '|--'+pci_device[0]
-            for i in range(1,len(pci_device)):
+            for i in range(1, len(pci_device)):
                 pci_device[i] = '|  '+pci_device[i]
             if (dev_id==0) and expand_pci_tree:
                 pci_device = ['|']+pci_device
diff --git a/hysop/backend/host/fortran/operator/diffusion.py b/hysop/backend/host/fortran/operator/diffusion.py
index a01cd44c1..036a2b77d 100644
--- a/hysop/backend/host/fortran/operator/diffusion.py
+++ b/hysop/backend/host/fortran/operator/diffusion.py
@@ -10,8 +10,13 @@ from hysop.core.graph.graph import op_apply
 class DiffusionFFTW(FortranFFTWOperator):
 
     @debug
-    def __init__(self, Fin, Fout,
-                       nu, variables, dt, **kargs):
+    def __new__(cls, Fin, Fout, nu, variables, dt, **kargs):
+        return super(DiffusionFFTW, cls).__new__(cls,
+                input_fields=None, output_fields=None,
+                input_params=None, **kargs)
+
+    @debug
+    def __init__(self, Fin, Fout, nu, variables, dt, **kargs):
         """Diffusion operator base.
 
         Parameters
diff --git a/hysop/backend/host/fortran/operator/fortran_fftw.py b/hysop/backend/host/fortran/operator/fortran_fftw.py
index 643edfdaa..03067d585 100644
--- a/hysop/backend/host/fortran/operator/fortran_fftw.py
+++ b/hysop/backend/host/fortran/operator/fortran_fftw.py
@@ -25,6 +25,11 @@ class FortranFFTWOperator(FortranOperator):
     fftw2py interface.
     """
 
+    @debug
+    def __new__(cls, input_fields, output_fields, **kwds):
+        return super(FortranFFTWOperator, cls).__new__(cls,
+            input_fields=input_fields, output_fields=output_fields, **kwds)
+
     @debug
     def __init__(self, input_fields, output_fields, **kwds):
         """
@@ -44,8 +49,8 @@ class FortranFFTWOperator(FortranOperator):
         check_instance(output_fields, dict, keys=Field,
                        values=CartesianTopologyDescriptors)
 
-        nb_components =input_fields.keys()[0].nb_components
-        tensor_fields = set(input_fields.keys()+output_fields.keys())
+        nb_components = next(iter(input_fields)).nb_components
+        tensor_fields = set(input_fields.keys()).union(output_fields.keys())
         for tf in tensor_fields:
             for fi in tf.fields:
                 if (fi.dtype != HYSOP_REAL):
@@ -66,7 +71,7 @@ class FortranFFTWOperator(FortranOperator):
         # Special case: 3D diffusion of a scalar
         self._scalar_3d = (nb_components == 1) and all(tf.nb_components==1 for tf in tensor_fields)
 
-        domain = self.input_fields.keys()[0].domain
+        domain = next(iter(self.input_fields)).domain
         self.dim      = domain.dim
         self.domain   = domain
 
@@ -92,7 +97,7 @@ class FortranFFTWOperator(FortranOperator):
     def handle_topologies(self, input_topology_states, output_topology_states):
         super(FortranFFTWOperator,self).handle_topologies(input_topology_states, output_topology_states)
 
-        topology = self.input_fields.values()[0].topology
+        topology = next(iter(self.input_fields.values())).topology
         for (field,topoview) in self.input_fields.items():
             assert all(topoview.topology.cart_shape == topology.cart_shape), 'topology mismatch'
         for (field,topoview) in self.output_fields.items():
@@ -104,7 +109,7 @@ class FortranFFTWOperator(FortranOperator):
         if self.discretized:
             return
         super(FortranFFTWOperator,self).discretize()
-        topo_view = self.input_discrete_fields.values()[0].topology
+        topo_view = next(iter(self.input_discrete_fields.values())).topology
         self._fftw_discretize(topo_view)
 
     @debug
diff --git a/hysop/backend/host/host_array_backend.py b/hysop/backend/host/host_array_backend.py
index af003b8df..ae05dcf9b 100644
--- a/hysop/backend/host/host_array_backend.py
+++ b/hysop/backend/host/host_array_backend.py
@@ -264,9 +264,9 @@ class HostArrayBackend(ArrayBackend):
         self._unsupported_argument('empty_like', 'subok', subok, True)
         if (order is None) or (order == MemoryOrdering.SAME_ORDER):
             try:
-                if a.flags['C_CONTIGUOUS']:
+                if a.flags.c_contiguous:
                     order = MemoryOrdering.C_CONTIGUOUS
-                elif a.flags['F_CONTIGUOUS']:
+                elif a.flags.f_contiguous:
                     order = MemoryOrdering.F_CONTIGUOUS
                 else:
                     order = default_order
diff --git a/hysop/backend/host/host_operator.py b/hysop/backend/host/host_operator.py
index f52750eee..b2bbf16ae 100644
--- a/hysop/backend/host/host_operator.py
+++ b/hysop/backend/host/host_operator.py
@@ -41,26 +41,29 @@ class HostOperator(ComputationalGraphOperator, metaclass=ABCMeta):
         return set([Backend.HOST])
 
 
+class OpenClMappedMemoryObjectGetter(object):
+    def __new__(cls, obj, evt, **kwds):
+        return super(OpenClMappedMemoryObjectGetter, cls).__new__(cls, **kwds)
+
+    def __init__(self, obj, evt, **kwds):
+        super(OpenClMappedMemoryObjectGetter, self).__init__(**kwds)
+        check_instance(obj, OpenClMappable)
+        self.__obj = obj
+        self.__evt = evt
+
+    def __getitem__(self, key):
+        return obj.get_mapped_object(key=key)
+
+    @property
+    def evt(self):
+        return self.__evt
+
 
 class OpenClMappable(object):
     """
     Extend host operator capabilities to work on mapped opencl buffers
     """
 
-    class OpenClMappedMemoryObjectGetter(object):
-        def __new__(cls, obj, evt):
-            return super(OpenClMappedMemoryObjectGetter, cls).__new__(cls)
-        def __init__(self, obj, evt):
-            super(OpenClMappedMemoryObjectGetter, self).__init__()
-            check_instance(obj, OpenClMappable)
-            self.__obj = obj
-            self.__evt = evt
-        def __getitem__(self, key):
-            return obj.get_mapped_object(key=key)
-        @property
-        def evt(self):
-            return self.__evt
-
     @classmethod
     def supported_backends(cls):
         sb = super(OpenClMappable, cls).supported_backends()
@@ -119,7 +122,7 @@ class OpenClMappable(object):
         self.__mapped_objects     = {}
 
     def __del__(self):
-        self.unmap_objects()
+        self.unmap_objects(force=True)
 
     @property
     def cl_env(self):
@@ -207,9 +210,9 @@ class OpenClMappable(object):
         self.__mapped = True
         return evt
 
-    def unmap_objects(self):
+    def unmap_objects(self, force=False):
         msg='Device memory objects have already been unmapped from host.'
-        assert self.__mapped, msg
+        assert force or self.__mapped, msg
         self.__mapped_objects.clear()
         self.__mapped = False
 
@@ -234,7 +237,7 @@ class OpenClMappable(object):
             queue = first_not_None(queue, self.cl_env.default_queue)
             try:
                 evt = self.map_objects(queue, is_blocking)
-                yield self.OpenClMappedMemoryObjectGetter(self, evt)
+                yield OpenClMappedMemoryObjectGetter(self, evt)
             except:
                 raise
             finally:
diff --git a/hysop/backend/host/python/operator/curl.py b/hysop/backend/host/python/operator/curl.py
index 14666ad2e..305017c3d 100644
--- a/hysop/backend/host/python/operator/curl.py
+++ b/hysop/backend/host/python/operator/curl.py
@@ -68,12 +68,12 @@ def filter_curl_3d__1(Fin, K, Fout):
             for k in range(0, Fin.shape[2]):
                 Fout[i,j,k] = -K[k]*Fin[i,j,k]
 
+
 class PythonSpectralCurl(SpectralCurlOperatorBase, HostOperator):
     """
     Compute the curl by using an python FFT backend.
     """
 
-
     def setup(self, work):
         super(PythonSpectralCurl, self).setup(work=work)
 
diff --git a/hysop/backend/host/python/operator/diffusion.py b/hysop/backend/host/python/operator/diffusion.py
index ec36df387..9f8e5d0b9 100644
--- a/hysop/backend/host/python/operator/diffusion.py
+++ b/hysop/backend/host/python/operator/diffusion.py
@@ -20,7 +20,7 @@ class PythonDiffusion(DiffusionOperatorBase, OpenClMappable, HostOperator):
     """
     Solves the implicit diffusion equation using numpy fft.
     """
-        
+
     @classmethod
     def build_diffusion_filter(cls, dim, *args, **kwds):
         target = kwds.get('target', __DEFAULT_NUMBA_TARGET__)
@@ -71,18 +71,18 @@ class PythonDiffusion(DiffusionOperatorBase, OpenClMappable, HostOperator):
             raise NotImplementedError(msg)
         return functools.partial(F, *args)
 
-    
+
     def setup(self, work):
         super(PythonDiffusion, self).setup(work=work)
         diffusion_filters = ()
-        for (Fo,Ft,Kd) in zip(self.dFout.dfields, 
-                              self.forward_transforms, 
-                              self.all_dkds): 
+        for (Fo,Ft,Kd) in zip(self.dFout.dfields,
+                              self.forward_transforms,
+                              self.all_dkds):
             args = (Ft.full_output_buffer,) + tuple(Kd)
             F = self.build_diffusion_filter(Fo.dim, *args)
             diffusion_filters += (F,)
         self.diffusion_filters = diffusion_filters
-   
+
 
     @op_apply
     def apply(self, simulation, **kwds):
diff --git a/hysop/backend/host/python/operator/directional/stretching_dir.py b/hysop/backend/host/python/operator/directional/stretching_dir.py
index 259eedb43..ead3d9c42 100644
--- a/hysop/backend/host/python/operator/directional/stretching_dir.py
+++ b/hysop/backend/host/python/operator/directional/stretching_dir.py
@@ -14,6 +14,11 @@ from hysop.numerics.odesolvers.runge_kutta import ExplicitRungeKutta, Euler, RK2
 
 class PythonDirectionalStretching(DirectionalStretchingBase, HostDirectionalOperator):
 
+    @debug
+    def __new__(cls, formulation, **kwds):
+        return super(PythonDirectionalStretching, cls).__new__(cls,
+                formulation=formulation, **kwds)
+
     @debug
     def __init__(self, formulation, **kwds):
         super(PythonDirectionalStretching, self).__init__(formulation=formulation, **kwds)
diff --git a/hysop/backend/host/python/operator/poisson_curl.py b/hysop/backend/host/python/operator/poisson_curl.py
index 5e724b48c..3e8182938 100644
--- a/hysop/backend/host/python/operator/poisson_curl.py
+++ b/hysop/backend/host/python/operator/poisson_curl.py
@@ -9,7 +9,7 @@ from hysop.tools.numba_utils import make_numba_signature, prange
 from hysop.core.graph.graph import op_apply
 from hysop.backend.host.host_operator import HostOperator, OpenClMappable
 from hysop.operator.base.poisson_curl import SpectralPoissonCurlOperatorBase
-        
+
 from hysop.backend.host.python.operator.diffusion import PythonDiffusion
 from hysop.backend.host.python.operator.poisson import PythonPoisson
 from hysop.backend.host.python.operator.solenoidal_projection import PythonSolenoidalProjection
@@ -19,7 +19,7 @@ class PythonPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClMappable, HostOpe
     """
     Solves the poisson rotational equation using numpy fftw.
     """
-        
+
     @classmethod
     def build_filter_curl_2d__0_m(cls, FIN, K, FOUT, target=__DEFAULT_NUMBA_TARGET__):
         args=(FIN,K,FOUT)
@@ -51,7 +51,7 @@ class PythonPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClMappable, HostOpe
     def build_filter_curl_3d__0_n(cls, FIN, K, FOUT, target=__DEFAULT_NUMBA_TARGET__):
         args=(FIN,K,FOUT)
         signature, _ = make_numba_signature(*args)
-        layout='(n,m,p),(n)->(n,m,p)' 
+        layout='(n,m,p),(n)->(n,m,p)'
         @nb.guvectorize([signature], layout,
             target=target, nopython=True, cache=True)
         def filter_curl_3d__0_n(Fin, K, Fout):
@@ -93,7 +93,7 @@ class PythonPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClMappable, HostOpe
     def build_filter_curl_3d__1_n(cls, FIN, K, FOUT, target=__DEFAULT_NUMBA_TARGET__):
         args=(FIN,K,FOUT)
         signature, _ = make_numba_signature(*args)
-        layout='(n,m,p),(n)->(n,m,p)' 
+        layout='(n,m,p),(n)->(n,m,p)'
         @nb.guvectorize([signature], layout,
             target=target, nopython=True, cache=True)
         def filter_curl_3d__1_n(Fin, K, Fout):
@@ -107,7 +107,7 @@ class PythonPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClMappable, HostOpe
     def build_filter_curl_3d__1_m(cls, FIN, K, FOUT, target=__DEFAULT_NUMBA_TARGET__):
         args=(FIN,K,FOUT)
         signature, _ = make_numba_signature(*args)
-        layout='(n,m,p),(m)->(n,m,p)' 
+        layout='(n,m,p),(m)->(n,m,p)'
         @nb.guvectorize([signature], layout,
             target=target, nopython=True, cache=True)
         def filter_curl_3d__1_m(Fin, K, Fout):
@@ -121,7 +121,7 @@ class PythonPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClMappable, HostOpe
     def build_filter_curl_3d__1_p(cls, FIN, K, FOUT, target=__DEFAULT_NUMBA_TARGET__):
         args=(FIN,K,FOUT)
         signature, _ = make_numba_signature(*args)
-        layout='(n,m,p),(p)->(n,m,p)' 
+        layout='(n,m,p),(p)->(n,m,p)'
         @nb.guvectorize([signature], layout,
             target=target, nopython=True, cache=True)
         def filter_curl_3d__1_p(Fin, K, Fout):
@@ -130,17 +130,17 @@ class PythonPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClMappable, HostOpe
                     for k in range(0, Fin.shape[2]):
                         Fout[i,j,k] += K[k]*Fin[i,j,k]
         return functools.partial(filter_curl_3d__1_p, *args)
-        
+
     def setup(self, work):
         super(PythonPoissonCurl, self).setup(work=work)
-        
+
         dim = self.dim
         WIN, UIN, UOUT = self.WIN, self.UIN, self.UOUT
         K, KK = self.K, self.KK
         UK = self.UK
         assert len(WIN)==len(KK), (len(WIN),len(KK))
         assert len(UIN)==len(UOUT)==len(UK), (len(UIN),len(UOUT),len(UK))
-        
+
         # diffusion filters
         if self.should_diffuse:
             diffusion_filters = ()
@@ -166,7 +166,7 @@ class PythonPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClMappable, HostOpe
             F = PythonPoisson.build_poisson_filter(dim, *((Wi,)+KKi+(Wi,)))
             poisson_filters += (F,)
         self.poisson_filters = poisson_filters
-        
+
         # curl filter
         if (dim==2):
             curl0 = self.build_filter_curl_2d__0_m(UIN[0], UK[0], UOUT[0])
@@ -194,7 +194,7 @@ class PythonPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClMappable, HostOpe
     def apply(self, simulation=None, **kwds):
         """Apply analytic formula."""
         super(PythonPoissonCurl, self).apply(**kwds)
-        
+
         diffuse=self.should_diffuse
         project=self.do_project(simu=simulation)
 
@@ -216,8 +216,8 @@ class PythonPoissonCurl(SpectralPoissonCurlOperatorBase, OpenClMappable, HostOpe
                 for curl_filter in curl_filters:
                     curl_filter()
                 Bt(simulation=simulation)
-        
+
         self.dU.exchange_ghosts()
         if (diffuse or project):
             self.dW.exchange_ghosts()
-        self.update_energy(simulation=simulation) 
+        self.update_energy(simulation=simulation)
diff --git a/hysop/core/checkpoints.py b/hysop/core/checkpoints.py
index efa9697bd..12e7fd33b 100644
--- a/hysop/core/checkpoints.py
+++ b/hysop/core/checkpoints.py
@@ -446,7 +446,7 @@ class CheckpointHandler(object):
                 nbytes += np.prod(shape, dtype=np.int64) * dtype.itemsize
 
         if (root is not None):
-            root.attrs['nbytes'] = nbytes
+            root.attrs['nbytes'] = int(nbytes)
             msg='>Maximum checkpoint size will be {}, without compression and metadata.'
             vprint(root.tree())
             vprint(msg.format(bytes2str(nbytes)))
diff --git a/hysop/core/graph/computational_graph.py b/hysop/core/graph/computational_graph.py
index df3fa8692..e3400d790 100644
--- a/hysop/core/graph/computational_graph.py
+++ b/hysop/core/graph/computational_graph.py
@@ -94,7 +94,7 @@ class ComputationalGraph(ComputationalGraphNode, metaclass=ABCMeta):
                     reqs.enforce_unique_topology_shape, reqs.enforce_unique_transposition_state,
                     reqs.enforce_unique_ghosts, reqs.enforce_unique_memory_order,
                     ) + optypes
-            vals = map(str, vals)
+            vals = tuple(map(str, vals))
             values.append(vals)
 
         template = '\n   {:<{name_size}}   {:^{topology_size}}      {:^{tstates_size}}      {:^{ghosts_size}}      {:^{order_size}}      {:<{type_size0}}      {:<{type_size1}}      {:<{type_size2}}'
@@ -182,15 +182,14 @@ class ComputationalGraph(ComputationalGraphNode, metaclass=ABCMeta):
 
         titles = [[('OPERATOR', 'FIELD', 'DISCRETIZATION', 'GHOSTS',
                     'MEMORY ORDER', 'CAN_SPLIT', 'TSTATES')]]
-        name_size = max(len(s[0]) for ss in sinputs.values()+soutputs.values()+titles for s in ss)
-        field_size = max(len(s[1]) for ss in sinputs.values()+soutputs.values()+titles for s in ss)
-        discr_size = max(len(s[2]) for ss in sinputs.values()+soutputs.values()+titles for s in ss)
-        ghosts_size = max(len(s[3]) for ss in sinputs.values()+soutputs.values()+titles for s in ss)
-        order_size = max(len(s[4]) for ss in sinputs.values()+soutputs.values()+titles for s in ss)
-        cansplit_size = max(len(s[5]) for ss in sinputs.values() +
-                            soutputs.values()+titles for s in ss)
-        tstates_size = max(len(s[6]) for ss in sinputs.values() +
-                           soutputs.values()+titles for s in ss)
+        vals = tuple(sinputs.values()) + tuple(soutputs.values()) + tuple(titles)
+        name_size     = max(len(s[0]) for ss in vals for s in ss)
+        field_size    = max(len(s[1]) for ss in vals for s in ss)
+        discr_size    = max(len(s[2]) for ss in vals for s in ss)
+        ghosts_size   = max(len(s[3]) for ss in vals for s in ss)
+        order_size    = max(len(s[4]) for ss in vals for s in ss)
+        cansplit_size = max(len(s[5]) for ss in vals for s in ss)
+        tstates_size  = max(len(s[6]) for ss in vals for s in ss)
 
         template = '\n   {:<{name_size}}   {:^{field_size}}     {:^{discr_size}}      {:^{ghosts_size}}      {:^{order_size}}      {:^{cansplit_size}}      {:^{tstates_size}}'
 
@@ -370,7 +369,7 @@ class ComputationalGraph(ComputationalGraphNode, metaclass=ABCMeta):
                     field_topologies.setdefault(topo, []).append(node)
             for topo in sorted(field_topologies.keys(), key=lambda x: x.tag):
                 pnames = set(node.pretty_name for node in field_topologies[topo])
-                pnames = sorted(pnames)
+                pnames = tuple(sorted(pnames))
                 nbyline = 4
                 nentries = len(pnames)//nbyline
                 n0 = len(str(topo.backend.kind).lower())
@@ -393,8 +392,9 @@ class ComputationalGraph(ComputationalGraphNode, metaclass=ABCMeta):
                     topologies.setdefault(field, []).append(entries)
 
         titles = [[('BACKEND', 'TOPOLOGY', 'OPERATORS')]]
-        backend_size = max(len(s[0]) for ss in topologies.values()+titles for s in ss)
-        topo_size = max(len(s[1]) for ss in topologies.values()+titles for s in ss)
+        vals = tuple(topologies.values()) + tuple(titles)
+        backend_size = max(len(s[0]) for ss in vals for s in ss)
+        topo_size    = max(len(s[1]) for ss in vals for s in ss)
         template = '\n   {:<{backend_size}}   {:<{topo_size}}   {}'
         sizes = {'backend_size': backend_size,
                  'topo_size': topo_size}
diff --git a/hysop/core/graph/computational_node.py b/hysop/core/graph/computational_node.py
index 3ef5a57c7..5222ab09c 100644
--- a/hysop/core/graph/computational_node.py
+++ b/hysop/core/graph/computational_node.py
@@ -556,7 +556,7 @@ class ComputationalGraphNode(OperatorBase, metaclass=ABCMeta):
         if fills the 'None' domain.
         """
         domains = {}
-        for field in set(self.input_fields.keys()+self.output_fields.keys()):
+        for field in set(self.input_fields.keys()).union(self.output_fields.keys()):
             domains.setdefault(field.domain, set()).add(self)
         if self.is_domainless:
             domains.setdefault(None, set()).add(self)
@@ -696,7 +696,7 @@ class ComputationalGraphNode(OperatorBase, metaclass=ABCMeta):
     def get_topo_descriptor(cls, variables, field):
         if (field in variables):
             return variables[field]
-        tfields = filter(lambda x: x.is_tensor, variables.keys())
+        tfields = tuple(filter(lambda x: x.is_tensor, variables.keys()))
         for tfield in tfields:
             if field in tfield:
                 return variables[tfield]
diff --git a/hysop/core/graph/graph_builder.py b/hysop/core/graph/graph_builder.py
index 2a1a7bfbf..d8ea7e0a9 100644
--- a/hysop/core/graph/graph_builder.py
+++ b/hysop/core/graph/graph_builder.py
@@ -374,7 +374,7 @@ class GraphBuilder(object):
             for (fields, io_params, op_kwds) in target_node._output_fields_to_dump:
                 if not fields:
                     fields = self.output_fields.keys()
-                fields = sorted(fields, key=lambda x: x.name)
+                fields = tuple(sorted(fields, key=lambda x: x.name))
                 for field in fields:
                     msg='{} is not an output field.'.format(field.name)
                     assert field in self.output_fields, msg
@@ -668,7 +668,7 @@ class GraphBuilder(object):
                 if not write_nodes:
                     return
 
-                src_topos = write_nodes.keys()
+                src_topos = tuple(write_nodes.keys())
                 if (src_topo is not None):
                     assert src_topo in src_topos
                     src_topos = (src_topo,)
@@ -892,7 +892,7 @@ class GraphBuilder(object):
                 if (target_topo.backend.kind is Backend.HOST) and write_nodes:
                     # give source topo priority according to topology_affinity
                     src_topos = write_nodes.keys()
-                    src_topos = sorted(src_topos, key=topology_affinity, reverse=True)
+                    src_topos = tuple(sorted(src_topos, key=topology_affinity, reverse=True))
                     src_topo = src_topos[0]
                     if (src_topo is not target_topo):
                         msg='   >Redistributing field {} from up to date topologies {} '
@@ -907,7 +907,7 @@ class GraphBuilder(object):
                 elif (target_topo.backend.kind is Backend.OPENCL) and write_nodes:
                     # give source topo priority according to topology_affinity
                     src_topos = write_nodes.keys()
-                    src_topos = sorted(src_topos, key=topology_affinity, reverse=True)
+                    src_topos = tuple(sorted(src_topos, key=topology_affinity, reverse=True))
                     src_topo = src_topos[0]
                     if (src_topo is not target_topo):
                         msg='   >Redistributing field {} from up to date topologies {} '
diff --git a/hysop/core/graph/node_requirements.py b/hysop/core/graph/node_requirements.py
index 243ec3217..273c087ec 100644
--- a/hysop/core/graph/node_requirements.py
+++ b/hysop/core/graph/node_requirements.py
@@ -21,12 +21,12 @@ class NodeRequirements(object):
 
 
 class OperatorRequirements(NodeRequirements):
-    __slots__ = ('_enforce_unique_transposition_state', 
+    __slots__ = ('_enforce_unique_transposition_state',
                  '_enforce_unique_topology_shape',
                  '_enforce_unique_memory_order',
                  '_enforce_unique_ghosts')
 
-    def __init__(self, operator, 
+    def __init__(self, operator,
             enforce_unique_transposition_state=None,
             enforce_unique_topology_shape=None,
             enforce_unique_memory_order=None,
@@ -85,7 +85,7 @@ class OperatorRequirements(NodeRequirements):
         return self._enforce_unique_ghosts
     enforce_unique_ghosts = property(_get_enforce_unique_ghosts,
                                      _set_enforce_unique_ghosts)
-    
+
     def check_and_update_reqs(self, field_requirements):
         """Check and possibly update field requirements against global node requirements."""
         if (field_requirements is None):
@@ -93,7 +93,7 @@ class OperatorRequirements(NodeRequirements):
         axes = set()
         memory_order, can_split, min_ghosts, max_ghosts = None, None, None, None
         field_names = ''
-        for (is_input, freqs) in field_requirements.iter_requirements(): 
+        for (is_input, freqs) in field_requirements.iter_requirements():
             if (freqs is None):
                 continue
             (field, td, req) = freqs
@@ -103,8 +103,8 @@ class OperatorRequirements(NodeRequirements):
                         msg ='::GLOBAL OPERATOR REQUIREMENTS ERROR::\n'
                         msg+='Previous axes: \n  {} required by fields {}\nare incompatible '
                         msg+='with axes requirements \n  {} enforced by {} field {}.\n'
-                        msg=msg.format(tuple(axes),     field_names, 
-                                       tuple(req.axes), 
+                        msg=msg.format(tuple(axes),     field_names,
+                                       tuple(req.axes),
                                       'input' if is_input else 'output',
                                        field.name)
                         raise RuntimeError(msg)
@@ -119,8 +119,8 @@ class OperatorRequirements(NodeRequirements):
                             msg ='::GLOBAL OPERATOR REQUIREMENTS ERROR::\n'
                             msg+='Previous memory order: \n  {} required by fields {}\nare incompatible '
                             msg+='with memory order requirements \n  {} enforced by {} field {}.\n'
-                            msg=msg.format(memory_order, field_names, 
-                                           req.memory_order, 
+                            msg=msg.format(memory_order, field_names,
+                                           req.memory_order,
                                           'input' if is_input else 'output',
                                            field.name)
                             raise RuntimeError(msg)
@@ -136,8 +136,8 @@ class OperatorRequirements(NodeRequirements):
                             msg ='::GLOBAL OPERATOR REQUIREMENTS ERROR::\n'
                             msg+='Previous cartesian split directions: \n  {} required by fields {}\nare incompatible '
                             msg+='with cartesian split directions requirements \n  {} enforced by {} field {}.'
-                            msg=msg.format(can_split, field_names, 
-                                           req.can_split, 
+                            msg=msg.format(can_split, field_names,
+                                           req.can_split,
                                           'input' if is_input else 'output',
                                            field.name)
                             msg+='\nDomain cannot be splitted accross multiple processes.\n'
@@ -154,14 +154,14 @@ class OperatorRequirements(NodeRequirements):
                     if any(maxg<min_ghosts):
                         msg+='Previous cartesian min_ghosts: \n  {} required by fields {}\nare incompatible '
                         msg+='with cartesian max_ghosts requirements \n  {} enforced by {} field {}.\n'
-                        msg=msg.format(min_ghosts, field_names, 
+                        msg=msg.format(min_ghosts, field_names,
                                        maxg, 'input' if is_input else 'output',
                                        field.name)
                         raise RuntimeError(msg)
                     elif any(ming>max_ghosts):
                         msg+='Previous cartesian max_ghosts: \n  {} required by fields {}\nare incompatible '
                         msg+='with cartesian min_ghosts requirements \n  {} enforced by {} field {}.\n'
-                        msg=msg.format(max_ghosts, field_names, 
+                        msg=msg.format(max_ghosts, field_names,
                                        ming, 'input' if is_input else 'output',
                                        field.name)
                         raise RuntimeError(msg)
@@ -175,19 +175,18 @@ class OperatorRequirements(NodeRequirements):
             field_names += field.name + ', '
 
         axes = tuple(axes)
-        
+
         has_single_input = (len(field_requirements._input_field_requirements)<=1)
-    
+
         # enforce global operator requirements onto field requirements
-        for (is_input, freqs) in field_requirements.iter_requirements(): 
+        for (is_input, freqs) in field_requirements.iter_requirements():
             if (freqs is None):
                 continue
             (field, td, req) = freqs
-            
+
             # enforce transposition axes
             if self.enforce_unique_transposition_state:
                 if axes:
-                    #axes = sorted(axes, key=lambda x: sum((xi-i)**2 for (i,xi) in zip(range(len(x)),x)))
                     req.axes = (axes[0],)
                 elif not has_single_input:
                     req.axes = (TranspositionState[field.dim].default_axes(),)
@@ -198,7 +197,7 @@ class OperatorRequirements(NodeRequirements):
                     req.memory_order = memory_order
                 elif not has_single_input:
                     req.memory_order = MemoryOrdering.C_CONTIGUOUS
-            
+
             # enforce topology shape (indirectly by enforcing split directions)
             if self.enforce_unique_topology_shape:
                 assert (can_split is not None)
@@ -213,4 +212,4 @@ class OperatorRequirements(NodeRequirements):
             if self.enforce_unique_ghosts:
                 req.min_ghosts = min_ghosts
                 req.max_ghosts = min_ghosts
-            
+
diff --git a/hysop/core/memory/memory_request.py b/hysop/core/memory/memory_request.py
index 9b80a06a9..98387c3d6 100644
--- a/hysop/core/memory/memory_request.py
+++ b/hysop/core/memory/memory_request.py
@@ -300,7 +300,7 @@ class MultipleOperatorMemoryRequests(object):
     def operators(self):
         ops = []
         for requests in self._all_requests_per_backend.values():
-            ops += requests.keys()
+            ops += list(requests.keys())
         return ops
 
     def __iadd__(self, other):
diff --git a/hysop/core/mpi/topo_tools.py b/hysop/core/mpi/topo_tools.py
index a109623b1..07cec6a77 100644
--- a/hysop/core/mpi/topo_tools.py
+++ b/hysop/core/mpi/topo_tools.py
@@ -156,7 +156,6 @@ class TopoTools(object):
         # Get the list of processes
         assert child is not None
         assert parent is not None
-        #child_ranks = [i for i in range(child.Get_size())]
         child_group = child.Get_group()
         parent_group = parent.Get_group()
         inter_group = MPI.Group.Intersect(child_group, parent_group)
diff --git a/hysop/fields/continuous_field.py b/hysop/fields/continuous_field.py
index c378e015b..a6339c78a 100644
--- a/hysop/fields/continuous_field.py
+++ b/hysop/fields/continuous_field.py
@@ -372,7 +372,7 @@ class FieldContainerI(TaggedObject):
                         return False
                 return True
             if isinstance(a, dict):
-                for k in set(a.keys()+b.keys()):
+                for k in set(a.keys()).union(b.keys()):
                     if (k not in a) or (k not in b):
                         return False
                     ak, bk = a[k], b[k]
diff --git a/hysop/fields/discrete_field.py b/hysop/fields/discrete_field.py
index de291ea94..c979bfd87 100644
--- a/hysop/fields/discrete_field.py
+++ b/hysop/fields/discrete_field.py
@@ -286,7 +286,7 @@ class DiscreteScalarFieldViewContainerI(object, metaclass=ABCMeta):
                         return False
                 return True
             if isinstance(a, dict):
-                for k in set(a.keys()+b.keys()):
+                for k in set(a.keys()).union(b.keys()):
                     if (k not in a) or (k not in b):
                         return False
                     ak, bk = a[k], b[k]
diff --git a/hysop/fields/ghost_exchangers.py b/hysop/fields/ghost_exchangers.py
index 1866d6325..99bcbfa87 100644
--- a/hysop/fields/ghost_exchangers.py
+++ b/hysop/fields/ghost_exchangers.py
@@ -27,7 +27,7 @@ def gprint_buffer(msg, buf, *args, **kwds):
     else:
         mpi_type=None
     gprint('{}: mpi_type={}, shape={}, dtype={}, c_contiguous={}, f_contiguous={}'.format(msg,
-        mpi_type, buf.shape, buf.dtype, buf.flags['C_CONTIGUOUS'], buf.flags['F_CONTIGUOUS']))
+        mpi_type, buf.shape, buf.dtype, buf.flags.c_contiguous, buf.flags.f_contiguous))
     if no_data:
         return
     gprint('   ' + '\n   '.join(str(buf).split('\n')))
diff --git a/hysop/fields/tests/test_cartesian.py b/hysop/fields/tests/test_cartesian.py
index a5a415458..baacac2e8 100644
--- a/hysop/fields/tests/test_cartesian.py
+++ b/hysop/fields/tests/test_cartesian.py
@@ -511,16 +511,16 @@ def test_mpi_ghost_exchange_periodic(comm=None):
             if rank==0:
                 print('    >DTYPE={}'.format(dtype))
 
-            F0 = Field(domain=domain, name='F0', nb_components=1, dtype=dtype, _register=False)
-            F1 = Field(domain=domain, name='F1', nb_components=2, dtype=dtype, _register=False)
-            F2 = Field(domain=domain, name='F2', shape=(2,2), dtype=dtype, _register=False)
+            F0 = Field(domain=domain, name='F0', nb_components=1, dtype=dtype)
+            F1 = Field(domain=domain, name='F1', nb_components=2, dtype=dtype)
+            F2 = Field(domain=domain, name='F2', shape=(2,2), dtype=dtype)
 
             for (backend, cl_env) in iter_backends():
                 if rank==0:
                     print('      >BACKEND.{}{}'.format(backend,
                             '' if (cl_env is None) else '::{}.{}'.format(
                                 cl_env.platform.name.strip(), cl_env.device.name.strip())))
-                for shape in it.product(range(0,size+1), repeat=dim):
+                for shape in it.product(range(0, size+1), repeat=dim):
                     if np.prod(shape, dtype=np.uint32)!=size:
                         continue
                     if rank==0:
@@ -610,14 +610,14 @@ def test_mpi_ghost_exchange_runtime(comm=None):
         print('*'*len(msg))
         print('test_mpi_ghost_exchange_runtime()'.format(size))
 
-    for dim in range(1,3+__ENABLE_LONG_TESTS__):
+    for dim in range(1, 3+__ENABLE_LONG_TESTS__):
         if rank==0:
             sys.stdout.write('>DIM={}\n'.format(dim))
 
         npts = (17,16,19)[:dim]
         nghosts = (2,1,3)[:dim]
 
-        for shape in it.product(range(0,size+1), repeat=dim):
+        for shape in it.product(range(0, size+1), repeat=dim):
             if np.prod(shape, dtype=np.uint32)!=size:
                 continue
             if rank==0:
@@ -638,11 +638,9 @@ def test_mpi_ghost_exchange_runtime(comm=None):
                 brk = False
                 try:
                     for (lbd, rbd) in domain_boundary_iterator(dim):
-                        domain = Box(dim=dim, lboundaries=lbd,
-                                              rboundaries=rbd)
+                        domain = Box(dim=dim, lboundaries=lbd, rboundaries=rbd)
 
-                        F = Field(domain=domain, name='F', nb_components=1,
-                                    dtype=dtype, _register=False)
+                        F = Field(domain=domain, name='F', nb_components=1, dtype=dtype)
 
                         discretization = CartesianDiscretization(npts, nghosts,
                                                 lboundaries=F.lboundaries,
@@ -708,9 +706,9 @@ def test_mpi_ghost_accumulate_periodic(comm=None):
             if rank==0:
                 print('    >DTYPE={}'.format(dtype))
 
-            F0 = Field(domain=domain, name='F0', nb_components=1, dtype=dtype, _register=False)
-            F1 = Field(domain=domain, name='F1', nb_components=2, dtype=dtype, _register=False)
-            F2 = Field(domain=domain, name='F2', shape=(2,2), dtype=dtype, _register=False)
+            F0 = Field(domain=domain, name='F0', nb_components=1, dtype=dtype)
+            F1 = Field(domain=domain, name='F1', nb_components=2, dtype=dtype)
+            F2 = Field(domain=domain, name='F2', shape=(2,2), dtype=dtype)
 
             for (backend, cl_env) in iter_backends():
                 if rank==0:
@@ -756,7 +754,7 @@ def test_mpi_ghost_accumulate_periodic(comm=None):
                                                 ExchangeMethod.NEIGHBOR_ALL_TO_ALL_V,
                                                 ExchangeMethod.NEIGHBOR_ALL_TO_ALL_W,):
                             if rank==0:
-                                print('             ExchangeMethod.{:<25} |'.format(exchange_method), end=' ')
+                                print('             ExchangeMethod.{:<25} |'.format(str(exchange_method)), end=' ')
                             sys.stdout.flush()
 
                             # test one direction at a time
@@ -804,6 +802,7 @@ def test_mpi_ghost_accumulate_periodic(comm=None):
 
                                             overlaping_neighbours = set(tuple((np.asarray(mask, dtype=np.int32)*displacements).tolist()) for mask in masks)
                                             overlaping_neighbours = filter(lambda x: 0<sum(_!=0 for _ in x)<=max_displacements, overlaping_neighbours) #diagonals
+                                            overlaping_neighbours = tuple(overlaping_neighbours)
 
                                             expected_values = (rank+1)*(i+1)
                                             for disp in overlaping_neighbours:
diff --git a/hysop/mesh/cartesian_mesh.py b/hysop/mesh/cartesian_mesh.py
index 9249db90e..8053befc8 100644
--- a/hysop/mesh/cartesian_mesh.py
+++ b/hysop/mesh/cartesian_mesh.py
@@ -703,7 +703,7 @@ class CartesianMeshView(MeshView):
             def __init__(self, dim, ghosts, compute_resolution, local_compute_indices):
                 iter_shape = tuple(compute_resolution[:axes])
                 gI = np.ix_(*local_compute_indices[axes:])
-                I = tuple(gI[i]-ghosts[j]  for i,j in enumerate(range(axes, dim)))
+                I = tuple(gI[i]-ghosts[j] for i,j in enumerate(range(axes, dim)))
 
                 self._iter_shape = iter_shape
                 self._data = (dim, ghosts, I, gI)
diff --git a/hysop/numerics/fft/gpyfft_fft.py b/hysop/numerics/fft/gpyfft_fft.py
index e5de67e09..abf665291 100644
--- a/hysop/numerics/fft/gpyfft_fft.py
+++ b/hysop/numerics/fft/gpyfft_fft.py
@@ -563,7 +563,7 @@ Post callback source code:
             else:
                 oip += ('if (b>={}) {{ return; }};'.format(fake_batchsize, fp),)
             oip += ('{K} -= {b}*{D};'.format(K=K, b=b, D=D),)
-        for i in range(ndim-1,-1,-1):
+        for i in range(ndim-1, -1, -1):
             Ki = idx.format('xyz'[i])
             Si = S[i]
             oip += ('const uint {Ki} = {K}/{Si};'.format(Ki=Ki, K=K, Si=Si),)
@@ -582,7 +582,7 @@ Post callback source code:
             real_offset = ('{base_offset}uL'.format(base_offset=base_offset),)
             if (fake_batchsize > 1):
                 real_offset += ('{b}*{D}'.format(b=b, D=real_distance),)
-            for i in range(ndim-1,0,-1):
+            for i in range(ndim-1, 0, -1):
                 Ki = idx.format('xyz'[i])
                 Si = real_strides[i]
                 real_offset += ('{Ki}*{Si}'.format(Ki, Si),)
diff --git a/hysop/numerics/fft/host_fft.py b/hysop/numerics/fft/host_fft.py
index ceb5c70e6..2e19a8e31 100644
--- a/hysop/numerics/fft/host_fft.py
+++ b/hysop/numerics/fft/host_fft.py
@@ -70,9 +70,9 @@ class HostFFTI(FFTI):
         super(HostFFTI, self).__init__(backend=backend, **kwds)
 
     @classmethod
-    def default_interface(cls, threads=None, 
+    def default_interface(cls, threads=None,
                        backend=None, allocator=None,
-                       planner_effort=None, 
+                       planner_effort=None,
                        planning_timelimit=None,
                        destroy_input=False,
                        warn_on_allocation=True,
@@ -96,8 +96,8 @@ class HostFFTI(FFTI):
             threads            = first_not_None(threads,            __FFTW_NUM_THREADS__)
             planner_effort     = first_not_None(planner_effort,     __FFTW_PLANNER_EFFORT__)
             planning_timelimit = first_not_None(planning_timelimit, __FFTW_PLANNER_TIMELIMIT__)
-            return FftwFFT(threads=threads,  
-                           planner_effort=planner_effort, 
+            return FftwFFT(threads=threads,
+                           planner_effort=planner_effort,
                            planning_timelimit=planning_timelimit,
                            backend=backend, allocator=allocator,
                            destroy_input=destroy_input,
@@ -105,19 +105,19 @@ class HostFFTI(FFTI):
                            warn_on_misalignment=warn_on_misalignment,
                            error_on_allocation=error_on_allocation,
                            **kwds)
-    
+
     def new_queue(self, tg, name):
         return HostFFTQueue(name=name)
-                       
+
     def plan_copy(self, tg, src, dst):
         src = self.ensure_callable(src)
         dst = self.ensure_callable(dst)
-        
+
         @static_vars(numba_copy=None)
         def exec_copy(src=src, dst=dst):
             src, dst = src(), dst()
             if can_exec_hptt(src, dst):
-                hptt.tensorTransposeAndUpdate(perm=range(src.ndim),
+                hptt.tensorTransposeAndUpdate(perm=tuple(range(src.ndim)),
                         alpha=1.0, A=src, beta=0.0, B=dst)
             elif HAS_NUMBA:
                 if (exec_copy.numba_copy is None):
@@ -126,11 +126,11 @@ class HostFFTI(FFTI):
             else:
                 dst[...] = src
         return exec_copy
-    
+
     def plan_accumulate(self, tg, src, dst):
         src = self.ensure_callable(src)
         dst = self.ensure_callable(dst)
-        
+
         @static_vars(numba_accumulate=None)
         def exec_accumulate(src=src, dst=dst):
             src, dst = src(), dst()
@@ -162,7 +162,7 @@ class HostFFTI(FFTI):
             else:
                 dst[...] = np.transpose(a=src, axes=axes)
         return exec_transpose
-    
+
     def plan_fill_zeros(self, tg, a, slices):
         assert slices
         a = self.ensure_callable(a)
@@ -171,21 +171,21 @@ class HostFFTI(FFTI):
             for slc in slices:
                 buf[slc] = 0
         return exec_fill_zeros
-    
-    def plan_compute_energy(self, tg, fshape, src, dst, transforms, 
+
+    def plan_compute_energy(self, tg, fshape, src, dst, transforms,
             method='round', target=None, **kwds):
         """
-        Plan to compute energy from src array to dst array using given transforms, 
+        Plan to compute energy from src array to dst array using given transforms,
         method (round or weighted) and numba target.
         """
         (N, NS2, C2C, R2C, S) = super(HostFFTI, self).plan_compute_energy(tg, fshape, src, dst, transforms, **kwds)
         dim = src.ndim
-        
+
         # We do not forget the hermitian symmetry:
-        #  normally: E = 1/2 |Fi|**2 
+        #  normally: E = 1/2 |Fi|**2
         #  but for a R2C transform:  E = 1/2 |Fi|**2 + 1/2 |Fi*|**2 = 1 |Fi|**2
         C = 2**(R2C-1)
-        
+
         target = first_not_None(target, __DEFAULT_NUMBA_TARGET__)
         args = (src,)+N+NS2+C2C+(C, S, dst)
         signature, layout  = make_numba_signature(*args)
diff --git a/hysop/numerics/interpolation/polynomial.py b/hysop/numerics/interpolation/polynomial.py
index c1e2e6277..3f2d446ff 100644
--- a/hysop/numerics/interpolation/polynomial.py
+++ b/hysop/numerics/interpolation/polynomial.py
@@ -74,7 +74,7 @@ class PolynomialInterpolator(object):
         else:
             msg='Unknown PolynomialInterpolation value {}.'.format(pi)
             raise NotImplementedError(msg)
-        for i in range(1,5):
+        for i in range(1, 5):
             if spi.endswith(str(2*i)):
                 fd=2*i
                 break
@@ -321,7 +321,7 @@ class PolynomialInterpolator(object):
             warnings.warn(msg, HysopCacheWarning)
 
         P0 = 0
-        for idx in it.product(*tuple(range(0,pi) for pi in p)):
+        for idx in it.product(*tuple(range(0, pi) for pi in p)):
             P0 += pvals[idx] * np.prod(np.power(xvals, idx))
         self.P0 = P0
 
@@ -341,7 +341,7 @@ class PolynomialInterpolator(object):
             print('\nBuilding system...')
 
         eqs = []
-        for dvec in it.product(*tuple(range(0,ki+1) for ki in k)):
+        for dvec in it.product(*tuple(range(0, ki+1) for ki in k)):
             if verbose:
                 print('  => derivative {}'.format(dvec))
 
@@ -354,7 +354,7 @@ class PolynomialInterpolator(object):
                 print('     stencil:')
                 print(stencil)
 
-            for idx in it.product(*tuple(range(gi,gi+2) for gi in ghosts)):
+            for idx in it.product(*tuple(range(gi, gi+2) for gi in ghosts)):
                 if verbose:
                     print('    -> point {}'.format(idx))
 
@@ -571,15 +571,15 @@ class PolynomialSubgridInterpolator(object):
             msg='Failed to load data from cache because:\n{}'.format(e)
             warnings.warn(msg, HysopCacheWarning)
 
-        X = tuple(np.asarray(tuple(sm.Rational(j,gr) for j in range(0,si)))
+        X = tuple(np.asarray(tuple(sm.Rational(j,gr) for j in range(0, si)))
                 for i,(gr,si) in enumerate(zip(gr, s)))
         V = np.vander(X[0], N=p[0], increasing=True)
         for i in range(1, dim):
             Vi = np.vander(X[i], N=p[i], increasing=True)
             V = np.multiply.outer(V,Vi)
 
-        even_axes = range(0,V.ndim,2)
-        odd_axes  = range(1,V.ndim,2)
+        even_axes = tuple(range(0, V.ndim,2))
+        odd_axes  = tuple(range(1, V.ndim,2))
         axes = even_axes + odd_axes
 
         V = np.transpose(V, axes=axes).copy()
diff --git a/hysop/numerics/odesolvers/runge_kutta.py b/hysop/numerics/odesolvers/runge_kutta.py
index 52c3d3fcb..ca555e42b 100644
--- a/hysop/numerics/odesolvers/runge_kutta.py
+++ b/hysop/numerics/odesolvers/runge_kutta.py
@@ -56,10 +56,10 @@ class ExplicitRungeKutta(RungeKutta):
         """Buffers = dict of (nb_stages+1) np.ndarray of size compatible with Xin"""
         check_instance(Xin, dict, keys=str, values=np.ndarray)
         varnames = Xin.keys()
-        views   = first_not_None(views,   { k:Ellipsis for (k,v) in Xin.items() })
-        Xout    = first_not_None(Xout,    { k:np.empty_like(v[views[k]])
+        views   = first_not_None(views,   { k: Ellipsis for (k,v) in Xin.items() })
+        Xout    = first_not_None(Xout,    { k: np.empty_like(v[views[k]])
                                                 for (k,v) in Xin.items() })
-        buffers = first_not_None(buffers, { k:tuple(np.empty_like(v)
+        buffers = first_not_None(buffers, { k: tuple(np.empty_like(v)
                                                 for i in range(self.stages+1))
                                                 for (k,v) in Xin.items()})
         check_instance(views,   dict, keys=str, values=(type(Ellipsis),slice,tuple))
@@ -67,7 +67,7 @@ class ExplicitRungeKutta(RungeKutta):
         check_instance(buffers, dict, keys=str, values=tuple)
         assert callable(RHS), type(RHS)
         if __debug__:
-            compute_shape = Xout.values()[0].shape
+            compute_shape = next(iter(Xout.values())).shape
             assert varnames == Xout.keys()
             assert varnames == views.keys()
             assert varnames == buffers.keys()
diff --git a/hysop/numerics/remesh/kernel_generator.py b/hysop/numerics/remesh/kernel_generator.py
index 3711b9042..1fd7be909 100644
--- a/hysop/numerics/remesh/kernel_generator.py
+++ b/hysop/numerics/remesh/kernel_generator.py
@@ -76,7 +76,7 @@ class Kernel(object):
             print()
 
         #split polynomials
-        X = np.arange(-Ms,+Ms+1)
+        X = np.arange(-Ms, +Ms+1)
         if split_polys:
             Pt_l = []
             Pt_r = []
@@ -335,10 +335,6 @@ class SymmetricKernelGenerator(object):
             Pp.append(sm.polys.polytools.poly(pexpr,x))
             Pm.append(sm.polys.polytools.poly(mexpr,x))
         P    = Pm[::-1] + Pp
-        Pbeg   = -Ms
-        Pend   = +Ms-1
-        Prange = lambda: range(-Ms,+Ms)
-        Penum  = lambda: enumerate(P,start=Pbeg)
 
         #precompute the r first polynomial derivatives
         dPs = []
@@ -354,7 +350,7 @@ class SymmetricKernelGenerator(object):
 
         ### Continuity equations (Gamma is Cr continuous)
         # Parity in x=0 -- Gamma is EVEN     -> (r+1)//2 eqs
-        for d in range(1,r+1,2):
+        for d in range(1, r+1,2):
             eq = coeffs[0][d] # =0
             eqs.append(eq)
         # Right-most point, zero derivatives -> (r+1)    eqs
@@ -374,9 +370,9 @@ class SymmetricKernelGenerator(object):
 
         ### Discrete moments
         s = sm.symbols('s')
-        for m in range(0,n+1):
+        for m in range(0, n+1):
             expr = f2q(0)
-            for l in range(-Ms+1,Ms+1):
+            for l in range(-Ms+1, Ms+1):
                 if m>0 and l==0: continue
                 i = Ms-l
                 e = P[i].xreplace({x:s-f2q(l)}).as_expr()
@@ -422,7 +418,7 @@ if __name__=='__main__':
     sg = StencilGenerator()
     sg.configure(dim=1, derivative=2)
 
-    for i in range(1,5):
+    for i in range(1, 5):
         p = 2*i
         H = sg.generate_exact_stencil(order=p-1, origin=i)
         print(H.coeffs)
diff --git a/hysop/numerics/stencil/stencil.py b/hysop/numerics/stencil/stencil.py
index 7135ac813..85dd6c1a6 100644
--- a/hysop/numerics/stencil/stencil.py
+++ b/hysop/numerics/stencil/stencil.py
@@ -140,7 +140,7 @@ class Stencil(object):
         assert sdim<=adim, 'Stencil dimension greater than array dimension.'
         assert set(symbols.keys())==self.variables(), 'Missing symbols {}.'.format(self.variables()-set(symbols.keys()))
         out  = first_not_None(out, np.empty_like(a[iview]))
-        axis = first_not_None(to_tuple(axis), range(adim)[-sdim:])
+        axis = first_not_None(to_tuple(axis), tuple(range(adim))[-sdim:])
         assert len(axis) == sdim
         assert out.ndim  == a.ndim
         assert out.shape == a[iview].shape
diff --git a/hysop/numerics/stencil/stencil_generator.py b/hysop/numerics/stencil/stencil_generator.py
index 7c8fe7391..2d78ac6dd 100644
--- a/hysop/numerics/stencil/stencil_generator.py
+++ b/hysop/numerics/stencil/stencil_generator.py
@@ -458,7 +458,7 @@ class StencilGenerator(object):
             def preficate(it):
                 return (it<=N).all() and sum(it)==n
 
-            nn = range(n+1)
+            nn  = range(n+1)
             itd = it.product(nn, repeat=dim)
             itd = filter(preficate, itd)
             expr = 0
diff --git a/hysop/numerics/tests/bench_fft.py b/hysop/numerics/tests/bench_fft.py
index 7b625c1bb..ee6d45fc0 100644
--- a/hysop/numerics/tests/bench_fft.py
+++ b/hysop/numerics/tests/bench_fft.py
@@ -141,7 +141,7 @@ class BenchFFT(object):
     def iter_shapes(self, offset=0):
         minj=12
         maxj=27
-        for j in range(minj,maxj):
+        for j in range(minj, maxj):
             shape = (2**j+offset,)
             cshape = list(shape)
             cshape[-1] = cshape[-1]//2 + 1
diff --git a/hysop/numerics/tests/test_fft.py b/hysop/numerics/tests/test_fft.py
index ca4dc8b83..d4923c894 100644
--- a/hysop/numerics/tests/test_fft.py
+++ b/hysop/numerics/tests/test_fft.py
@@ -636,7 +636,7 @@ class TestFFT(object):
             base = 2+i
             print('  '+msg[i])
             for ghosts in ((0,0,0),): #(2,0,0),(0,1,0),(0,0,3)):
-                for j1 in range(minj[i],maxj[i]):
+                for j1 in range(minj[i], maxj[i]):
                     shape = (3,2,base**j1,)
                     cshape = list(shape)
                     cshape[-1] = cshape[-1]//2 + 1
diff --git a/hysop/operator/analytic.py b/hysop/operator/analytic.py
index 7d243477d..8f6da9093 100644
--- a/hysop/operator/analytic.py
+++ b/hysop/operator/analytic.py
@@ -21,7 +21,7 @@ class AnalyticField(ComputationalGraphNodeGenerator):
         """
         AnalyticField operator frontend.
 
-        Apply a user-defined formula onto a field, possibly 
+        Apply a user-defined formula onto a field, possibly
         dependent on space variables and external fields/parameters.
 
         Parameters
@@ -50,23 +50,23 @@ class AnalyticField(ComputationalGraphNodeGenerator):
         check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
         extra_input_kwds = first_not_None(extra_input_kwds, {})
         base_kwds = first_not_None(base_kwds, {})
-            
+
         assert 'extra_input_kwds' not in kwds
         assert 'component' not in kwds
         assert 'coords' not in kwds
 
         if (implementation is Implementation.PYTHON) and (extra_input_kwds is not None):
-            candidate_input_tensors = filter(lambda f: isinstance(f, Field), extra_input_kwds.values())
+            candidate_input_tensors = tuple(filter(lambda f: isinstance(f, Field), extra_input_kwds.values()))
         else:
             extra_input_kwds = {}
-            candidate_input_tensors = ()
+            candidate_input_tensors = tuple()
         candidate_output_tensors = (field,)
-        
+
         formula = to_tuple(formula)
         if len(formula) == 1:
             formula = formula*field.nb_components
         check_instance(formula, tuple, size=field.nb_components)
-        
+
         super(AnalyticField, self).__init__(
                 candidate_input_tensors=candidate_input_tensors,
                 candidate_output_tensors=candidate_output_tensors,
@@ -78,7 +78,7 @@ class AnalyticField(ComputationalGraphNodeGenerator):
         self._extra_input_kwds = extra_input_kwds
         self._implementation = implementation
         self._kwds = kwds
-    
+
     @debug
     def _generate(self):
         nodes = []
@@ -95,7 +95,7 @@ class AnalyticField(ComputationalGraphNodeGenerator):
                     field=field,
                     formula=formula,
                     variables=variables,
-                    implementation=impl, 
+                    implementation=impl,
                     extra_input_kwds=extra_input_kwds,
                     **kwds)
             nodes.append(node)
@@ -118,7 +118,7 @@ class AnalyticScalarField(ComputationalGraphNodeFrontend):
                 Implementation.OPENCL: OpenClAnalyticField
         }
         return __implementations
-    
+
     @classmethod
     def default_implementation(cls):
         return Implementation.PYTHON
@@ -129,7 +129,7 @@ class AnalyticScalarField(ComputationalGraphNodeFrontend):
         """
         AnalyticField operator frontend.
 
-        Apply a user-defined formula onto a field, possibly 
+        Apply a user-defined formula onto a field, possibly
         dependent on space variables and external fields/parameters.
 
         Parameters
@@ -158,6 +158,6 @@ class AnalyticScalarField(ComputationalGraphNodeFrontend):
         if (implementation is Implementation.PYTHON):
             kwds['extra_input_kwds'] = extra_input_kwds
 
-        super(AnalyticScalarField, self).__init__(field=field, formula=formula, 
-                variables=variables, implementation=implementation, 
+        super(AnalyticScalarField, self).__init__(field=field, formula=formula,
+                variables=variables, implementation=implementation,
                 base_kwds=base_kwds, **kwds)
diff --git a/hysop/operator/base/convergence.py b/hysop/operator/base/convergence.py
index eeca70420..7d3cfcdbc 100644
--- a/hysop/operator/base/convergence.py
+++ b/hysop/operator/base/convergence.py
@@ -49,7 +49,7 @@ class ConvergenceBase(object):
         check_instance(convergence, TensorParameter, allow_none=True)
         check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
 
-        field = variables.keys()[0]
+        field = next(iter(variables))
         if convergence is None:
             convergence = TensorParameter(name="|residual|", dtype=field.dtype,
                                           shape=(field.nb_components, ),
diff --git a/hysop/operator/base/curl.py b/hysop/operator/base/curl.py
index f6651ca53..e586e32fb 100644
--- a/hysop/operator/base/curl.py
+++ b/hysop/operator/base/curl.py
@@ -19,14 +19,19 @@ class CurlOperatorBase(SpectralOperatorBase):
     """
     Compute the curl using a specific implementation.
     """
-    
+
+    @debug
+    def __new__(cls, Fin, Fout, variables, **kwds):
+        return super(CurlOperatorBase, cls).__new__(cls,
+                input_fields=None, output_fields=None, **kwds)
+
     @debug
     def __init__(self, Fin, Fout, variables, **kwds):
         """
         Create an operator that computes the curl of an input field Fin.
 
         Given Fin, a 2D ScalarField, a 2D VectorField or a 3D VectorField, compute Fout = curl(Fin).
-        
+
         Only the following configurations are supported:
                  dim   nb_components  |   dim   nb_components
         Input:    2        (1,2)      |    3          3
@@ -45,7 +50,7 @@ class CurlOperatorBase(SpectralOperatorBase):
         kwds: dict, optional
             Extra parameters passed towards base class (MultiSpaceDerivatives).
         """
-        
+
         check_instance(Fin,  Field)
         check_instance(Fout, Field)
         check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
@@ -83,18 +88,18 @@ class CurlOperatorBase(SpectralOperatorBase):
         else:
             msg='Unsupported dimension {}.'.format(dim)
             raise RuntimeError(msg)
-        
+
         # input and output fields
         input_fields  = { Fin:  variables[Fin] }
         output_fields = { Fout: variables[Fout] }
 
-        super(CurlOperatorBase, self).__init__(input_fields=input_fields, 
+        super(CurlOperatorBase, self).__init__(input_fields=input_fields,
                 output_fields=output_fields, **kwds)
 
         self.Fin  = Fin
         self.Fout = Fout
         self.dim  = dim
-        
+
     @debug
     def discretize(self):
         if self.discretized:
@@ -109,19 +114,22 @@ class SpectralCurlOperatorBase(CurlOperatorBase):
     """
     Compute the curl using a specific spectral implementation.
     """
+    def __new__(cls, **kwds):
+        return super(SpectralCurlOperatorBase, cls).__new__(cls, **kwds)
+
     def __init__(self, **kwds):
         super(SpectralCurlOperatorBase, self).__init__(**kwds)
-        
+
         dim = self.dim
         Fin, Fout = self.Fin, self.Fout
         assert Fin is not Fout, 'Cannot compute curl inplace!'
-        
+
         if (dim==2):
             tg0 = self.new_transform_group()
             tg1 = self.new_transform_group()
             if (Fin.nb_components==1):
                 assert (Fout.nb_components==2)
-                F0 = tg0.require_forward_transform(Fin, axes=0, 
+                F0 = tg0.require_forward_transform(Fin, axes=0,
                         custom_output_buffer='auto')
                 B0 = tg0.require_backward_transform(Fout[0], axes=0,
                         custom_input_buffer='auto', matching_forward_transform=F0)
@@ -130,7 +138,7 @@ class SpectralCurlOperatorBase(CurlOperatorBase):
                         custom_output_buffer='auto')
                 B1 = tg1.require_backward_transform(Fout[1], axes=1,
                         custom_input_buffer='auto', matching_forward_transform=F1)
-            
+
                 expr0 = Assignment(B0.s, +F0.s.diff(F0.s.frame.coords[1]))
                 expr1 = Assignment(B1.s, -F1.s.diff(F1.s.frame.coords[0]))
             elif (Fin.nb_components==2):
@@ -141,10 +149,10 @@ class SpectralCurlOperatorBase(CurlOperatorBase):
 
                 F1 = tg1.require_forward_transform(Fin[0], axes=0,
                         custom_output_buffer='auto')
-                B1 = tg1.require_backward_transform(Fout, axes=0, 
+                B1 = tg1.require_backward_transform(Fout, axes=0,
                         custom_input_buffer='auto', matching_forward_transform=F1,
                         action=SpectralTransformAction.ACCUMULATE)
-                
+
                 expr0 = Assignment(B0.s, +F0.s.diff(F0.s.frame.coords[0]))
                 expr1 = Assignment(B1.s, -F1.s.diff(F1.s.frame.coords[1]))
             else:
@@ -212,11 +220,11 @@ class SpectralCurlOperatorBase(CurlOperatorBase):
             backward_transforms = (B0, B1, B2, B3, B4, B5)
         else:
             raise NotImplementedError
-       
+
         self.forward_transforms  = forward_transforms
         self.backward_transforms = backward_transforms
         self.K = K
-    
+
     @debug
     def discretize(self):
         if self.discretized:
@@ -224,13 +232,13 @@ class SpectralCurlOperatorBase(CurlOperatorBase):
         super(SpectralCurlOperatorBase, self).discretize()
         dFin, dFout = self.dFin, self.dFout
         K = self.K
-        
+
         dK = ()
         for (tg,Ki) in K:
             _, dKi, _ = tg.discrete_wave_numbers[Ki]
             dK += (dKi,)
         self.dK = dK
-    
+
     def setup(self, work):
         super(SpectralCurlOperatorBase, self).setup(work)
 
diff --git a/hysop/operator/base/custom_symbolic_operator.py b/hysop/operator/base/custom_symbolic_operator.py
index f83231f7c..3a3df48c8 100644
--- a/hysop/operator/base/custom_symbolic_operator.py
+++ b/hysop/operator/base/custom_symbolic_operator.py
@@ -351,7 +351,7 @@ class SymbolicExpressionInfo(object):
                 if ((k in self.input_dfields) and (self.input_dfields[k].dfield is v.dfield))}
         self.stencils = SortedDict()
 
-        dfields  = tuple(input_dfields.values()) + tuple(output_dfields.values())
+        dfields = tuple(input_dfields.values()) + tuple(output_dfields.values())
         if (force_symbolic_axes is not None):
             if isinstance(force_symbolic_axes, tuple):
                 axes = force_symbolic_axes
@@ -393,7 +393,7 @@ class SymbolicExpressionInfo(object):
 
     def check_arrays(self):
         compute_resolution = self.compute_resolution
-        arrays = set(a for a in (self.input_arrays.values() + self.output_arrays.values()))
+        arrays = set(self.input_arrays.values()).union(self.output_arrays.values())
         for a in arrays:
             if not a.is_bound:
                 msg='FATAL ERROR: {}::{} has not been bound to any memory '
@@ -419,7 +419,7 @@ class SymbolicExpressionInfo(object):
         self.compute_resolution = tuple(compute_resolution)
 
     def check_buffers(self):
-        buffers = set(b for b in (self.input_buffers.values() + self.output_buffers.values()))
+        buffers = set(self.input_buffers.values()).union(self.output_buffers.values())
         for b in buffers:
             if not b.is_bound:
                 msg='FATAL ERROR: {}::{} has not been bound to any memory '
@@ -531,7 +531,7 @@ class SymbolicExpressionParser(object):
         kind = None
         fields = SortedDict()
         arrays = SortedDict()
-        exprs = filter(lambda e: isinstance(e, ValidExpressions), cls.extract_expressions(exprs))
+        exprs = tuple(filter(lambda e: isinstance(e, ValidExpressions), cls.extract_expressions(exprs)))
         for expr in exprs:
             check_instance(expr, ValidExpressions)
             lhs = expr.args[0]
diff --git a/hysop/operator/base/diffusion.py b/hysop/operator/base/diffusion.py
index 227a43257..817562ead 100644
--- a/hysop/operator/base/diffusion.py
+++ b/hysop/operator/base/diffusion.py
@@ -13,9 +13,17 @@ class DiffusionOperatorBase(PoissonOperatorBase):
     Common base for spectral diffusion operator.
     """
 
+    @debug
+    def __new__(cls, Fin, Fout, variables,
+                 nu, dt, name=None, pretty_name=None, **kwds):
+        return super(DiffusionOperatorBase, cls).__new__(cls,
+                Fin=Fin, Fout=Fout, variables=variables,
+                name=name, pretty_name=pretty_name,
+                input_params=None, **kwds)
+
     @debug
     def __init__(self, Fin, Fout, variables,
-                       nu, dt, 
+                       nu, dt,
                        name=None, pretty_name=None,
                        **kwds):
         """
@@ -41,7 +49,7 @@ class DiffusionOperatorBase(PoissonOperatorBase):
                 dF/dt = nu*Laplacian(F)
                 in  = Win
                 out = Wout
-                
+
             *Implicit resolution in spectral space:
                 F_hat(tn+1) = 1/(1-nu*dt*sum(Ki**2)) * F_hat(tn)
         """
@@ -50,7 +58,7 @@ class DiffusionOperatorBase(PoissonOperatorBase):
 
         input_params   = { dt.name: dt,
                            nu.name: nu }
-        
+
         default_name = 'Diffusion_{}_{}'.format(Fin.name, Fout.name)
         default_pretty_name = 'Diffusion_{}_{}'.format(Fin.pretty_name, Fout.pretty_name)
         name = first_not_None(name, default_name)
diff --git a/hysop/operator/base/memory_reordering.py b/hysop/operator/base/memory_reordering.py
index 66ec7b3dd..85813c916 100644
--- a/hysop/operator/base/memory_reordering.py
+++ b/hysop/operator/base/memory_reordering.py
@@ -98,7 +98,7 @@ class MemoryReorderingBase(object, metaclass=ABCMeta):
                         output_field=output_field,
                         input_topology_states=input_topology_states)
         assert len(input_topology_states)==1
-        istate = input_topology_states.values()[0]
+        istate = next(iter(input_topology_states.values()))
         ostate.memory_order = self.target_memory_order
         return ostate
 
diff --git a/hysop/operator/base/poisson.py b/hysop/operator/base/poisson.py
index 99c3380b1..7940f6ef6 100644
--- a/hysop/operator/base/poisson.py
+++ b/hysop/operator/base/poisson.py
@@ -13,6 +13,16 @@ class PoissonOperatorBase(SpectralOperatorBase):
     Solves the poisson equation using a specific implementation.
     """
 
+    @debug
+    def __new__(cls, Fin, Fout, variables,
+            name=None, pretty_name=None,
+            dump_energy=None, dump_input_energy=None, dump_output_energy=None,
+            plot_energy=None, plot_input_energy=None, plot_output_energy=None, plot_inout_energy=None,
+            **kwds):
+        return super(PoissonOperatorBase, cls).__new__(cls,
+                name=name, pretty_name=pretty_name,
+                input_fields=None, output_fields=None, **kwds)
+
     @debug
     def __init__(self, Fin, Fout, variables,
             name=None, pretty_name=None,
diff --git a/hysop/operator/base/poisson_curl.py b/hysop/operator/base/poisson_curl.py
index c7c6e4e2c..ee049b87a 100644
--- a/hysop/operator/base/poisson_curl.py
+++ b/hysop/operator/base/poisson_curl.py
@@ -21,6 +21,19 @@ class PoissonCurlOperatorBase(object):
     Solves the poisson-rotational equation using a specific implementation.
     """
 
+    @debug
+    def __new__(cls, vorticity, velocity, variables,
+            diffusion=None, dt=None, projection=None,
+            dump_energy=None, dump_velocity_energy=None,
+            dump_input_vorticity_energy=None, dump_output_vorticity_energy=None,
+            plot_energy=None, plot_velocity_energy=None,
+            plot_input_vorticity_energy=None, plot_output_vorticity_energy=None, plot_inout_vorticity_energy=None,
+            **kwds):
+        return super(PoissonCurlOperatorBase, cls).__new__(cls,
+                input_fields=None, output_fields=None,
+                input_params=None, output_params=None, **kwds)
+
+
     @debug
     def __init__(self, vorticity, velocity, variables,
             diffusion=None, dt=None, projection=None,
@@ -218,6 +231,14 @@ class PoissonCurlOperatorBase(object):
 
 
 class SpectralPoissonCurlOperatorBase(PoissonCurlOperatorBase, SpectralOperatorBase):
+
+    @debug
+    def __new__(cls, vorticity, velocity, variables,
+            diffusion=None, dt=None, projection=None, **kwds):
+        return super(SpectralPoissonCurlOperatorBase, cls).__new__(cls,
+                vorticity=vorticity, velocity=velocity, variables=variables,
+                diffusion=diffusion, dt=dt, projection=projection, **kwds)
+
     @debug
     def __init__(self, vorticity, velocity, variables,
             diffusion=None, dt=None, projection=None, **kwds):
@@ -256,14 +277,14 @@ class SpectralPoissonCurlOperatorBase(PoissonCurlOperatorBase, SpectralOperatorB
         kd1s = ()
         for Wi in W_Fts:
             expr1 = grad(Wi, Wi.frame)
-            kd1 = sorted(tg.push_expressions(*to_tuple(expr1)), key=lambda k: k.axis)
+            kd1 = tuple(sorted(tg.push_expressions(*to_tuple(expr1)), key=lambda k: k.axis))
             kd1s += (kd1,)
 
         # laplacian
         kd2s = ()
         for Wi in W_Fts:
             expr2 = laplacian(Wi, Wi.frame)
-            kd2 = sorted(tg.push_expressions(*to_tuple(expr2)), key=lambda k: k.axis)
+            kd2 = tuple(sorted(tg.push_expressions(*to_tuple(expr2)), key=lambda k: k.axis))
             kd2s += (kd2,)
 
         # curl
diff --git a/hysop/operator/base/solenoidal_projection.py b/hysop/operator/base/solenoidal_projection.py
index 1e6880b88..8181ea633 100644
--- a/hysop/operator/base/solenoidal_projection.py
+++ b/hysop/operator/base/solenoidal_projection.py
@@ -15,14 +15,20 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
     """
     Solves solenoidal projection (project a 3d field F such that div(F)=0)
     """
-    
+
+    @debug
+    def __new__(cls, input_field, output_field, variables,
+            input_field_div=None, output_field_div=None, **kwds):
+        return super(SolenoidalProjectionOperatorBase, cls).__new__(cls,
+                input_fields=input_field, output_fields=output_field, **kwds)
+
     @debug
-    def __init__(self, input_field, output_field, variables, 
-            input_field_div=None, output_field_div=None, **kwds): 
+    def __init__(self, input_field, output_field, variables,
+            input_field_div=None, output_field_div=None, **kwds):
         """
-        SolenoidalProjection projects a 3D vector field 
+        SolenoidalProjection projects a 3D vector field
         onto the space of divergence free fields.
-        
+
         Parameters
         ----------
         input_field:  :class:`~hysop.fields.continuous_field.Field`
@@ -35,10 +41,10 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
             Optionally compute output field divergence.
         variables: dict
             dictionary of fields as keys and topologies as values.
-        kwds : 
+        kwds :
             base class parameters.
         """
-        
+
         check_instance(output_field, Field)
         check_instance(input_field,  Field)
         check_instance(input_field_div,  Field, allow_none=True)
@@ -56,11 +62,11 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
             output_fields[input_field_div] = variables[input_field_div]
         if (output_field_div is not None):
             output_fields[output_field_div] = variables[output_field_div]
-        
+
         dim   = input_field.domain.dim
         icomp = input_field.nb_components
         ocomp = output_field.nb_components
-        
+
         if (input_field.domain != output_field.domain):
             msg = 'input_field and output_field do not share the same domain.'
             raise ValueError(msg)
@@ -81,10 +87,10 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
             msg='input_field_div component mistmach, got {} components but expected 1.'
             msg=msg.format(output_field_div.nb_components)
             raise RuntimeError(msg)
-        
-        super(SolenoidalProjectionOperatorBase, self).__init__(input_fields=input_fields, 
+
+        super(SolenoidalProjectionOperatorBase, self).__init__(input_fields=input_fields,
                 output_fields=output_fields, **kwds)
-        
+
         Fin  = input_field
         Fout = output_field
         compute_divFin  = (input_field_div is not None)
@@ -100,9 +106,9 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
         kd1s, kd2s = (), ()
         for Wi in Fts:
             expr = grad(Wi, Wi.frame)
-            kd1 = sorted(tg.push_expressions(*to_tuple(expr)), key=lambda k: k.axis)
+            kd1 = tuple(sorted(tg.push_expressions(*to_tuple(expr)), key=lambda k: k.axis))
             expr = laplacian(Wi, Wi.frame)
-            kd2 = sorted(tg.push_expressions(expr), key=lambda k: k.axis)
+            kd2 = tuple(sorted(tg.push_expressions(expr), key=lambda k: k.axis))
             kd1s += (kd1,)
             kd2s += (kd2,)
 
@@ -122,16 +128,16 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
         self.Fout    = output_field
         self.divFin  = input_field_div
         self.divFout = output_field_div
-        
+
         self.compute_divFin  = compute_divFout
         self.compute_divFout = compute_divFout
-        
+
         self.tg = tg
         self.forward_transforms  = forward_transforms
         self.backward_transforms = backward_transforms
         self.backward_divFin_transform  = backward_divFin_transform
         self.backward_divFout_transform = backward_divFout_transform
-        
+
         self.Bts   = Bts
         self.Fts   = Fts
         self.kd1s  = kd1s
@@ -144,17 +150,17 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
         super(SolenoidalProjectionOperatorBase, self).discretize()
         dFin  = self.get_input_discrete_field(self.Fin)
         dFout = self.get_output_discrete_field(self.Fout)
-        
+
         if self.compute_divFin:
             ddivFin  = self.output_discrete_fields[self.divFin]
         else:
             ddivFin = None
-        
+
         if self.compute_divFout:
             ddivFout  = self.output_discrete_fields[self.divFout]
         else:
             ddivFout = None
-        
+
         kd1s, kd2s = self.kd1s, self.kd2s
         dkd1s = ()
         for kd1 in kd1s:
@@ -163,7 +169,7 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
                 idx, dwi, _ = self.tg.discrete_wave_numbers[wi]
                 dkd1[idx] = dwi
             dkd1s += (tuple(dkd1),)
-        
+
         dkd2s = ()
         for kd2 in kd2s:
             dkd2 = [None,]*len(kd1)
@@ -171,7 +177,7 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
                 idx, dwi, _ = self.tg.discrete_wave_numbers[wi]
                 dkd2[idx] = dwi
             dkd2s += (tuple(dkd2),)
-        
+
         self.dFin     = dFin
         self.dFout    = dFout
         self.ddivFin  = ddivFin
@@ -179,17 +185,17 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
         self.dkd1s = tuple(dkd1s)
         self.dkd2s = tuple(dkd2s)
 
-    
+
     def get_work_properties(self):
         requests = super(SolenoidalProjectionOperatorBase, self).get_work_properties()
-        for (i,(Ft,Bt)) in enumerate(zip(self.forward_transforms, 
+        for (i,(Ft,Bt)) in enumerate(zip(self.forward_transforms,
                                          self.backward_transforms)):
             assert (Ft.backend == Bt.backend)
             assert (Ft.output_dtype == Bt.input_dtype), (Ft.output_dtype, Bt.input_dtype)
             assert (Ft.output_shape == Bt.input_shape), (Ft.output_shape, Bt.input_shape)
             shape = Ft.output_shape
             dtype = Ft.output_dtype
-            request = MemoryRequest(backend=self.tg.backend, dtype=dtype, 
+            request = MemoryRequest(backend=self.tg.backend, dtype=dtype,
                                     shape=shape, nb_components=1,
                                     alignment=self.min_fft_alignment)
             requests.push_mem_request('fft_buffer_{}'.format(i), request)
@@ -197,7 +203,7 @@ class SolenoidalProjectionOperatorBase(SpectralOperatorBase):
 
     def setup(self, work):
         dkd1s, dkd2s = self.dkd1s, self.dkd2s
-        
+
         output_axes = self.forward_transforms[0].output_axes
         for (i,(Ft,Bt)) in enumerate(zip(self.forward_transforms, self.backward_transforms)):
             dtmp, = work.get_buffer(self, 'fft_buffer_{}'.format(i))
diff --git a/hysop/operator/base/spectral_operator.py b/hysop/operator/base/spectral_operator.py
index e5adba837..700508771 100644
--- a/hysop/operator/base/spectral_operator.py
+++ b/hysop/operator/base/spectral_operator.py
@@ -33,13 +33,17 @@ from hysop.numerics.fft.fft import FFTI, simd_alignment, is_byte_aligned, HysopF
 
 class SpectralComputationalGraphNodeFrontend(ComputationalGraphNodeFrontend):
 
-    def __init__(self, implementation, **kwds):
-        impl, extra_kwds = self.get_actual_implementation(implementation=implementation, **kwds)
+    def __new__(cls, implementation, enforce_implementation=True, **kwds):
+        return super(SpectralComputationalGraphNodeFrontend, cls).__new__(cls,
+            implementation=implementation, **kwds)
+
+    def __init__(self, implementation, enforce_implementation=True, **kwds):
+        impl, extra_kwds = self.get_actual_implementation(implementation=implementation,
+                enforce_implementation=enforce_implementation, **kwds)
         for k in extra_kwds.keys():
             assert k not in kwds
         kwds.update(extra_kwds)
-        super(SpectralComputationalGraphNodeFrontend, self).__init__(
-            implementation=impl, **kwds)
+        super(SpectralComputationalGraphNodeFrontend, self).__init__(implementation=impl, **kwds)
 
     @classmethod
     def get_actual_implementation(cls, implementation,
@@ -81,7 +85,7 @@ class SpectralComputationalGraphNodeFrontend(ComputationalGraphNodeFrontend):
         """
         implementation = first_not_None(implementation, cls.default_implementation())
         assert implementation in cls.implementations()
-        extra_kwds = {'enable_opencl_host_buffer_mapping': False}
+        extra_kwds = {}
         if (enforce_implementation):
             return (implementation, extra_kwds)
         if (implementation == Implementation.OPENCL):
@@ -92,7 +96,6 @@ class SpectralComputationalGraphNodeFrontend(ComputationalGraphNodeFrontend):
                 raise RuntimeError(msg)
             from hysop.backend.device.opencl import cl
             if (cl_env.device.type == cl.device_type.CPU):
-                extra_kwds['enable_opencl_host_buffer_mapping'] = True
                 if (Implementation.PYTHON in cls.implementations()):
                     from hysop.backend.host.host_operator import HostOperator, OpenClMappable
                     op_cls = cls.implementations()[Implementation.PYTHON]
@@ -106,6 +109,7 @@ class SpectralComputationalGraphNodeFrontend(ComputationalGraphNodeFrontend):
                         raise TypeError(msg)
                     assert Backend.HOST in op_cls.supported_backends()
                     assert Backend.OPENCL in op_cls.supported_backends()
+                    extra_kwds['enable_opencl_host_buffer_mapping'] = True
                     return (Implementation.PYTHON, extra_kwds)
         return (implementation, extra_kwds)
 
@@ -118,8 +122,11 @@ class SpectralOperatorBase(object):
     min_fft_alignment = simd_alignment  # FFTW SIMD.
 
     @debug
-    def __init__(self, fft_interface=None, fft_interface_kwds=None,
-                 **kwds):
+    def __new__(cls, fft_interface=None, fft_interface_kwds=None, **kwds):
+        return super(SpectralOperatorBase, cls).__new__(cls, **kwds)
+
+    @debug
+    def __init__(self, fft_interface=None, fft_interface_kwds=None, **kwds):
         """
         Initialize a spectral operator base.
         kwds: dict
@@ -217,7 +224,7 @@ class SpectralOperatorBase(object):
         for tg in self.transform_groups.values():
             for (k, v) in tg.get_mem_requests(**kwds).items():
                 check_instance(k, str)  # temporary buffer name
-                check_instance(v, int)  # nbytes
+                check_instance(v, (int,np.integer))  # nbytes
                 K = (k, tg.backend)
                 if K in memory_requests:
                     memory_requests[K] = max(memory_requests[K], v)
@@ -229,7 +236,7 @@ class SpectralOperatorBase(object):
         requests = super(SpectralOperatorBase, self).get_work_properties(**kwds)
         for ((k, backend), v) in self.get_mem_requests(**kwds).items():
             check_instance(k, str)
-            check_instance(v, int)
+            check_instance(v, (int,np.integer))
             if (v > 0):
                 mrequest = MemoryRequest(backend=backend, size=v,
                                          alignment=self.min_fft_alignment)
@@ -262,6 +269,9 @@ class SpectralTransformGroup(object):
     """
     DEBUG = False
 
+    def __new__(cls, op, tag, mem_tag, **kwds):
+        return super(SpectralTransformGroup, cls).__new__(cls, **kwds)
+
     def __init__(self, op, tag, mem_tag, **kwds):
         """
         Parameters
@@ -288,6 +298,7 @@ class SpectralTransformGroup(object):
         All forward_fields and backward_fields have to live on the same domain and
         their boundary conditions should comply with given expressions.
         """
+        super(SpectralTransformGroup, self).__init__(**kwds)
         mem_tag = first_not_None(mem_tag, 'fft_pool')
         check_instance(op, SpectralOperatorBase)
         check_instance(tag, str)
@@ -339,11 +350,11 @@ class SpectralTransformGroup(object):
 
     @property
     def forward_fields(self):
-        return map(lambda x: x[0], self._forward_transforms.keys())
+        return tuple(map(lambda x: x[0], self._forward_transforms.keys()))
 
     @property
     def backward_fields(self):
-        return map(lambda x: x[0], self._backward_transforms.keys())
+        return tuple(map(lambda x: x[0], self._backward_transforms.keys()))
 
     @property
     def forward_transforms(self):
@@ -549,7 +560,7 @@ class SpectralTransformGroup(object):
         memory_requests = {}
         for fwd in self.forward_transforms.values():
             mem_requests = fwd.get_mem_requests(**kwds)
-            check_instance(mem_requests, dict, keys=str, values=int)
+            check_instance(mem_requests, dict, keys=str, values=(int,np.integer))
             for (k, v) in mem_requests.items():
                 if k in memory_requests:
                     memory_requests[k] = max(memory_requests[k], v)
@@ -557,7 +568,7 @@ class SpectralTransformGroup(object):
                     memory_requests[k] = v
         for bwd in self.backward_transforms.values():
             mem_requests = bwd.get_mem_requests(**kwds)
-            check_instance(mem_requests, dict, keys=str, values=int)
+            check_instance(mem_requests, dict, keys=str, values=(int,np.integer))
             for (k, v) in mem_requests.items():
                 if k in memory_requests:
                     memory_requests[k] = max(memory_requests[k], v)
@@ -828,7 +839,7 @@ class SpectralTransformGroup(object):
     @property
     def output_parameters(self):
         parameters = set()
-        for pt in self._forward_transforms.values() + self._backward_transforms.values():
+        for pt in tuple(self._forward_transforms.values()) + tuple(self._backward_transforms.values()):
             parameters.update(pt.output_parameters)
         return parameters
 
@@ -889,10 +900,20 @@ class PlannedSpectralTransform(object):
     """
     DEBUG = False
 
+    def __new__(cls, transform_group, tag, symbolic_transform, action,
+                 custom_input_buffer=None, custom_output_buffer=None,
+                 matching_forward_transform=None,
+                 dump_energy=None, plot_energy=None, compute_energy_frequencies=None,
+                 **kwds):
+        return super(PlannedSpectralTransform, cls).__new__(cls, **kwds)
+
     def __init__(self, transform_group, tag, symbolic_transform, action,
                  custom_input_buffer=None, custom_output_buffer=None,
                  matching_forward_transform=None,
-                 dump_energy=None, plot_energy=None, compute_energy_frequencies=None):
+                 dump_energy=None, plot_energy=None, compute_energy_frequencies=None,
+                 **kwds):
+
+        super(PlannedSpectralTransform, self).__init__(**kwds)
 
         check_instance(transform_group, SpectralTransformGroup)
         check_instance(transform_group.op, SpectralOperatorBase)
@@ -1707,10 +1728,11 @@ class PlannedSpectralTransform(object):
                     B1_tag:   nbytes,
                     TMP_tag:  tmp_nbytes}
 
-        if (self._energy_nbytes > 0):
-            requests[ENERGY_tag] = self._energy_nbytes
-        if (self._mutexes_nbytes > 0):
-            requests[MUTEXES_tag] = self._mutexes_nbytes
+        if self._do_compute_energy:
+            if (self._energy_nbytes > 0):
+                requests[ENERGY_tag] = self._energy_nbytes
+            if (self._mutexes_nbytes > 0):
+                requests[MUTEXES_tag] = self._mutexes_nbytes
 
         return requests
 
@@ -1747,7 +1769,7 @@ class PlannedSpectralTransform(object):
         except ValueError:
             TMP = None
 
-        if (self._energy_nbytes > 0):
+        if (self._energy_nbytes is not None) and (self._energy_nbytes > 0):
             ENERGY, = work.get_buffer(op, ENERGY_tag, handle=True)
             energy_buffer = ENERGY[:self._energy_nbytes].view(dtype=self.dfield.dtype)
             assert energy_buffer.size == self._max_wavenumber+1
@@ -1755,7 +1777,7 @@ class PlannedSpectralTransform(object):
             ENERGY = None
             energy_buffer = None
 
-        if (self._mutexes_nbytes > 0):
+        if (self._mutexes_nbytes is not None) and (self._mutexes_nbytes > 0):
             MUTEXES, = work.get_buffer(op, MUTEXES_tag, handle=True)
             mutexes_buffer = MUTEXES[:self._mutexes_nbytes].view(dtype=np.int32)
             assert mutexes_buffer.size == self._max_wavenumber+1
diff --git a/hysop/operator/base/stretching_dir.py b/hysop/operator/base/stretching_dir.py
index 9ad570e3a..2939654f0 100644
--- a/hysop/operator/base/stretching_dir.py
+++ b/hysop/operator/base/stretching_dir.py
@@ -38,6 +38,15 @@ class DirectionalStretchingBase(object):
         am.update(cls.__available_methods)
         return am
 
+    @debug
+    def __new__(cls, formulation, velocity, vorticity, variables, dt,
+                 C=None, A=None, name=None,
+                 implementation=None, base_kwds=None, **kwds):
+        return super(DirectionalStretchingBase, cls).__new__(cls,
+            input_fields=None, output_fields=None,
+            input_params=None, output_params=None,
+            **kwds)
+
     @debug
     def __init__(self, formulation, velocity, vorticity, variables, dt,
                  C=None, A=None, name=None,
@@ -91,8 +100,8 @@ class DirectionalStretchingBase(object):
             Should only be given for MIXED_GRAD_UW formulations.
             ValueError will be raised on other formulations.
             The linear combination coefficients A is a scalar.
-            Contained value should be a numerical coefficient, a parameter 
-            (or a generic symbolic expression) such that C*(A*grad(U).W + (1-A)*grad^T(U).W) 
+            Contained value should be a numerical coefficient, a parameter
+            (or a generic symbolic expression) such that C*(A*grad(U).W + (1-A)*grad^T(U).W)
             is directionally splittable.
             Here * is the classical elementwise multiplication.
             Default value is 0.5.
diff --git a/hysop/operator/directional/diffusion_dir.py b/hysop/operator/directional/diffusion_dir.py
index 993dc4ee4..27222c4cb 100644
--- a/hysop/operator/directional/diffusion_dir.py
+++ b/hysop/operator/directional/diffusion_dir.py
@@ -21,6 +21,19 @@ class DirectionalDiffusion(DirectionalSymbolic):
     Directional diffusion using the symbolic code generation framework.
     """
 
+    @debug
+    def __new__(cls, fields, coeffs, variables, dt,
+                laplacian_formulation=True,
+                name=None, implementation=None, base_kwds=None, **kwds):
+        base_kwds  = first_not_None(base_kwds, {})
+        return super(DirectionalDiffusion, cls).__new__(cls,
+                name=name, variables=variables, dt=dt,
+                base_kwds=base_kwds, implementation=implementation, exprs=None,
+                candidate_input_tensors=None,
+                candidate_output_tensors=None,
+                **kwds)
+
+
     @debug
     def __init__(self, fields, coeffs, variables, dt,
                 laplacian_formulation=True,
diff --git a/hysop/operator/directional/stretching_dir.py b/hysop/operator/directional/stretching_dir.py
index 8a1d84d6a..905b31bb7 100644
--- a/hysop/operator/directional/stretching_dir.py
+++ b/hysop/operator/directional/stretching_dir.py
@@ -22,6 +22,17 @@ class DirectionalStretching(DirectionalSymbolic):
     Directional stretching using the symbolic code generation framework.
     """
 
+    @debug
+    def __new__(cls, formulation, velocity, vorticity, variables, dt,
+                 C=None, A=None, name=None, implementation=None, base_kwds=None, **kwds):
+        base_kwds = first_not_None(base_kwds, {})
+        return super(DirectionalStretching, cls).__new__(cls,
+                name=name, variables=variables, dt=dt,
+                candidate_input_tensors=None,
+                candidate_output_tensors=None,
+                base_kwds=base_kwds, implementation=implementation,
+                exprs=None, **kwds)
+
     @debug
     def __init__(self, formulation, velocity, vorticity, variables, dt,
                  C=None, A=None, name=None, implementation=None, base_kwds=None, **kwds):
@@ -181,6 +192,7 @@ class DirectionalStretching(DirectionalSymbolic):
         return exprs
 
 
+
 class StaticDirectionalStretching(DirectionalOperatorFrontend):
     """
     Directional stretching using the symbolic code generation framework.
@@ -196,6 +208,19 @@ class StaticDirectionalStretching(DirectionalOperatorFrontend):
     def default_implementation(cls):
         return Implementation.PYTHON
 
+    @debug
+    def __new__(cls, formulation, velocity, vorticity, variables, dt,
+                 C=None, A=None, name=None, implementation=None, base_kwds=None, **kwds):
+        base_kwds = first_not_None(base_kwds, {})
+        return super(StaticDirectionalStretching, cls).__new__(cls,
+            name=name, variables=variables, dt=dt,
+            formulation=formulation,
+            velocity=velocity, vorticity=vorticity,
+            C=C, A=A,
+            candidate_input_tensors=None,
+            candidate_output_tensors=None,
+            base_kwds=base_kwds, implementation=implementation, **kwds)
+
     @debug
     def __init__(self, formulation, velocity, vorticity, variables, dt,
                  C=None, A=None, name=None, implementation=None, base_kwds=None, **kwds):
@@ -322,3 +347,4 @@ class StaticDirectionalStretching(DirectionalOperatorFrontend):
             candidate_input_tensors=(velocity, vorticity,),
             candidate_output_tensors=(vorticity,),
             base_kwds=base_kwds, implementation=implementation, **kwds)
+
diff --git a/hysop/operator/hdf_io.py b/hysop/operator/hdf_io.py
index 3cc7e687c..76dcab47f 100755
--- a/hysop/operator/hdf_io.py
+++ b/hysop/operator/hdf_io.py
@@ -41,6 +41,11 @@ class HDF_IO(ComputationalGraphOperator, metaclass=ABCMeta):
         """
         return Backend.all
 
+    def __new__(cls, var_names=None,
+                name_prefix='', name_postfix='',
+                force_backend=None, **kwds):
+        return super(HDF_IO, cls).__new__(cls, **kwds)
+
     def __init__(self, var_names=None,
                 name_prefix='', name_postfix='',
                 force_backend=None, **kwds):
@@ -184,7 +189,7 @@ class HDF_IO(ComputationalGraphOperator, metaclass=ABCMeta):
 
     def discretize(self):
         super(HDF_IO, self).discretize()
-        topo = self.input_fields.values()[0]
+        topo = next(iter(self.input_fields.values()))
         use_local_hdf5 = (topo.cart_size == 1)
         use_local_hdf5 |= (topo.proc_shape[0] == topo.cart_size) and (topo.cart_size <= 16) and (not self.io_params.hdf5_disable_slicing)
         # XDMF JOIN do not support more than 16 arguments
@@ -216,7 +221,7 @@ class HDF_IO(ComputationalGraphOperator, metaclass=ABCMeta):
                 local_compute_slices[field] = mesh.local_compute_slices
                 global_compute_slices[field] = mesh.global_compute_slices
             else:
-                local_compute_slices[field] = tuple(slice(0, 0) for _ in range(self.domain.dim))
+                local_compute_slices[field]  = tuple(slice(0, 0) for _ in range(self.domain.dim))
                 global_compute_slices[field] = tuple(slice(0, 0) for _ in range(self.domain.dim))
         self._local_compute_slices = local_compute_slices
         self._global_compute_slices = global_compute_slices
@@ -294,6 +299,12 @@ class HDF_Writer(HDF_IO):
 </Xdmf>
 """
 
+    def __new__(cls, variables,
+            name=None, pretty_name=None, **kwds):
+        return super(HDF_Writer, cls).__new__(cls,
+                input_fields=None, output_fields=None,
+                name=name, pretty_name=pretty_name, **kwds)
+
     def __init__(self, variables,
             name=None, pretty_name=None, **kwds):
         """
@@ -436,7 +447,7 @@ class HDF_Writer(HDF_IO):
         ds_names = self.dataset.keys()
         joinrkfiles = None
         if self.use_local_hdf5 and (self.topology.cart_size > 1):
-            joinrkfiles = range(self.topology.cart_size)
+            joinrkfiles = tuple(range(self.topology.cart_size))
         grid_attributes = XMF.prepare_grid_attributes(
             ds_names,
             resolution, origin, step, joinrkfiles=joinrkfiles)
@@ -593,6 +604,11 @@ class HDF_Reader(HDF_IO):
     Parallel reading of hdf/xdmf files to fill some fields in.
     """
 
+    def __new__(cls, variables, restart=None, name=None, **kwds):
+        return super(HDF_Reader, cls).__new__(cls,
+                input_fields=None, output_fields=None,
+                name=name, **kwds)
+
     def __init__(self, variables, restart=None, name=None, **kwds):
         """Read some fields data from hdf/xmdf files.
         Parallel readings.
diff --git a/hysop/operator/memory_reordering.py b/hysop/operator/memory_reordering.py
index b75230f04..6d9f8dbc2 100644
--- a/hysop/operator/memory_reordering.py
+++ b/hysop/operator/memory_reordering.py
@@ -21,8 +21,8 @@ class MemoryReorderingNotImplementedError(NotImplementedError):
 class MemoryReordering(ComputationalGraphNodeGenerator):
     """
     Operator generator for inplace and out of place field memory reordering.
-    
-    Available implementations are: 
+
+    Available implementations are:
         *python (numpy based memory reordering, n-dimensional)
 
     This is currently used to convert C_CONTIGUOUS fields to F_CONTIGUOUS fields.
@@ -31,10 +31,10 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
     graph node generator generates an operator per supplied field,
     possibly on different implementation backends.
 
-    See hysop.operator.base.reorder.MemoryReorderingBase for operator backend 
+    See hysop.operator.base.reorder.MemoryReorderingBase for operator backend
     implementation interface.
     """
-    
+
     @classmethod
     def implementations(cls):
         from hysop.backend.host.python.operator.memory_reordering import PythonMemoryReordering
@@ -44,25 +44,25 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
                 Implementation.OPENCL: OpenClMemoryReordering,
         }
         return _implementations
-    
+
     @classmethod
     def default_implementation(cls):
         msg='MemoryReordering has no default implementation, '
         msg+='implementation should match the discrete field topology backend.'
         raise RuntimeError(msg)
-    
+
     @debug
-    def __init__(self, fields, variables, 
+    def __init__(self, fields, variables,
                 target_memory_order,
                 output_fields=None,
-                implementation=None, 
+                implementation=None,
                 name=None,
-                base_kwds=None, 
+                base_kwds=None,
                 **kwds):
         """
         Initialize a MemoryReordering operator generator operating on CartesianTopology topologies.
         MemoryReordering is deduced from topology requirements.
-        
+
         Parameters
         ----------
         fields: Field, list or tuple of Fields
@@ -71,7 +71,7 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
         output_fields: Field, list or tuple of Fields, optional
             Output continuous fields where the results are stored.
             Reordered shapes should match the input shapes.
-            By default output_fields are the same as input_fields 
+            By default output_fields are the same as input_fields
             resulting in inplace memory reordering.
             Input and output are matched by order int list/tuple.
         variables: dict
@@ -80,10 +80,10 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
             Target memory order to achieve.
         implementation: Implementation, optional, defaults to None
             target implementation, should be contained in available_implementations().
-            If implementation is set and topology does not match backend, 
+            If implementation is set and topology does not match backend,
                 RuntimeError will be raised on _generate.
             If None, implementation will be set according to topologies backend,
-            different implementations may be choosen for different Fields if 
+            different implementations may be choosen for different Fields if
             defined on different backends.
         name: string
             prefix for generated operator names
@@ -91,9 +91,9 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
             Base class keywords arguments.
             If None, an empty dict will be passed.
         kwds:
-            Keywords arguments that will be passed towards implementation 
+            Keywords arguments that will be passed towards implementation
             memory reordering operator __init__.
-        
+
         Notes
         -----
         Out of place reordering will always be faster to process.
@@ -105,10 +105,10 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
         A MemoryReordering operator implementation should support the MemoryReorderingBase
         interface (see hysop.operator.base.memory_ordering.MemoryOrdering).
 
-        This ComputationalGraphNodeFrontend will generate a operator for each 
+        This ComputationalGraphNodeFrontend will generate a operator for each
         input and output ScalarField pair.
 
-        All implementations should raise MemoryReorderingNotImplementedError if the user supplied 
+        All implementations should raise MemoryReorderingNotImplementedError if the user supplied
         parameters leads to unimplemented or unsupported memory reordering features.
         """
         input_fields = to_tuple(fields)
@@ -119,23 +119,23 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
                                     for (ifield,ofield) in zip(input_fields, output_fields) )
         else:
             output_fields = tuple( ifield for ifield in input_fields )
-        
+
         check_instance(input_fields,  tuple, values=Field)
         check_instance(output_fields, tuple, values=Field, size=len(input_fields))
         check_instance(target_memory_order, MemoryOrdering)
         assert target_memory_order in (MemoryOrdering.C_CONTIGUOUS,
                                        MemoryOrdering.F_CONTIGUOUS)
 
-        
-        candidate_input_tensors  = filter(lambda x: x.is_tensor, input_fields)
-        candidate_output_tensors = filter(lambda x: x.is_tensor, output_fields)
-        
+
+        candidate_input_tensors  = tuple(filter(lambda x: x.is_tensor, input_fields))
+        candidate_output_tensors = tuple(filter(lambda x: x.is_tensor, output_fields))
+
         base_kwds = first_not_None(base_kwds, {})
-        super(MemoryReordering, self).__init__(name=name, 
+        super(MemoryReordering, self).__init__(name=name,
                                                candidate_input_tensors=candidate_input_tensors,
                                                candidate_output_tensors=candidate_output_tensors,
                                                **base_kwds)
-        
+
         # expand tensors
         ifields, ofields = (), ()
         for (ifield, ofield) in zip(input_fields, output_fields):
@@ -143,7 +143,7 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
             msg+='and field {} of shape {}.'
             if (ifield.is_tensor ^ ofield.is_tensor):
                 if ifield.is_tensor:
-                    msg=msg.format(ifield.short_description(), ifield.shape, 
+                    msg=msg.format(ifield.short_description(), ifield.shape,
                                    ofield.short_description(), '(1,)')
                 else:
                     msg=msg.format(ifield.short_description(), '(1,)',
@@ -155,13 +155,13 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
                 raise RuntimeError(msg)
             ifields += ifield.fields
             ofields += ofield.fields
-        
+
         input_fields  = ifields
         output_fields = ofields
 
         check_instance(input_fields,  tuple, values=ScalarField)
         check_instance(output_fields, tuple, values=ScalarField, size=len(input_fields))
-        
+
         check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
         check_instance(base_kwds, dict, keys=str)
 
@@ -204,14 +204,14 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
                 msg += 'input_field {} dimension {}.'
                 msg.format(input_field, idim, input_fields[0].name, dim)
                 raise ValueError(msg)
-        
+
         self.input_fields   = input_fields
         self.output_fields  = output_fields
         self.variables      = variables
         self.implementation = implementation
         self.target_memory_order = target_memory_order
         self.kwds           = kwds
-    
+
     def _get_op_and_check_implementation(self, src_topo, dst_topo):
         backend = getattr(src_topo,'backend', None) or getattr(dst_topo, 'backend', None)
 
@@ -259,23 +259,23 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
             raise TypeError(msg)
 
         return op_cls
-    
+
     @debug
     def _generate(self):
         nodes = []
         for (ifield, ofield) in zip(self.input_fields, self.output_fields):
             src_topo = ComputationalGraphNode.get_topo_descriptor(self.variables, ifield)
             dst_topo = ComputationalGraphNode.get_topo_descriptor(self.variables, ofield)
-            
+
             MemoryReorderingOp = self._get_op_and_check_implementation(src_topo, dst_topo)
             kwds = self.kwds.copy()
 
             variables = { ifield: src_topo }
             variables[ofield] = dst_topo
-                
+
             # instantiate operator
-            node = MemoryReorderingOp(input_field=ifield, output_field=ofield, 
-                                    variables=variables, 
+            node = MemoryReorderingOp(input_field=ifield, output_field=ofield,
+                                    variables=variables,
                                     target_memory_order=self.target_memory_order,
                                     **kwds)
             nodes.append(node)
@@ -284,4 +284,4 @@ class MemoryReordering(ComputationalGraphNodeGenerator):
             raise RuntimeError(msg)
         return nodes
 
-    
+
diff --git a/hysop/operator/tests/test_bilevel_advection.py b/hysop/operator/tests/test_bilevel_advection.py
index 516011b9d..3b2f64eec 100644
--- a/hysop/operator/tests/test_bilevel_advection.py
+++ b/hysop/operator/tests/test_bilevel_advection.py
@@ -111,9 +111,8 @@ class TestBilevelAdvectionOperator(object):
         dt.value = (0.99 * velocity_cfl) / (max(shape)-1)
 
         ref_impl = Implementation.FORTRAN
-        implementations = DirectionalAdvection.implementations().keys()
-        implementations += Advection.implementations().keys()
-        implementations = list(set(implementations))
+        implementations = set(DirectionalAdvection.implementations().keys()).union(Advection.implementations().keys())
+        implementations = list(implementations)
         assert (ref_impl in implementations)
         implementations.remove(ref_impl)
         implementations = [ref_impl] + implementations
diff --git a/hysop/operator/tests/test_custom_symbolic.py b/hysop/operator/tests/test_custom_symbolic.py
index 2bf80f1aa..935b43d91 100644
--- a/hysop/operator/tests/test_custom_symbolic.py
+++ b/hysop/operator/tests/test_custom_symbolic.py
@@ -186,11 +186,11 @@ class TestCustomSymbolic(object):
         for dtype in dtypes:
             print('  DTYPE {}'.format(dtype.__name__))
             A = Field(domain=domain, name='A', dtype=dtype,
-                    nb_components=1, register_object=False)
+                    nb_components=1)
             B = Field(domain=domain, name='B', dtype=dtype,
-                    nb_components=1, register_object=False)
+                    nb_components=1)
             C = Field(domain=domain, name='C', dtype=dtype,
-                    nb_components=3, register_object=False)
+                    nb_components=3)
 
             P0 = ScalarParameter('P0', dtype=np.float32, initial_value=1.0)
             P1 = ScalarParameter('P1', dtype=np.float64, initial_value=2.0, const=True)
@@ -390,11 +390,11 @@ class TestCustomSymbolic(object):
 
         print(' DISCRETIZATION {}'.format(discretization))
         A = Field(domain=domain, name='A', dtype=np.float32,
-                nb_components=1, register_object=False)
+                nb_components=1)
         B = Field(domain=domain, name='B', dtype=np.float32,
-                nb_components=2, register_object=False)
+                nb_components=2)
         C = Field(domain=domain, name='C', dtype=np.float64,
-                nb_components=3, register_object=False)
+                nb_components=3)
 
         P0 = ScalarParameter('P0', dtype=np.float32, initial_value=3.14)
         T0 = TensorParameter('T',  shape=(3,), dtype=np.int32, initial_value=[4,3,2])
@@ -518,11 +518,11 @@ class TestCustomSymbolic(object):
 
         print(' DISCRETIZATION {}'.format(discretization))
         A = Field(domain=domain, name='A', dtype=np.float32,
-                nb_components=1, register_object=False)
+                nb_components=1)
         B = Field(domain=domain, name='B', dtype=np.float32,
-                nb_components=2, register_object=False)
+                nb_components=2)
         C = Field(domain=domain, name='C', dtype=np.float32,
-                nb_components=3, register_object=False)
+                nb_components=3)
 
         P0 = ScalarParameter('P0', dtype=np.float32, initial_value=3.14)
         T0 = TensorParameter('T',  shape=(3,), dtype=np.int32, initial_value=[4,3,2])
diff --git a/hysop/operator/tests/test_diffusion.py b/hysop/operator/tests/test_diffusion.py
index ce25f3ab4..9dad287c7 100644
--- a/hysop/operator/tests/test_diffusion.py
+++ b/hysop/operator/tests/test_diffusion.py
@@ -52,9 +52,9 @@ class TestDiffusionOperator(object):
             nu = ScalarParameter('nu', dtype=dtype, initial_value=random.random(), const=True)
             nb_components = 5 if (dim==2) else 6
             Fin  = Field(domain=domain, name='Fin', dtype=dtype,
-                    nb_components=nb_components, register_object=False)
+                    nb_components=nb_components)
             Fout = Field(domain=domain, name='Fout', dtype=dtype,
-                   nb_components=nb_components, register_object=False)
+                   nb_components=nb_components)
             self._test_one(shape=shape, dim=dim, dtype=dtype,
                     is_inplace=is_inplace, domain=domain,
                     Fin=Fin, Fout=Fout, nu=nu)
diff --git a/hysop/operator/tests/test_directional_diffusion.py b/hysop/operator/tests/test_directional_diffusion.py
index ba4c87ef8..e6ae504d6 100644
--- a/hysop/operator/tests/test_directional_diffusion.py
+++ b/hysop/operator/tests/test_directional_diffusion.py
@@ -64,9 +64,9 @@ class TestDirectionalDiffusionOperator(object):
         domain = Box(length=(2*npw.pi,)*dim)
         for dtype in flt_types:
             Fin  = Field(domain=domain, name='Fin', dtype=dtype,
-                    nb_components=3, register_object=False)
+                    nb_components=3)
             Fout = Field(domain=domain, name='Fout', dtype=dtype,
-                   nb_components=3, register_object=False)
+                   nb_components=3)
             coeffs = npw.random.rand(Fin.nb_components, dim).astype(dtype)
             for order in orders:
                 for time_integrator in time_integrators:
@@ -241,7 +241,7 @@ class TestDirectionalDiffusionOperator(object):
         csg.configure(dtype=MPQ, dim=1)
         D2 = csg.generate_exact_stencil(derivative=2, order=order)
 
-        directions = range(dim) + range(dim)[::-1]
+        directions = tuple(range(dim)) + tuple(range(dim))[::-1]
         dt_coeffs = (0.5,)*(2*dim)
         Xin   = { 'F{}'.format(i): dfin[i] for i in range(len(dfin)) }
         views = { 'F{}'.format(i): view    for i in range(len(dfin)) }
diff --git a/hysop/operator/tests/test_directional_stretching.py b/hysop/operator/tests/test_directional_stretching.py
index dd795a428..811caabb0 100644
--- a/hysop/operator/tests/test_directional_stretching.py
+++ b/hysop/operator/tests/test_directional_stretching.py
@@ -68,11 +68,11 @@ class TestDirectionalStretchingOperator(object):
         domain = Box(length=(2*npw.pi,)*dim)
         for dtype in flt_types:
             Vin  = Field(domain=domain, name='Vin', dtype=dtype,
-                   nb_components=dim, register_object=False)
+                   nb_components=dim)
             Win  = Field(domain=domain, name='Win', dtype=dtype,
-                   nb_components=dim, register_object=False)
+                   nb_components=dim)
             Wout = Field(domain=domain, name='Wout', dtype=dtype,
-                   nb_components=dim, register_object=False)
+                   nb_components=dim)
             C = npw.random.rand(dim).astype(dtype)
             for A in As:
                 for order in orders:
@@ -265,7 +265,7 @@ class TestDirectionalStretchingOperator(object):
         Vin = tuple(d.get().handle.copy() for d in dvin.data)
         Win = tuple(d.get().handle.copy() for d in dwin.data)
 
-        directions = range(dim) + range(dim)[::-1]
+        directions = tuple(range(dim)) + tuple(range(dim))[::-1]
         dt_coeffs = (0.5,)*(2*dim)
 
         wis   = tuple('W{}'.format(i) for i in range(3))
diff --git a/hysop/operator/tests/test_poisson.py b/hysop/operator/tests/test_poisson.py
index 3de240114..151f0477f 100644
--- a/hysop/operator/tests/test_poisson.py
+++ b/hysop/operator/tests/test_poisson.py
@@ -121,9 +121,9 @@ class TestPoissonOperator(object):
                          rboundaries=rboundaries)
 
             Psi = Field(domain=domain, name='Psi', dtype=dtype,
-                        nb_components=nb_components, register_object=False)
+                        nb_components=nb_components)
             W = Field(domain=domain, name='W', dtype=dtype,
-                      nb_components=nb_components, register_object=False)
+                      nb_components=nb_components)
 
             self._test_one(shape=shape, dim=dim, dtype=dtype,
                            domain=domain, Psi=Psi, W=W,
diff --git a/hysop/operator/tests/test_solenoidal_projection.py b/hysop/operator/tests/test_solenoidal_projection.py
index 966a8cf56..79ad2fe2d 100644
--- a/hysop/operator/tests/test_solenoidal_projection.py
+++ b/hysop/operator/tests/test_solenoidal_projection.py
@@ -129,9 +129,9 @@ class TestSolenoidalProjectionOperator(object):
             U     = VelocityField(domain=domain, name='U', dtype=dtype)
             U0    = VelocityField(domain=domain, name='U0', dtype=dtype)
             U1    = VelocityField(domain=domain, name='U1', dtype=dtype)
-            divU  = U.div(name='divU', register_object=False)
-            divU0 = U0.div(name='divU0', register_object=False)
-            divU1 = U1.div(name='divU1', register_object=False)
+            divU  = U.div(name='divU')
+            divU0 = U0.div(name='divU0')
+            divU1 = U1.div(name='divU1')
 
             self._test_one(shape=shape, dtype=dtype, polynomial=polynomial,
                     domain=domain, U=U, U0=U0, U1=U1,
diff --git a/hysop/operator/transpose.py b/hysop/operator/transpose.py
index cf8228cb8..b80b8e825 100644
--- a/hysop/operator/transpose.py
+++ b/hysop/operator/transpose.py
@@ -159,8 +159,8 @@ class Transpose(ComputationalGraphNodeGenerator):
         check_instance(input_fields,  tuple, values=Field)
         check_instance(output_fields, tuple, values=Field, size=len(input_fields))
 
-        candidate_input_tensors  = filter(lambda x: x.is_tensor, input_fields)
-        candidate_output_tensors = filter(lambda x: x.is_tensor, output_fields)
+        candidate_input_tensors  = tuple(filter(lambda x: x.is_tensor, input_fields))
+        candidate_output_tensors = tuple(filter(lambda x: x.is_tensor, output_fields))
 
         base_kwds = {} if (base_kwds is None) else {}
         super(Transpose,self).__init__(name=name,
diff --git a/hysop/problem.py b/hysop/problem.py
index 236f2d9cf..b7938075a 100644
--- a/hysop/problem.py
+++ b/hysop/problem.py
@@ -87,7 +87,7 @@ class Problem(ComputationalGraph):
     def check_unique_clenv(self):
         cl_env, first_op = None, None
         for op in self.nodes:
-            for topo in set(op.input_fields.values() + op.output_fields.values()):
+            for topo in set(op.input_fields.values()).union(op.output_fields.values()):
                 if (topo.backend.kind == Backend.OPENCL):
                     if (cl_env is None):
                         first_op = op
diff --git a/hysop/symbolic/spectral.py b/hysop/symbolic/spectral.py
index 4813b28cc..8ba27630e 100644
--- a/hysop/symbolic/spectral.py
+++ b/hysop/symbolic/spectral.py
@@ -22,6 +22,9 @@ class WaveNumberIndex(sm.Symbol):
         obj._real_index = None
         return obj
 
+    def __init__(self, axis):
+        super(WaveNumberIndex, self).__init__()
+
     def bind_axes(self, axes):
         assert (self._axes is None) or (axes == self._axes)
         dim = len(axes)
@@ -108,6 +111,9 @@ class WaveNumber(Dummy):
 
         return obj
 
+    def __init__(self, axis, transform, exponent, **kwds):
+        super(WaveNumber, self).__init__(name=None, pretty_name=None, **kwds)
+
     @property
     def axis(self):
         return self._axis
diff --git a/hysop/tools/callback.py b/hysop/tools/callback.py
index 2d51b476b..f5d1f605b 100644
--- a/hysop/tools/callback.py
+++ b/hysop/tools/callback.py
@@ -271,7 +271,7 @@ class CallbackGroup(CallbackTask):
     def report(self,offset=0):
         s = ''
         s += '{}{}'.format(self.offset_str(offset), self.name)
-        tasks = sorted(self.tasks, key=lambda x:x.total(), reverse=True)
+        tasks = tuple(sorted(self.tasks, key=lambda x:x.total(), reverse=True))
         for task in tasks:
             s += '\n'+task.report(offset+1)
         s += '\n{}{:15s} {}'.format(self.offset_str(offset+1),'total:',
@@ -432,7 +432,7 @@ class CallbackProfiler(object):
             self.groups[groupname] = group
 
     def registered_targets(self):
-        return self.groups.keys() + self.tasks.keys()
+        return tuple(self.groups.keys()) + tuple(self.tasks.keys())
 
     def register_callbacks(self,target,tic_callbacks=[],tac_callbacks=[]):
         dic = self._check_registered(target)
@@ -466,12 +466,12 @@ class CallbackProfiler(object):
                         s += '\n'+task.report(1)
         elif mode =='recursive':
             if self.has_groups():
-                groups = sorted(self.groups.values(), key=lambda x:x.total(), reverse=True)
+                groups = tuple(sorted(self.groups.values(), key=lambda x:x.total(), reverse=True))
                 for group in groups:
                     s += '\n'+group.report(1)
             if self.has_tasks():
                 individual_tasknames = set(self.tasks.keys()).difference(self._tasks_in_group)
-                tasknames = sorted(individual_tasknames, key=lambda x:self.tasks[x].total(), reverse=True)
+                tasknames = tuple(sorted(individual_tasknames, key=lambda x:self.tasks[x].total(), reverse=True))
                 for taskname in tasknames:
                     task = self.tasks[taskname]
                     s += '\n'+task.report(1)
diff --git a/hysop/tools/debug_dumper.py b/hysop/tools/debug_dumper.py
index 3592cf40d..bf94fd838 100644
--- a/hysop/tools/debug_dumper.py
+++ b/hysop/tools/debug_dumper.py
@@ -60,10 +60,10 @@ class DebugDumper(object):
     def lformat(cls, id_, iteration, time, tag, min_, max_, mean, variance, dtype, shape, description='', dump_precision=None):
         try:
             return '{:<4}  {:<10}  {:<20.{p}f}  {:<40}  {:<+20.{p}f}  {:<+20.{p}f}  {:<+20.{p}f}  {:<+20.{p}f}  {:<20}  {:<20}  {}'.format(
-                id_, iteration, time, tag, min_, max_, mean, variance, dtype, shape, description, p=dump_precision)
+                id_, iteration, time, tag, min_, max_, mean, variance, str(dtype), str(shape), description, p=dump_precision)
         except:
             return '{:<4}  {:<10}  {:<20}  {:<40}  {:<20}  {:<20}  {:<20}  {:<20}  {:<20}  {:<20}  {}'.format(
-                id_, iteration, time, tag, min_, max_, mean, variance, dtype, shape, description)
+                id_, iteration, time, tag, min_, max_, mean, variance, str(dtype), str(shape), description)
 
     def print_header(self, with_datetime=False):
         now = datetime.datetime.now()
diff --git a/hysop/tools/decorators.py b/hysop/tools/decorators.py
index 6b57a9144..23a859037 100644
--- a/hysop/tools/decorators.py
+++ b/hysop/tools/decorators.py
@@ -179,7 +179,7 @@ def not_implemented(f):
     fname = f.__name__
     @wraps(f)
     def wrapper(*args, **kargs):
-        args = list(args) + kargs.values()
+        args = list(args) + list(kargs.values())
         if len(args)>0:
             arg0  = args[0]
             if isinstance(arg0,type): # @classmethod
diff --git a/hysop/tools/enum.py b/hysop/tools/enum.py
index ccf4eb5d1..e5aed8ad1 100644
--- a/hysop/tools/enum.py
+++ b/hysop/tools/enum.py
@@ -206,7 +206,7 @@ class EnumFactory(object):
         mcls = type(name+'MetaEnum', (EnumFactory.MetaEnum,), mcls_dic)
 
         class Enum(base_cls, metaclass=mcls):
-            def __init__(self, field=sorted(fields.keys())[0]):
+            def __init__(self, field=tuple(sorted(fields.keys()))[0]):
                 assert isinstance(field, str) and len(field)>0
                 self._field = field
                 self._value = self.__class__.dtype(self.__class__._fields[field])
diff --git a/hysop/tools/field_utils.py b/hysop/tools/field_utils.py
index 10090f484..20e335d20 100644
--- a/hysop/tools/field_utils.py
+++ b/hysop/tools/field_utils.py
@@ -411,7 +411,7 @@ if __name__ == '__main__':
     print
     bxvars, pxvars, vxvars, lxvars = (('x',)*5,)*4
     varpows = (1,)*5
-    xvars_components = range(5)
+    xvars_components = tuple(range(5))
     var_components=(0,4,3,2)
     _print(DifferentialStringFormatter.format_pd(bvar, pvar, vvar, lvar,
                                               bxvars, pxvars, vxvars, lxvars,
diff --git a/hysop/tools/hptt_utils.py b/hysop/tools/hptt_utils.py
index ec7af8190..2594ac704 100644
--- a/hysop/tools/hptt_utils.py
+++ b/hysop/tools/hptt_utils.py
@@ -26,11 +26,11 @@ if HAS_HPTT:
             return False
         if src.dtype not in (np.float32, np.float64, np.complex64, np.complex128):
             return False
-        if src.flags['C_CONTIGUOUS'] != dst.flags['C_CONTIGUOUS']:
+        if src.flags.c_contiguous != dst.flags.c_contiguous:
             return False
-        if src.flags['F_CONTIGUOUS'] != dst.flags['F_CONTIGUOUS']:
+        if src.flags.f_contiguous != dst.flags.f_contiguous:
             return False
-        if not (src.flags['C_CONTIGUOUS'] ^ src.flags['F_CONTIGUOUS']):
+        if not (src.flags.c_contiguous ^ src.flags.f_contiguous):
             return False
         return not array_share_data(src, dst)
 else:
diff --git a/hysop/tools/misc.py b/hysop/tools/misc.py
index 22f3b47b4..98df7fdf7 100644
--- a/hysop/tools/misc.py
+++ b/hysop/tools/misc.py
@@ -122,7 +122,7 @@ class Utils(object):
     """
     @staticmethod
     def argsort(seq):
-        return sorted(range(len(seq)), key=seq.__getitem__)
+        return tuple(sorted(range(len(seq)), key=seq.__getitem__))
 
     @staticmethod
     def upper_pow2(x):
diff --git a/hysop/tools/mpi_utils.py b/hysop/tools/mpi_utils.py
index 0fab45c3d..3a211ba67 100644
--- a/hysop/tools/mpi_utils.py
+++ b/hysop/tools/mpi_utils.py
@@ -50,9 +50,9 @@ def dtype_to_mpi_type(dtype):
 
 def order_to_mpi_order(order):
     from hysop.constants import MemoryOrdering
-    if (order in 'cC') or (order==MemoryOrdering.C_CONTIGUOUS) or (order==MPI.ORDER_C): 
+    if (order in 'cC') or (order==MemoryOrdering.C_CONTIGUOUS) or (order==MPI.ORDER_C):
         return MPI.ORDER_C
-    elif (order in 'fF') or (order==MemoryOrdering.F_CONTIGUOUS) or (order==MPI.ORDER_F): 
+    elif (order in 'fF') or (order==MemoryOrdering.F_CONTIGUOUS) or (order==MPI.ORDER_F):
         return MPI.ORDER_F
     else:
         msg='Unknown value of type {}.'.format(type(order))
@@ -65,8 +65,8 @@ def get_mpi_order(data):
         is_f_contiguous = data.is_fortran_contiguous
     else:
         # assume numpy like interface
-        is_c_contiguous = data.flags['C_CONTIGUOUS']
-        is_f_contiguous = data.flags['F_CONTIGUOUS']
+        is_c_contiguous = data.flags.c_contiguous
+        is_f_contiguous = data.flags.f_contiguous
     if is_c_contiguous:
         return MPI.ORDER_C
     elif is_f_contiguous:
diff --git a/hysop/tools/numba_utils.py b/hysop/tools/numba_utils.py
index 7c5bc6bbb..62ad05cb7 100644
--- a/hysop/tools/numba_utils.py
+++ b/hysop/tools/numba_utils.py
@@ -76,9 +76,9 @@ def make_numba_signature(*args, **kwds):
         elif isinstance(a, np.ndarray):
             assert a.dtype.type in dtype_to_ntype, a.dtype.type
             dtype = dtype_to_ntype[a.dtype.type]
-            if a.flags['C_CONTIGUOUS']:
+            if a.flags.c_contiguous:
                 na = nb.types.Array(dtype=dtype, ndim=a.ndim, layout='C')
-            elif a.flags['F_CONTIGUOUS']:
+            elif a.flags.f_contiguous:
                 na = nb.types.Array(dtype=dtype, ndim=a.ndim, layout='F')
             else:
                 na = nb.types.Array(dtype=dtype, ndim=a.ndim, layout='A')
diff --git a/hysop/tools/profiler.py b/hysop/tools/profiler.py
index dd2330c0f..57d9e49be 100644
--- a/hysop/tools/profiler.py
+++ b/hysop/tools/profiler.py
@@ -123,8 +123,8 @@ class Profiler(object):
         summary = self.summary
         if len(summary) > 0:
             if (self._l>1) and (len(summary) == 1) and \
-                    isinstance(summary.values()[0], FProfiler):
-                s = '>{}::{}'.format(self.get_name(), summary.values()[0])
+                    isinstance(next(iter(summary.values())), FProfiler):
+                s = '>{}::{}'.format(self.get_name(), next(iter(summary.values())))
             else:
                 s = '{}[{}]>{}{}'.format(
                     '\n' if (self._l==1) else '',
diff --git a/hysop/tools/spectral_utils.py b/hysop/tools/spectral_utils.py
index 045b24eec..c7c3dbfec 100644
--- a/hysop/tools/spectral_utils.py
+++ b/hysop/tools/spectral_utils.py
@@ -823,7 +823,7 @@ class EnergyPlotter(object):
 
             if len(energy_parameters)==1:
                 default_fig_title = 'Energy parameter {}'.format(
-                    energy_parameters.values()[0].pretty_name)
+                    next(iter(energy_parameters.values())).pretty_name)
             else:
                 default_fig_title = 'Energy parameters {}'.format(
                         ' | '.join(p.pretty_name for p in energy_parameters.values()))
diff --git a/hysop/tools/transposition_states.py b/hysop/tools/transposition_states.py
index 79e8c5ab5..a3d9ed47e 100644
--- a/hysop/tools/transposition_states.py
+++ b/hysop/tools/transposition_states.py
@@ -55,7 +55,7 @@ class TranspositionStateType(type):
             return it.permutations(range(dim), dim)
         def __filter_axes(cls, predicate):
             """Return a filtered iterator on all possible permutations."""
-            return filter(predicate, cls.all_axes())
+            return tuple(filter(predicate, cls.all_axes()))
 
         cls_methods['dimension']        = classmethod(__dimension)
         cls_methods['default_axes']     = classmethod(__default_axes)
diff --git a/hysop/topology/cartesian_descriptor.py b/hysop/topology/cartesian_descriptor.py
index 7e7ee3482..747109168 100644
--- a/hysop/topology/cartesian_descriptor.py
+++ b/hysop/topology/cartesian_descriptor.py
@@ -178,8 +178,7 @@ class CartesianTopologyDescriptor(TopologyDescriptor):
         If None is returned, create_topology will be called instead.
         """
         if known_topologies:
-            ordered_topologies = sorted(known_topologies,
-                                            key=lambda topo: sum(topo.ghosts))
+            ordered_topologies = tuple(sorted(known_topologies, key=lambda topo: sum(topo.ghosts)))
             return ordered_topologies[0]
         else:
             return None
diff --git a/hysop_examples/example_utils.py b/hysop_examples/example_utils.py
index cf67e6a98..41dc24d53 100644
--- a/hysop_examples/example_utils.py
+++ b/hysop_examples/example_utils.py
@@ -195,6 +195,7 @@ class HysopArgParser(argparse.ArgumentParser):
             dump_dirs = map(lambda ddir: os.path.abspath(ddir), dump_dirs)
             dump_dirs = filter(lambda ddir: os.path.isdir(ddir) and \
                     (ddir not in ('/','/home','~',os.path.expanduser('~'))), dump_dirs)
+            dump_dirs = tuple(dump_dirs)
             if args.no_interactive:
                 confirm_deletion = True
             else:
@@ -1298,7 +1299,7 @@ class HysopArgParser(argparse.ArgumentParser):
         return graphical_io
 
     def _check_graphical_io_args(self, args):
-        if (args.visu_rank<0):
+        if (args.visu_rank is not None) and (args.visu_rank<0):
             args.visu_rank = None
         self._check_default(args, 'visu_rank', int, allow_none=True)
         self._check_positive(args, 'visu_rank', strict=False, allow_none=True)
diff --git a/hysop_examples/examples/analytic/analytic.py b/hysop_examples/examples/analytic/analytic.py
index 4c27e7253..011865c2a 100755
--- a/hysop_examples/examples/analytic/analytic.py
+++ b/hysop_examples/examples/analytic/analytic.py
@@ -1,8 +1,6 @@
-#!/usr/bin/env python2
 import numpy as np
 import sympy as sm
 
-
 def compute(args):
     '''
     HySoP Analytic Example: Initialize a field with a space and time dependent analytic formula.
diff --git a/hysop_examples/examples/particles_above_salt/particles_above_salt_bc.py b/hysop_examples/examples/particles_above_salt/particles_above_salt_bc.py
index 1054d6ace..d0b941991 100644
--- a/hysop_examples/examples/particles_above_salt/particles_above_salt_bc.py
+++ b/hysop_examples/examples/particles_above_salt/particles_above_salt_bc.py
@@ -1,6 +1,6 @@
 ## See Meiburg 2012 & 2014
 ## Sediment-laden fresh water above salt water.
-    
+
 import numpy as np
 import scipy as sp
 import sympy as sm
@@ -21,7 +21,7 @@ def delta(Ys, l0):
     for Yi in Ys:
         Y0 = Y0*Yi
     return 0.1*l0*(np.random.rand(*Y0.shape)-0.5)
-    
+
 def init_concentration(data, coords, l0, component):
     assert (component==0)
     X = coords[0]
@@ -55,7 +55,7 @@ def compute(args):
 
     from hysop.methods import SpaceDiscretization, Remesh, TimeIntegrator, \
                               ComputeGranularity, Interpolation
-    
+
     from hysop.symbolic import sm, space_symbols, local_indices_symbols
     from hysop.symbolic.base import SymbolicTensor
     from hysop.symbolic.field import curl
@@ -82,7 +82,7 @@ def compute(args):
     npts=N[::-1]
     Xo=Xo[::-1]
     Xn=Xn[::-1]
-    
+
     lboundaries = (BoxBoundaryCondition.PERIODIC,)*(dim-1)+(BoxBoundaryCondition.SYMMETRIC,)
     rboundaries = (BoxBoundaryCondition.PERIODIC,)*(dim-1)+(BoxBoundaryCondition.SYMMETRIC,)
 
@@ -93,7 +93,7 @@ def compute(args):
 
     box = Box(origin=Xo, length=np.subtract(Xn,Xo),
                 lboundaries=lboundaries, rboundaries=rboundaries)
-    
+
     # Get default MPI Parameters from domain (even for serial jobs)
     mpi_params = MPIParams(comm=box.task_comm,
                            task_id=box.current_task())
@@ -107,18 +107,18 @@ def compute(args):
     elif (impl is Implementation.OPENCL):
         # For the OpenCL implementation we need to setup the compute device
         # and configure how the code is generated and compiled at runtime.
-                
+
         # Create an explicit OpenCL context from user parameters
         from hysop.backend.device.opencl.opencl_tools import get_or_create_opencl_env, get_device_number
         cl_env = get_or_create_opencl_env(
             mpi_params=mpi_params,
-            platform_id=args.cl_platform_id, 
+            platform_id=args.cl_platform_id,
             device_id=box.machine_rank%get_device_number() if args.cl_device_id is None else None)
-        
+
         # Configure OpenCL kernel generation and tuning (already done by HysopArgParser)
         from hysop.methods import OpenClKernelConfig
         method = { OpenClKernelConfig: args.opencl_kernel_config }
-        
+
         # Setup opencl specific extra operator keyword arguments
         extra_op_kwds['cl_env'] = cl_env
     else:
@@ -131,7 +131,7 @@ def compute(args):
     vorti = VorticityField(velocity=velo)
     C = Field(domain=box, name='C', dtype=args.dtype, lboundaries=C_lboundaries, rboundaries=C_rboundaries)
     S = Field(domain=box, name='S', dtype=args.dtype, lboundaries=S_lboundaries, rboundaries=S_rboundaries)
-    
+
     # Symbolic fields
     frame = velo.domain.frame
     Us    = velo.s(*frame.vars)
@@ -141,43 +141,43 @@ def compute(args):
     dts   = dt.s
 
     ### Build the directional operators
-    #> Directional advection 
+    #> Directional advection
     advec = DirectionalAdvection(implementation=impl,
             name='advec',
-            velocity = velo,       
+            velocity = velo,
             advected_fields = (vorti,S),
             velocity_cfl = args.cfl,
             variables = {velo: npts, vorti: npts, S: npts},
             dt=dt, **extra_op_kwds)
-   
+
     V0 = [0]*dim
     VP = [0]*dim
     VP[0] = Vp
     advec_C = DirectionalAdvection(implementation=impl,
             name='advec_C',
-            velocity = velo,       
+            velocity = velo,
             advected_fields = (C,),
             relative_velocity = VP,
             velocity_cfl = args.cfl,
             variables = {velo: npts, C: npts},
             dt=dt, **extra_op_kwds)
-    
+
     #> Stretch vorticity
     if (dim==3):
         stretch = DirectionalStretching(implementation=impl,
                  name='stretch',
                  pretty_name='stretch',
                  formulation = args.stretching_formulation,
-                 velocity  = velo,       
+                 velocity  = velo,
                  vorticity = vorti,
                  variables = {velo: npts, vorti: npts},
                  dt=dt, **extra_op_kwds)
     elif (dim==2):
-        stretch = None 
+        stretch = None
     else:
         msg='Unsupported dimension {}.'.format(dim)
         raise RuntimeError(msg)
-    
+
     #> Diffusion of vorticity, S and C
     diffuse_S = Diffusion(implementation=impl,
              enforce_implementation=enforce_implementation,
@@ -186,7 +186,7 @@ def compute(args):
              nu = nu_S,
              Fin = S,
              variables = {S: npts},
-             dt=dt, 
+             dt=dt,
              **extra_op_kwds)
     diffuse_C = Diffusion(implementation=impl,
              enforce_implementation=enforce_implementation,
@@ -204,27 +204,27 @@ def compute(args):
     lhs = Ws.diff(frame.time)
     rhs = curl(Fext, frame)
     exprs = Assignment.assign(lhs, rhs)
-    external_force = DirectionalSymbolic(name='Fext', 
+    external_force = DirectionalSymbolic(name='Fext',
                                     implementation=impl,
                                     exprs=exprs, dt=dt,
                                     variables={vorti: npts,
                                                S: npts,
                                                C: npts},
                                     **extra_op_kwds)
-    
-    splitting = StrangSplitting(splitting_dim=dim, 
+
+    splitting = StrangSplitting(splitting_dim=dim,
                     order=args.strang_order)
     splitting.push_operators(advec, advec_C, stretch, external_force)
-    
+
     ### Build standard operators
     #> Poisson operator to recover the velocity from the vorticity
-    poisson = PoissonCurl(name='poisson', velocity=velo, vorticity=vorti, 
+    poisson = PoissonCurl(name='poisson', velocity=velo, vorticity=vorti,
                             variables={velo:npts, vorti: npts},
                             diffusion=nu_W, dt=dt,
                             implementation=impl,
                             enforce_implementation=enforce_implementation,
                             **extra_op_kwds)
-    
+
     #> Operator to compute the infinite norm of the velocity
     min_max_U = MinMaxFieldStatistics(name='min_max_U', field=velo,
             Finf=True, implementation=impl, variables={velo:npts},
@@ -233,14 +233,14 @@ def compute(args):
     min_max_W = MinMaxFieldStatistics(field=vorti,
             Finf=True, implementation=impl, variables={vorti:npts},
             **extra_op_kwds)
-    
+
     #> Operators to dump all fields
     dump_fields = HDF_Writer(name='dump',
                              io_params=args.io_params.clone(filename='fields'),
                              force_backend=Backend.OPENCL,
-                             variables={velo: npts, 
+                             variables={velo: npts,
                                         vorti: npts,
-                                        C: npts, 
+                                        C: npts,
                                         S: npts},
                              **extra_op_kwds)
 
@@ -250,7 +250,7 @@ def compute(args):
     view[-1] = (-200.0,+200.0)
     view = tuple(view)
     io_params = IOParams(filename='horizontally_averaged_profiles', frequency=0)
-    compute_mean_fields = ComputeMeanField(name='mean', 
+    compute_mean_fields = ComputeMeanField(name='mean',
             fields={C: (view, axes), S: (view, axes)},
             variables={C: npts, S: npts},
             io_params=io_params)
@@ -258,23 +258,23 @@ def compute(args):
     ### Adaptive timestep operator
     adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True,
                                     name='merge_dt', pretty_name='dt', )
-    dt_cfl = adapt_dt.push_cfl_criteria(cfl=args.cfl, 
+    dt_cfl = adapt_dt.push_cfl_criteria(cfl=args.cfl,
                                         Fmin=min_max_U.Fmin,
                                         Fmax=min_max_U.Fmax,
-                                        equivalent_CFL=True, 
+                                        equivalent_CFL=True,
                                         relative_velocities=[V0, VP],
                                         name='dt_cfl', pretty_name='CFL')
     dt_advec = adapt_dt.push_advection_criteria(lcfl=args.lcfl, Finf=min_max_W.Finf,
                                                  criteria=AdvectionCriteria.W_INF,
                                                  name='dt_lcfl', pretty_name='LCFL')
 
-    
-    ## Create the problem we want to solve and insert our 
+
+    ## Create the problem we want to solve and insert our
     # directional splitting subgraph and the standard operators.
     # The method dictionnary passed to this graph will be dispatched
     # accross all operators contained in the graph.
     method.update(
-            { 
+            {
                ComputeGranularity:    args.compute_granularity,
                SpaceDiscretization:   args.fd_order,
                TimeIntegrator:        args.time_integrator,
@@ -283,42 +283,42 @@ def compute(args):
         )
 
     problem = Problem(method=method)
-    problem.insert(poisson, 
+    problem.insert(poisson,
                    diffuse_S, diffuse_C,
-                   splitting, 
+                   splitting,
                    dump_fields,
                    compute_mean_fields,
-                   min_max_U, min_max_W, 
+                   min_max_U, min_max_W,
                    adapt_dt)
     problem.build(args)
-    
+
     # If a visu_rank was provided, and show_graph was set,
     # display the graph on the given process rank.
     if args.display_graph:
         problem.display(args.visu_rank)
-    
+
     # Create a simulation
     # (do not forget to specify the t and dt parameters here)
-    simu = Simulation(start=args.tstart, end=args.tend, 
+    simu = Simulation(start=args.tstart, end=args.tend,
                       nb_iter=args.nb_iter,
                       max_iter=args.max_iter,
                       dt0=args.dt, times_of_interest=args.times_of_interest,
                       t=t, dt=dt)
-    simu.write_parameters(t, dt_cfl, dt_advec, dt, 
+    simu.write_parameters(t, dt_cfl, dt_advec, dt,
             min_max_U.Finf, min_max_W.Finf, adapt_dt.equivalent_CFL,
             filename='parameters.txt', precision=8)
-    
+
     # Initialize vorticity, velocity, S and C on all topologies
     problem.initialize_field(field=velo,  formula=init_velocity)
     problem.initialize_field(field=vorti, formula=init_vorticity)
     problem.initialize_field(field=C,     formula=init_concentration, l0=l0)
     problem.initialize_field(field=S,     formula=init_salinity, l0=l0)
-    
+
     # Finally solve the problem
-    problem.solve(simu, dry_run=args.dry_run, 
+    problem.solve(simu, dry_run=args.dry_run,
             debug_dumper=args.debug_dumper,
             checkpoint_handler=args.checkpoint_handler)
-    
+
     # Finalize
     problem.finalize()
 
@@ -329,18 +329,18 @@ if __name__=='__main__':
     class ParticleAboveSaltArgParser(HysopArgParser):
         def __init__(self):
             prog_name = 'particle_above_salt_bc'
-            default_dump_dir = '{}/hysop_examples/{}'.format(HysopArgParser.tmp_dir(), 
+            default_dump_dir = '{}/hysop_examples/{}'.format(HysopArgParser.tmp_dir(),
                     prog_name)
 
             description=colors.color('HySoP Particles Above Salt Example: ', fg='blue',
                     style='bold')
             description+=colors.color('[Meiburg 2014]', fg='yellow', style='bold')
-            description+=colors.color('\nSediment-laden fresh water above salt water.', 
+            description+=colors.color('\nSediment-laden fresh water above salt water.',
                     fg='yellow')
             description+='\n'
             description+='\nThis example focuses on a validation study for the '
             description+='hybrid particle-mesh vortex method in the Boussinesq approximation.'
-    
+
             super(ParticleAboveSaltArgParser, self).__init__(
                  prog_name=prog_name,
                  description=description,
@@ -357,10 +357,10 @@ if __name__=='__main__':
     parser = ParticleAboveSaltArgParser()
 
     parser.set_defaults(impl='cl', ndim=2, npts=(64,),
-                        box_origin=(0.0,), box_length=(1.0,), 
-                        tstart=0.0, tend=500.0, 
+                        box_origin=(0.0,), box_length=(1.0,),
+                        tstart=0.0, tend=500.0,
                         dt=1e-6, cfl=4.00, lcfl=0.95,
-                        dump_times=tuple(float(x) for x in range(0,500,10)),
+                        dump_times=tuple(float(x) for x in range(0, 500, 10)),
                         dump_freq=0)
 
     parser.run(compute)
diff --git a/hysop_examples/examples/particles_above_salt/particles_above_salt_periodic.py b/hysop_examples/examples/particles_above_salt/particles_above_salt_periodic.py
index 1e9905192..7965443c5 100644
--- a/hysop_examples/examples/particles_above_salt/particles_above_salt_periodic.py
+++ b/hysop_examples/examples/particles_above_salt/particles_above_salt_periodic.py
@@ -372,7 +372,7 @@ if __name__=='__main__':
                         box_origin=(0.0,), box_length=(1.0,),
                         tstart=0.0, tend=500.0,
                         dt=1e-6, cfl=0.5, lcfl=0.125,
-                        dump_times=tuple(float(x) for x in range(0,500,5)),
+                        dump_times=tuple(float(x) for x in range(0, 500, 5)),
                         dump_freq=0)
 
     parser.run(compute)
diff --git a/hysop_examples/examples/particles_above_salt/particles_above_salt_symmetrized.py b/hysop_examples/examples/particles_above_salt/particles_above_salt_symmetrized.py
index cf3aa0b13..36e8009be 100644
--- a/hysop_examples/examples/particles_above_salt/particles_above_salt_symmetrized.py
+++ b/hysop_examples/examples/particles_above_salt/particles_above_salt_symmetrized.py
@@ -359,7 +359,7 @@ if __name__=='__main__':
                         box_origin=(0.0,), box_length=(1.0,),
                         tstart=0.0, tend=500.0,
                         dt=1e-6, cfl=0.5, lcfl=0.125,
-                        dump_times=tuple(float(x) for x in range(0,500,5)),
+                        dump_times=tuple(float(x) for x in range(0, 500, 5)),
                         dump_freq=0)
 
     parser.run(compute)
diff --git a/setup.py.in b/setup.py.in
index da7b1c75a..c00757c87 100644
--- a/setup.py.in
+++ b/setup.py.in
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """setup.py file for @PYPACKAGE_NAME@
 
-- 
GitLab