From 79b1ebe44b09ef8a896c6c42ce44798c26be1e53 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Keck <Jean-Baptiste.Keck@imag.fr>
Date: Wed, 4 Nov 2020 13:32:21 +0100
Subject: [PATCH] fix min max mpi reduction bug

---
 ci/docker_images/ubuntu/groovy/Dockerfile     |   6 +-
 ci/scripts/test.sh                            |  55 +++---
 .../backend/host/fortran/operator/poisson.py  |   5 +
 .../host/fortran/operator/scales_advection.py |  10 ++
 hysop/backend/host/host_operator.py           |  19 +-
 hysop/core/graph/computational_graph.py       |  10 +-
 .../core/graph/computational_node_frontend.py |   5 +-
 hysop/core/graph/continuous.py                |  16 +-
 hysop/core/tests/test_checkpoint.sh           | 168 +++++++++---------
 hysop/fields/tests/test_cartesian.sh          |   2 +-
 hysop/operator/adapt_timestep.py              |  58 +++---
 hysop/operator/advection.py                   |   2 +-
 hysop/operator/base/external_force.py         |   2 +
 hysop/operator/base/memory_reordering.py      |   6 +-
 hysop/operator/base/min_max.py                |  23 ++-
 hysop/operator/directional/directional.py     |   7 +-
 hysop/operator/hdf_io.py                      |   6 +-
 hysop/operator/mean_field.py                  |   4 +-
 hysop/operator/misc.py                        |   2 +-
 hysop/operator/parameter_plotter.py           |   6 +-
 .../operator/tests/test_bilevel_advection.py  |   6 +-
 hysop/operator/tests/test_scales_advection.py |   6 +-
 hysop/tools/debug_dumper.py                   |  23 +--
 hysop/tools/handle.py                         |  25 ++-
 hysop_examples/argparser.py                   |   2 +-
 .../examples/cylinder/oscillating_cylinder.py |   4 +-
 .../flow_around_sphere/flow_around_sphere.py  |   2 +-
 .../particles_above_salt_bc.py                |   4 +-
 .../particles_above_salt_bc_3d.py             |   4 +-
 .../particles_above_salt_symmetrized.py       |   4 +-
 .../turbulent_scalar_advection.py             |   2 +-
 .../sediment_deposit/sediment_deposit.py      |   4 +-
 .../sediment_deposit_levelset.py              |   4 +-
 .../examples/taylor_green/taylor_green.py     |   5 +-
 .../taylor_green/taylor_green_cpuFortran.py   |   2 +-
 35 files changed, 282 insertions(+), 227 deletions(-)

diff --git a/ci/docker_images/ubuntu/groovy/Dockerfile b/ci/docker_images/ubuntu/groovy/Dockerfile
index 20099148b..8f4cc1513 100644
--- a/ci/docker_images/ubuntu/groovy/Dockerfile
+++ b/ci/docker_images/ubuntu/groovy/Dockerfile
@@ -50,9 +50,9 @@ RUN cd /tmp && \
  cd /tmp && \
  rm -Rf /tmp/hptt
 
-# HDF5 + h5py (v1.10.6 is currently the last supported by h5py)
+# HDF5 1.12.0 + h5py 3.0.0
 RUN cd /tmp && \
- wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.6/src/hdf5-1.10.6.tar.gz && \
+ wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.12/hdf5-1.12.0/src/hdf5-1.12.0.tar.gz && \
  tar -xvzf hdf5-*.tar.gz && \
  rm -f hdf5-*.tar.gz && \
  cd hdf5-* && \
@@ -60,7 +60,7 @@ RUN cd /tmp && \
  make && \
  make install && \
  rm -rf /tmp/hdf5-*
-RUN CC="${MPICC}" HDF5_MPI="ON" HDF5_VERSION="1.10.6" HDF5_DIR="${MPI_ROOT}" pip3.8 install --upgrade --no-binary=h5py h5py
+RUN CC="${MPICC}" HDF5_MPI="ON" HDF5_VERSION="1.12.0" HDF5_DIR="${MPI_ROOT}" pip3.8 install --upgrade --no-binary=h5py h5py
 
 # other python packages (standard primefac package does not support python3)
 RUN pip3.8 install --upgrade scipy sympy matplotlib gmpy2 psutil py-cpuinfo Mako editdistance portalocker colors.py tee pycairo argparse_color_formatter networkx pyvis zarr numcodecs jsonpickle memory-tempfile
diff --git a/ci/scripts/test.sh b/ci/scripts/test.sh
index 19c9a1029..7bd8ca02f 100755
--- a/ci/scripts/test.sh
+++ b/ci/scripts/test.sh
@@ -78,7 +78,8 @@ export PYOPENCL_COMPILER_OUTPUT=0
 # OpenMPI specific variables
 export OMPI_ALLOW_RUN_AS_ROOT=1
 export OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1
-export OMPI_MCA_rmaps_base_oversubscribe=1
+export OMPI_MCA_rmaps_base_oversubscribe=1 
+export OMPI_MCA_btl_vader_single_copy_mechanism=none  # see https://github.com/open-mpi/ompi/issues/4948
 
 echo "Trying to load hysop module:"
 ${PYTHON_EXECUTABLE} -c 'import hysop; print(hysop)'
@@ -93,9 +94,9 @@ RUN_EXAMPLES=${RUN_EXAMPLES:-true}
 RUN_LONG_TESTS=${RUN_LONG_TESTS:-false}
 
 COMMON_TEST_OPTIONS=''
-TEST_DIR="$HYSOP_DIR"
+TEST_DIR="${HYSOP_DIR}"
 COMMON_EXAMPLE_OPTIONS='-VNC -d16 -cp float -maxit 2 --autotuner-max-candidates 1 --save-checkpoint --checkpoint-dump-freq 0 --checkpoint-dump-period 0 --checkpoint-dump-last --checkpoint-dump-times'
