From ec4f05d490c28272a4d7f63011396fec25f2edb5 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Sun, 1 Nov 2020 02:20:24 +0100
Subject: [PATCH] fix many operators

---
 ci/scripts/test.sh                            | 14 ++---
 .../backend/device/codegen/base/variables.py  |  2 +-
 .../device/codegen/functions/runge_kutta.py   |  2 +-
 .../device/codegen/kernels/custom_symbolic.py |  2 +-
 hysop/backend/device/kernel_autotuner.py      |  4 +-
 .../device/opencl/opencl_kernel_autotuner.py  |  3 +-
 .../backend/device/opencl/opencl_operator.py  | 10 ++--
 .../host/python/operator/derivative.py        |  5 ++
 .../operator/directional/advection_dir.py     | 10 ++--
 .../python/operator/flowrate_correction.py    |  8 +++
 .../host/python/operator/penalization.py      |  8 +++
 .../host/python/operator/spatial_filtering.py |  3 ++
 .../backend/host/python/operator/transpose.py |  4 ++
 .../python/operator/vorticity_absorption.py   |  9 ++++
 hysop/core/graph/computational_node.py        |  2 +-
 .../core/graph/computational_node_frontend.py | 15 ++++++
 hysop/core/graph/continuous.py                | 16 +++++-
 hysop/core/graph/graph_builder.py             |  8 +--
 hysop/core/graph/node_generator.py            | 10 ++--
 hysop/fields/cartesian_discrete_field.py      |  5 +-
 hysop/fields/ghost_exchangers.py              |  2 +-
 hysop/numerics/remesh/kernel_generator.py     | 53 ++++++++++---------
 .../splitting/directional_splitting.py        |  5 ++
 hysop/numerics/splitting/strang.py            | 10 +++-
 hysop/numerics/stencil/stencil.py             |  2 +-
 hysop/numerics/stencil/stencil_generator.py   |  2 +-
 hysop/operator/base/advection_dir.py          | 17 ++++--
 .../operator/base/custom_symbolic_operator.py |  7 ++-
 hysop/operator/base/derivative.py             | 18 ++++++-
 hysop/operator/base/spatial_filtering.py      | 15 ++++++
 hysop/operator/base/transpose_operator.py     | 15 ++++--
 hysop/operator/derivative.py                  | 23 +++++++-
 hysop/operator/directional/advection_dir.py   | 26 ++++++---
 hysop/operator/directional/directional.py     |  6 ++-
 hysop/operator/gradient.py                    | 29 ++++++++++
 hysop/operator/misc.py                        | 12 ++++-
 hysop/operator/spatial_filtering.py           | 23 ++++++++
 hysop/operator/tests/test_absorption.py       |  4 +-
 .../tests/test_directional_advection.py       | 10 ++--
 hysop/operator/tests/test_fd_derivative.py    |  8 +--
 hysop/operator/tests/test_penalization.py     |  6 +--
 .../operator/tests/test_restriction_filter.py |  2 +-
 hysop/operator/tests/test_transpose.py        |  4 +-
 .../tests/test_velocity_correction.py         |  4 +-
 hysop/operator/transpose.py                   | 15 +++++-
 hysop/symbolic/parameter.py                   |  4 +-
 46 files changed, 357 insertions(+), 105 deletions(-)

diff --git a/ci/scripts/test.sh b/ci/scripts/test.sh
index a2a1f0382..1d67f1da2 100755
--- a/ci/scripts/test.sh
+++ b/ci/scripts/test.sh
@@ -119,13 +119,13 @@ if [ "$RUN_TESTS" = true ]; then
     #hysop_test "core/graph/tests/test_graph.py"
     #hysop_test "fields/tests/test_fields.py"
     #hysop_test "numerics/tests/test_fft.py"
-    hysop_test "operator/tests/test_analytic.py"
-    hysop_test "operator/tests/test_transpose.py"
-    hysop_test "operator/tests/test_fd_derivative.py"
-    hysop_test "operator/tests/test_absorption.py"
-    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_analytic.py"
+    #hysop_test "operator/tests/test_transpose.py"
+    #hysop_test "operator/tests/test_fd_derivative.py"
+    #hysop_test "operator/tests/test_absorption.py"
+    #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"
diff --git a/hysop/backend/device/codegen/base/variables.py b/hysop/backend/device/codegen/base/variables.py
index 55489215e..910fc0cd3 100644
--- a/hysop/backend/device/codegen/base/variables.py
+++ b/hysop/backend/device/codegen/base/variables.py
@@ -853,7 +853,7 @@ class CodegenVectorClBuiltin(CodegenVector):
                 value = [self.svalue[i] for i in key]
                 return '({})({})'.format(ctype, ','.join(value))
             return access
-        elif isinstance(key, int) :
+        elif isinstance(key, (int, np.integer)) :
             if key<0:
                 key += dim
             if key<0 or key>=dim:
diff --git a/hysop/backend/device/codegen/functions/runge_kutta.py b/hysop/backend/device/codegen/functions/runge_kutta.py
index b581e6ebb..325be7b1f 100644
--- a/hysop/backend/device/codegen/functions/runge_kutta.py
+++ b/hysop/backend/device/codegen/functions/runge_kutta.py
@@ -65,7 +65,7 @@ class RungeKuttaFunction(OpenClFunctionCodeGenerator):
                 args[arg] = rhs.args[arg]
                 rhs_args[arg] = args[arg]
 
-        rhs_name = hashlib.md5(rhs.fname).hexdigest()[0:8]
+        rhs_name = hashlib.md5(rhs.fname.encode('utf-8')).hexdigest()[0:8]
         ctype    = args[y].ctype
         basename = 'apply_{}_{}_{}'.format(method.name(),ctype,rhs_name)
 
diff --git a/hysop/backend/device/codegen/kernels/custom_symbolic.py b/hysop/backend/device/codegen/kernels/custom_symbolic.py
index a7e5c9a7f..30af9a64f 100644
--- a/hysop/backend/device/codegen/kernels/custom_symbolic.py
+++ b/hysop/backend/device/codegen/kernels/custom_symbolic.py
@@ -1281,7 +1281,7 @@ class CustomSymbolicKernelGenerator(KernelCodeGenerator, metaclass=ABCMeta):
         for fcall in self.fcalls:
             fcall.fn_kwds['offset'] = self.csc.local_offset
             if 'dx' in fcall.fn_kwds:
-                mesh_info_0 = self.field_mesh_infos.values()[0]
+                mesh_info_0 = next(iter(self.field_mesh_infos.values()))
                 dx = mesh_info_0['dx'][0]
                 fcall.fn_kwds['dx'] = dx
 
diff --git a/hysop/backend/device/kernel_autotuner.py b/hysop/backend/device/kernel_autotuner.py
index f1ad2f8e0..41a1c12ea 100644
--- a/hysop/backend/device/kernel_autotuner.py
+++ b/hysop/backend/device/kernel_autotuner.py
@@ -196,7 +196,7 @@ class KernelAutotuner(object, metaclass=ABCMeta):
                         tuning_mode=False, dry_run=False)
 
         hasher = self._hash_func()
-        hasher.update(kernel_src)
+        hasher.update(kernel_src.encode('utf-8'))
         src_hash = hasher.hexdigest()
 
         if (kernel_name != cached_kernel_name):
@@ -529,7 +529,7 @@ class KernelAutotuner(object, metaclass=ABCMeta):
                         tuning_mode=False, dry_run=False)
 
         hasher = self._hash_func()
-        hasher.update(kernel_src)
+        hasher.update(kernel_src.encode('utf-8'))
         src_hash = hasher.hexdigest()
 
         (prg, kernel) = self.build_from_source(kernel_name=kernel_name,
diff --git a/hysop/backend/device/opencl/opencl_kernel_autotuner.py b/hysop/backend/device/opencl/opencl_kernel_autotuner.py
index 4a948228e..fe974165b 100644
--- a/hysop/backend/device/opencl/opencl_kernel_autotuner.py
+++ b/hysop/backend/device/opencl/opencl_kernel_autotuner.py
@@ -81,7 +81,8 @@ class OpenClKernelAutotuner(KernelAutotuner):
                     bytes2str(kernel_local_mem_size))
             raise RuntimeError(msg)
 
-        if cl_env.cl_version >= (1,2) and (device.type == cl.device_type.CUSTOM):
+        cl_version = tuple(map(int, cl_env.cl_version))
+        if cl_version >= (1,2) and (device.type == cl.device_type.CUSTOM):
             max_kernel_global_size = kernel.get_work_group_info(kwgi.GLOBAL_WORK_SIZE, device)
             if np.any(np.greater(global_work_size, max_kernel_global_size)):
                 msg='Global size {} exceeded for kernel {} which allows only {} for '