-EXAMPLE_DIR="$HYSOP_DIR/../hysop_examples/examples"
+EXAMPLE_DIR="${HYSOP_DIR}/../hysop_examples/examples"
 
 hysop_test() {
      test=$1
@@ -115,29 +116,31 @@ example_test() {
 }
 
 if [ "$RUN_TESTS" = true ]; then
-    hysop_test "core/arrays/tests/test_array.py"
-    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_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
+    #hysop_test "core/arrays/tests/test_array.py"
+    #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_scales_advection.py"
+    #hysop_test "operator/tests/test_bilevel_advection.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
 
 if [ "${RUN_LONG_TESTS}" = true ]; then
diff --git a/hysop/backend/host/fortran/operator/poisson.py b/hysop/backend/host/fortran/operator/poisson.py
index 6f5638ff9..c44b8608b 100644
--- a/hysop/backend/host/fortran/operator/poisson.py
+++ b/hysop/backend/host/fortran/operator/poisson.py
@@ -9,6 +9,11 @@ from hysop.constants import HYSOP_REAL
 
 class PoissonFFTW(FortranFFTWOperator):
 
+    def __new__(cls, Fin, Fout, variables,
+                 extra_input_kwds=None, **kwds):
+        return super(PoissonFFTW, cls).__new__(cls,
+            input_fields=None, output_fields=None, **kwds)
+
     def __init__(self, Fin, Fout, variables,
                  extra_input_kwds=None, **kwds):
         """Operator to solve Poisson equation using FFTW in Fortran.
diff --git a/hysop/backend/host/fortran/operator/scales_advection.py b/hysop/backend/host/fortran/operator/scales_advection.py
index 6d9206b67..1b930d04e 100644
--- a/hysop/backend/host/fortran/operator/scales_advection.py
+++ b/hysop/backend/host/fortran/operator/scales_advection.py
@@ -78,6 +78,16 @@ class ScalesAdvection(FortranOperator):
         am.update(cls.__available_methods)
         return am
 
+    @debug
+    def __new__(cls, velocity,
+                 advected_fields_in, advected_fields_out,
+                 variables, dt, **kwds):
+        return super(ScalesAdvection, cls).__new__(cls,
+            input_fields=None, output_fields=None,
+            input_params=None, output_params=None,
+            **kwds)
+
+
     @debug
     def __init__(self, velocity,
                  advected_fields_in, advected_fields_out,
diff --git a/hysop/backend/host/host_operator.py b/hysop/backend/host/host_operator.py
index b2bbf16ae..368cb9c3e 100644
--- a/hysop/backend/host/host_operator.py
+++ b/hysop/backend/host/host_operator.py
@@ -14,25 +14,32 @@ from hysop.core.graph.computational_operator import ComputationalGraphOperator
 from hysop.topology.topology_descriptor import TopologyDescriptor
 
 
-class HostOperator(ComputationalGraphOperator, metaclass=ABCMeta):
+class HostOperatorBase(ComputationalGraphOperator, metaclass=ABCMeta):
     """
-    Abstract class for discrete operators working on OpenCL backends.
+    Abstract class for discrete operators working on cpu.
+    HostOperatorBase ignore the extra cl_env keyword parameter.
     """
 
     @debug
-    def __new__(cls, **kwds):
-        return super(HostOperator, cls).__new__(cls, **kwds)
+    def __new__(cls, cl_env=None, **kwds):
+        return super(HostOperatorBase, cls).__new__(cls, **kwds)
 
     @debug
-    def __init__(self, **kwds):
+    def __init__(self, cl_env=None, **kwds):
         """
         Create the common attributes of all host operators.
 
         All input and output variable topologies should be of kind
         Backend.HOST and share the same HostEnvironment.
         """
-        super(HostOperator, self).__init__(**kwds)
+        super(HostOperatorBase, self).__init__(**kwds)
+
 
+class HostOperator(ComputationalGraphOperator, metaclass=ABCMeta):
+    """
+    Abstract class for discrete operators working on cpu.
+    HostOperator extra cl_env keyword parameter and enforces HOST backend.
+    """
     @classmethod
     def supported_backends(cls):
         """
diff --git a/hysop/core/graph/computational_graph.py b/hysop/core/graph/computational_graph.py
index 5535f39c6..e20bc9402 100644
--- a/hysop/core/graph/computational_graph.py
+++ b/hysop/core/graph/computational_graph.py
@@ -195,7 +195,7 @@ class ComputationalGraph(ComputationalGraphNode, metaclass=ABCMeta):
 
         ss = '>INPUTS:'
         if sinputs:
-            for (td, sreqs) in sinputs.items():
+            for (td, sreqs) in sorted(sinputs.items()):
                 if isinstance(td, Topology):
                     ss += '\n {}'.format(td.short_description())
                 else:
@@ -205,7 +205,7 @@ class ComputationalGraph(ComputationalGraphNode, metaclass=ABCMeta):
                                       discr_size=discr_size, ghosts_size=ghosts_size,
                                       order_size=order_size, cansplit_size=cansplit_size,
                                       tstates_size=tstates_size)
-                for (opname, fname, discr, ghosts, order, can_split, tstates) in sreqs:
+                for (opname, fname, discr, ghosts, order, can_split, tstates) in sorted(sreqs):
                     ss += template.format(
                         opname, fname, discr, ghosts, order, can_split, tstates,
                         name_size=name_size, field_size=field_size,
@@ -216,7 +216,7 @@ class ComputationalGraph(ComputationalGraphNode, metaclass=ABCMeta):
             ss += ' None'
         ss += '\n>OUTPUTS:'
         if soutputs:
-            for (td, sreqs) in soutputs.items():
+            for (td, sreqs) in sorted(soutputs.items()):
                 if isinstance(td, Topology):
                     ss += '\n {}'.format(td.short_description())
                 else:
@@ -226,7 +226,7 @@ class ComputationalGraph(ComputationalGraphNode, metaclass=ABCMeta):
                                       discr_size=discr_size, ghosts_size=ghosts_size,
                                       order_size=order_size, cansplit_size=cansplit_size,
                                       tstates_size=tstates_size)
-                for (opname, fname, discr, ghosts, order, can_split, tstates) in sreqs:
+                for (opname, fname, discr, ghosts, order, can_split, tstates) in sorted(sreqs):
                     ss += template.format(
                         opname, fname, discr, ghosts, order, can_split, tstates,
                         name_size=name_size, field_size=field_size,
@@ -716,7 +716,7 @@ class ComputationalGraph(ComputationalGraphNode, metaclass=ABCMeta):
         net = self.to_pyvis()
 
         import tempfile
-        with tempfile.NamedTemporaryFile(suffix='.html') as f:
+        with tempfile.NamedTemporaryFile(suffix='.html', delete=False) as f:
             net.show(f.name)
 
     def to_html(self, path, io_rank=0, show_buttons=False):
diff --git a/hysop/core/graph/computational_node_frontend.py b/hysop/core/graph/computational_node_frontend.py
index 35d830954..4595a5348 100644
--- a/hysop/core/graph/computational_node_frontend.py
+++ b/hysop/core/graph/computational_node_frontend.py
@@ -122,10 +122,11 @@ class ComputationalGraphNodeFrontend(ComputationalGraphNodeGenerator):
             sargs = ['*{} = {}'.format(k,v.__class__)
                         for (k,v) in self.impl_kwds.items()]
             msg  = 'FATAL ERROR during {}.generate():\n'
-            msg += ' => failed to call {}.__init__()\n    with the following keywords:'
+            msg += ' => failed to initialize an instance of type {}'
+            msg += '\n    by using the following keyword arguments:'
             msg += '\n     '+'\n     '.join(sargs)
             msg = msg.format(self.__class__, self.impl)
-            print('\n{}\n'.format(msg))
+            print('\n{}'.format(msg))
             raise
 
         for kwds in self._input_fields_to_dump:
diff --git a/hysop/core/graph/continuous.py b/hysop/core/graph/continuous.py
index e5e545d54..5374fd349 100755
--- a/hysop/core/graph/continuous.py
+++ b/hysop/core/graph/continuous.py
@@ -25,13 +25,7 @@ class OperatorBase(TaggedObject, metaclass=ABCMeta):
     @debug
     def __new__(cls, name, fields, tensor_fields, parameters,
                        mpi_params=None, io_params=False, **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)
+        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):
@@ -61,13 +55,7 @@ class OperatorBase(TaggedObject, metaclass=ABCMeta):
         mpi_params: MPIParams
             File i/o config (filename, format ...)
         """
-        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)
+        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/tests/test_checkpoint.sh b/hysop/core/tests/test_checkpoint.sh
index e64700965..f748ad47c 100755
--- a/hysop/core/tests/test_checkpoint.sh
+++ b/hysop/core/tests/test_checkpoint.sh
@@ -1,7 +1,8 @@
 #!/usr/bin/env bash
 set -feu -o pipefail
+set -x
 PYTHON_EXECUTABLE=${PYTHON_EXECUTABLE:-python3.8}
-MPIRUN_EXECUTABLE=${MPIRUN_EXECUTABLE:-mpirun --allow-run-as-root}
+MPIRUN_EXECUTABLE=${MPIRUN_EXECUTABLE:-mpirun}
 
 SCRIPT_DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null 2>&1 && pwd )"
 EXAMPLE_DIR="$(realpath ${SCRIPT_DIR}/../../../hysop_examples/examples)"
@@ -30,94 +31,93 @@ function compare_files {
 #
 # Basic test with 2D diffusion (serial)
 #
-EXAMPLE_FILE="${EXAMPLE_DIR}/scalar_diffusion/scalar_diffusion.py"
-TEST_DIR='/tmp/hysop_tests/checkpoints/scalar_diffusion'
-COMMON_OPTIONS="-NC -impl opencl -cp fp32 -d64 --debug-dump-target dump -nu 0.02 -niter 20 -te 0.1 --dump-tstart 0.05 --dump-freq 1 "
-
-echo
-echo "TEST SCALAR DIFFUSION CHECKPOINT (SERIAL)"
-if [[ ! -f "${EXAMPLE_FILE}" ]]; then
-    echo "Cannot find example file '${EXAMPLE_FILE}'."
-    exit 1
-fi
-
-echo ' Running simulations...'
-"${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -S "${TEST_DIR}/checkpoint0.tar" --dump-dir "${TEST_DIR}/run0" --checkpoint-dump-time 0.05 --checkpoint-dump-freq 0
-"${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -S "${TEST_DIR}/checkpoint1.tar" --dump-dir "${TEST_DIR}/run1" --checkpoint-dump-time 0.05 --checkpoint-dump-freq 0
-
-echo ' Comparing solutions...'
-echo "  >debug dumps match"
-compare_files "${TEST_DIR}/run0/dump/run.txt" "${TEST_DIR}/run1/dump/run.txt"
-for f0 in $(find "${TEST_DIR}/run0" -name '*.h5' | sort -n); do
-    f1=$(echo "${f0}" | sed 's/run0/run1/')
-    compare_files "${f0}" "${f1}"
-    echo "  >$(basename ${f0}) match"
-done
-
-echo
-echo ' Running simulations from checkpoints...'
-"${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -L "${TEST_DIR}/checkpoint0.tar" --dump-dir "${TEST_DIR}/run2"
-"${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -L "${TEST_DIR}/checkpoint1.tar" --dump-dir "${TEST_DIR}/run3"
-
-echo ' Comparing solutions...'
-compare_files "${TEST_DIR}/run2/dump/run.txt" "${TEST_DIR}/run3/dump/run.txt"
-echo "  >debug dumps match"
-for f0 in $(find "${TEST_DIR}/run2" -name '*.h5' | sort -n); do
-    f1=$(echo "${f0}" | sed 's/run2/run3/')
-    f2=$(echo "${f0}" | sed 's/run2/run0/')
-    f3=$(echo "${f0}" | sed 's/run2/run1/')
-    compare_files "${f0}" "${f1}"
-    compare_files "${f0}" "${f2}"
-    compare_files "${f0}" "${f3}"
-    echo "  >$(basename ${f0}) match"
-done
+#EXAMPLE_FILE="${EXAMPLE_DIR}/scalar_diffusion/scalar_diffusion.py"
+#TEST_DIR='/tmp/hysop_tests/checkpoints/scalar_diffusion'
+#COMMON_OPTIONS="-NC -impl opencl -cp fp32 -d64 --debug-dump-target dump -nu 0.02 -niter 20 -te 0.1 --dump-tstart 0.05 --dump-freq 1 "
+
+#echo
+#echo "TEST SCALAR DIFFUSION CHECKPOINT (SERIAL)"
+#if [[ ! -f "${EXAMPLE_FILE}" ]]; then
+    #echo "Cannot find example file '${EXAMPLE_FILE}'."
+    #exit 1
+#fi
+
+#echo ' Running simulations...'
+#"${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -S "${TEST_DIR}/checkpoint0.tar" --dump-dir "${TEST_DIR}/run0" --checkpoint-dump-time 0.05 --checkpoint-dump-freq 0
+#"${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -S "${TEST_DIR}/checkpoint1.tar" --dump-dir "${TEST_DIR}/run1" --checkpoint-dump-time 0.05 --checkpoint-dump-freq 0
+
+#echo ' Comparing solutions...'
+#echo "  >debug dumps match"
+#compare_files "${TEST_DIR}/run0/dump/run.txt" "${TEST_DIR}/run1/dump/run.txt"
+#for f0 in $(find "${TEST_DIR}/run0" -name '*.h5' | sort -n); do
+    #f1=$(echo "${f0}" | sed 's/run0/run1/')
+    #compare_files "${f0}" "${f1}"
+    #echo "  >$(basename ${f0}) match"
+#done
+
+#echo
+#echo ' Running simulations from checkpoints...'
+#"${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -L "${TEST_DIR}/checkpoint0.tar" --dump-dir "${TEST_DIR}/run2"
+#"${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -L "${TEST_DIR}/checkpoint1.tar" --dump-dir "${TEST_DIR}/run3"
+
+#echo ' Comparing solutions...'
+#compare_files "${TEST_DIR}/run2/dump/run.txt" "${TEST_DIR}/run3/dump/run.txt"
+#echo "  >debug dumps match"
+#for f0 in $(find "${TEST_DIR}/run2" -name '*.h5' | sort -n); do
+    #f1=$(echo "${f0}" | sed 's/run2/run3/')
+    #f2=$(echo "${f0}" | sed 's/run2/run0/')
+    #f3=$(echo "${f0}" | sed 's/run2/run1/')
+    #compare_files "${f0}" "${f1}"
+    #compare_files "${f0}" "${f2}"
+    #compare_files "${f0}" "${f3}"
+    #echo "  >$(basename ${f0}) match"
+#done
 
 
 #
 # Basic test with 2D diffusion (MPI)
 #
-
-EXAMPLE_FILE="${EXAMPLE_DIR}/scalar_diffusion/scalar_diffusion.py"
-TEST_DIR='/tmp/hysop_tests/checkpoints/scalar_diffusion_mpi'
-COMMON_OPTIONS="-NC -impl opencl -cp fp32 -d64 --debug-dump-target dump -nu 0.02 -niter 20 -te 0.1 --dump-tstart 0.05 --dump-freq 1 "
-
-echo
-echo "TEST SCALAR DIFFUSION CHECKPOINT (MPI)"
-if [[ ! -f "${EXAMPLE_FILE}" ]]; then
-    echo "Cannot find example file '${EXAMPLE_FILE}'."
-    exit 1
-fi
-
-echo ' Running simulations...'
-${MPIRUN_EXECUTABLE} -np 4 "${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -S "${TEST_DIR}/checkpoint0.tar" --dump-dir "${TEST_DIR}/run0" --checkpoint-dump-time 0.05 --checkpoint-dump-freq 0
-${MPIRUN_EXECUTABLE} -np 4 "${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -S "${TEST_DIR}/checkpoint1.tar" --dump-dir "${TEST_DIR}/run1" --checkpoint-dump-time 0.05 --checkpoint-dump-freq 0
-
-echo ' Comparing solutions...'
-echo "  >debug dumps match"
-compare_files "${TEST_DIR}/run0/dump/run.txt" "${TEST_DIR}/run1/dump/run.txt"
-for f0 in $(find "${TEST_DIR}/run0" -name '*.h5' | sort -n); do
-    f1=$(echo "${f0}" | sed 's/run0/run1/')
-    compare_files "${f0}" "${f1}"
-    echo "  >$(basename ${f0}) match"
-done
-
-echo
-echo ' Running simulations from checkpoints...'
-${MPIRUN_EXECUTABLE} -np 4 "${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -L "${TEST_DIR}/checkpoint0.tar" --dump-dir "${TEST_DIR}/run2"
-${MPIRUN_EXECUTABLE} -np 4 "${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -L "${TEST_DIR}/checkpoint1.tar" --dump-dir "${TEST_DIR}/run3"
-
-echo ' Comparing solutions...'
-compare_files "${TEST_DIR}/run2/dump/run.txt" "${TEST_DIR}/run3/dump/run.txt"
-echo "  >debug dumps match"
-for f0 in $(find "${TEST_DIR}/run2" -name '*.h5' | sort -n); do
-    f1=$(echo "${f0}" | sed 's/run2/run3/')
-    f2=$(echo "${f0}" | sed 's/run2/run0/')
-    f3=$(echo "${f0}" | sed 's/run2/run1/')
-    compare_files "${f0}" "${f1}"
-    compare_files "${f0}" "${f2}"
-    compare_files "${f0}" "${f3}"
-    echo "  >$(basename ${f0}) match"
-done
+#EXAMPLE_FILE="${EXAMPLE_DIR}/scalar_diffusion/scalar_diffusion.py"
+#TEST_DIR='/tmp/hysop_tests/checkpoints/scalar_diffusion_mpi'
+#COMMON_OPTIONS="-NC -impl opencl -cp fp32 -d64 --debug-dump-target dump -nu 0.02 -niter 20 -te 0.1 --dump-tstart 0.05 --dump-freq 1 "
+
+#echo
+#echo "TEST SCALAR DIFFUSION CHECKPOINT (MPI)"
+#if [[ ! -f "${EXAMPLE_FILE}" ]]; then
+    #echo "Cannot find example file '${EXAMPLE_FILE}'."
+    #exit 1
+#fi
+
+#echo ' Running simulations...'
+#${MPIRUN_EXECUTABLE} -np 4 "${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -S "${TEST_DIR}/checkpoint0.tar" --dump-dir "${TEST_DIR}/run0" --checkpoint-dump-time 0.05 --checkpoint-dump-freq 0
+#${MPIRUN_EXECUTABLE} -np 4 "${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -S "${TEST_DIR}/checkpoint1.tar" --dump-dir "${TEST_DIR}/run1" --checkpoint-dump-time 0.05 --checkpoint-dump-freq 0
+
+#echo ' Comparing solutions...'
+#echo "  >debug dumps match"
+#compare_files "${TEST_DIR}/run0/dump/run.txt" "${TEST_DIR}/run1/dump/run.txt"
+#for f0 in $(find "${TEST_DIR}/run0" -name '*.h5' | sort -n); do
+    #f1=$(echo "${f0}" | sed 's/run0/run1/')
+    #compare_files "${f0}" "${f1}"
+    #echo "  >$(basename ${f0}) match"
+#done
+
+#echo
+#echo ' Running simulations from checkpoints...'
+#${MPIRUN_EXECUTABLE} -np 4 "${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -L "${TEST_DIR}/checkpoint0.tar" --dump-dir "${TEST_DIR}/run2"
+#${MPIRUN_EXECUTABLE} -np 4 "${PYTHON_EXECUTABLE}" "${EXAMPLE_FILE}" ${COMMON_OPTIONS} -L "${TEST_DIR}/checkpoint1.tar" --dump-dir "${TEST_DIR}/run3"
+
+#echo ' Comparing solutions...'
+#compare_files "${TEST_DIR}/run2/dump/run.txt" "${TEST_DIR}/run3/dump/run.txt"
+#echo "  >debug dumps match"
+#for f0 in $(find "${TEST_DIR}/run2" -name '*.h5' | sort -n); do
+    #f1=$(echo "${f0}" | sed 's/run2/run3/')
+    #f2=$(echo "${f0}" | sed 's/run2/run0/')
+    #f3=$(echo "${f0}" | sed 's/run2/run1/')
+    #compare_files "${f0}" "${f1}"
+    #compare_files "${f0}" "${f2}"
+    #compare_files "${f0}" "${f3}"
+    #echo "  >$(basename ${f0}) match"
+#done
 
 
 #
diff --git a/hysop/fields/tests/test_cartesian.sh b/hysop/fields/tests/test_cartesian.sh
index 46bad9664..6c840f628 100755
--- a/hysop/fields/tests/test_cartesian.sh
+++ b/hysop/fields/tests/test_cartesian.sh
@@ -1,7 +1,7 @@
 #!/usr/bin/env bash
 set -feu -o pipefail
 PYTHON_EXECUTABLE=${PYTHON_EXECUTABLE:-python3.8}
-MPIRUN_EXECUTABLE=${MPIRUN_EXECUTABLE:-mpirun --allow-run-as-root}
+MPIRUN_EXECUTABLE=${MPIRUN_EXECUTABLE:-mpirun}
 
 SCRIPT_DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null 2>&1 && pwd )"
 TEST_FILE=${SCRIPT_DIR}/test_cartesian.py
diff --git a/hysop/operator/adapt_timestep.py b/hysop/operator/adapt_timestep.py
index 46e5769ad..465cb51aa 100755
--- a/hysop/operator/adapt_timestep.py
+++ b/hysop/operator/adapt_timestep.py
@@ -2,20 +2,21 @@
 Update time-step, depending on the flow field.
 See :ref:`adaptive time_step` for details.
 """