diff --git a/hysop/backend/device/opencl/opencl_operator.py b/hysop/backend/device/opencl/opencl_operator.py
index 773b6c879..4c942f37f 100644
--- a/hysop/backend/device/opencl/opencl_operator.py
+++ b/hysop/backend/device/opencl/opencl_operator.py
@@ -76,7 +76,8 @@ class OpenClOperator(ComputationalGraphOperator, metaclass=ABCMeta):
 
         msg='mpi_params was {} and cl_env was {}.'
 
-        for topo in kwds.get('input_fields',{}).values()+kwds.get('output_fields',{}).values():
+        for topo in set(kwds.get('input_fields', []).values()).union(
+                        kwds.get('output_fields',[]).values()):
             if isinstance(topo, Topology):
                 if (cl_env is None):
                     cl_env = topo.backend.cl_env
@@ -269,12 +270,15 @@ class OpenClOperator(ComputationalGraphOperator, metaclass=ABCMeta):
         Check if all topologies are on OpenCL backend and check that all opencl environments
         match.
         """
-        topo = (self.input_fields.values()+self.output_fields.values())[0]
+        topologies = set(self.input_fields.values()).union(self.output_fields.values())
+        if not topologies:
+            return
+        topo = next(iter(topologies))
         assert isinstance(topo, TopologyView)
         assert topo.backend.kind == Backend.OPENCL
         ref_env = self.cl_env
 
-        for topo in set(self.input_fields.values()+self.output_fields.values()):
+        for topo in topologies:
             assert isinstance(topo, TopologyView)
             assert topo.backend.kind == Backend.OPENCL
             assert topo.backend.cl_env == ref_env
diff --git a/hysop/backend/host/python/operator/derivative.py b/hysop/backend/host/python/operator/derivative.py
index e4a16f48b..2625c15a5 100644
--- a/hysop/backend/host/python/operator/derivative.py
+++ b/hysop/backend/host/python/operator/derivative.py
@@ -56,6 +56,10 @@ class PythonFiniteDifferencesSpaceDerivative(FiniteDifferencesSpaceDerivativeBas
     using explicit finite differences.
     """
 
+    @debug
+    def __new__(cls, **kwds):
+        return super(PythonFiniteDifferencesSpaceDerivative, cls).__new__(cls, **kwds)
+
     @debug
     def __init__(self, **kwds):
         """
@@ -154,3 +158,4 @@ class PythonFiniteDifferencesSpaceDerivative(FiniteDifferencesSpaceDerivativeBas
             out[...] *= scale
 
         self.dFout.exchange_ghosts()
+
diff --git a/hysop/backend/host/python/operator/directional/advection_dir.py b/hysop/backend/host/python/operator/directional/advection_dir.py
index fb80b950f..28477a5cd 100644
--- a/hysop/backend/host/python/operator/directional/advection_dir.py
+++ b/hysop/backend/host/python/operator/directional/advection_dir.py
@@ -17,6 +17,10 @@ DEBUG = False
 class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperator):
     counter = 0
 
+    @debug
+    def __kwds__(cls, **kwds):
+        return super(PythonDirectionalAdvection, cls).__new__(cls, **kwds)
+
     @debug
     def __init__(self, **kwds):
         super(PythonDirectionalAdvection, self).__init__(**kwds)
@@ -92,7 +96,7 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
         velo_compute_view = velo_mesh.local_compute_slices
         self._velocity_mesh_attributes = (velo_mesh_iterator, dx, inv_dx, velo_compute_view, X0)
 
-        dsinputs0 = dsinputs.values()[0]
+        dsinputs0 = next(iter(dsinputs.values()))
         scalar_mesh = dsinputs0.mesh
         scalar_mesh_iterator = scalar_mesh.build_compute_mesh_iterator(cr)
         N0          = scalar_mesh.global_start[-1]
@@ -140,8 +144,8 @@ class PythonDirectionalAdvection(DirectionalAdvectionBase, HostDirectionalOperat
             else:
                 def dump(*args, **kwds):
                     pass
-            Sin  = self.dadvected_fields_in.values()[0]
-            Sout = self.dadvected_fields_out.values()[0]
+            Sin  = next(iter(self.dadvected_fields_in.values()))
+            Sout = next(iter(self.dadvected_fields_out.values()))
             P    = self.dposition
             print('DT= {}'.format(dt))
             self._compute_advection(dt)
diff --git a/hysop/backend/host/python/operator/flowrate_correction.py b/hysop/backend/host/python/operator/flowrate_correction.py
index b8382624a..f380575a3 100644
--- a/hysop/backend/host/python/operator/flowrate_correction.py
+++ b/hysop/backend/host/python/operator/flowrate_correction.py
@@ -21,6 +21,14 @@ class PythonFlowRateCorrection(HostOperator):
     Velocity is corrected in Y and Z-direction from mean vorticity
     """
 