+import numpy as np
 
 from abc import ABCMeta, abstractmethod
 from hysop.constants import HYSOP_REAL, StretchingCriteria, AdvectionCriteria
 from hysop.tools.types import check_instance, first_not_None
 from hysop.tools.decorators import debug
-from hysop.tools.numpywrappers import npw
 from hysop.core.graph.node_generator import ComputationalGraphNodeGenerator
 from hysop.core.graph.computational_operator import ComputationalGraphOperator
 from hysop.core.graph.graph import op_apply
 from hysop.fields.continuous_field import Field
 from hysop.parameters import ScalarParameter, TensorParameter
 from hysop.core.mpi import MPI
+from hysop.backend.host.host_operator import HostOperatorBase
 
-class TimestepCriteria(ComputationalGraphOperator, metaclass=ABCMeta):
+class TimestepCriteria(HostOperatorBase, metaclass=ABCMeta):
 
     @debug
     def __new__(cls, parameter, input_params, output_params,
@@ -63,7 +64,9 @@ class TimestepCriteria(ComputationalGraphOperator, metaclass=ABCMeta):
         self.min_dt   = 0.0 if (min_dt   is None) else min_dt
         self.max_dt   = 1e8 if (max_dt   is None) else max_dt
         self.dt_coeff = 1.0 if (dt_coeff is None) else dt_coeff
-        self.dt       = parameter
+
+        dt = parameter
+        self.dt = dt
 
         # Collect values from all MPI process
         if self.mpi_params.size == 1:
@@ -71,8 +74,8 @@ class TimestepCriteria(ComputationalGraphOperator, metaclass=ABCMeta):
             self._collect_max = lambda e: e
         else:
             comm = self.mpi_params.comm
-            self._sendbuff = npw.zeros((1, ))
-            self._recvbuff = npw.zeros((1, ))
+            self._sendbuff = np.zeros((1, ), dtype=dt.dtype)
+            self._recvbuff = np.zeros((1, ), dtype=dt.dtype)
             def _collect_max(val):
                 self._sendbuff[0] = val
                 comm.Allreduce(self._sendbuff, self._recvbuff, op=MPI.MAX)
@@ -88,8 +91,8 @@ class TimestepCriteria(ComputationalGraphOperator, metaclass=ABCMeta):
     def apply(self, **kwds):
         dt = self.compute_criteria(**kwds)
         dt *= self.dt_coeff
-        dt = self._collect_max(npw.maximum(dt, self.min_dt))
-        dt = self._collect_min(npw.minimum(dt, self.max_dt))
+        dt = self._collect_max(np.maximum(dt, self.min_dt, dtype=dt.dtype))
+        dt = self._collect_min(np.minimum(dt, self.max_dt, dtype=dt.dtype))
         assert (dt > 0.0), 'negative or zero timestep encountered.'
         self.dt.set_value(dt)
 
@@ -137,7 +140,7 @@ class ConstantTimestepCriteria(TimestepCriteria):
         if isinstance(cst, int):
             cst = float(cst)
         assert (cst > 0.0), 'negative cst factor.'
-        check_instance(cst, (float, npw.ndarray, list, tuple))
+        check_instance(cst, (float, np.ndarray, list, tuple))
         check_instance(Finf, (ScalarParameter, TensorParameter))
         if isinstance(Finf, ScalarParameter):
             assert isinstance(cst, float)
@@ -145,10 +148,10 @@ class ConstantTimestepCriteria(TimestepCriteria):
         else:
             is_scalar = False
             if isinstance(cst, float):
-                cst = npw.full(shape=Finf.shape, dtype=Finf.dtype, fill_value=cst)
+                cst = np.full(shape=Finf.shape, dtype=Finf.dtype, fill_value=cst)
             if isinstance(cst, (list, tuple)):
                 assert Finf.ndim == 1
-                cst = npw.asarray(cst)
+                cst = np.asarray(cst)
             msg='Shape mismatch between parameter {} and cst {}.'
             msg=msg.format(Finf.shape, cst.shape)
             assert Finf.shape == cst.shape, msg
@@ -170,13 +173,13 @@ class ConstantTimestepCriteria(TimestepCriteria):
         if self.is_scalar:
             assert Finf >= 0
             if (Finf == 0):
-                return npw.inf
+                return np.inf
             else:
                 return cst/Finf
         else:
             assert Finf.min() >= 0
             mask = (Finf!=0)
-            dt = npw.full_like(cst, fill_value=npw.inf)
+            dt = np.full_like(cst, fill_value=np.inf)
             dt[mask] = cst[mask] / Finf[mask]
             return dt.min()
 
@@ -267,7 +270,7 @@ class CflTimestepCriteria(TimestepCriteria):
 
         rv = ()
         for Vr in relative_velocities:
-            Vr = npw.asarray(Vr, dtype=dtype)
+            Vr = np.asarray(Vr, dtype=dtype)
             assert Vr.shape == shape
             rv += (Vr,)
         relative_velocities = rv
@@ -305,15 +308,15 @@ class CflTimestepCriteria(TimestepCriteria):
         assert len(dx) == Fmin.size == Fmax.size
         assert len(self.relative_velocities)>=1
 
-        dt = npw.inf
+        dt = np.inf
         for Vr in self.relative_velocities:
             Vmin = Fmin - Vr
             Vmax = Fmax - Vr
-            Vinf = npw.maximum(npw.abs(Vmin), npw.abs(Vmax))
-            if npw.all(npw.divide(Vinf, dx)==0):
-                cdt = cfl*npw.inf
+            Vinf = np.maximum(np.abs(Vmin), np.abs(Vmax))
+            if np.all(np.divide(Vinf, dx)==0):
+                cdt = cfl*np.inf
             else:
-                cdt = cfl / npw.max(npw.divide(Vinf, dx))
+                cdt = cfl / np.max(np.divide(Vinf, dx))
             dt = min(dt, cdt)
         return dt
 
@@ -403,19 +406,19 @@ class AdvectionTimestepCriteria(TimestepCriteria):
         lcfl = self.lcfl
         if (criteria is AdvectionCriteria.W_INF):
             Finf = self.Finf()
-            if npw.max(Finf)==0:
-                return lcfl*npw.inf
+            if np.max(Finf)==0:
+                return lcfl*np.inf
             else:
-                return lcfl / npw.max(Finf)
+                return lcfl / np.max(Finf)
         elif (criteria is AdvectionCriteria.GRAD_U):
             gradFinf = self.gradFinf()
             if (gradFinf.ndim == 2):
-                gradFinf = npw.diag(gradFinf) #extract diagonal
-            return lcfl / npw.max(gradFinf)
+                gradFinf = np.diag(gradFinf) #extract diagonal
+            return lcfl / np.max(gradFinf)
         elif (criteria is AdvectionCriteria.DEFORMATION):
             gradFinf = self.gradFinf()
             gradFinf = (gradFinf + gradFinf.T)/2.0
-            return lcfl / npw.max( gradFinf.sum(axis=0) )
+            return lcfl / np.max( gradFinf.sum(axis=0) )
         else:
             msg='Unsupported stretching criteria {}.'.format(criteria)
             raise RuntimeError(msg)
@@ -445,7 +448,7 @@ class StretchingTimestepCriteria(TimestepCriteria):
 
         where |dFi/dXj| = |gradF|_inf_ij
 
-        ie. dt = cst / npw.max( gradFinf.sum(axis=1) )
+        ie. dt = cst / np.max( gradFinf.sum(axis=1) )
 
         Parameters
         ----------
@@ -482,11 +485,12 @@ class StretchingTimestepCriteria(TimestepCriteria):
         criteria = self.criteria
         if (criteria is StretchingCriteria.GRAD_U):
             gradFinf = self.gradFinf()
-            return self.cst / npw.max( gradFinf.sum(axis=1) )
+            return self.cst / np.max( gradFinf.sum(axis=1) )
         else:
             msg='Unsupported stretching criteria {}.'.format(criteria)
             raise RuntimeError(msg)
 
+
 class MergeTimeStepCriterias(TimestepCriteria):
 
     @debug
@@ -529,7 +533,7 @@ class MergeTimeStepCriterias(TimestepCriteria):
     @debug
     def apply(self, simulation, **kwds):
         assert simulation.dt is self.dt, 'Parameter mismatch between Simulation and AdaptiveTimeStep.'
-        if self._start_time is None or simulation.t() > self._start_time:
+        if (self._start_time is None) or (simulation.t() > self._start_time):
             super(MergeTimeStepCriterias, self).apply(simulation=simulation, **kwds)
 
 
diff --git a/hysop/operator/advection.py b/hysop/operator/advection.py
index f165a25c0..e2c29b5b6 100644
--- a/hysop/operator/advection.py
+++ b/hysop/operator/advection.py
@@ -3,7 +3,7 @@
 Advection operator generator.
 """
 from hysop.constants import Implementation
-from hysop.tools.types import check_instance, to_list
+from hysop.tools.types import check_instance, to_list, first_not_None
 from hysop.tools.decorators import debug
 from hysop.fields.continuous_field import Field
 from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
diff --git a/hysop/operator/base/external_force.py b/hysop/operator/base/external_force.py
index 0fef2bb95..5435e8659 100644
--- a/hysop/operator/base/external_force.py
+++ b/hysop/operator/base/external_force.py
@@ -11,6 +11,7 @@ from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
 from hysop.fields.continuous_field import Field
 from hysop.symbolic.relational import Assignment
 from hysop.core.memory.memory_request import MemoryRequest
+from hysop.core.graph.graph import op_apply
 from hysop.parameters.tensor_parameter import TensorParameter
 from hysop.parameters.scalar_parameter import ScalarParameter
 from hysop.tools.interface import NamedObjectI
@@ -267,6 +268,7 @@ class SpectralExternalForceOperatorBase(SpectralOperatorBase):
         super(SpectralExternalForceOperatorBase, self).setup(work)
         self.Fext.post_setup(self, work)
 
+    @op_apply
     def apply(self, **kwds):
         self.Fext.apply(self, **kwds)
 
diff --git a/hysop/operator/base/memory_reordering.py b/hysop/operator/base/memory_reordering.py
index cbcf5d9cf..bdba3fb55 100644
--- a/hysop/operator/base/memory_reordering.py
+++ b/hysop/operator/base/memory_reordering.py
@@ -18,9 +18,9 @@ class MemoryReorderingBase(object, metaclass=ABCMeta):
                     target_memory_order, name=None, pretty_name=None,
                     **kwds):
         return super(MemoryReorderingBase, cls).__new__(cls,
-                input_fields=input_fields,
-                output_fields=output_fields,
-                name=name, pretty_name=pname,
+                input_fields=None,
+                output_fields=None,
+                name=name, pretty_name=pretty_name,
                 **kwds)
 
     @debug
diff --git a/hysop/operator/base/min_max.py b/hysop/operator/base/min_max.py
index 7a1faeedd..ff8bc693e 100644
--- a/hysop/operator/base/min_max.py
+++ b/hysop/operator/base/min_max.py
@@ -3,11 +3,11 @@
 MinMaxFieldStatisticsBase: compute min(F), max(F) and/or max(|F|) for a given field
 MinMaxGradientStatisticsBase: compute min(gradF), max(gradF) and/or max(|gradF|) for a given field, direction-wise.
 """
+import numpy as np
 from abc import abstractmethod
 from hysop.tools.types       import check_instance, first_not_None, to_tuple
 from hysop.tools.enum        import EnumFactory
 from hysop.tools.decorators  import debug
-from hysop.tools.numpywrappers import npw
 from hysop.fields.continuous_field import Field
 from hysop.constants import Backend
 from hysop.parameters.tensor_parameter import TensorParameter
@@ -77,7 +77,7 @@ class MinMaxFieldStatisticsBase(object):
             else:
                 param = None
             if (param is not None):
-                assert npw.prod(param.shape) == nb_components
+                assert np.prod(param.shape) == nb_components
             parameters[k] = param
         return parameters
 
@@ -174,7 +174,7 @@ class MinMaxFieldStatisticsBase(object):
         check_instance(field, Field)
         check_instance(components, tuple, values=int,
                 allow_none=True, minval=0, maxval=field.nb_components-1)
-        check_instance(coeffs, dict, keys=str, values=(int, float, npw.number), allow_none=True)
+        check_instance(coeffs, dict, keys=str, values=(int, float, np.number), allow_none=True)
         check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors,
                 allow_none=True)
         check_instance(name, str, allow_none=True)
@@ -241,8 +241,8 @@ class MinMaxFieldStatisticsBase(object):
             comm = self.mpi_params.comm
             Fmin, Fmax = self.Fmin, self.Fmax
             if (self.Fmax is not None):
-                sendbuff = npw.zeros_like(Fmax.value)
-                recvbuff = npw.zeros_like(Fmax.value)
+                sendbuff = np.zeros_like(Fmax.tensor_value)
+                recvbuff = np.zeros_like(Fmax.tensor_value)
                 def _collect_max(val, sendbuff=sendbuff, recvbuff=recvbuff):
                     sendbuff[...] = val
                     comm.Allreduce(sendbuff, recvbuff, op=MPI.MAX)
@@ -250,8 +250,8 @@ class MinMaxFieldStatisticsBase(object):
             else:
                 _collect_max = None
             if (self.Fmin is not None):
-                sendbuff = npw.zeros_like(Fmin.value)
-                recvbuff = npw.zeros_like(Fmin.value)
+                sendbuff = np.zeros_like(Fmin.tensor_value)
+                recvbuff = np.zeros_like(Fmin.tensor_value)
                 def _collect_min(val, sendbuff=sendbuff, recvbuff=recvbuff):
                     sendbuff[...] = val
                     comm.Allreduce(sendbuff, recvbuff, op=MPI.MIN)