+    @debug
+    def __new__(cls, velocity, vorticity,
+                 dt, flowrate, variables, absorption_start=None,
+                 implementation=None, **kwds):
+        return super(PythonFlowRateCorrection, cls).__new__(cls,
+            input_fields=None, output_fields=None,
+            input_params=None, **kwds)
+
     @debug
     def __init__(self, velocity, vorticity,
                  dt, flowrate, variables, absorption_start=None,
diff --git a/hysop/backend/host/python/operator/penalization.py b/hysop/backend/host/python/operator/penalization.py
index 4fb5ae297..d1aa890e3 100644
--- a/hysop/backend/host/python/operator/penalization.py
+++ b/hysop/backend/host/python/operator/penalization.py
@@ -37,6 +37,14 @@ class PythonPenalizeVorticity(HostOperator):
         am.update(cls.__available_methods)
         return am
 
+    @debug
+    def __new__(cls, obstacles, variables,
+                 velocity, vorticity,
+                 dt, coeff=None, ubar=None, formulation=None, **kwds):
+        return super(PythonPenalizeVorticity, cls).__new__(cls,
+            input_fields=None, output_fields=None,
+            input_params=None, **kwds)
+
     @debug
     def __init__(self, obstacles, variables,
                  velocity, vorticity,
diff --git a/hysop/backend/host/python/operator/spatial_filtering.py b/hysop/backend/host/python/operator/spatial_filtering.py
index eabfbf5a2..77d23cef7 100644
--- a/hysop/backend/host/python/operator/spatial_filtering.py
+++ b/hysop/backend/host/python/operator/spatial_filtering.py
@@ -16,6 +16,7 @@ from hysop.operator.base.spatial_filtering import (
 
 
 class PythonPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, HostOperator):
+
     def discretize(self, **kwds):
         if self.discretized:
             return
@@ -41,6 +42,7 @@ class PythonPolynomialInterpolationFilter(PolynomialInterpolationFilterBase, Hos
 
 
 class PythonPolynomialRestrictionFilter(PolynomialRestrictionFilterBase, HostOperator):
+
     def discretize(self, **kwds):
         if self.discretized:
             return
@@ -64,6 +66,7 @@ class PythonPolynomialRestrictionFilter(PolynomialRestrictionFilterBase, HostOpe
             fout[idx] = (Rr*fin[islc]).sum()
         self.dFout.exchange_ghosts()
 
+
 class PythonRemeshRestrictionFilter(RemeshRestrictionFilterBase, HostOperator):
     """
     Python implementation for lowpass spatial filtering: small grid -> coarse grid
diff --git a/hysop/backend/host/python/operator/transpose.py b/hysop/backend/host/python/operator/transpose.py
index fab4aee04..de9da3702 100644
--- a/hysop/backend/host/python/operator/transpose.py
+++ b/hysop/backend/host/python/operator/transpose.py
@@ -22,6 +22,10 @@ class PythonTranspose(TransposeOperatorBase, HostOperator):
     Uses basic numpy transposition facility.
     """
 
+    @debug
+    def __new__(cls, **kwds):
+        return super(PythonTranspose, cls).__new__(cls, **kwds)
+
     @debug
     def __init__(self, **kwds):
         """Initialize a Transpose operator on the python backend."""
diff --git a/hysop/backend/host/python/operator/vorticity_absorption.py b/hysop/backend/host/python/operator/vorticity_absorption.py
index d5f317291..5f43942e6 100644
--- a/hysop/backend/host/python/operator/vorticity_absorption.py
+++ b/hysop/backend/host/python/operator/vorticity_absorption.py
@@ -18,6 +18,15 @@ class PythonVorticityAbsorption(HostOperator):
     at domain outlet.
     """
 
+    @debug
+    def __new__(cls, velocity, vorticity,
+                 dt, flowrate, start_coord, variables,
+                 custom_filter=None,
+                 implementation=None, **kwds):
+        return super(PythonVorticityAbsorption, cls).__new__(cls,
+            input_fields=None, output_fields=None,
+            input_params=None, **kwds)
+
     @debug
     def __init__(self, velocity, vorticity,
                  dt, flowrate, start_coord, variables,
diff --git a/hysop/core/graph/computational_node.py b/hysop/core/graph/computational_node.py
index e1cfedcd8..3ef5a57c7 100644
--- a/hysop/core/graph/computational_node.py
+++ b/hysop/core/graph/computational_node.py
@@ -380,7 +380,7 @@ class ComputationalGraphNode(OperatorBase, metaclass=ABCMeta):
 
         if ('mpi_params' in self.__kwds) and ('ComputationalGraph' not in map(lambda c: c.__name__, self.__class__.__mro__)):
             mpi_params = self.__kwds['mpi_params']
-            for topo in set(self.input_fields.values() + self.output_fields.values()):
+            for topo in set(self.input_fields.values()).union(self.output_fields.values()):
                 if isinstance(topo, Topology) and (topo.mpi_params != mpi_params):
                     msg = 'MPI parameters mismatch between already specified topology mpi_params '
                     msg += 'and operator MPI paramaters in operator {}.'.format(self.name)
diff --git a/hysop/core/graph/computational_node_frontend.py b/hysop/core/graph/computational_node_frontend.py
index c67bf0626..35d830954 100644
--- a/hysop/core/graph/computational_node_frontend.py
+++ b/hysop/core/graph/computational_node_frontend.py
@@ -11,6 +11,16 @@ from hysop.topology.topology import Topology
 
 class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
 
+    @debug
+    def __new__(cls, implementation=None, base_kwds=None,
+            candidate_input_tensors=None, candidate_output_tensors=None,
+            **impl_kwds):
+        base_kwds = {} if (base_kwds is None) else base_kwds
+        return super(ComputationalGraphNodeFrontend, cls).__new__(cls,
+                candidate_input_tensors=candidate_input_tensors,
+                candidate_output_tensors=candidate_output_tensors,
+                **base_kwds)
+
     @debug
     def __init__(self, implementation=None, base_kwds=None,
             candidate_input_tensors=None, candidate_output_tensors=None,
@@ -191,6 +201,11 @@ class MultiComputationalGraphNodeFrontend(ComputationalGraphNodeFrontend):
     (ComputationalGraphNodeFrontend only generates OPERATOR from IMPLEMENTATION).
     """
 
+    @debug
+    def __new__(cls, implementation_key, implementation=None, **kwds):
+        return super(MultiComputationalGraphNodeFrontend, cls).__new__(cls,
+                implementation=implementation, **kwds)
+
     @debug
     def __init__(self, implementation_key, implementation=None, **kwds):
         """
diff --git a/hysop/core/graph/continuous.py b/hysop/core/graph/continuous.py
index 5374fd349..e5e545d54 100755
--- a/hysop/core/graph/continuous.py
+++ b/hysop/core/graph/continuous.py
@@ -25,7 +25,13 @@ class OperatorBase(TaggedObject, metaclass=ABCMeta):
     @debug
     def __new__(cls, name, fields, tensor_fields, parameters,
                        mpi_params=None, io_params=False, **kwds):
-        return super(OperatorBase, cls).__new__(cls, tagged_cls=OperatorBase, tag_prefix='node', **kwds)
+        try:
+            return super(OperatorBase, cls).__new__(cls, tagged_cls=OperatorBase, tag_prefix='node', **kwds)
+        except TypeError:
+            if ('cl_env' not in kwds):
+                raise
+            kwds.pop('cl_env')
+            return super(OperatorBase, cls).__new__(cls, tagged_cls=OperatorBase, tag_prefix='node', **kwds)
 
     @debug
     def __init__(self, name, fields, tensor_fields, parameters, mpi_params=None, io_params=False, **kwds):
@@ -55,7 +61,13 @@ class OperatorBase(TaggedObject, metaclass=ABCMeta):
         mpi_params: MPIParams
             File i/o config (filename, format ...)
         """
-        super(OperatorBase,self).__init__(tagged_cls=OperatorBase, tag_prefix='node', **kwds)
+        try:
+            super(OperatorBase,self).__init__(tagged_cls=OperatorBase, tag_prefix='node', **kwds)
+        except TypeError:
+            if ('cl_env' not in kwds):
+                raise
+            kwds.pop('cl_env')
+            super(OperatorBase,self).__init__(tagged_cls=OperatorBase, tag_prefix='node', **kwds)
 
         check_instance(fields, tuple, values=ScalarField)
         check_instance(tensor_fields, tuple, values=TensorField)
diff --git a/hysop/core/graph/graph_builder.py b/hysop/core/graph/graph_builder.py
index 4527f390d..2a1a7bfbf 100644
--- a/hysop/core/graph/graph_builder.py
+++ b/hysop/core/graph/graph_builder.py
@@ -619,10 +619,10 @@ class GraphBuilder(object):
 
                 assert len(op.input_fields)  == 1
                 assert len(op.output_fields) == 1
-                assert op.input_fields.keys()[0]   == field
-                assert op.output_fields.keys()[0]  == field
-                assert op.input_fields.values()[0] == src_topo
-                dst_topo = op.output_fields.values()[0]
+                assert next(iter(op.input_fields))   == field
+                assert next(iter(op.output_fields))  == field
+                assert next(iter(op.input_fields.values())) == src_topo
+                dst_topo = next(iter(op.output_fields.values()))
 
                 op_node = self.add_vertex(graph, op)
 
diff --git a/hysop/core/graph/node_generator.py b/hysop/core/graph/node_generator.py
index 402f61fac..2ded034a2 100644
--- a/hysop/core/graph/node_generator.py
+++ b/hysop/core/graph/node_generator.py
@@ -9,12 +9,16 @@ class ComputationalGraphNodeGenerator(object, metaclass=ABCMeta):
     A class that can generate multiple hysop.core.graph.ComputationalGraphNode.
     """
 
+    @debug
+    def __new__(cls, candidate_input_tensors, candidate_output_tensors,
+            name=None, pretty_name=None, **kwds):
+        return super(ComputationalGraphNodeGenerator, cls).__new__(cls, **kwds)
+
     @debug
     def __init__(self, candidate_input_tensors, candidate_output_tensors,
-            name=None, pretty_name=None,
-            **kwds):
+            name=None, pretty_name=None, **kwds):
         super(ComputationalGraphNodeGenerator, self).__init__(**kwds)
-        candidate_input_tensors = first_not_None(candidate_input_tensors, ())
+        candidate_input_tensors  = first_not_None(candidate_input_tensors, ())
         candidate_output_tensors = first_not_None(candidate_output_tensors, ())
 
         self.name = name or self.__class__.__name__
diff --git a/hysop/fields/cartesian_discrete_field.py b/hysop/fields/cartesian_discrete_field.py
index 259e13760..4cb537849 100644
--- a/hysop/fields/cartesian_discrete_field.py
+++ b/hysop/fields/cartesian_discrete_field.py
@@ -1114,12 +1114,13 @@ class CartesianDiscreteScalarFieldView(CartesianDiscreteScalarFieldViewContainer
             tstate._is_read_only = is_read_only
 
         bfield = self._dfield._field
-        btopo = self._dfield._topology
+        btopo  = self._dfield._topology
 
         field = bfield.field_like(name=name, pretty_name=pretty_name,
                                   latex_name=latex_name, var_name=var_name,
                                   initial_values=initial_values, dtype=dtype,
-                                  lboundaries=lboundaries, rboundaries=rboundaries)
+                                  lboundaries=lboundaries, rboundaries=rboundaries,
+                                  nb_components=kwds.pop('nb_components', None))
 
         topology = btopo.topology_like(backend=backend,
                                        grid_resolution=grid_resolution, ghosts=ghosts,
diff --git a/hysop/fields/ghost_exchangers.py b/hysop/fields/ghost_exchangers.py
index f22ffbecf..1866d6325 100644
--- a/hysop/fields/ghost_exchangers.py
+++ b/hysop/fields/ghost_exchangers.py
@@ -211,7 +211,7 @@ class GhostExchanger(GhostExchangerI):
         self.mesh = topology.mesh
         self.backend = data[0].backend if hasattr(data[0], 'backend') else topology.backend
         self.host_backend = self.backend.host_array_backend
-        self.base_tag = int(hashlib.sha1(name).hexdigest(), 16) % (104729)
+        self.base_tag = int(hashlib.sha1(name.encode('utf-8')).hexdigest(), 16) % (104729)
         self.exchange_method = exchange_method
         self.ghost_op = ghost_op
         self.ghost_mask = ghost_mask
diff --git a/hysop/numerics/remesh/kernel_generator.py b/hysop/numerics/remesh/kernel_generator.py
index 34666fe26..3711b9042 100644
--- a/hysop/numerics/remesh/kernel_generator.py
+++ b/hysop/numerics/remesh/kernel_generator.py
@@ -4,7 +4,7 @@ from hysop.tools.numerics   import mpq,mpfr,mpqize,f2q
 from hysop.tools.cache      import load_data_from_cache, update_cache
 
 class Kernel(object):
-    def __init__(self, register=False, verbose=False, split_polys=False, **kargs):
+    def __init__(self, register=False, verbose=True, split_polys=False, **kargs):
         """
         Use SymmetricKernelGenerator or KernelGenerator to generate a kernel.
         Do not call directly.
@@ -28,7 +28,7 @@ class Kernel(object):
     @classmethod
     def hash_kernel_key(cls,n,r,deg,Ms,H,remesh):
         s = '{}_{}_{}_{}_{}_{}'.format(n,r,deg,Ms,H,int(remesh))
-        return hashlib.sha256(s).hexdigest()
+        return hashlib.sha256(s.encode('utf-8')).hexdigest()
 
     @classmethod
     def cache_file(cls):
@@ -82,19 +82,19 @@ class Kernel(object):
             Pt_r = []
             Pt_L = []
             Pt_R = []
-            Cl = np.empty(shape=(deg+1,2*Ms),dtype=float)
-            Cr = np.empty(shape=(deg+1,2*Ms),dtype=float)
-            CL = np.empty(shape=(deg+1,2*Ms),dtype=float)
-            CR = np.empty(shape=(deg+1,2*Ms),dtype=float)
+            Cl = np.empty(shape=(deg+1,2*Ms), dtype=np.float64)
+            Cr = np.empty(shape=(deg+1,2*Ms), dtype=np.float64)
+            CL = np.empty(shape=(deg+1,2*Ms), dtype=np.float64)
+            CR = np.empty(shape=(deg+1,2*Ms), dtype=np.float64)
             for i,Pix in enumerate(P):
-                Pit_l = Pix.xreplace({x:t+X[i]})
-                Pit_r = Pix.xreplace({x:X[i]+1-t})
-                Pit_L = Pix.xreplace({x:+1-t+X[i]})
-                Pit_R = Pix.xreplace({x:X[i]-t})
-                Cl[:,i] = np.asarray(Pit_l.all_coeffs(), dtype=float)
-                Cr[:,i] = np.asarray(Pit_r.all_coeffs(), dtype=float)
-                CL[:,i] = np.asarray(Pit_L.all_coeffs(), dtype=float)
-                CR[:,i] = np.asarray(Pit_R.all_coeffs(), dtype=float)
+                Pit_l = sm.polys.polytools.poly(Pix.as_expr().xreplace({x:t+X[i]}), t, domain='QQ')
+                Pit_r = sm.polys.polytools.poly(Pix.as_expr().xreplace({x:X[i]+1-t}), t, domain='QQ')
+                Pit_L = sm.polys.polytools.poly(Pix.as_expr().xreplace({x:+1-t+X[i]}), t, domain='QQ')
+                Pit_R = sm.polys.polytools.poly(Pix.as_expr().xreplace({x:X[i]-t}), t, domain='QQ')
+                Cl[:,i] = np.asarray(Pit_l.all_coeffs(), dtype=np.float64)
+                Cr[:,i] = np.asarray(Pit_r.all_coeffs(), dtype=np.float64)
+                CL[:,i] = np.asarray(Pit_L.all_coeffs(), dtype=np.float64)
+                CR[:,i] = np.asarray(Pit_R.all_coeffs(), dtype=np.float64)
                 Pt_l.append(Pit_l)
                 Pt_r.append(Pit_r)
                 Pt_L.append(Pit_L)
@@ -125,13 +125,13 @@ class Kernel(object):
         else:
             Pt   = []
             Pt_L = []
-            C = np.empty(shape=(deg+1,2*Ms),dtype=float)
-            CL = np.empty(shape=(deg+1,2*Ms),dtype=float)
+            C  = np.empty(shape=(deg+1,2*Ms), dtype=np.float64)
+            CL = np.empty(shape=(deg+1,2*Ms), dtype=np.float64)
             for i,Pix in enumerate(P):
-                Pit = Pix.xreplace({x:t+X[i]})
-                Pit_L = Pix.xreplace({x:1-t+X[i]})
-                C[:,i] = np.asarray(Pit.all_coeffs(), dtype=float)
-                CL[:,i] = np.asarray(Pit_L.all_coeffs(), dtype=float)
+                Pit   = sm.polys.polytools.poly(Pix.as_expr().xreplace({x:t+X[i]}), t, domain='QQ')
+                Pit_L = sm.polys.polytools.poly(Pix.as_expr().xreplace({x:1-t+X[i]}), t, domain='QQ')
+                C[:,i]  = np.asarray(Pit.all_coeffs(), dtype=np.float64)
+                CL[:,i] = np.asarray(Pit_L.all_coeffs(), dtype=np.float64)
                 Pt.append(Pit)
                 Pt_L.append(Pit_L)
             self.Pt = Pt
@@ -289,7 +289,7 @@ class SymmetricKernelGenerator(object):
                 print('  Cache overwrite requested.')
         else:
             if verbose:
-                print('  Checking if kernel was already computed... ',)
+                print('  Checking if kernel was already computed... ', end=' ')
             data = load_data_from_cache(Kernel.cache_file(), key)
             if (data is not None):
                 if verbose:
@@ -359,17 +359,17 @@ class SymmetricKernelGenerator(object):
             eqs.append(eq)
         # Right-most point, zero derivatives -> (r+1)    eqs
         for d in range(r+1):
-            eq = dPs[-1][d].xreplace({x:Ms}) # =0
+            eq = dPs[-1][d].xreplace({x:sm.Integer(Ms)}) # =0
             eqs.append(eq.as_expr())
         # Cr-continuity on inner points      -> (Ms-1)*(r+1) eqs
         for d in range(r+1):
             for i in range(Ms-1):
-                eq = dPs[i][d].xreplace({x:i+1}) - dPs[i+1][d].xreplace({x:i+1}) # = 0
+                eq = dPs[i][d].xreplace({x:sm.Integer(i+1)}) - dPs[i+1][d].xreplace({x:sm.Integer(i+1)}) # = 0
                 eqs.append(eq.as_expr())
 
         ### Interpolation condition on the left -> Ms equations
         for i in range(Ms):
-            eq = Pp[i].xreplace({x:i}) - H[Ms+i] # = 0
+            eq = Pp[i].xreplace({x:sm.Integer(i)}) - H[Ms+i] # = 0
             eqs.append(eq.as_expr())
 
         ### Discrete moments
@@ -379,7 +379,7 @@ class SymmetricKernelGenerator(object):
             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)})
+                e = P[i].xreplace({x:s-f2q(l)}).as_expr()
                 if m>0: e *= f2q(l**m)
                 expr += e
             Pm = sm.polys.polytools.poly(expr,s)
@@ -438,13 +438,14 @@ if __name__=='__main__':
 
         if len(kernels)==0:
             continue
+        continue
         k0 = kernels[0]
 
         fig = plt.figure()
         plt.xlabel(r'$x$')
         plt.ylabel(r'$\Lambda_{'+'{},{}'.format(p,'r')+'}$')
         X = np.linspace(-k0.Ms-1,+k0.Ms+1,1000)
-        s = plt.subplot(1,1,1)
+        s = plt.subplot(1,1,1, label=i)
         for i,k in enumerate(kernels):
             s.plot(X,k(X),label=r'$\Lambda_{'+'{},{}'.format(p,k.r)+'}$')
         s.plot(k0.I,k0.H,'or')
diff --git a/hysop/numerics/splitting/directional_splitting.py b/hysop/numerics/splitting/directional_splitting.py
index 023705755..3e928c74f 100644
--- a/hysop/numerics/splitting/directional_splitting.py
+++ b/hysop/numerics/splitting/directional_splitting.py
@@ -6,6 +6,11 @@ from hysop.tools.types import first_not_None
 
 class DirectionalSplitting(ComputationalGraphNodeGenerator):
 
+    @debug
+    def __new__(cls, splitting_dim, extra_kwds, **kargs):
+        return super(DirectionalSplitting, cls).__new__(cls,
+                candidate_input_tensors=None, candidate_output_tensors=None, **kargs)
+
     @debug
     def __init__(self, splitting_dim, extra_kwds, **kargs):
         super(DirectionalSplitting,self).__init__(candidate_input_tensors=None,
diff --git a/hysop/numerics/splitting/strang.py b/hysop/numerics/splitting/strang.py
index 3d9fbec97..ac842f22e 100644
--- a/hysop/numerics/splitting/strang.py
+++ b/hysop/numerics/splitting/strang.py
@@ -8,6 +8,12 @@ StrangOrder = EnumFactory.create('StrangOrder',
 
 class StrangSplitting(DirectionalSplitting):
 
+    @debug
+    def __new__(cls, splitting_dim, order, extra_kwds=None, **kargs):
+        return super(StrangSplitting, cls).__new__(cls,
+                splitting_dim=splitting_dim,
+                extra_kwds=extra_kwds, **kargs)
+
     @debug
     def __init__(self, splitting_dim, order, extra_kwds=None, **kargs):
         super(StrangSplitting,self).__init__(splitting_dim=splitting_dim,
@@ -24,12 +30,12 @@ class StrangSplitting(DirectionalSplitting):
             msg=msg.format(order)
             raise ValueError(msg)
 
-        directions = range(splitting_dim)
+        directions = tuple(range(splitting_dim))
         if order==1:
             dt_coeff = (1.0,)*splitting_dim
         elif order==2:
             dt_coeff = (0.5,)*(2*splitting_dim)
-            directions += range(splitting_dim-1,-1,-1)
+            directions += tuple(range(splitting_dim-1,-1,-1))
         self.directions = directions
         self.dt_coeff = dt_coeff
 
diff --git a/hysop/numerics/stencil/stencil.py b/hysop/numerics/stencil/stencil.py
index cdf3d4efb..7135ac813 100644
--- a/hysop/numerics/stencil/stencil.py
+++ b/hysop/numerics/stencil/stencil.py
@@ -276,7 +276,7 @@ class Stencil(object):
             value = factor*self.coeffs[x]
             return (offset,value)
         iterator = np.ndindex(self.shape)
-        iterator = it.imap(mapfun, iterator)
+        iterator = map(mapfun, iterator)
         iterator = filter(lambda x: x[1]!=0, iterator)
         return iterator
 
diff --git a/hysop/numerics/stencil/stencil_generator.py b/hysop/numerics/stencil/stencil_generator.py
index f621e998c..7c8fe7391 100644
--- a/hysop/numerics/stencil/stencil_generator.py
+++ b/hysop/numerics/stencil/stencil_generator.py
@@ -299,7 +299,7 @@ class StencilGenerator(object):
     def _hash_equations(cls, eqs):
         hasher = cls._hash_algorithm()
         for e in eqs:
-            hasher.update(str(e))
+            hasher.update(str(e).encode('utf-8'))
         return hasher.hexdigest()
 
     def get_config(self):
diff --git a/hysop/operator/base/advection_dir.py b/hysop/operator/base/advection_dir.py
index 7c6331e0d..addedf7b3 100644
--- a/hysop/operator/base/advection_dir.py
+++ b/hysop/operator/base/advection_dir.py
@@ -37,6 +37,17 @@ class DirectionalAdvectionBase(object):
         am.update(cls.__available_methods)
         return am
 
+    @debug
+    def __new__(cls, velocity, relative_velocity, splitting_direction,
+                    advected_fields_in, advected_fields_out,
+                    variables, dt, velocity_cfl,
+                    remesh_criteria_eps=None,
+                    **kwds):
+        return super(DirectionalAdvectionBase, cls).__new__(cls,
+                input_fields=None, output_fields=None,
+                input_params=None, output_params=None,
+                splitting_direction=splitting_direction, **kwds)
+
     @debug
     def __init__(self, velocity, relative_velocity, splitting_direction,
                     advected_fields_in, advected_fields_out,
@@ -279,17 +290,17 @@ class DirectionalAdvectionBase(object):
 
         self.is_bilevel = None
         if any(self.dvelocity.compute_resolution != \
-               self.dadvected_fields_out.values()[0].compute_resolution):
+               next(iter(self.dadvected_fields_out.values())).compute_resolution):
             self.is_bilevel = self.dvelocity.compute_resolution
 
     @debug
     def get_work_properties(self):
         requests  = super(DirectionalAdvectionBase, self).get_work_properties()
-        f = self.dadvected_fields_in.values()[0]
+        f = next(iter(self.dadvected_fields_in.values()))
         (pos, request, request_id) = MemoryRequest.cartesian_dfield_like(name='position', dfield=f,
                 ghosts=0, nb_components=1, is_read_only=False)
         requests.push_mem_request(request_id, request)
-        assert all(self.dadvected_fields_in.values()[0].compute_resolution == pos.resolution)
+        assert all(f.compute_resolution == pos.resolution)
         self.dposition = pos
         return requests
 
diff --git a/hysop/operator/base/custom_symbolic_operator.py b/hysop/operator/base/custom_symbolic_operator.py
index 5b03ebb86..f83231f7c 100644
--- a/hysop/operator/base/custom_symbolic_operator.py
+++ b/hysop/operator/base/custom_symbolic_operator.py
@@ -1,5 +1,6 @@
 from abc import ABCMeta
 
+import numpy as np
 from hysop.deps import sm
 from hysop.tools.numpywrappers import npw
 from hysop.tools.types       import check_instance, to_tuple, InstanceOf, first_not_None, to_set
@@ -134,7 +135,8 @@ class ExprDiscretizationInfo(object):
             self.parameters[p.name] = param
 
     def __iadd__(self, rhs):
-        check_instance(rhs, int)
+        check_instance(rhs, (np.integer, int))
+        rhs = int(rhs)
         for (obj, counts) in self.read_counter.items():
             if isinstance(obj, self.IndexedCounterTypes):
                 counts[counts>0] += rhs
@@ -154,7 +156,8 @@ class ExprDiscretizationInfo(object):
         return self
 
     def __imul__(self, rhs):
-        check_instance(rhs, int)
+        check_instance(rhs, (np.integer, int))
+        rhs = int(rhs)
         for (obj, counts) in self.read_counter.items():
             if isinstance(obj, self.IndexedCounterTypes):
                 counts[...] = rhs*counts
diff --git a/hysop/operator/base/derivative.py b/hysop/operator/base/derivative.py
index fd6446bb2..259a996a6 100644
--- a/hysop/operator/base/derivative.py
+++ b/hysop/operator/base/derivative.py
@@ -22,6 +22,14 @@ class SpaceDerivativeBase(object, metaclass=ABCMeta):
     Common implementation interface for derivative operators.
     """
 
+    @debug
+    def __new__(cls, F, dF, A=None,
+            derivative=None, direction=None,
+            variables=None, name=None, pretty_name=None, require_tmp=None, **kwds):
+        return super(SpaceDerivativeBase, cls).__new__(cls, name=name, pretty_name=pretty_name,
+                input_fields=None, output_fields=None,
+                input_params=None, **kwds)
+
     @debug
     def __init__(self, F, dF, A=None,
             derivative=None, direction=None,
@@ -74,7 +82,7 @@ class SpaceDerivativeBase(object, metaclass=ABCMeta):
                  => derivative=(k0,...,kn)
         """
         assert (derivative is not None)
-        A          = first_not_None(A, 1)
+        A = first_not_None(A, 1)
 
         check_instance(F,  ScalarField)
         check_instance(dF, ScalarField, allow_none=True)
@@ -193,6 +201,9 @@ class SpaceDerivativeBase(object, metaclass=ABCMeta):
 
 class FiniteDifferencesSpaceDerivativeBase(SpaceDerivativeBase):
 
+    def __new__(cls, **kwds):
+        return super(FiniteDifferencesSpaceDerivativeBase, cls).__new__(cls, **kwds)
+
     def __init__(self, **kwds):
         super(FiniteDifferencesSpaceDerivativeBase, self).__init__(**kwds)
 
@@ -226,6 +237,11 @@ class FiniteDifferencesSpaceDerivativeBase(SpaceDerivativeBase):
 
 
 class SpectralSpaceDerivativeBase(SpectralOperatorBase, SpaceDerivativeBase):
+    @debug
+    def __new__(cls, testing=False, **kwds):
+        kwds['require_tmp'] = False
+        return super(SpectralSpaceDerivativeBase, cls).__new__(cls, **kwds)
+
     @debug
     def __init__(self, testing=False, **kwds):
         kwds['require_tmp'] = False
diff --git a/hysop/operator/base/spatial_filtering.py b/hysop/operator/base/spatial_filtering.py
index 077bbe60a..1d2fad979 100644
--- a/hysop/operator/base/spatial_filtering.py
+++ b/hysop/operator/base/spatial_filtering.py
@@ -29,6 +29,14 @@ class SpatialFilterBase(object):
     Common base implementation for lowpass spatial filtering: small grid -> coarse grid
     """
 
+    def __new__(cls, input_field, output_field,
+                     input_topo,  output_topo,
+                     **kwds):
+        return super(SpatialFilterBase, cls).__new__(cls,
+                input_fields=None,
+                output_fields=None,
+                **kwds)
+
     def __init__(self, input_field, output_field,
                        input_topo,  output_topo,
                        **kwds):
@@ -145,6 +153,13 @@ class SpectralRestrictionFilterBase(RestrictionFilterBase, SpectralOperatorBase)
     Base implementation for lowpass spatial filtering: small grid -> coarse grid
     using the spectral method.
     """
+
+    @debug
+    def __new__(cls, plot_input_energy=None,
+                     plot_output_energy=None,
+                     **kwds):
+        return super(SpectralRestrictionFilterBase, cls).__new__(cls, **kwds)
+
     @debug
     def __init__(self, plot_input_energy=None,
                        plot_output_energy=None,
diff --git a/hysop/operator/base/transpose_operator.py b/hysop/operator/base/transpose_operator.py
index 867b0cb01..0bf258ce7 100644
--- a/hysop/operator/base/transpose_operator.py
+++ b/hysop/operator/base/transpose_operator.py
@@ -14,6 +14,15 @@ class TransposeOperatorBase(object, metaclass=ABCMeta):
     Common implementation interface for transposition operators.
     """
 
+    @debug
+    def __new__(cls, input_field, output_field, variables, axes,
+                   name=None, pretty_name=None, **kwds):
+        return super(TransposeOperatorBase, cls).__new__(cls,
+                input_fields=None,
+                output_fields=None,
+                name=name, pretty_name=pretty_name,
+                **kwds)
+
     @debug
     def __init__(self, input_field, output_field, variables, axes,
                    name=None, pretty_name=None, **kwds):
@@ -87,7 +96,7 @@ class TransposeOperatorBase(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()))
         axes = self.axes
         ostate.axes = tuple( istate.axes[i] for i in axes )
         return ostate
@@ -140,7 +149,7 @@ class TransposeOperatorBase(object, metaclass=ABCMeta):
         """
         from hysop.tools.transposition_states import TranspositionState
         assert candidate_axes, 'candidate axes is None or empty.'
-        dim = len(candidate_axes.keys()[0])
+        dim = len(next(iter(candidate_axes)))
         tstates = TranspositionState[dim]
         check_instance(candidate_axes, dict, keys=tuple,
                 values=(tstates, type(None)))
@@ -150,5 +159,5 @@ class TransposeOperatorBase(object, metaclass=ABCMeta):
         else:
             idx = 0
 
-        axes = candidate_axes.keys()[idx]
+        axes = tuple(candidate_axes.keys())[idx]
         return axes
diff --git a/hysop/operator/derivative.py b/hysop/operator/derivative.py
index 35e03a2bd..b53ce8b84 100644
--- a/hysop/operator/derivative.py
+++ b/hysop/operator/derivative.py
@@ -34,6 +34,18 @@ class SpaceDerivative(ComputationalGraphNodeFrontend):
         """SpaceDerivative.fd(...) <=> FiniteDifferencesSpaceDerivative(...)"""
         return FiniteDifferencesSpaceDerivative(*args, **kwds)
 
+    @debug
+    def __new__(cls, F, dF, A=None,
+            derivative=None, direction=None,
+            name=None, pretty_name=None,
+            variables=None, implementation=None,
+            base_kwds=None, **kwds):
+        return super(SpaceDerivative, cls).__new__(cls,
+                F=F, dF=dF, A=A,
+                direction=direction, derivative=derivative,
+                variables=variables, implementation=implementation,
+                base_kwds=base_kwds, name=name, pretty_name=pretty_name, **kwds)
+
     @debug
     def __init__(self, F, dF, A=None,
             derivative=None, direction=None,
@@ -175,6 +187,15 @@ class FiniteDifferencesSpaceDerivative(SpaceDerivative):
 class MultiSpaceDerivatives(DirectionalOperatorGeneratorI, ComputationalGraphNodeGenerator):
     """Generate multiple SpaceDerivative operators at once."""
 
+    @debug
+    def __new__(mcls, Fs, dFs, As=None,
+             cls=FiniteDifferencesSpaceDerivative,
+             directions=None, derivatives=None,
+             extra_params=None, base_kwds=None,
+             variables=None, **op_kwds):
+        base_kwds = {} if (base_kwds is None) else base_kwds
+        return super(MultiSpaceDerivatives, mcls).__new__(mcls, **base_kwds)
+
     @debug
     def __init__(self, Fs, dFs, As=None,
              cls=FiniteDifferencesSpaceDerivative,
@@ -195,7 +216,7 @@ class MultiSpaceDerivatives(DirectionalOperatorGeneratorI, ComputationalGraphNod
          """
          from hysop.operator.min_max import MinMaxDerivativeStatistics
          if not issubclass(cls, (SpaceDerivative, MinMaxDerivativeStatistics)) or \
-                (cls in (SpaceDerivative, MinMaxDerivativeStatistics)):
+                        (cls in (SpaceDerivative, MinMaxDerivativeStatistics)):
              msg="cls should be a subclass of SpaceDerivative or MinMaxSpaceDerivativeStatistics, got {}."
              msg+='\ncls MRO is:\n  '
              msg+='\n  '.join(str(t) for t in cls.__mro__)
diff --git a/hysop/operator/directional/advection_dir.py b/hysop/operator/directional/advection_dir.py
index f3d057bfd..8966b1e64 100644
--- a/hysop/operator/directional/advection_dir.py
+++ b/hysop/operator/directional/advection_dir.py
@@ -38,6 +38,16 @@ class DirectionalAdvection(DirectionalOperatorFrontend):
     def default_implementation(cls):
         return Implementation.PYTHON
 
+    @debug
+    def __new__(cls, velocity, advected_fields, variables, dt,
+                advected_fields_out=None, relative_velocity=None,
+                implementation=None, base_kwds=None, **kwds):
+        return super(DirectionalAdvection, cls).__new__(cls,
+                velocity=velocity, relative_velocity=relative_velocity,
+                dt=dt, advected_fields_in=advected_fields, advected_fields_out=advected_fields_out,
+                variables=variables, implementation=implementation, base_kwds=base_kwds,
+                candidate_input_tensors=None, candidate_output_tensors=None, **kwds)
+
     @debug
     def __init__(self, velocity, advected_fields, variables, dt,
                                  advected_fields_out=None,
@@ -55,8 +65,8 @@ class DirectionalAdvection(DirectionalOperatorFrontend):
         advected_fields: Field or array like of Fields
             instance or list of continuous fields to be advected.
         advected_fields_out: Field or array like of Field, optional, defaults to None
-            advection output, if set to None, advection is done inplace 
-            on a per variable basis. 
+            advection output, if set to None, advection is done inplace
+            on a per variable basis.
         relative_velocity: array_like of constants representing relative velocity, optional
             Relative velocity that has to be taken into account.
             Vi = Ui - Urel[i] for each components.
@@ -85,13 +95,13 @@ class DirectionalAdvection(DirectionalOperatorFrontend):
         check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors)
         check_instance(base_kwds, dict, keys=str)
         check_instance(dt, ScalarParameter)
-    
+
         advected_fields = to_list(advected_fields)
         nb_fields = len(advected_fields)
         assert len(set(advected_fields))==len(advected_fields)
 
         relative_velocity = to_list(first_not_None(relative_velocity, [0]))
-        relative_velocity = map(lambda x: x if isinstance(x, str) else float(x), relative_velocity)
+        relative_velocity = tuple(map(lambda x: x if isinstance(x, str) else float(x), relative_velocity))
         if len(relative_velocity)==1:
             relative_velocity *= velocity.nb_components
         relative_velocity = tuple(relative_velocity)
@@ -107,9 +117,9 @@ class DirectionalAdvection(DirectionalOperatorFrontend):
                     advected_fields_out[i] = advected_fields[i]
         advected_fields = tuple(advected_fields)
         advected_fields_out = tuple(advected_fields_out)
-        
+
         check_instance(advected_fields,  tuple, values=Field)
-        check_instance(advected_fields_out, tuple, values=Field, 
+        check_instance(advected_fields_out, tuple, values=Field,
                 size=len(advected_fields))
         check_instance(velocity, Field)
         check_instance(relative_velocity, tuple, values=(str,float), size=velocity.nb_components)
@@ -119,10 +129,10 @@ class DirectionalAdvection(DirectionalOperatorFrontend):
 
         candidate_input_tensors  = advected_fields + (velocity,)
         candidate_output_tensors = advected_fields_out
-        
+
         super(DirectionalAdvection, self).__init__(velocity=velocity, relative_velocity=relative_velocity,
                 dt=dt, advected_fields_in=advected_fields, advected_fields_out=advected_fields_out,
-                variables=variables, implementation=implementation, base_kwds=base_kwds, 
+                variables=variables, implementation=implementation, base_kwds=base_kwds,
                 candidate_input_tensors=candidate_input_tensors,
                 candidate_output_tensors=candidate_output_tensors,
                 **kwds)
diff --git a/hysop/operator/directional/directional.py b/hysop/operator/directional/directional.py
index 56e753349..d74fa8c83 100644
--- a/hysop/operator/directional/directional.py
+++ b/hysop/operator/directional/directional.py
@@ -79,9 +79,10 @@ class DirectionalOperatorBase(object):
 
         return requirements
 
+
 class DirectionalOperatorGeneratorI(object):
     def __new__(cls, **kwds):
-        super(DirectionalOperatorGeneratorI, cls).__new__(cls, **kwds)
+        return super(DirectionalOperatorGeneratorI, cls).__new__(cls, **kwds)
 
     def __init__(self, **kwds):
         super(DirectionalOperatorGeneratorI, self).__init__(**kwds)
@@ -307,6 +308,7 @@ class DirectionalOperatorFrontend(DirectionalOperatorGenerator, metaclass=ABCMet
 
     @debug
     def __new__(cls, implementation=None, base_kwds=None, **op_kwds):
+        base_kwds = {} if (base_kwds is None) else base_kwds
         return super(DirectionalOperatorFrontend, cls).__new__(cls, operator=None,
                 base_kwds=base_kwds, **op_kwds)
 
@@ -341,7 +343,7 @@ class DirectionalOperatorFrontend(DirectionalOperatorGenerator, metaclass=ABCMet
             the implementation implementation
         """
 
-        base_kwds = base_kwds or dict()
+        base_kwds = {} if (base_kwds is None) else base_kwds
 
         check_instance(implementation, Implementation, allow_none=True)
         check_instance(base_kwds, dict, keys=str)
diff --git a/hysop/operator/gradient.py b/hysop/operator/gradient.py
index 9eb9c1bbe..b93b9264d 100644
--- a/hysop/operator/gradient.py
+++ b/hysop/operator/gradient.py
@@ -33,6 +33,22 @@ class Gradient(MultiSpaceDerivatives):
     def implementations(cls):
         return SpaceDerivative.implementations()
 
+    @debug
+    def __new__(mcls, F, gradF, directions=None, implementation=None,
+                        cls=FiniteDifferencesSpaceDerivative,
+                        base_kwds=None, **kwds):
+        base_kwds = {} if (base_kwds is None) else base_kwds
+        base_kwds.update(dict(
+                candidate_input_tensors=(F,),
+                candidate_output_tensors=(gradF,)))
+        return super(Gradient, mcls).__new__(mcls,
+                Fs=None, dFs=None, cls=cls,
+                candidate_input_tensors=(F,),
+                candidate_output_tensors=(gradF,),
+                derivatives=None, directions=directions,
+                implementation=implementation,
+                base_kwds=base_kwds, **kwds)
+
     @debug
     def __init__(self, F, gradF, directions=None, implementation=None,
                         cls=FiniteDifferencesSpaceDerivative,
@@ -147,6 +163,19 @@ class MinMaxGradientStatistics(Gradient):
     This will generate multiple MinMaxDerivativeStatistics operators.
     """
 
+    @debug
+    def __new__(mcls, F, gradF=None, directions=None, coeffs=None,
+            Fmin=None, Fmax=None, Finf=None,
+            all_quiet=True, print_tensors=True,
+            name=None, pretty_name=None, pbasename=None, ppbasename=None,
+            variables=None, implementation=None, base_kwds=None,
+            cls=MinMaxFiniteDifferencesDerivativeStatistics,
+            **kwds):
+        return super(MinMaxGradientStatistics, mcls).__new__(mcls,
+                F=F, gradF=gradF,
+                directions=directions, extra_params=None,
+                cls=cls, variables=variables, **kwds)
+
     @debug
     def __init__(self, F, gradF=None, directions=None, coeffs=None,
             Fmin=None, Fmax=None, Finf=None,
diff --git a/hysop/operator/misc.py b/hysop/operator/misc.py
index a7ccd4c63..628c4bdb7 100644
--- a/hysop/operator/misc.py
+++ b/hysop/operator/misc.py
@@ -53,6 +53,16 @@ class ForceTopologyState(Noop):
     If backend is not given, all input fields will be imposed to live on Backend.HOST.
     """
 
+    @debug
+    def __new__(cls, fields, variables,
+            tstate=None, memorder=None,
+            backend=None, extra_kwds=None,
+            mpi_params=None, cl_env=None,
+            **kwds):
+        kwds.setdefault('mpi_params', None)
+        return super(ForceTopologyState, cls).__new__(cls,
+                input_fields=None, output_fields=None, **kwds)
+
     @debug
     def __init__(self, fields, variables,
             tstate=None, memorder=None,
@@ -82,8 +92,8 @@ class ForceTopologyState(Noop):
         cl_env = first_not_None(cl_env, extra_kwds.get('cl_env', None))
         mpi_params = first_not_None(mpi_params, extra_kwds.get('mpi_params', None), getattr(cl_env, 'mpi_params', None))
 
-        extra_kwds.setdefault('cl_env',  cl_env)
         extra_kwds.setdefault('mpi_params', mpi_params)
+        extra_kwds.setdefault('cl_env',  cl_env)
         kwds.setdefault('cl_env',  cl_env)
         kwds.setdefault('mpi_params', mpi_params)
 
diff --git a/hysop/operator/spatial_filtering.py b/hysop/operator/spatial_filtering.py
index 2ec022e1d..2f3746a3e 100644
--- a/hysop/operator/spatial_filtering.py
+++ b/hysop/operator/spatial_filtering.py
@@ -17,6 +17,18 @@ FilteringMethod = EnumFactory.create('FilteringMethod',
 ['SPECTRAL', 'REMESH', 'POLYNOMIAL', 'SUBGRID'])
 
 class SpatialFilterFrontend(MultiComputationalGraphNodeFrontend):
+
+    def __new__(cls, input_variable, output_variable,
+                 filtering_method, implementation=None,
+                 base_kwds=None,
+                 **kwds):
+        return super(SpatialFilterFrontend, cls).__new__(cls,
+                input_field=None, input_topo=None,
+                output_field=None, output_topo=None,
+                implementation_key=filtering_method,
+                implementation=implementation,
+                base_kwds=base_kwds, **kwds)
+
     def __init__(self, input_variable, output_variable,
                  filtering_method, implementation=None,
                  base_kwds=None,
@@ -145,6 +157,17 @@ class SpatialFilter(ComputationalGraphNodeGenerator):
     """
     Graphnode generator to generate interpolation or restriction filter for multiple fields at once.
     """
+
+    @debug
+    def __new__(cls, input_variables, output_variables,
+                 filtering_method, implementation=None,
+                 base_kwds=None, **kwds):
+        base_kwds = first_not_None(base_kwds, {})
+        return super(SpatialFilter, cls).__new__(cls,
+                candidate_input_tensors=None,
+                candidate_output_tensors=None,
+                **base_kwds)
+
     @debug
     def __init__(self, input_variables, output_variables,
                  filtering_method, implementation=None,
diff --git a/hysop/operator/tests/test_absorption.py b/hysop/operator/tests/test_absorption.py
index e790ac7dd..cd6b4430a 100644
--- a/hysop/operator/tests/test_absorption.py
+++ b/hysop/operator/tests/test_absorption.py
@@ -72,9 +72,9 @@ class TestVorticityAbsorption(object):
 
         domain = Box(length=(1,)*dim)
         velo = Field(domain=domain, name='velo', dtype=dtype,
-                     nb_components=3, register_object=False)
+                     nb_components=3)
         vorti = Field(domain=domain, name='vorti', dtype=dtype,
-                      nb_components=3, register_object=False)
+                      nb_components=3)
 
         self._test_one(shape=shape, dim=dim, dtype=dtype,
                        domain=domain, velo=velo, vorti=vorti,
diff --git a/hysop/operator/tests/test_directional_advection.py b/hysop/operator/tests/test_directional_advection.py
index 6e87321a5..c860edb21 100644
--- a/hysop/operator/tests/test_directional_advection.py
+++ b/hysop/operator/tests/test_directional_advection.py
@@ -67,11 +67,11 @@ class TestDirectionalAdvectionOperator(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)
             Sin  = Field(domain=domain, name='Sin', dtype=dtype,
-                    nb_components=2, register_object=False)
+                    nb_components=2)
             Sout = Field(domain=domain, name='Sout', dtype=dtype,
-                    nb_components=2, register_object=False)
+                    nb_components=2)
             for time_integrator in time_integrators:
                 for remesh_kernel in remesh_kernels:
                     for velocity_cfl in velocity_cfls:
@@ -115,10 +115,10 @@ class TestDirectionalAdvectionOperator(object):
 
         # Use optimal timestep, ||Vi||_inf is 1 on a per-axis basis
         dt = velocity_cfl * np.divide(domain.length, shape).min()
-        dt = ScalarParameter('dt', initial_value=dt, constant=True)
+        dt = ScalarParameter('dt', initial_value=dt, const=True)
 
         ref_impl = Implementation.PYTHON
-        implementations = DirectionalAdvection.implementations().keys()
+        implementations = list(DirectionalAdvection.implementations().keys())
         assert ref_impl in implementations
         implementations.remove(ref_impl)
         implementations = [ref_impl] + implementations
diff --git a/hysop/operator/tests/test_fd_derivative.py b/hysop/operator/tests/test_fd_derivative.py
index 02f26f68a..776d803a9 100644
--- a/hysop/operator/tests/test_fd_derivative.py
+++ b/hysop/operator/tests/test_fd_derivative.py
@@ -116,8 +116,8 @@ class TestFiniteDifferencesDerivative(object):
 
         domain = Box(length=(1,)*dim)
         F = Field(domain=domain, name='F', dtype=dtype,
-                nb_components=3, register_object=False)
-        gradF = F.gradient(register_object=False)
+                nb_components=3)
+        gradF = F.gradient()
 
         self._test_one(shape=shape, dim=dim, dtype=dtype,
                 domain=domain, F=F, gradF=gradF)
@@ -126,8 +126,8 @@ class TestFiniteDifferencesDerivative(object):
     def _test_one(self, shape, dim, dtype,
             domain, F, gradF):
 
-        print('\nTesting {}D Gradient: dtype={} shape={}'.format()
-                dim, dtype.__name__, shape)
+        print('\nTesting {}D Gradient: dtype={} shape={}'.format(
+                dim, dtype.__name__, shape))
         print(' >Input analytic functions:')
         for (i,Fi) in enumerate(self.analytic_expressions[dim]['F']):
             print('   *F{}(x,t) = {}'.format(i, Fi))
diff --git a/hysop/operator/tests/test_penalization.py b/hysop/operator/tests/test_penalization.py
index 3d060497b..fe33fc94b 100644
--- a/hysop/operator/tests/test_penalization.py
+++ b/hysop/operator/tests/test_penalization.py
@@ -82,11 +82,11 @@ class TestPenalizeVorticity(object):
 
         domain = Box(length=(1,)*dim)
         velo = Field(domain=domain, name='velo', dtype=dtype,
-                     nb_components=3, register_object=False)
+                     nb_components=3)
         vorti = Field(domain=domain, name='vorti', dtype=dtype,
-                      nb_components=3, register_object=False)
+                      nb_components=3)
         obstacle = Field(domain=domain, name='sphere', dtype=dtype,
-                         nb_components=1, register_object=False)
+                         nb_components=1)
 
         self._test_one(shape=shape, dim=dim, dtype=dtype,
                        domain=domain, velo=velo, vorti=vorti,
diff --git a/hysop/operator/tests/test_restriction_filter.py b/hysop/operator/tests/test_restriction_filter.py
index b522834fd..45228d948 100755
--- a/hysop/operator/tests/test_restriction_filter.py
+++ b/hysop/operator/tests/test_restriction_filter.py
@@ -64,7 +64,7 @@ class TestMultiresolutionFilter(object):
 
         domain = Box(length=(1,)*dim)
         f = Field(domain=domain, name='velo', dtype=dtype,
-                  nb_components=1, register_object=False)
+                  nb_components=1)
 
         self._test_one(shape_f=shape_f, shape_c=shape_c,
                        dim=dim, dtype=dtype,
diff --git a/hysop/operator/tests/test_transpose.py b/hysop/operator/tests/test_transpose.py
index df2b808d5..9008a3880 100644
--- a/hysop/operator/tests/test_transpose.py
+++ b/hysop/operator/tests/test_transpose.py
@@ -66,9 +66,9 @@ class TestTransposeOperator(object):
         domain = Box(length=(1.0,)*dim)
         for nb_components in (2,):
             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)
             for axes in all_axes:
                 for shape in shapes:
                     self._test_one(shape=shape, axes=axes,
diff --git a/hysop/operator/tests/test_velocity_correction.py b/hysop/operator/tests/test_velocity_correction.py
index 0d32f7e14..125c41042 100644
--- a/hysop/operator/tests/test_velocity_correction.py
+++ b/hysop/operator/tests/test_velocity_correction.py
@@ -90,9 +90,9 @@ class TestFlowRateCorrection(object):
 
         domain = Box(length=(1,)*dim)
         velo = Field(domain=domain, name='velo', dtype=dtype,
-                     nb_components=3, register_object=False)
+                     nb_components=3)
         vorti = Field(domain=domain, name='vorti', dtype=dtype,
-                      nb_components=3, register_object=False)
+                      nb_components=3)
 
         self._test_one(shape=shape, dim=dim, dtype=dtype,
                        domain=domain, velo=velo, vorti=vorti)
diff --git a/hysop/operator/transpose.py b/hysop/operator/transpose.py
index 26a2232b4..cf8228cb8 100644
--- a/hysop/operator/transpose.py
+++ b/hysop/operator/transpose.py
@@ -53,6 +53,19 @@ class Transpose(ComputationalGraphNodeGenerator):
         msg+='implementation should match the discrete field topology backend.'
         raise RuntimeError(msg)
 
+    @debug
+    def __new__(cls, fields, variables, axes,
+                output_fields=None,
+                implementation=None,
+                name=None,
+                base_kwds=None,
+                **kwds):
+        base_kwds = {} if (base_kwds is None) else {}
+        return super(Transpose, cls).__new__(cls, name=name,
+                candidate_input_tensors=None,
+                candidate_output_tensors=None,
+                **base_kwds)
+
     @debug
     def __init__(self, fields, variables, axes,
                 output_fields=None,
@@ -149,7 +162,7 @@ class Transpose(ComputationalGraphNodeGenerator):
         candidate_input_tensors  = filter(lambda x: x.is_tensor, input_fields)
         candidate_output_tensors = filter(lambda x: x.is_tensor, output_fields)
 
-        base_kwds = base_kwds or dict()
+        base_kwds = {} if (base_kwds is None) else {}
         super(Transpose,self).__init__(name=name,
                 candidate_input_tensors=candidate_input_tensors,
                 candidate_output_tensors=candidate_output_tensors,
diff --git a/hysop/symbolic/parameter.py b/hysop/symbolic/parameter.py
index dbc9f8db9..c38e93d0d 100644
--- a/hysop/symbolic/parameter.py
+++ b/hysop/symbolic/parameter.py
@@ -22,10 +22,10 @@ class SymbolicTensorParameter(SymbolicTensor):
                 make_scalar_kwds=make_scalar_kwds,
                 scalar_kwds=scalar_kwds, **kwds)
 
-    def __init__(cls, parameter, name=None, pretty_name=None,
+    def __init__(self, parameter, name=None, pretty_name=None,
             value=None, shape=None, scalar_cls=None, **kwds):
         super(SymbolicTensorParameter, self).__init__(
-                name=name, pretty_name=pname,
+                name=name, pretty_name=pretty_name,
                 shape=shape, value=value, scalar_cls=scalar_cls,
                 make_scalar_kwds=None,
                 scalar_kwds=None, **kwds)
-- 
GitLab