@@ -265,6 +265,13 @@ class MinMaxFieldStatisticsBase(object):
         """Backend agnostic computation of min and max parameters."""
         dfield, components, coeffs = self._dfield, self._components, self._coeffs
         Fmin, Fmax, Finf = self.Fmin, self.Fmax, self.Finf
+
+        # For now field.min/field.max take into account ghosts...
+        # This is because pyopencl.reduction does not support reductions on array views.
+        # Here we force synchronize ghosts in every direction, including diagonals,
+        #  in order to overwrite potential bad boundary values.
+        dfield.exchange_ghosts()
+
         if (Fmin is not None):
             fmin = Fmin.tensor_value
             for i in components:
@@ -276,7 +283,7 @@ class MinMaxFieldStatisticsBase(object):
                 fmax[i] = dfield.data[i].max().get()
             Fmax.value = self._collect_max(fmax * coeffs['Fmax'])
         if (Finf is not None):
-            self.Finf.value = npw.maximum(npw.abs(Fmin()), npw.abs(Fmax())) * coeffs['Finf']
+            self.Finf.value = np.maximum(np.abs(Fmin()), np.abs(Fmax())) * coeffs['Finf']
 
 
 class MinMaxDerivativeStatisticsBase(MinMaxFieldStatisticsBase):
diff --git a/hysop/operator/directional/directional.py b/hysop/operator/directional/directional.py
index d74fa8c83..5c09b05d3 100644
--- a/hysop/operator/directional/directional.py
+++ b/hysop/operator/directional/directional.py
@@ -267,10 +267,11 @@ class DirectionalOperatorGenerator(DirectionalOperatorGeneratorI, metaclass=ABCM
             sargs = ['*{} = {}'.format(k,v.__class__)
                         for (k,v) in kargs.items()]
             msg  = 'FATAL ERROR during {}.generate():\n'
-            msg += ' => failed to call {}.__init__()\n    with the following keywords:'
+            msg += ' => failed to initialize an instance of type {}'
+            msg += '\n    by using the following keyword arguments:'
             msg += '\n     '+'\n     '.join(sargs)
-            msg = msg.format(self.__class__, self._operator)
-            print('\n{}\n'.format(msg))
+            msg = msg.format(self.__class__, self.impl)
+            print('\n{}'.format(msg))
             raise
 
         if (self._input_fields_to_dump is not None):
diff --git a/hysop/operator/hdf_io.py b/hysop/operator/hdf_io.py
index ac0494264..c04a0b82d 100755
--- a/hysop/operator/hdf_io.py
+++ b/hysop/operator/hdf_io.py
@@ -27,14 +27,14 @@ from hysop.tools.types import check_instance, first_not_None
 from hysop.tools.numpywrappers import npw
 from hysop.tools.io_utils import IO, IOParams, XMF
 from hysop.core.graph.graph import op_apply
-from hysop.core.graph.computational_graph import ComputationalGraphOperator
+from hysop.backend.host.host_operator import HostOperatorBase
 from hysop.fields.continuous_field import Field
 from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
 from hysop.core.memory.memory_request import MemoryRequest
 from hysop.topology.topology_descriptor import TopologyDescriptor
 
 
-class HDF_IO(ComputationalGraphOperator, metaclass=ABCMeta):
+class HDF_IO(HostOperatorBase, metaclass=ABCMeta):
     """
     Abstract interface to read/write from/to hdf files, for
     hysop fields.
@@ -434,9 +434,11 @@ class HDF_Writer(HDF_IO):
         topo = self.topology
         dim = topo.domain.dim
         dx = list(topo.mesh.space_step)
+
         mesh = self.refmesh
         res = list(mesh.grid_resolution)
         orig = list(topo.domain.origin)
+
         resolution = [1, ]*3
         origin = [0.0, ]*3
         step = [0.0, ]*3
diff --git a/hysop/operator/mean_field.py b/hysop/operator/mean_field.py
index fcf0105ce..c71c55c6b 100644
--- a/hysop/operator/mean_field.py
+++ b/hysop/operator/mean_field.py
@@ -18,9 +18,9 @@ from hysop.core.graph.graph import op_apply
 from hysop.core.graph.computational_graph import ComputationalGraphOperator
 from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
 from hysop.fields.continuous_field import Field
+from hysop.backend.host.host_operator import HostOperatorBase
 
-
-class ComputeMeanField(ComputationalGraphOperator):
+class ComputeMeanField(HostOperatorBase):
     """
     Interface to compute the mean of a field in chosen axes.
     """
diff --git a/hysop/operator/misc.py b/hysop/operator/misc.py
index 628c4bdb7..73a2158ff 100644
--- a/hysop/operator/misc.py
+++ b/hysop/operator/misc.py
@@ -94,7 +94,7 @@ class ForceTopologyState(Noop):
 
         extra_kwds.setdefault('mpi_params', mpi_params)
         extra_kwds.setdefault('cl_env',  cl_env)
-        kwds.setdefault('cl_env',  cl_env)
+        #kwds.setdefault('cl_env',  cl_env)
         kwds.setdefault('mpi_params', mpi_params)
 
         super(ForceTopologyState, self).__init__(input_fields=input_fields,
diff --git a/hysop/operator/parameter_plotter.py b/hysop/operator/parameter_plotter.py
index 53b33b6c4..1b510d2fd 100644
--- a/hysop/operator/parameter_plotter.py
+++ b/hysop/operator/parameter_plotter.py
@@ -2,12 +2,13 @@ from abc import abstractmethod
 from hysop.tools.types import to_tuple, check_instance, first_not_None
 from hysop.tools.numpywrappers import npw
 from hysop.tools.io_utils import IO
+from hysop.core.graph.graph import op_apply
 from hysop.core.graph.computational_graph import ComputationalGraphOperator
 from hysop.parameters.scalar_parameter import ScalarParameter
 from hysop.parameters.tensor_parameter import TensorParameter
+from hysop.backend.host.host_operator import HostOperatorBase
 
-
-class PlottingOperator(ComputationalGraphOperator):
+class PlottingOperator(HostOperatorBase):
     """
     Base operator for plotting.
     """
@@ -77,6 +78,7 @@ class PlottingOperator(ComputationalGraphOperator):
         self.fig.show()
         self.plt.pause(0.001)
 
+    @op_apply
     def apply(self, **kwds):
         self._update(**kwds)
         self._save(**kwds)
diff --git a/hysop/operator/tests/test_bilevel_advection.py b/hysop/operator/tests/test_bilevel_advection.py
index 45777b9e2..9854135c3 100644
--- a/hysop/operator/tests/test_bilevel_advection.py
+++ b/hysop/operator/tests/test_bilevel_advection.py
@@ -61,11 +61,11 @@ class TestBilevelAdvectionOperator(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=1, register_object=False)
+                    nb_components=1)
             Sout = Field(domain=domain, name='Sout', dtype=dtype,
-                    nb_components=1, register_object=False)
+                    nb_components=1)
             for time_integrator in time_integrators:
                 for remesh_kernel in remesh_kernels:
                     for velocity_cfl in velocity_cfls:
diff --git a/hysop/operator/tests/test_scales_advection.py b/hysop/operator/tests/test_scales_advection.py
index e3430a8e6..62e6d083f 100644
--- a/hysop/operator/tests/test_scales_advection.py
+++ b/hysop/operator/tests/test_scales_advection.py
@@ -68,11 +68,11 @@ class TestScalesAdvectionOperator(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=5, register_object=False)
+                        nb_components=5)
             Sout = Field(domain=domain, name='Sout', dtype=dtype,
-                         nb_components=5, register_object=False)
+                         nb_components=5)
             for time_integrator in time_integrators:
                 for remesh_kernel in remesh_kernels:
                     for velocity_cfl in velocity_cfls:
diff --git a/hysop/tools/debug_dumper.py b/hysop/tools/debug_dumper.py
index bf94fd838..e9ff4e467 100644
--- a/hysop/tools/debug_dumper.py
+++ b/hysop/tools/debug_dumper.py
@@ -14,17 +14,21 @@ from hysop.fields.discrete_field import DiscreteScalarFieldView
 class DebugDumper(object):
     def __init__(self, name, path, force_overwrite=False,
                  enable_on_op_apply=False, dump_precision=10,
-                 comm=MPI.COMM_WORLD, io_leader=0):
+                 comm=MPI.COMM_WORLD, io_leader=0, dump_data=False):
         assert isinstance(name, str), name
         assert isinstance(path, str), path
+
         directory = os.path.join(path, name)
+
         blobs_directory = os.path.join(directory, 'data')
-        if not os.path.isdir(blobs_directory) and comm.rank==0:
-            os.makedirs(blobs_directory)
+        if dump_data:
+            if not os.path.isdir(blobs_directory) and comm.rank==0:
+                os.makedirs(blobs_directory)
+        self.blobs_directory = blobs_directory
+        self.dump_data = dump_data
 
         self.name = name
         self.directory = directory
-        self.blobs_directory = blobs_directory
         self.dump_id = 0
         self.enable_on_op_apply = enable_on_op_apply
         self.dump_precision = dump_precision
@@ -113,18 +117,17 @@ class DebugDumper(object):
                 dtype = d.dtype
                 shape = None
                 id_ = self.dump_id
-                _min = comm.allreduce(float(np.nanmin(d)),  op=MPI.MIN)
-                _max = comm.allreduce(float(np.nanmax(d)),  op=MPI.MAX)
-                mean = comm.allreduce(float(np.nanmean(d))) / comm_size
-                variance = comm.allreduce(float(np.nansum((d-mean)**2))) / \
-                    float(comm.allreduce(int(d.size)))
+                _min = comm.allreduce(np.nanmin(d),  op=MPI.MIN)
+                _max = comm.allreduce(np.nanmax(d),  op=MPI.MAX)
+                mean = comm.allreduce(np.nanmean(d), op=MPI.SUM) / comm_size
+                variance = comm.allreduce(np.nansum((d-mean)**2)) / comm.allreduce(d.size)
                 entry = '\n'+self.lformat(id_, iteration, time, tag_, _min, _max,
                                           mean, variance, dtype, shape, description_, self.dump_precision)
 
                 if self.is_io_leader:
                     self.runfile.write(entry)
 
-                if (comm_size == 1):
+                if self.dump_data and (comm_size == 1):
                     dst = '{}/{}'.format(self.blobs_directory, self.dump_id)
                     np.savez_compressed(dst, data=d)
 
diff --git a/hysop/tools/handle.py b/hysop/tools/handle.py
index ee553ede7..48cedf0ec 100644
--- a/hysop/tools/handle.py
+++ b/hysop/tools/handle.py
@@ -96,7 +96,18 @@ class TaggedObject(object, metaclass=ABCMeta):
         assert (tag_formatter is None) or callable(tag_formatter)
         tagged_cls = first_not_None(tagged_cls, cls)
 
-        obj = super(TaggedObject, cls).__new__(cls, **kwds)
+        try:
+            obj = super(TaggedObject, cls).__new__(cls, **kwds)
+        except TypeError:
+            msg = '\nFATAL ERROR during {}.__new__(cls, **kwds).'.format(cls.__name__)
+            msg+= '\nThis may be due to extra keyword arguments:\n  *'
+            msg+= '\n  *'.join(kwds.keys())
+            msg+= '\nIf you believe that those arguments are valid, check the following types: \n  *'
+            msg+= '\n  *'.join(map(str, cls.__mro__[:-1]))
+            msg+= '\n'
+            print(msg)
+            raise
+
         if tagged_cls in TaggedObject.__ids:
             obj.__tag_id = TaggedObject.__ids[tagged_cls]
             TaggedObject.__ids[tagged_cls] += 1
@@ -119,7 +130,17 @@ class TaggedObject(object, metaclass=ABCMeta):
         assert (tag_postfix is None) or isinstance(tag_postfix, str)
         assert (tag_formatter is None) or callable(tag_formatter)
 
-        super(TaggedObject, self).__init__(**kwds)
+        try:
+            super(TaggedObject, self).__init__(**kwds)
+        except TypeError:
+            cls = self.__class__
+            msg = '\nFATAL ERROR during {}.__init__(**kwds).'.format(cls.__name__)
+            msg+= '\nThis may be due to extra keyword arguments:\n  *'
+            msg+= '\n  *'.join(kwds.keys())
+            msg+= '\nIf you believe that those arguments are valid, check the following types: \n  *'
+            msg+= '\n  *'.join(map(str, cls.__mro__[:-1]))
+            print(msg)
+            raise
 
         # reaffect attributes (some classes use only __init__ for simplicity)
         self.__tag_prefix    = first_not_None(self.__tag_prefix, tag_prefix, '')
diff --git a/hysop_examples/argparser.py b/hysop_examples/argparser.py
index 41dc24d53..92913a4cb 100644
--- a/hysop_examples/argparser.py
+++ b/hysop_examples/argparser.py
@@ -1247,7 +1247,7 @@ class HysopArgParser(argparse.ArgumentParser):
                 ndumps = int(np.floor(T/dt)) + 1
                 toi = tstart + np.arange(ndumps)*dt
                 dump_times.update(toi)
-            dump_times = filter(lambda t: (t>=tstart) & (t<=tend), dump_times)
+            dump_times = tuple(filter(lambda t: (t>=tstart) & (t<=tend), dump_times))
 
             setattr(args, '{}_times'.format(bname), tuple(sorted(dump_times)))
             times_of_interest.update(dump_times)
diff --git a/hysop_examples/examples/cylinder/oscillating_cylinder.py b/hysop_examples/examples/cylinder/oscillating_cylinder.py
index 38392623f..37eb902e4 100644
--- a/hysop_examples/examples/cylinder/oscillating_cylinder.py
+++ b/hysop_examples/examples/cylinder/oscillating_cylinder.py
@@ -149,8 +149,8 @@ def compute(args):
     #> Directional stretching + diffusion
     if (dim==3):
         stretch = DirectionalStretching(implementation=impl,
-                 name='stretch',
-                 pretty_name='stretch',
+                 name='S',
+                 pretty_name='S',
                  formulation = args.stretching_formulation,
                  velocity  = velo,
                  vorticity = vorti,
diff --git a/hysop_examples/examples/flow_around_sphere/flow_around_sphere.py b/hysop_examples/examples/flow_around_sphere/flow_around_sphere.py
index b68cd5681..f406b5657 100644
--- a/hysop_examples/examples/flow_around_sphere/flow_around_sphere.py
+++ b/hysop_examples/examples/flow_around_sphere/flow_around_sphere.py
@@ -183,7 +183,7 @@ def compute(args):
         StretchOp = StaticDirectionalStretching
     stretch = StretchOp(
         implementation=Implementation.PYTHON if implIsFortran else impl,
-        name='stretch',
+        name='S',
         formulation=StretchingFormulation.CONSERVATIVE,
         velocity=velo,
         vorticity=vorti,
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 5ae967ba2..1837f8924 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
@@ -165,8 +165,8 @@ def compute(args):
     #> Stretch vorticity
     if (dim==3):
         stretch = DirectionalStretching(implementation=impl,
-                 name='stretch',
-                 pretty_name='stretch',
+                 name='S',
+                 pretty_name='S',
                  formulation = args.stretching_formulation,
                  velocity  = velo,
                  vorticity = vorti,
diff --git a/hysop_examples/examples/particles_above_salt/particles_above_salt_bc_3d.py b/hysop_examples/examples/particles_above_salt/particles_above_salt_bc_3d.py
index 8fad5c0be..cdd314b74 100644
--- a/hysop_examples/examples/particles_above_salt/particles_above_salt_bc_3d.py
+++ b/hysop_examples/examples/particles_above_salt/particles_above_salt_bc_3d.py
@@ -179,8 +179,8 @@ def compute(args):
     #> Stretch vorticity
     if (dim==3):
         stretch = DirectionalStretching(implementation=impl,
-                 name='stretch',
-                 pretty_name='stretch',
+                 name='S',
+                 pretty_name='S',
                  formulation = args.stretching_formulation,
                  velocity  = velo,
                  vorticity = vorti,
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 7c3e8942c..b481f4154 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
@@ -161,8 +161,8 @@ def compute(args):
     #> Stretch vorticity
     if (dim==3):
         stretch = DirectionalStretching(implementation=impl,
-                 name='stretch',
-                 pretty_name='stretch',
+                 name='S',
+                 pretty_name='S',
                  formulation = args.stretching_formulation,
                  velocity  = velo,
                  vorticity = vorti,
diff --git a/hysop_examples/examples/scalar_advection/turbulent_scalar_advection.py b/hysop_examples/examples/scalar_advection/turbulent_scalar_advection.py
index 22f4ffb7c..3c1f91519 100644
--- a/hysop_examples/examples/scalar_advection/turbulent_scalar_advection.py
+++ b/hysop_examples/examples/scalar_advection/turbulent_scalar_advection.py
@@ -143,7 +143,7 @@ advec = Advection(
 #> Directional stretching
 stretch = StaticDirectionalStretching(
     implementation=Implementation.PYTHON,
-    name='stretch',
+    name='S',
     formulation=StretchingFormulation.CONSERVATIVE,
     velocity=velo,
     vorticity=vorti,
diff --git a/hysop_examples/examples/sediment_deposit/sediment_deposit.py b/hysop_examples/examples/sediment_deposit/sediment_deposit.py
index 3379c68a5..4e88960a8 100644
--- a/hysop_examples/examples/sediment_deposit/sediment_deposit.py
+++ b/hysop_examples/examples/sediment_deposit/sediment_deposit.py
@@ -200,8 +200,8 @@ def compute(args):
     #> Stretch vorticity
     if (dim==3):
         stretch = DirectionalStretching(implementation=impl,
-                 name='stretch',
-                 pretty_name='stretch',
+                 name='S',
+                 pretty_name='S',
                  formulation = args.stretching_formulation,
                  velocity  = velo,
                  vorticity = vorti,
diff --git a/hysop_examples/examples/sediment_deposit/sediment_deposit_levelset.py b/hysop_examples/examples/sediment_deposit/sediment_deposit_levelset.py
index e0f376654..099577c90 100644
--- a/hysop_examples/examples/sediment_deposit/sediment_deposit_levelset.py
+++ b/hysop_examples/examples/sediment_deposit/sediment_deposit_levelset.py
@@ -254,8 +254,8 @@ def compute(args):
     #> Stretch vorticity
     if (dim==3):
         stretch = DirectionalStretching(implementation=impl,
-                 name='stretch',
-                 pretty_name='stretch',
+                 name='S',
+                 pretty_name='S',
                  formulation = args.stretching_formulation,
                  velocity  = velo,
                  vorticity = vorti,
diff --git a/hysop_examples/examples/taylor_green/taylor_green.py b/hysop_examples/examples/taylor_green/taylor_green.py
index 0de08e7fe..39e4a016b 100644
--- a/hysop_examples/examples/taylor_green/taylor_green.py
+++ b/hysop_examples/examples/taylor_green/taylor_green.py
@@ -95,7 +95,6 @@ def compute(args):
                 name='advec',
                 velocity = velo,
                 advected_fields = (vorti,),
-                velocity_cfl = args.cfl,
                 variables = {velo: npts, vorti: npts},
                 dt=dt, **extra_op_kwds)
         advec_dir = None
@@ -112,7 +111,7 @@ def compute(args):
     #> Directional stretching
     if (impl is Implementation.PYTHON) or (impl is Implementation.FORTRAN):
         stretch_dir = StaticDirectionalStretching(implementation=Implementation.PYTHON,
-                 name='stretch',
+                 name='S',
                  formulation = args.stretching_formulation,
                  velocity  = velo,
                  vorticity = vorti,
@@ -120,7 +119,7 @@ def compute(args):
                  dt=dt, **extra_op_kwds)
     else:
         stretch_dir = DirectionalStretching(implementation=impl,
-                 name='stretch',
+                 name='S',
                  formulation = args.stretching_formulation,
                  velocity  = velo,
                  vorticity = vorti,
diff --git a/hysop_examples/examples/taylor_green/taylor_green_cpuFortran.py b/hysop_examples/examples/taylor_green/taylor_green_cpuFortran.py
index 4d8de7b01..7f4b07c5c 100644
--- a/hysop_examples/examples/taylor_green/taylor_green_cpuFortran.py
+++ b/hysop_examples/examples/taylor_green/taylor_green_cpuFortran.py
@@ -76,7 +76,7 @@ def compute(args):
             dt=dt, **extra_op_kwds)
     #> Directional stretching
     stretch = StaticDirectionalStretching(implementation=impl,
-             name='stretch',
+             name='S',
              formulation = args.stretching_formulation,
              velocity  = velo,
              vorticity = vorti,
-- 
GitLab