From c0e6d37e991eaa5cead9709e831d9e651cea6bed Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Tue, 5 Apr 2016 16:53:38 +0200
Subject: [PATCH] fix/add some tests.

---
 HySoP/hysop/__init__.py.in                    |   5 +-
 HySoP/hysop/constants.py                      |   7 -
 HySoP/hysop/domain/tests/test_subset.py       |  16 +-
 HySoP/hysop/f2py/fftw2py.pyf                  |   4 +-
 HySoP/hysop/fields/continuous.py              | 164 +++++++----
 HySoP/hysop/fields/discrete.py                |  24 +-
 HySoP/hysop/fields/tests/func_for_tests.py    |  39 ++-
 HySoP/hysop/fields/tests/test_field.py        | 193 +++++++++----
 .../tests/test_gpu_multiresolution_filter.py  |   8 +-
 .../gpu/tests/test_multiphase_baroclinic.py   |  12 +-
 HySoP/hysop/numerics/integrators/odesolver.py |   1 +
 .../hysop/numerics/tests/test_integrators.py  |   0
 .../numerics/tests/test_update_ghosts.py      |   2 +-
 HySoP/hysop/operator/computational.py         |  39 +--
 HySoP/hysop/operator/continuous.py            |  48 ++--
 HySoP/hysop/operator/discrete/discrete.py     |  10 +-
 HySoP/hysop/operator/discrete/spectrum.py     |   0
 HySoP/hysop/operator/drag_and_lift.py         |  12 +-
 HySoP/hysop/operator/hdf_io.py                | 152 ++++++-----
 HySoP/hysop/operator/tests/test_analytic.py   |  32 +--
 .../operator/tests/test_drag_and_lift.py      |  21 +-
 HySoP/hysop/operator/tests/test_hdf5_io.py    |  68 +++--
 .../operator/tests/test_multiphase_gradp.py   |   6 +-
 .../tests/test_multiresolution_filter.py      |   6 +-
 .../operator/tests/test_particle_advection.py |  13 +-
 .../hysop/operator/tests/test_penalization.py |   3 +-
 HySoP/hysop/operator/tests/test_spectrum.py   |   5 +-
 HySoP/hysop/problem/tests/test_simulation.py  |   9 +-
 HySoP/hysop/tools/io_utils.py                 | 258 +++++++++++++-----
 HySoP/hysop/tools/parameters.py               |  69 +----
 HySoP/hysop/tools/tests/test_parameters.py    |   6 +-
 HySoP/src/fftw/fft3d.f90                      |   4 +-
 32 files changed, 787 insertions(+), 449 deletions(-)
 mode change 100644 => 100755 HySoP/hysop/constants.py
 mode change 100644 => 100755 HySoP/hysop/fields/continuous.py
 mode change 100644 => 100755 HySoP/hysop/fields/discrete.py
 mode change 100644 => 100755 HySoP/hysop/fields/tests/func_for_tests.py
 mode change 100644 => 100755 HySoP/hysop/numerics/integrators/odesolver.py
 mode change 100644 => 100755 HySoP/hysop/numerics/tests/test_integrators.py
 mode change 100644 => 100755 HySoP/hysop/operator/computational.py
 mode change 100644 => 100755 HySoP/hysop/operator/continuous.py
 mode change 100644 => 100755 HySoP/hysop/operator/discrete/discrete.py
 mode change 100644 => 100755 HySoP/hysop/operator/discrete/spectrum.py
 mode change 100644 => 100755 HySoP/hysop/operator/drag_and_lift.py
 mode change 100644 => 100755 HySoP/hysop/operator/hdf_io.py
 mode change 100644 => 100755 HySoP/hysop/operator/tests/test_spectrum.py
 mode change 100644 => 100755 HySoP/hysop/problem/tests/test_simulation.py
 mode change 100644 => 100755 HySoP/hysop/tools/io_utils.py
 mode change 100644 => 100755 HySoP/hysop/tools/parameters.py

diff --git a/HySoP/hysop/__init__.py.in b/HySoP/hysop/__init__.py.in
index e97188c94..4e2d010be 100755
--- a/HySoP/hysop/__init__.py.in
+++ b/HySoP/hysop/__init__.py.in
@@ -13,6 +13,7 @@ __DEBUG__ = "@DEBUG@" in ["2", "3"]
 __PROFILE__ = "@PROFILE@" in ["0", "1"]
 __OPTIMIZE__ = "@OPTIM@" is "ON"
 
+import os
 from hysop.tools.sys_utils import SysUtils
 # Box-type physical domain
 from hysop.domain.box import Box
@@ -23,8 +24,8 @@ from hysop.fields.variable_parameter import VariableParameter
 # Simulation parameters
 from hysop.problem.simulation import Simulation
 # Tools (io, mpi ...)
-from hysop.tools.io_utils import IO
-from hysop.tools.parameters import IOParams, MPIParams, Discretization
+from hysop.tools.io_utils import IO, IOParams
+from hysop.tools.parameters import MPIParams, Discretization
 import hysop.mpi
 # Problem
 from hysop.problem.problem import Problem
diff --git a/HySoP/hysop/constants.py b/HySoP/hysop/constants.py
old mode 100644
new mode 100755
index a849cb686..0dd0fe4cb
--- a/HySoP/hysop/constants.py
+++ b/HySoP/hysop/constants.py
@@ -56,13 +56,6 @@ WITH_GUESS = 1
 ## y is different from result arg.
 NOALIAS = 2
 
-# File format types for output
-HDF5 = 998
-"""HDF5 format id"""
-
-ASCII = 997
-"""ascii format id"""
-
 ## Default value for task id (mpi task)
 DEFAULT_TASK_ID = 999
 
diff --git a/HySoP/hysop/domain/tests/test_subset.py b/HySoP/hysop/domain/tests/test_subset.py
index 5faf5e19e..00bb7bcc7 100755
--- a/HySoP/hysop/domain/tests/test_subset.py
+++ b/HySoP/hysop/domain/tests/test_subset.py
@@ -1,8 +1,12 @@
+"""Tests computation of subsets of
+a given geometry/domain
+"""
 from hysop.domain.subsets import Sphere, HemiSphere
 from hysop.domain.subsets import Subset, SubBox
 from hysop.domain.subsets import Cylinder, HemiCylinder
 from hysop.operator.hdf_io import HDF_Reader
-from hysop.tools.parameters import Discretization, IOParams
+from hysop.tools.parameters import Discretization
+from hysop.tools.io_utils import IOParams
 from hysop import Field, Box
 from hysop.mpi.topology import Cartesian
 import hysop.tools.numpywrappers as npw
@@ -51,7 +55,7 @@ def check_subset(discr, fileref, class_name):
     subs.discretize(topo)
     assert subs.ind.keys()[0] == topo
     assert len(subs.ind[topo][0]) == dom.dimension
-    scal = Field(domain=topo.domain, is_vector=False)
+    scal = Field(domain=topo.domain, is_vector=False, name='s0')
     sd = scal.discretize(topo)
     sd[0][subs.ind[topo][0]] = 2.
     icompute = topo.mesh.iCompute
@@ -96,7 +100,7 @@ def init_multi(discr, fileref):
     reader.setup()
     reader.apply()
     reader.finalize()
-    scal = Field(domain=topo.domain, is_vector=False)
+    scal = Field(domain=topo.domain, is_vector=False, name='s1')
     sd = scal.discretize(topo)
     return topo, sd, vd
 
@@ -113,7 +117,7 @@ def init_obst_list(discr, fileref):
     obs_list = [s1, s2, s3]
     for obs in obs_list:
         obs.discretize(topo)
-    sd = Field(domain=topo.domain, is_vector=False).discretize(topo)
+    sd = Field(domain=topo.domain, is_vector=False, name='s2').discretize(topo)
     indices = Subset.union(obs_list, topo)
     sd[0][indices] = 2.
     icompute = topo.mesh.iCompute
@@ -138,7 +142,7 @@ def init_subtract(discr, fileref):
     box = SubBox(origin=pos, length=ll, parent=topo.domain)
     s1.discretize(topo)
     box.discretize(topo)
-    sd = Field(domain=topo.domain, is_vector=False).discretize(topo)
+    sd = Field(domain=topo.domain, is_vector=False, name='s2').discretize(topo)
     indices = Subset.subtract(box, s1, topo)
     sd[0][indices] = 2.
     icompute = topo.mesh.iCompute
@@ -176,7 +180,7 @@ def init_subtract_lists(discr, fileref):
         ob.discretize(topo)
     for ob in box:
         ob.discretize(topo)
-    sd = Field(domain=topo.domain, is_vector=False).discretize(topo)
+    sd = Field(domain=topo.domain, is_vector=False, name='s2').discretize(topo)
     indices = Subset.subtract_list_of_sets(box, obs, topo)
     sd[0][indices] = 2.
     return np.allclose(sd[0][topo.mesh.iCompute], vd[0][topo.mesh.iCompute])
diff --git a/HySoP/hysop/f2py/fftw2py.pyf b/HySoP/hysop/f2py/fftw2py.pyf
index a9a833e88..860933a2d 100644
--- a/HySoP/hysop/f2py/fftw2py.pyf
+++ b/HySoP/hysop/f2py/fftw2py.pyf
@@ -119,8 +119,8 @@ module fftw2py ! in fftw2py.f90
     end subroutine solve_curl_2d
     subroutine spectrum_3d(field,spectrum,wavelengths,ghosts,length) ! in fftw2py.f90:fftw2py
         real(kind=pk) dimension(:,:,:),intent(in) :: field
-        real(kind=pk) dimension(:),intent(inout) :: spectrum
-        real(kind=pk) dimension(:),intent(inout) :: wavelengths
+        real(kind=pk) dimension(:),intent(in,out) :: spectrum
+        real(kind=pk) dimension(:),intent(in,out) :: wavelengths
         integer dimension(3),intent(in) :: ghosts
         real(kind=pk) intent(in) :: length
     end subroutine spectrum_3d
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
old mode 100644
new mode 100755
index 4e497228d..38f713105
--- a/HySoP/hysop/fields/continuous.py
+++ b/HySoP/hysop/fields/continuous.py
@@ -17,6 +17,9 @@ from hysop.constants import debug
 from hysop.fields.discrete import DiscreteField
 from hysop.mpi import main_rank
 from hysop.tools.profiler import Profiler
+from hysop.tools.io_utils import IOParams, IO
+from hysop.operator.hdf_io import HDF_Writer, HDF_Reader
+from hysop.mpi.topology import Cartesian
 
 
 class Field(object):
@@ -29,8 +32,8 @@ class Field(object):
     discretize it.
 
     Example :
-    if topo1 and topo2 are two hysop.mpi.topology.Cartesian objects:
-    \code
+    if topo1 and topo2 are two hysop.mpi.topology.Cartesian objects::
+
     scal = Field(domain=dom, name='Scalar')
     # Discretize scal on two different topologies
     scal.discretize(topo1)
@@ -40,7 +43,6 @@ class Field(object):
     scal_discr2 = scal.discreteFields[topo2]
     # or in a more  straightforward way:
     scal_discr2 = scal.discretize(topo2)
-    \endcode
 
     """
 
@@ -49,7 +51,7 @@ class Field(object):
         return object.__new__(cls, *args, **kw)
 
     @debug
-    def __init__(self, domain, formula=None, name=None, is_vector=False,
+    def __init__(self, domain, name, formula=None, is_vector=False,
                  nb_components=None, vectorize_formula=False):
         """Defines a continuous field (scalar or vector) on a domain.
 
@@ -58,10 +60,10 @@ class Field(object):
 
         domain : domain.Domain
             physical domain where this field is defined.
-        formula : python function, optional
-            a user-defined python function. See notes for details.
         name : string
             a name for the variable. Required because of h5py i/o.
+         formula : python function, optional
+            a user-defined python function. See notes for details.
         is_vector : bool, optional
             true if the field is a vector.
         nb_components : int, optional
@@ -85,24 +87,22 @@ class Field(object):
 
         """
 
-        ## Physical domain of definition of the variable.
+        # Physical domain of definition of the variable.
         self.domain = domain
-        ## dimension (1, 2 or 3D), equal to domain dimension.
+        # dimension (1, 2 or 3D), equal to domain dimension.
         self.dimension = self.domain.dimension
-        ## Id (optional)
-        if name is None:
-            name = 'unnamed'
+        # Id
         self.name = name
-        ## Dictionnary of all the discretizations of this field.
-        ## Key = hysop.mpi.topology.Cartesian,
-        ## value =hysop.fields.discrete.DiscreteField.
-        ## Example : vel = discreteFields[topo]
+        # Dictionnary of all the discretizations of this field.
+        # Key = hysop.mpi.topology.Cartesian,
+        # value =hysop.fields.discrete.DiscreteField.
+        # Example : vel = discreteFields[topo]
         self.discreteFields = {}
-        ## Analytic formula.
+        # Analytic formula.
         self.formula = formula
-        ## A list of extra parameters, used in _formula
+        # A list of extra parameters, used in _formula
         self.extraParameters = tuple([])
-        ## Number of components of the field
+        # Number of components of the field
         if nb_components is None:
             self.nb_components = self.dimension if is_vector else 1
         else:
@@ -111,10 +111,16 @@ class Field(object):
         # True if the given formula and must be vectorized
         self._vectorize = vectorize_formula
 
-        ## Time profiling
+        # Time profiling
         self.profiler = Profiler(self, self.domain.comm_task)
 
+        # i/o : a dictionnary of hdf writers.
+        # None by default, if needed must be
+        # initialized with init_writer function.
+        self.writer = {}
+
     def get_profiling_info(self):
+        """Update profiler"""
         for d in self.discreteFields.values():
             self.profiler += d.profiler
 
@@ -136,8 +142,8 @@ class Field(object):
         simply returns the discretized field.
 
         """
-
-        if topo in self.discreteFields.keys():
+        assert isinstance(topo, Cartesian)
+        if topo in self.discreteFields:
             return self.discreteFields[topo]
         else:
             self.discreteFields[topo] = DiscreteField(
@@ -191,6 +197,32 @@ class Field(object):
             df.initialize(self.formula, self._vectorize, time,
                           *self.extraParameters)
 
+    def randomize(self, topo):
+        """Initialize a field on topo with random values.
+
+        Parameters
+        -----------
+        topo : :class:`mpi.topology.Cartesian`
+            the topology on which initialization is performed.
+        """
+        df = self.discretize(topo)
+        df.randomize()
+        return self.discreteFields[topo]
+
+    def copy(self, field_in, topo):
+        """Initialize a field on topo with values from another field.
+
+        Parameters
+        -----------
+        field_in : :class:`fields.continuous.Continuous`
+            field to be copied
+        topo : :class:`mpi.topology.Cartesian`
+            the topology on which initialization is performed.
+        """
+        df = self.discretize(topo)
+        field_in_d = field_in.discretize(topo)
+        df.copy(field_in_d)
+
     def value(self, *pos):
         """Evaluation of the field at a given position, according
         to the analytic formula given during construction.
@@ -297,37 +329,81 @@ class Field(object):
         """
         self.extraParameters = args
 
-    def hdf_dump(self, discretization, io_params=None):
-        from hysop.operator.hdf_io import HDF_Writer
-        from hysop.problem.simulation import Simulation
-        simu = Simulation(nbIter=1)
+    def init_writer(self, discretization, io_params=None):
+        """Create an hdf writer for a given discretization
+
+        Parameters
+        ----------
+        discretization : :class:`~hysop.tools.params.Discretization` or
+         :class:`~hysop.mpi.topology.Cartesian`
+            set a data distribution description for the output, i.e.
+            choose the space discretisation used for hdf file.
+        io_params : :class:`~hysop.tools.io_utils.IOParams`, optional
+            output file description
+        """
         if io_params is None:
-            from hysop.tools.parameters import IOParams
-            from hysop.constants import HDF5
-            io_params = IOParams(self.name, fileformat=HDF5)
-        wr = HDF_Writer(variables={self: discretization},
-                        io_params=io_params)
+            io_params = IOParams(self.name, fileformat=IO.HDF5)
+        wr = HDF_Writer(variables={self: discretization}, io_params=io_params)
         wr.discretize()
         wr.setup()
-        wr.apply(simu)
-        wr.finalize()
+        self.writer[self.discreteFields[wr.topology]] = wr
 
-    def hdf_load(self, discretization, io_params=None, restart=None):
-        from hysop.operator.hdf_io import HDF_Reader
+    def hdf_dump(self, topo, simu):
+        """Save field data into hdf file
 
-        if io_params is None:
-            from hysop.tools.parameters import IOParams
-            from hysop.constants import HDF5
-            io_params = IOParams(self.name, fileformat=HDF5)
-        read = HDF_Reader(variables={self: discretization},
-                          io_params=io_params, restart=restart)
+        Parameters
+        ----------
+        topo : :class:`~hysop.mpi.topology.Cartesian`
+            data distribution description
+        simu : :class:`~hysop.problem.simulation.Simulation`
+            time discretisation
+        """
+        df = self.discretize(topo)
+        if df not in self.writer:
+            self.init_writer(topo)
+        self.writer[df].apply(simu)
+        self.writer[df].finalize()
+
+    def hdf_load(self, topo, io_params=None,
+                 restart=None, dataset_name=None):
+        """Load field from hdf fileformat
+
+        Parameters
+        ----------
+        topo : :class:`~hysop.mpi.topology.Cartesian`
+            data distribution description
+        io_params : :class:`~hysop.tools.io_utils.IOParams`, optional
+            output file description
+        restart : int, optional
+            iteration number. Default=None, i.e. read first iteration.
+        dataset_name : string, optional
+            name of the dataset in h5 file used to fill this field in.
+            If None, look for self.name as dataset name.
+
+        Notes
+        -----
+
+        * a call to hdf_load(topo, ...) will discretize the
+        field on topo, if not already done
+        * to read the dataset 'field1' from file ff_00001.h5::
+
+            myfield = Field(box, name='f1')
+            myfield.hdf_load(topo, IOParams('ff', restart=1,
+                             dataset_name='field1')
+
+        """
+        if dataset_name is None:
+            var_names = None
+        else:
+            var_names = {self: dataset_name}
+        read = HDF_Reader(variables={self: topo},
+                          io_params=io_params, restart=restart,
+                          var_names=var_names)
         read.discretize()
         read.setup()
         read.apply()
         read.finalize()
 
-        return 0
-
     def zero(self, topo=None):
         """Reset to 0.0 all components of this field.
 
@@ -339,7 +415,3 @@ class Field(object):
         else:
             for dfield in self.discreteFields.values():
                 dfield.zero()
-
-    def finalize(self):
-        for dfield in self.discreteFields.values():
-            dfield.finalize()
diff --git a/HySoP/hysop/fields/discrete.py b/HySoP/hysop/fields/discrete.py
old mode 100644
new mode 100755
index eef86310a..3f941c4f9
--- a/HySoP/hysop/fields/discrete.py
+++ b/HySoP/hysop/fields/discrete.py
@@ -5,7 +5,7 @@
 Documentation and examples can be found in :ref:`fields`.
 
 """
-from hysop.constants import debug
+from hysop.constants import debug, ORDER
 from itertools import count
 from hysop.tools.profiler import Profiler
 import numpy as np
@@ -160,6 +160,24 @@ class DiscreteField(object):
             "WARNING: The shape of " + self.name + " has changed during"\
             " field initialisation."
 
+    def randomize(self):
+        """Initialize a the with random values.
+
+        """
+        for d in xrange(self.nb_components):
+            self.data[d][...] = np.random.random(self.topology.mesh.resolution)
+
+    def copy(self, field_in):
+        """Initialize a field on topo with values from another field.
+
+        Parameters
+        -----------
+        field_in : :class:`fields.discrete.DiscreteField`
+            field to be copied
+        """
+        for d in xrange(self.nb_components):
+            self.data[d][...] = field_in.data[d].copy(order=ORDER)
+
     def norm(self):
         """Compute the euclidian norm of the discrete field
         p-norm = (sum data[d](i,...)**p)**1/p for d = 1..dim
@@ -200,8 +218,6 @@ class DiscreteField(object):
         for dim in xrange(self.nb_components):
             self.data[dim][...] = 0.0
 
-    def finalize(self):
-        pass
-
     def get_profiling_info(self):
+        """update profiler"""
         pass
diff --git a/HySoP/hysop/fields/tests/func_for_tests.py b/HySoP/hysop/fields/tests/func_for_tests.py
old mode 100644
new mode 100755
index ce3ea9d16..3e2eaf562
--- a/HySoP/hysop/fields/tests/func_for_tests.py
+++ b/HySoP/hysop/fields/tests/func_for_tests.py
@@ -1,10 +1,10 @@
-"""
-@file func_for_tests.py a set of functions, used in tests
+"""A set of functions, used in tests
 to initialize the fields
 """
 import numpy as np
-sin = np.sin
 import math
+sin = np.sin
+cos = np.cos
 pi = math.pi
 
 
@@ -88,3 +88,36 @@ def v2d(res, x, y, t):
     res[0][...] = sin(x) + t
     res[1][...] = y
     return res
+
+
+def v_time_2d(res, x, y, t):
+    """For 2d vector field, all components depend on time
+    """
+    res[0][...] = sin(x) + t
+    res[1][...] = y * t
+    return res
+
+
+def v_time_3d(res, x, y, z, t):
+    """For 3d vector field, all components depend on time
+    """
+    res[0][...] = sin(x) + t
+    res[1][...] = y * t
+    res[2][...] = z * t ** 2
+    return res
+
+
+def v_TG(res, x, y, z, t):
+    """Taylor Green velocity"""
+    res[0][...] = sin(x) * cos(y) * cos(z)
+    res[1][...] = - cos(x) * sin(y) * cos(z)
+    res[2][...] = 0.
+    return res
+
+
+def w_TG(res, x, y, z, t): 
+    """Taylor Green vorticity"""
+    res[0][...] = - cos(x) * sin(y) * sin(z)
+    res[1][...] = - sin(x) * cos(y) * sin(z)
+    res[2][...] = 2. * sin(x) * sin(y) * cos(z)
+    return res
diff --git a/HySoP/hysop/fields/tests/test_field.py b/HySoP/hysop/fields/tests/test_field.py
index b5181bbd2..6247b5dee 100755
--- a/HySoP/hysop/fields/tests/test_field.py
+++ b/HySoP/hysop/fields/tests/test_field.py
@@ -1,29 +1,37 @@
+"""Test fields build and usage (discrete and continuous)
 """
-Testing hysop.field.continuous.Field
-"""
-import hysop
-print hysop.__file__
-from hysop.fields.continuous import Field
-from hysop.domain.box import Box
-from hysop.tools.parameters import Discretization
+from hysop import Simulation, Box, Field, Discretization, IOParams, IO
 import numpy as np
+import os
 from hysop.fields.tests.func_for_tests import func_scal_1, func_scal_2, \
-    func_vec_1, func_vec_2, func_vec_3, func_vec_4, func_vec_5, func_vec_6
+    func_vec_1, func_vec_2, func_vec_3, func_vec_4, func_vec_5, func_vec_6,\
+    v_time_3d, v_time_2d
 from numpy import allclose
+import shutil
+from hysop.mpi import main_rank
+
 
-d3D = Discretization([33, 33, 33])
-d2D = Discretization([33, 33])
+d3d = Discretization([33, 33, 33])
+d2d = Discretization([33, 33])
 nbc = 4
 
 
+def formula(x, y, z):
+    return (x * y * z, 1.0, -x * y * z)
+
+
+def form2(x, y, z, dt):
+    return (x * y * z, 1.0, -x * y * z * dt)
+
+
 def test_continuous():
     """ Basic continuous field construction """
 
     dom = Box()
-    cf = Field(dom)
+    cf = Field(dom, name='f1')
     assert cf.nb_components == 1
     assert cf.domain is dom
-    cf = Field(dom, is_vector=True)
+    cf = Field(dom, name='f2', is_vector=True)
     assert cf.nb_components == dom.dimension
     assert cf.domain is dom
 
@@ -31,8 +39,7 @@ def test_continuous():
 def test_analytical():
     """Test formula"""
     dom = Box()
-    formula = lambda x, y, z: (x * y * z, 1.0, -x * y * z)
-    caf = Field(dom, formula, vectorize_formula=True)
+    caf = Field(dom, name='f1', formula=formula, vectorize_formula=True)
     assert caf.value(0., 0., 0.) == (0., 1., 0.)
     assert caf.value(1., 1., 1.) == (1., 1., -1.)
     assert caf.value(1., 2., 3.) == (6., 1., -6.)
@@ -41,12 +48,10 @@ def test_analytical():
 def test_analytical_reset():
     """Test formula"""
     dom = Box()
-    formula = lambda x, y, z: (x * y * z, 1.0, -x * y * z)
-    caf = Field(dom, formula, vectorize_formula=True)
+    caf = Field(dom, name='f1', formula=formula, vectorize_formula=True)
     assert caf.value(0., 0., 0.) == (0., 1., 0.)
     assert caf.value(1., 1., 1.) == (1., 1., -1.)
     assert caf.value(1., 2., 3.) == (6., 1., -6.)
-    form2 = lambda x, y, z, dt: (x * y * z, 1.0, -x * y * z * dt)
     caf.set_formula(form2, vectorize_formula=True)
     assert caf.value(0., 0., 0., 0.) == (0., 1., 0.)
     assert caf.value(1., 1., 1., 1.) == (1., 1., -1.)
@@ -56,8 +61,8 @@ def test_analytical_reset():
 def test_discretization():
     """Test multiple discretizations"""
     dom = Box()
-    csf = Field(dom)
-    cvf = Field(dom, is_vector=True)
+    csf = Field(dom, name='f1')
+    cvf = Field(dom, name='f2', is_vector=True)
     r3d = Discretization([33, 33, 17])
     topo = dom.create_topology(r3d)
     csf.discretize(topo)
@@ -75,10 +80,10 @@ def test_discretization():
 # Non-Vectorized formula for a scalar
 def test_analytical_field_1():
     box = Box()
-    topo = box.create_topology(discretization=d3D)
+    topo = box.create_topology(discretization=d3d)
     coords = topo.mesh.coords
-    caf = Field(box, formula=func_scal_1)
-    ref = Field(box)
+    caf = Field(box, name='f1', formula=func_scal_1)
+    ref = Field(box, name='f2')
     refd = ref.discretize(topo)
     cafd = caf.discretize(topo)
     refd = ref.discretize(topo)
@@ -97,10 +102,10 @@ def test_analytical_field_1():
 # Vectorized formula
 def test_analytical_field_2():
     box = Box()
-    topo = box.create_topology(d3D)
+    topo = box.create_topology(d3d)
     coords = topo.mesh.coords
-    caf = Field(box, formula=func_scal_2, vectorize_formula=True)
-    ref = Field(box)
+    caf = Field(box, name='f1', formula=func_scal_2, vectorize_formula=True)
+    ref = Field(box, name='f2')
     cafd = caf.discretize(topo)
     ids = id(cafd.data[0])
     refd = ref.discretize(topo)
@@ -118,10 +123,10 @@ def test_analytical_field_2():
 # Non-Vectorized formula for a vector
 def test_analytical_field_3():
     box = Box()
-    topo = box.create_topology(d3D)
+    topo = box.create_topology(d3d)
     coords = topo.mesh.coords
-    caf = Field(box, formula=func_vec_1, is_vector=True)
-    ref = Field(box, is_vector=True)
+    caf = Field(box, name='f1', formula=func_vec_1, is_vector=True)
+    ref = Field(box, name='f2', is_vector=True)
     refd = ref.discretize(topo)
     cafd = caf.discretize(topo)
     refd = ref.discretize(topo)
@@ -144,11 +149,11 @@ def test_analytical_field_3():
 # Vectorized formula for a vector
 def test_analytical_field_4():
     box = Box()
-    topo = box.create_topology(d3D)
+    topo = box.create_topology(d3d)
     coords = topo.mesh.coords
-    caf = Field(box, formula=func_vec_2, is_vector=True,
+    caf = Field(box, name='f1', formula=func_vec_2, is_vector=True,
                 vectorize_formula=True)
-    ref = Field(box, is_vector=True)
+    ref = Field(box, name='f2', is_vector=True)
     refd = ref.discretize(topo)
     cafd = caf.discretize(topo)
     refd = ref.discretize(topo)
@@ -171,12 +176,12 @@ def test_analytical_field_4():
 # Non-Vectorized formula for a vector, with extra-arguments
 def test_analytical_field_5():
     box = Box()
-    topo = box.create_topology(d3D)
+    topo = box.create_topology(d3d)
     coords = topo.mesh.coords
-    caf = Field(box, formula=func_vec_3, is_vector=True)
+    caf = Field(box, name='f1', formula=func_vec_3, is_vector=True)
     theta = 0.3
     caf.set_formula_parameters(theta)
-    ref = Field(box, is_vector=True)
+    ref = Field(box, name='f2', is_vector=True)
     refd = ref.discretize(topo)
     cafd = caf.discretize(topo)
     refd = ref.discretize(topo)
@@ -199,13 +204,13 @@ def test_analytical_field_5():
 # Vectorized formula for a vector, with extra-arguments
 def test_analytical_field_6():
     box = Box()
-    topo = box.create_topology(d3D)
+    topo = box.create_topology(d3d)
     coords = topo.mesh.coords
-    caf = Field(box, formula=func_vec_4, is_vector=True,
+    caf = Field(box, name='f1', formula=func_vec_4, is_vector=True,
                 vectorize_formula=True)
     theta = 0.3
     caf.set_formula_parameters(theta)
-    ref = Field(box, is_vector=True)
+    ref = Field(box, name='f2', is_vector=True)
     refd = ref.discretize(topo)
     cafd = caf.discretize(topo)
     refd = ref.discretize(topo)
@@ -229,12 +234,12 @@ def test_analytical_field_6():
 # different from domain dim and  with extra-arguments
 def test_analytical_field_7():
     box = Box()
-    topo = box.create_topology(d3D)
+    topo = box.create_topology(d3d)
     coords = topo.mesh.coords
-    caf = Field(box, formula=func_vec_5, nb_components=nbc)
+    caf = Field(box, name='f1', formula=func_vec_5, nb_components=nbc)
     theta = 0.3
     caf.set_formula_parameters(theta)
-    ref = Field(box, nb_components=nbc)
+    ref = Field(box, name='f2', nb_components=nbc)
     refd = ref.discretize(topo)
     cafd = caf.discretize(topo)
     refd = ref.discretize(topo)
@@ -259,12 +264,12 @@ def test_analytical_field_7():
 # different from domain dim and  with extra-arguments
 def test_analytical_field_8():
     box = Box(dimension=2, length=[1., 1.], origin=[0., 0.])
-    topo = box.create_topology(d2D)
+    topo = box.create_topology(d2d)
     coords = topo.mesh.coords
-    caf = Field(box, formula=func_vec_6, nb_components=nbc)
+    caf = Field(box, name='f1', formula=func_vec_6, nb_components=nbc)
     theta = 0.3
     caf.set_formula_parameters(theta)
-    ref = Field(box, nb_components=nbc)
+    ref = Field(box, name='f2', nb_components=nbc)
     refd = ref.discretize(topo)
     cafd = caf.discretize(topo)
     refd = ref.discretize(topo)
@@ -288,13 +293,13 @@ def test_analytical_field_8():
 # topologies.
 def test_analytical_field_9():
     box = Box()
-    topo = box.create_topology(d3D)
+    topo = box.create_topology(d3d)
     res2 = Discretization([65, 33, 65], [1, 1, 1])
     topo2 = box.create_topology(res2, dim=2)
     coords = topo.mesh.coords
     coords2 = topo2.mesh.coords
-    caf = Field(box, formula=func_vec_1, is_vector=True)
-    ref = Field(box, is_vector=True)
+    caf = Field(box, name='f1', formula=func_vec_1, is_vector=True)
+    ref = Field(box, name='f2', is_vector=True)
     refd = ref.discretize(topo)
     cafd = caf.discretize(topo)
     cafd2 = caf.discretize(topo2)
@@ -319,13 +324,13 @@ def test_analytical_field_9():
 # topologies.
 def test_analytical_field_10():
     box = Box()
-    topo = box.create_topology(d3D)
+    topo = box.create_topology(d3d)
     res2 = Discretization([65, 33, 65], [1, 1, 1])
     topo2 = box.create_topology(res2, dim=2)
     coords = topo.mesh.coords
     coords2 = topo2.mesh.coords
-    caf = Field(box, formula=func_vec_1, is_vector=True)
-    ref = Field(box, is_vector=True)
+    caf = Field(box, name='f1', formula=func_vec_1, is_vector=True)
+    ref = Field(box, name='f2', is_vector=True)
     refd = ref.discretize(topo)
     cafd = caf.discretize(topo)
     cafd2 = caf.discretize(topo2)
@@ -341,3 +346,93 @@ def test_analytical_field_10():
         assert allclose(cafd2[i], refd2.data[i])
         assert allclose(cafd[i], refd.data[i])
         assert id(cafd2.data[i]) == ids[i]
+
+
+def test_copy():
+    box = Box()
+    topo = box.create_topology(d3d)
+    f1 = Field(box, name='source', is_vector=True, formula=func_vec_1)
+    f2 = Field(box, name='target', is_vector=True)
+    f1d = f1.discretize(topo)
+    f1.initialize(topo=topo)
+    f2.copy(f1, topo)
+    f2d = f2.discretize(topo)
+    for d in xrange(f2.nb_components):
+        assert np.allclose(f2d[d], f1d[d])
+        assert f2d[d].flags.f_contiguous == f1d[d].flags.f_contiguous
+
+
+def test_randomize():
+    box = Box()
+    topo = box.create_topology(d3d)
+    f2 = Field(box, name='target', is_vector=True)
+    f2.randomize(topo)
+    f2d = f2.discretize(topo)
+    for d in xrange(f2.nb_components):
+        assert not np.allclose(f2d[d], 0.)
+
+
+def hdf_dump_load(discretisation, formula):
+    dimension = len(discretisation.resolution)
+    box = Box([1., ] * dimension)
+    topo = box.create_topology(discretisation)
+    fname = 'f1_' + str(dimension)
+    ff = Field(box, name=fname, is_vector=True, formula=formula)
+    iop = IOParams(ff.name)
+    simu = Simulation(nbIter=4)
+    simu.initialize()
+    ff.initialize(time=simu.time, topo=topo)
+    ff.hdf_dump(topo, simu)
+    simu.advance()
+    ff.initialize(time=simu.time, topo=topo)
+    # Write ff for current time and topo
+    ff.hdf_dump(topo, simu)
+    assert iop.filepath == IO.default_path()
+    assert os.path.isfile(iop.filename + '_00000.h5')
+    assert os.path.isfile(iop.filename + '_00001.h5')
+    assert os.path.isfile(iop.filename + '.xmf')
+    fname2 = 'f2_' + str(dimension)
+    gg = Field(box, name=fname2, is_vector=True)
+    iop_in = IOParams(iop.filename + '_00001.h5')
+    dsname = ff.name + '_' + str(topo.get_id())
+    # Load gg from ff values at current time, on topo
+    for i in xrange(gg.nb_components):
+        assert not np.allclose(gg.discretize(topo)[i], ff.discretize(topo)[i])
+    gg.hdf_load(topo, iop_in, dataset_name=dsname)
+    for i in xrange(gg.nb_components):
+        assert np.allclose(gg.discretize(topo)[i], ff.discretize(topo)[i])
+
+    # reset ff with its values at time 'O', on topo
+    ff.hdf_load(topo, restart=0)
+    for i in xrange(gg.nb_components):
+        assert not np.allclose(gg.discretize(topo)[i], ff.discretize(topo)[i])
+
+
+def test_hdf_dump_load_2d():
+    hdf_dump_load(d2d, v_time_2d)
+
+
+def test_hdf_dump_load_3d():
+    hdf_dump_load(d3d, v_time_3d)
+
+
+# This may be useful to run mpi tests
+if __name__ == "__main__":
+    test_continuous()
+    test_analytical()
+    test_analytical_reset()
+    test_discretization()
+    test_analytical_field_1()
+    test_analytical_field_2()
+    test_analytical_field_3()
+    test_analytical_field_4()
+    test_analytical_field_5()
+    test_analytical_field_6()
+    test_analytical_field_7()
+    test_analytical_field_8()
+    test_copy()
+    test_randomize()
+    test_hdf_dump_load_2d()
+    test_hdf_dump_load_3d()
+    if main_rank == 0:
+        shutil.rmtree(IO.default_path())
diff --git a/HySoP/hysop/gpu/tests/test_gpu_multiresolution_filter.py b/HySoP/hysop/gpu/tests/test_gpu_multiresolution_filter.py
index 4c87351fb..b44201bbd 100755
--- a/HySoP/hysop/gpu/tests/test_gpu_multiresolution_filter.py
+++ b/HySoP/hysop/gpu/tests/test_gpu_multiresolution_filter.py
@@ -33,7 +33,7 @@ def test_filter_linear():
     """This test compares the GPU linear filter with python implementation"""
     box = Box(length=L, origin=O, proc_tasks=PROC_TASKS)
     mpi_p = MPIParams(comm=box.comm_task, task_id=1)
-    f = Field(box, formula=func, is_vector=False)
+    f = Field(box, formula=func, is_vector=False, name='f1')
     d_fine = Discretization([513, 513, 513])
     d_coarse = Discretization([257, 257, 257], ghosts=[1, 1, 1])
     op = MultiresolutionFilter(d_in=d_fine, d_out=d_coarse,
@@ -62,7 +62,7 @@ def test_filter_linear():
             np.max(np.abs(valid[0][topo_coarse.mesh.iCompute] -
                           f_out[0][topo_coarse.mesh.iCompute]))
         if PY_COMPARE:
-            f_py = Field(box, formula=func)
+            f_py = Field(box, formula=func, name='fpy')
             op_py = MultiresolutionFilter(d_in=d_fine, d_out=d_coarse,
                                           variables={f_py: d_coarse},
                                           method={Remesh: Rmsh_Linear, },
@@ -85,7 +85,7 @@ def test_filter_L2_1():
     """
     box = Box(length=L, origin=O, proc_tasks=PROC_TASKS)
     mpi_p = MPIParams(comm=box.comm_task, task_id=1)
-    f = Field(box, formula=func)
+    f = Field(box, formula=func, name='f1')
     d_fine = Discretization([513, 513, 513])
     d_coarse = Discretization([257, 257, 257], ghosts=[2, 2, 2])
     if box.is_on_task(1):
@@ -113,7 +113,7 @@ def test_filter_L2_1():
             np.max(np.abs(valid[0][topo_coarse.mesh.iCompute] -
                           f_out[0][topo_coarse.mesh.iCompute]))
         if PY_COMPARE:
-            f_py = Field(box, formula=func)
+            f_py = Field(box, formula=func, name='fpy')
             op_py = MultiresolutionFilter(d_in=d_fine, d_out=d_coarse,
                                           variables={f_py: d_coarse},
                                           method={Remesh: L2_1, })
diff --git a/HySoP/hysop/gpu/tests/test_multiphase_baroclinic.py b/HySoP/hysop/gpu/tests/test_multiphase_baroclinic.py
index fbfd41963..742ba3629 100755
--- a/HySoP/hysop/gpu/tests/test_multiphase_baroclinic.py
+++ b/HySoP/hysop/gpu/tests/test_multiphase_baroclinic.py
@@ -78,12 +78,12 @@ def call_operator(func, grad_func, vfunc):
     """Call the baroclinic rhs operator from given initialization functions"""
     simu = Simulation(tinit=0.0, tend=1.0, timeStep=0.1, iterMax=1)
     box = Box()
-    rhs = Field(box, is_vector=True)
-    gradp = Field(box, is_vector=True, formula=vfunc)
-    rho = Field(box, is_vector=False, formula=func)
-    gradrho = Field(box, is_vector=True, formula=grad_func)
-    gradp_fine = Field(box, is_vector=True, formula=vfunc)
-    true_rhs = Field(box, is_vector=True)
+    rhs = Field(box, is_vector=True, name='rhs')
+    gradp = Field(box, is_vector=True, formula=vfunc, name='gradp')
+    rho = Field(box, is_vector=False, formula=func, name='rho')
+    gradrho = Field(box, is_vector=True, formula=grad_func, name='gradrho')
+    gradp_fine = Field(box, is_vector=True, formula=vfunc, name='fine')
+    true_rhs = Field(box, is_vector=True, name='ref')
     d_fine = Discretization([257, 257, 257])
     d_coarse = Discretization([129, 129, 129], ghosts=[2, 2, 2])
     op = MultiphaseBaroclinicRHS(rhs, rho, gradp,
diff --git a/HySoP/hysop/numerics/integrators/odesolver.py b/HySoP/hysop/numerics/integrators/odesolver.py
old mode 100644
new mode 100755
index ab12ac76e..6b6cf41f4
--- a/HySoP/hysop/numerics/integrators/odesolver.py
+++ b/HySoP/hysop/numerics/integrators/odesolver.py
@@ -51,6 +51,7 @@ class ODESolver(NumMethod):
         lwork = self.getWorkLengths(nb_components)
         assert len(self.work) == lwork, 'Wrong length for work arrays list.'
         self._memshape = tuple(topo.mesh.resolution)
+        print topo.mesh
         for wk in self.work:
             assert wk.shape == tuple(self._memshape)
         # Allocate buffers for ghost points synchronization.
diff --git a/HySoP/hysop/numerics/tests/test_integrators.py b/HySoP/hysop/numerics/tests/test_integrators.py
old mode 100644
new mode 100755
diff --git a/HySoP/hysop/numerics/tests/test_update_ghosts.py b/HySoP/hysop/numerics/tests/test_update_ghosts.py
index e1970b2f3..8ec8f843f 100755
--- a/HySoP/hysop/numerics/tests/test_update_ghosts.py
+++ b/HySoP/hysop/numerics/tests/test_update_ghosts.py
@@ -8,7 +8,7 @@ import numpy as np
 def _setup(ghosts, topo_dim):
     dom = Box()
 
-    f = Field(dom, is_vector=True)
+    f = Field(dom, is_vector=True, name='f')
     resolTopo = Discretization([33, 33, 33], ghosts=ghosts)
     topo = dom.create_topology(resolTopo, dim=topo_dim)
     f.discretize(topo)
diff --git a/HySoP/hysop/operator/computational.py b/HySoP/hysop/operator/computational.py
old mode 100644
new mode 100755
index 31ecca3eb..251e7e678
--- a/HySoP/hysop/operator/computational.py
+++ b/HySoP/hysop/operator/computational.py
@@ -1,10 +1,9 @@
 """ Abstract interface for computational operators
 
-
 """
 from abc import ABCMeta, abstractmethod
 from hysop.constants import debug
-from hysop.operator.continuous import Operator
+from hysop.operator.continuous import Operator, opapply
 from hysop.mpi.topology import Cartesian
 from hysop.tools.parameters import Discretization
 from hysop.tools.profiler import profile
@@ -44,7 +43,7 @@ class Computational(Operator):
         method : :class:`hysop.method`
             specific parameters for the operator
             (space discretisation ...). See methods.py for authorized values.
-        **kwds : arguments passed to operator.continuous.Operator base class.
+        **kwds : arguments passed to base class.
 
         Notes:
         ------
@@ -54,20 +53,20 @@ class Computational(Operator):
         # Set mpi related stuff
         self._set_domain_and_tasks()
 
-        ## A dictionnary of parameters, to define numerical methods used
-        ## used to discretize the operator.
-        ## When method is None, each operator must provide a default
-        ## behavior.
+        # A dictionnary of parameters, to define numerical methods used
+        # used to discretize the operator.
+        # When method is None, each operator must provide a default
+        # behavior.
         self.method = method
 
-        ## The discretization of this operator.
+        # The discretization of this operator.
         self.discrete_op = None
 
-        ## A dictionnary of discreteFields associated with this operator
-        ## key = continuous variable \n
-        ## Example : discrete_variable = discreteFields[velocity]
-        ## returns the discrete fields that corresponds to the continuous
-        ## variable velocity.
+        # A dictionnary of discreteFields associated with this operator
+        # key = continuous variable \n
+        # Example : discrete_variable = discreteFields[velocity]
+        # returns the discrete fields that corresponds to the continuous
+        # variable velocity.
         self.discreteFields = {}
 
         if not self._varsfromList:
@@ -311,14 +310,16 @@ class Computational(Operator):
 
     @debug
     @profile
+    @opapply
     def apply(self, simulation=None):
+        """Apply this operator to its variables.
+
+        Parameters
+        ----------
+        simulation : `:class::~hysop.problem.simulation.Simulation`
+
         """
-        Apply this operator to its variables.
-        @param simulation : object that describes the simulation
-        parameters (time, time step, iteration number ...), see
-        hysop.problem.simulation.Simulation for details.
-        """
-        super(Computational, self).apply(simulation)
+        #super(Computational, self).apply(simulation)
         if self.discrete_op is not None:
             self.discrete_op.apply(simulation)
 
diff --git a/HySoP/hysop/operator/continuous.py b/HySoP/hysop/operator/continuous.py
old mode 100644
new mode 100755
index 994ae33e0..a8ee8ab44
--- a/HySoP/hysop/operator/continuous.py
+++ b/HySoP/hysop/operator/continuous.py
@@ -4,8 +4,10 @@
 from abc import ABCMeta, abstractmethod
 from hysop.constants import debug
 from hysop.tools.profiler import Profiler
-from hysop.tools.parameters import MPIParams, IOParams
+from hysop.tools.io_utils import IOParams, IO
+from hysop.tools.parameters import MPIParams
 import hysop.tools.io_utils as io
+import inspect
 
 
 class Operator(object):
@@ -23,11 +25,11 @@ class Operator(object):
                  io_params=None):
         """
         Parameters
-        ----------------
+        ----------
         variables : list or dictionnary. See Notes for details.
         mpi_params : :class:`hysop.tools.parameters.MPIParams`
             mpi config for the operator (comm, task ...)
-        io_params : :class:`hysop.tools.parameters.IOParams`
+        io_params : :class:`hysop.tools.io_utils.IOParams`
             file i/o config (filename, format ...)
 
         Notes
@@ -95,20 +97,20 @@ class Operator(object):
         # _set_domain_and_tasks, called in derived class, since it may
         # require some specific initialization (check domain ...)
 
-        ## Input variables.
+        # Input variables.
         self.input = []
-        ## Output variables.
+        # Output variables.
         self.output = []
-        ## bool to check if the setup function has been called for
-        ## this operator
+        # bool to check if the setup function has been called for
+        # this operator
         self._is_uptodate = False
 
         self.name = self.__class__.__name__
-        ## List of operators that must be waited for.
+        # List of operators that must be waited for.
         self._wait_list = []
-        ## time monitoring
+        # time monitoring
         self.time_info = None
-        ## Dictionnary of optional parameters for output
+        # Dictionnary of optional parameters for output
         self.io_params = io_params
         # Object that deals with output file writing.
         # Optional.
@@ -140,7 +142,10 @@ class Operator(object):
         # Set profiler
         self.profiler = Profiler(self, self.domain.comm_task)
 
-    def _error_(self):
+    @staticmethod
+    def _error_():
+        """internal error message
+        """
         raise RuntimeError("This operator is not defined for the current task")
 
     def wait_for(self, op):
@@ -252,9 +257,19 @@ class Operator(object):
 
         """
         iopar = self.io_params
+        # if iopar is not None (i.e. set in operator init)
+        # and True or set as an IOParams , then
+        # build a writer
         if iopar:
             if isinstance(iopar, bool):
-                self.io_params = IOParams(filename)
+                # default values for iop
+                self.io_params = IOParams(filename, fileformat=IO.ASCII)
+            else:
+                if self.io_params.fileformat is not IO.ASCII:
+                    self.io_params.fileformat = IO.ASCII
+                    msg = 'Warning, wrong file format for operator output.'
+                    msg += 'This will be reset to ASCII.'
+                    print msg
             self._writer = io.Writer(io_params=self.io_params,
                                      mpi_params=self._mpis,
                                      buffshape=buffshape)
@@ -266,9 +281,6 @@ class Operator(object):
         return self._mpis.task_id
 
 
-import inspect
-
-
 def opsetup(f):
     """
     Setup decorator: what must be done by all operators
@@ -277,6 +289,8 @@ def opsetup(f):
     """
 
     def decorator(*args, **kwargs):
+        """Decorate 'setup' method
+        """
         # Job before setup of the function ...
         # nothing for the moment
         name = inspect.getmro(args[0].setup.im_class)
@@ -306,8 +320,10 @@ def opapply(f):
     Usage : add @opapply before apply class method
     """
     def decorator(*args, **kwargs):
+        """decorate 'apply' method"""
+        # get 'continuous' base class and run its apply function
         name = inspect.getmro(args[0].apply.im_class)
-        name[-3].setup(args[0])
+        name[-2].apply(args[0])
         #super(args[0].__class__, args[0]).apply()
         for op in args[0].wait_list():
             op.wait()
diff --git a/HySoP/hysop/operator/discrete/discrete.py b/HySoP/hysop/operator/discrete/discrete.py
old mode 100644
new mode 100755
index baff4c8f8..2069227a8
--- a/HySoP/hysop/operator/discrete/discrete.py
+++ b/HySoP/hysop/operator/discrete/discrete.py
@@ -143,11 +143,11 @@ class DiscreteOperator(object):
     @debug
     @abstractmethod
     def apply(self, simulation=None):
-        """
-        Apply this operator to its variables.
-        @param simulation : object that describes the simulation
-        parameters (time, time step, iteration number ...), see
-        hysop.problem.simulation.Simulation for details.
+        """Execute the operator for the current simulation state
+
+        Parameters
+        ----------
+        simulation : :class:`~hysop.problem.simulation.Simulation`
         """
 
     @debug
diff --git a/HySoP/hysop/operator/discrete/spectrum.py b/HySoP/hysop/operator/discrete/spectrum.py
old mode 100644
new mode 100755
diff --git a/HySoP/hysop/operator/drag_and_lift.py b/HySoP/hysop/operator/drag_and_lift.py
old mode 100644
new mode 100755
index a93587f24..b20ff7754
--- a/HySoP/hysop/operator/drag_and_lift.py
+++ b/HySoP/hysop/operator/drag_and_lift.py
@@ -3,11 +3,19 @@
 
 .. currentmodule:: hysop.operator
 
+* :class:`~drag_and_lift.MomentumForces` : Momentum Formula
+* :class:`~drag_and_lift.NocaForces` : Noca formulation
+ (formulation = 1, 2 or 3)
+
+See :ref:`forces`.
+
+
 """
 from hysop.operator.computational import Computational
 from hysop.operator.continuous import opsetup
 from abc import ABCMeta, abstractmethod
 from hysop.domain.control_box import ControlBox
+import numpy as np
 
 
 class Forces(Computational):
@@ -91,8 +99,6 @@ class MomentumForces(Forces):
 
     """
 
-    __metaclass__ = ABCMeta
-
     def __init__(self, velocity, penalisation_coeff, **kwds):
         """
         Parameters
@@ -141,8 +147,6 @@ class NocaForces(Forces):
     (See Noca99 or Plouhmans, 2002, Journal of Computational Physics)
     """
 
-    __metaclass__ = ABCMeta
-
     def __init__(self, velocity, vorticity, nu, formulation=1,
                  volume_of_control=None, surfdir=None, **kwds):
         """
diff --git a/HySoP/hysop/operator/hdf_io.py b/HySoP/hysop/operator/hdf_io.py
old mode 100644
new mode 100755
index e6eb02f95..e0771b289
--- a/HySoP/hysop/operator/hdf_io.py
+++ b/HySoP/hysop/operator/hdf_io.py
@@ -1,14 +1,17 @@
-"""
-@file printer.py
+"""I/O operators
+
+.. currentmodule:: hysop.operator.hdf_io
+
+* :class:`~HDF_Writer` : operator to write fields into an hdf file
+* :class:`~HDF_Reader` : operator to read fields from an hdf file
+* :class:`~HDF_IO` abstract interface for hdf io classes
 
-File output for field(s) value on a grid.
 """
-from hysop.constants import S_DIR, debug, HDF5, HYSOP_REAL
+from hysop.constants import S_DIR, debug, HYSOP_REAL
 from hysop.operator.computational import Computational
 from hysop.operator.continuous import opapply, opsetup
 import hysop.tools.numpywrappers as npw
-import hysop.tools.io_utils as io
-from hysop.tools.parameters import IOParams
+from hysop.tools.io_utils import IO, IOParams, XMF
 from abc import ABCMeta, abstractmethod
 
 try:
@@ -23,33 +26,37 @@ from hysop.tools.profiler import profile
 
 
 class HDF_IO(Computational):
-    """
-    Abstract interface for read/write from/to hdf files, for
+    """Abstract interface to read/write from/to hdf files, for
     hysop fields.
     """
 
     __metaclass__ = ABCMeta
 
     def __init__(self, var_names=None, subset=None, **kwds):
-        """
-        Read/write some fields data from/into hdf/xmdf files.
+        """Read/write some fields data from/into hdf/xmdf files.
         Parallel io.
-        @param var_names : a dictionnary of names to connect fields
-        to the dataset in the hdf file. See example below.
-        @param subset : a subset of the domain, on which data are read.
-        It must be a hysop.domain.subset.
 
-        Names paramater example:
-        if variables=[velo, vorti], and if hdf file contains
+
+        Parameters
+        ----------
+        var_names : a dictionnary, optional
+            keys = :class:`~hysop.fields.continuous.Field`,
+            values = string, field name. See notes below.
+        subset : :class:`~hysop.domain.subset.Subset`, optional
+            a subset of the domain, on which data are read or written,
+            default=the whole domain.
+
+        Notes
+        -----
+        Dataset in hdf files are identified with names. In hysop, when
+        writing a file, default dataset name is
+        'continuous_field.name + topo.id + component direction', but
+        this might be changed thanks to var_names argument.
+        For example, if variables=[velo, vorti], and if hdf file contains
         'vel_1_X, vel_1_Y, vel_1_Z, dat_2_X, dat_2_Y, dat_2_Z' keys, then
         use :
-        names = {velo: 'vel', vorti:'dat'} if you want to read vel/dat
+        var_names = {velo: 'vel', vorti:'dat'} if you want to read vel/dat
         into velo/vorti.
-        Mind that when writing an hdf file, dataset names are set as :
-        field.name + topo.id + component direction.
-
-        If names=None, field names will be used as keys to search for
-        dataset.
 
         """
         super(HDF_IO, self).__init__(**kwds)
@@ -65,29 +72,31 @@ class HDF_IO(Computational):
         # If no filename is given, set it to
         # the concatenation of variables'names.
         if self.io_params is None:
+            # if name is not set, name = concatenation of fields names,
+            # like v1_v2.
+            # WARNING FP: we must sort names (alph. order), else
+            # it seems the order may change from one mpi process
+            # to another.
             name = ''
+            names = []
             for var in self.input:
-                if var.name == 'unnamed':
-                    msg = 'Warning : you try to write an unnamed variable'
-                    msg += ' into hdf file. This may result in'
-                    msg += 'unexpected side effects.'
-                    print msg
-                name += var.name
-
-                name += '_'
-            self.io_params = IOParams(name, fileformat=HDF5)
+                names.append(var.name)
+            names.sort()
+            for nn in names:
+                name += nn + '_'
+            name = name[:-1]
+            self.io_params = IOParams(name, fileformat=IO.HDF5)
         else:
-            assert self.io_params.fileformat is HDF5
-
-        ## Set a subset of the original domain
+            assert self.io_params.fileformat is IO.HDF5
+        # Set a subset of the original domain
         self.subset = subset
 
-        ## Dictionnary of names to search in hdf file. May be None.
-        ## It will be checked during setup.
+        # Dictionnary of names to search in hdf file. May be None.
+        # It will be checked during setup.
         self.var_names = var_names
 
         # Local topology, that MUST be common to all variables.
-        self._topology = None
+        self.topology = None
         self._slices = None
         self._global_resolution = None
         self._sl = None
@@ -103,14 +112,14 @@ class HDF_IO(Computational):
     def discretize(self):
         super(HDF_IO, self)._standard_discretize()
         assert self._single_topo, 'Multi-resolution case is not allowed.'
-        self._topology = self.variables.values()[0]
+        self.topology = self.variables.values()[0]
 
         # Discretize the subset, if required
         if self.subset is not None:
-            self.subset.discretize(self._topology)
-            refmesh = self.subset.mesh[self._topology]
+            self.subset.discretize(self.topology)
+            refmesh = self.subset.mesh[self.topology]
         else:
-            refmesh = self._topology.mesh
+            refmesh = self.topology.mesh
         # Global resolution for hdf5 output (warning : this must
         # be the whole domain resolution, not the subset resolution)
         self._global_resolution = list(refmesh.global_resolution())
@@ -136,19 +145,19 @@ class HDF_IO(Computational):
         else:
             for var in self.var_names:
                 # Discrete field associated to var
-                var_d = var.discretize(self._topology)
+                var_d = var.discretize(self.topology)
                 for d in xrange(var_d.nb_components):
                     name = self.var_names[var] + S_DIR[d]
                     self.dataset[name] = var_d.data[d]
 
     def open_hdf(self, count, mode):
         filename = self._get_filename(count)
-        if self._topology.size == 1:
+        if self.topology.size == 1:
             self._hdf_file = h5py.File(filename, mode)
             compression = 'gzip'
         else:
             self._hdf_file = h5py.File(filename, mode, driver='mpio',
-                                       comm=self._topology.comm)
+                                       comm=self.topology.comm)
             compression = None
 
         return compression
@@ -169,31 +178,41 @@ class HDF_Writer(HDF_IO):
         """
         Write some fields data into hdf/xmdf files.
         Parallel writings.
-        @param xmfalways : true if xmf output must be done every time
-        an hdf5 file is created
+
+        Parameters
+        ----------
+        xmfalways : boolean, optional
+            true if xmf output must be done every time
+            an hdf5 file is created (i.e. at each time step),
+            default=True
+        kwds : base class arguments
         """
         super(HDF_Writer, self).__init__(**kwds)
 
         # count the number of calls
         self._count = 0
 
-        self.step = self._step_HDF5
         if xmfalways:
             self.step = self._step_HDF5_XMF
             self.finalize = lambda: 1
         else:
+            self.step = self._step_HDF5
             self.finalize = self.createXMFFile
         self._xdmf_data_files = []
-
         # filename = prefix_N, N = counter value
-        self._get_filename = lambda i: self.io_params.filename + \
-            "_{0:05d}".format(i) + '.h5'
+        self._get_filename = self._input_fname
         # In xdmf file, it is not possible to have two times
         # the same time value. So we save the last value of
         # time where xdmf has been written, to raise an exception
         # if that happens.
         self._last_written_time = None
 
+    def _input_fname(self, i):
+        """Set output file name for current iteration"""
+        msg = 'count < 0, simu must be initialized.'
+        assert i >= 0, msg
+        return self.io_params.filename + "_{0:05d}".format(i) + '.h5'
+
     @debug
     @profile
     @opapply
@@ -206,8 +225,7 @@ class HDF_Writer(HDF_IO):
             self._count += 1
 
     def createXMFFile(self):
-        """
-        Create and fill the xdmf file
+        """Create and fill the xdmf file
         """
 
         if self._mpis.rank == self.io_params.io_leader:
@@ -220,8 +238,8 @@ class HDF_Writer(HDF_IO):
             f.write("CollectionType=\"Temporal\">\n")
             ds_names = self.dataset.keys()
             for i, t in self._xdmf_data_files:
-                f.write(io.XMF.write_grid_attributes(
-                    self._topology, ds_names, i, t, self._get_filename(i),
+                f.write(XMF.write_grid_attributes(
+                    self.topology, ds_names, i, t, self._get_filename(i),
                     self.subset))
             f.write("  </Grid>\n")
             f.write(" </Domain>\n")
@@ -232,17 +250,17 @@ class HDF_Writer(HDF_IO):
 #        self.createXMFFile()
 
     def _step_HDF5(self, simu):
-        """
-        Write a h5 file with data on each mpi process.
+        """Write an h5 file with data on each mpi process.
         """
         # Remarks:
-        # - force np.float64, ParaView seems to not be able to read float32
+        # - force np.float64, ParaView seems unable to read float32
         # - writing compressed hdf5 files (gzip compression seems the best)
         # - gzip compression does not work in parallel.
 
         # Get 'current' filename. It depends on the number
         # of the current output (count) and on the current process
         # rank.
+        self._count = simu.currentIteration
         compression = self.open_hdf(self._count, mode='w')
         # Get the names of output variables and create the corresponding
         # datasets
@@ -275,11 +293,19 @@ class HDF_Reader(HDF_IO):
     Parallel reading of hdf/xdmf files to fill some fields in.
     """
     def __init__(self, restart=None, **kwds):
-        """
-        Read some fields data into hdf/xmdf files.
+        """Read some fields data from hdf/xmdf files.
         Parallel readings.
-        @param restart : set this to read a file corresponding to
-        a specific iteration.
+
+        Parameters
+        ----------
+        restart : int, optional
+            number of a specific iteration to be read, default=None,
+            i.e. read first iteration.
+        kwds : base class arguments
+
+        Notes: restart corresponds to the number which appears in
+        the hdf file name, corresponding to the number of the
+        iteration where writing occured.
         See examples in tests_hdf_io.py
         """
         super(HDF_Reader, self).__init__(**kwds)
@@ -315,10 +341,10 @@ class HDF_Reader(HDF_IO):
         # Do we need a CPU->GPU transfer here?
 
     def dataset_names(self):
-        """
-        Return the list of available names for datasets in
+        """Return the list of available names for datasets in
         the required file.
         """
+        assert self._is_discretized, 'a call to discretize() is required.'
         if self._hdf_file is None:
             self.open_hdf(count=self.restart, mode='r')
         return self._hdf_file.keys()
diff --git a/HySoP/hysop/operator/tests/test_analytic.py b/HySoP/hysop/operator/tests/test_analytic.py
index 3e2878ff9..b9d8a032c 100755
--- a/HySoP/hysop/operator/tests/test_analytic.py
+++ b/HySoP/hysop/operator/tests/test_analytic.py
@@ -1,6 +1,4 @@
-"""
-@file hysop.operator.tests.test_analytic
-Test initialization of fields with analytic formula
+"""Test initialization of fields with analytic formula
 """
 from numpy import allclose
 from hysop.domain.box import Box
@@ -21,8 +19,8 @@ simu = Simulation(tinit=0., tend=0.1, nbIter=1)
 # Non-Vectorized and vectorized formulas for a scalar
 def test_analytical_op_1():
     box = Box()
-    caf = Field(box, formula=func_scal_1)
-    caf2 = Field(box, formula=func_scal_2, vectorize_formula=True)
+    caf = Field(box, formula=func_scal_1, name='caf')
+    caf2 = Field(box, name='caf2', formula=func_scal_2, vectorize_formula=True)
     op = Analytic(variables={caf: d3D})
     op2 = Analytic(variables={caf2: d3D})
     op.discretize()
@@ -31,7 +29,7 @@ def test_analytical_op_1():
     op2.setup()
     topo = op.discreteFields[caf].topology
     coords = topo.mesh.coords
-    ref = Field(box)
+    ref = Field(box, name='ref')
     refd = ref.discretize(topo)
     cafd = caf.discreteFields[topo]
     cafd2 = caf2.discreteFields[topo]
@@ -49,8 +47,9 @@ def test_analytical_op_1():
 # Non-Vectorized and vectorized formulas for a vector
 def test_analytical_op_3():
     box = Box()
-    caf = Field(box, formula=func_vec_1, is_vector=True)
-    caf2 = Field(box, formula=func_vec_2, vectorize_formula=True, is_vector=True)
+    caf = Field(box, name='caf', formula=func_vec_1, is_vector=True)
+    caf2 = Field(box, name='caf', formula=func_vec_2,
+                 vectorize_formula=True, is_vector=True)
     op = Analytic(variables={caf: d3D})
     op2 = Analytic(variables={caf2: d3D})
     op.discretize()
@@ -59,7 +58,7 @@ def test_analytical_op_3():
     op2.setup()
     topo = op.discreteFields[caf].topology
     coords = topo.mesh.coords
-    ref = Field(box, is_vector=True)
+    ref = Field(box, is_vector=True, name='ref')
     refd = ref.discretize(topo)
     cafd = caf.discreteFields[topo]
     cafd2 = caf2.discreteFields[topo]
@@ -81,8 +80,9 @@ def test_analytical_op_3():
 # Non-Vectorized and vectorized formulas for a vector with extra-args
 def test_analytical_op_4():
     box = Box()
-    caf = Field(box, formula=func_vec_3, is_vector=True)
-    caf2 = Field(box, formula=func_vec_4, vectorize_formula=True, is_vector=True)
+    caf = Field(box, formula=func_vec_3, is_vector=True, name='caf')
+    caf2 = Field(box, formula=func_vec_4, vectorize_formula=True,
+                 name='caf2', is_vector=True)
     op = Analytic(variables={caf: d3D})
     op2 = Analytic(variables={caf2: d3D})
     op.discretize()
@@ -91,7 +91,7 @@ def test_analytical_op_4():
     op2.setup()
     topo = op.discreteFields[caf].topology
     coords = topo.mesh.coords
-    ref = Field(box, is_vector=True)
+    ref = Field(box, name='ref', is_vector=True)
     refd = ref.discretize(topo)
     cafd = caf.discreteFields[topo]
     cafd2 = caf2.discreteFields[topo]
@@ -116,13 +116,13 @@ def test_analytical_op_4():
 # Non-Vectorized formula for a nbc components field with extra-args
 def test_analytical_op_5():
     box = Box()
-    caf = Field(box, formula=func_vec_5, nb_components=nbc)
+    caf = Field(box, formula=func_vec_5, nb_components=nbc, name='caf')
     op = Analytic(variables={caf: d3D})
     op.discretize()
     op.setup()
     topo = op.discreteFields[caf].topology
     coords = topo.mesh.coords
-    ref = Field(box, nb_components=nbc)
+    ref = Field(box, nb_components=nbc, name='ref')
     refd = ref.discretize(topo)
     cafd = caf.discreteFields[topo]
     ids = [0, ] * nbc
@@ -140,13 +140,13 @@ def test_analytical_op_5():
 # Non-Vectorized formula for a nbc components field in 2D, with extra-args
 def test_analytical_op_6():
     box = Box(dimension=2, length=L2D, origin=origin2D)
-    caf = Field(box, formula=func_vec_6, nb_components=nbc)
+    caf = Field(box, formula=func_vec_6, nb_components=nbc, name='caf')
     op = Analytic(variables={caf: d2D})
     op.discretize()
     op.setup()
     topo = op.discreteFields[caf].topology
     coords = topo.mesh.coords
-    ref = Field(box, nb_components=nbc)
+    ref = Field(box, nb_components=nbc, name='ref')
     refd = ref.discretize(topo)
     cafd = caf.discreteFields[topo]
     ids = [0, ] * nbc
diff --git a/HySoP/hysop/operator/tests/test_drag_and_lift.py b/HySoP/hysop/operator/tests/test_drag_and_lift.py
index ddda768fc..06fd9fef8 100755
--- a/HySoP/hysop/operator/tests/test_drag_and_lift.py
+++ b/HySoP/hysop/operator/tests/test_drag_and_lift.py
@@ -2,7 +2,8 @@
 from hysop.domain.subsets import Sphere
 from hysop.operator.penalization import PenalizeVorticity
 from hysop.problem.simulation import Simulation
-from hysop.tools.parameters import Discretization, IOParams
+from hysop.tools.parameters import Discretization
+from hysop.tools.io_utils import IOParams, IO
 from hysop.mpi.topology import Cartesian
 import numpy as np
 import os
@@ -54,9 +55,9 @@ def s3d(res, x, y, z, t):
     return res
 
 
-Nx = 128
-Ny = 134
-Nz = 128
+Nx = 32
+Ny = 106
+Nz = 64
 g = 2
 
 
@@ -157,20 +158,24 @@ def test_all_drags():
     dg = {}
     dg['mom'] = MomentumForces(velocity=velo, discretization=topo,
                                obstacles=obst, penalisation_coeff=[1e8],
-                               io_params=IOParams('Heloise'))
+                               io_params=IOParams('Heloise',
+                                                  fileformat=IO.ASCII))
     sdir = [2]
     dg['noca_1'] = NocaForces(velocity=velo, vorticity=vorti,
                               discretization=topo, surfdir=sdir,
                               nu=0.3, obstacles=obst, formulation=1,
-                              io_params=IOParams('Noca1'))
+                              io_params=IOParams('Noca1',
+                                                 fileformat=IO.ASCII))
     dg['noca_2'] = NocaForces(velocity=velo, vorticity=vorti,
                               discretization=topo, surfdir=sdir,
                               nu=0.3, obstacles=obst, formulation=2,
-                              io_params=IOParams('Noca2'))
+                              io_params=IOParams('Noca2',
+                                                 fileformat=IO.ASCII))
     dg['noca_3'] = NocaForces(velocity=velo, vorticity=vorti,
                               discretization=topo, surfdir=sdir,
                               nu=0.3, obstacles=obst, formulation=3,
-                              io_params=IOParams('Noca3'))
+                              io_params=IOParams('Noca3',
+                                                 fileformat=IO.ASCII))
     for drag in dg.values():
         drag.discretize()
         drag.setup()
diff --git a/HySoP/hysop/operator/tests/test_hdf5_io.py b/HySoP/hysop/operator/tests/test_hdf5_io.py
index 0dee976d0..368836fd2 100755
--- a/HySoP/hysop/operator/tests/test_hdf5_io.py
+++ b/HySoP/hysop/operator/tests/test_hdf5_io.py
@@ -1,19 +1,17 @@
 # -*- coding: utf-8 -*-
-"""
-Tests for reader/writer of fields in hdf5 format.
+"""Tests for reader/writer of fields in hdf5 format.
 """
 
 from hysop import Box, Field
 import numpy as np
 import os
-from hysop.constants import HDF5
 from hysop.problem.simulation import Simulation
-import shutil
-from hysop.tools.parameters import Discretization, IOParams
+from hysop.tools.parameters import Discretization
 from hysop.operator.hdf_io import HDF_Writer, HDF_Reader
-from hysop.tools.io_utils import IO
+from hysop.tools.io_utils import IO, IOParams
 from hysop.mpi import main_rank, main_size
 from hysop.domain.subsets import SubBox
+from hysop.testsenv import postclean
 
 Lx = 2.
 nb = 65
@@ -63,27 +61,23 @@ def vort3D(res, x, y, z, t):
     return res
 
 
-def purgeFiles():
-    if main_rank == 0:
-        shutil.rmtree(working_dir)
-
-
+@postclean(working_dir)
 def test_write_read_scalar_3D():
     dom, topo = init1(3)
     scal3D = Field(domain=dom, name='Scal3D')
     scalRef = Field(domain=dom, formula=func3D, name='ScalRef3D')
 
     filename = working_dir + '/testIO_scal'
-    iop = IOParams(filename, fileformat=HDF5)
+    iop = IOParams(filename, fileformat=IO.HDF5)
     op = HDF_Writer(variables={scalRef: topo}, io_params=iop)
     simu = Simulation(nbIter=10)
     op.discretize()
     op.setup()
+    simu.initialize()
 
     scalRef.initialize(simu.time, topo=topo)
     op.apply(simu)
 
-    simu.initialize()
     simu.advance()
     simu.advance()
     # Print scalRef for other iterations
@@ -92,11 +86,11 @@ def test_write_read_scalar_3D():
     fullpath = iop.filename
     assert os.path.exists(fullpath + '.xmf')
     assert os.path.exists(fullpath + '_00000.h5')
-    assert os.path.exists(fullpath + '_00001.h5')
+    assert os.path.exists(fullpath + '_00002.h5')
 
     # Reader
-    iop_read = IOParams(working_dir + '/testIO_scal_00001.h5',
-                         fileformat=HDF5)
+    iop_read = IOParams(working_dir + '/testIO_scal_00002.h5',
+                        fileformat=IO.HDF5)
     reader = HDF_Reader(variables=[scal3D], discretization=topo,
                         io_params=iop_read,
                         var_names={scal3D: 'ScalRef3D_' + str(topo.get_id())})
@@ -115,6 +109,7 @@ def test_write_read_scalar_3D():
         assert np.allclose(scref.data[d][ind], sc3d.data[d][ind])
 
 
+@postclean(working_dir)
 def test_write_read_scalar_3D_defaults():
     dom, topo = init1(3)
     scal3D = Field(domain=dom, name='Scal3D')
@@ -137,9 +132,9 @@ def test_write_read_scalar_3D_defaults():
     filename = scalRef.name
     fullpath = os.path.join(IO.default_path(), filename)
 
-    assert os.path.exists(fullpath + '_.xmf')
-    assert os.path.exists(fullpath + '__00000.h5')
-    assert os.path.exists(fullpath + '__00001.h5')
+    assert os.path.exists(fullpath + '.xmf')
+    assert os.path.exists(fullpath + '_00000.h5')
+    assert os.path.exists(fullpath + '_00001.h5')
 
     sc3d = scal3D.discretize(topo)
     scref = scalRef.discretize(topo)
@@ -163,7 +158,7 @@ def test_write_read_scalar_3D_defaults():
     for d in xrange(scal3D.nb_components):
         assert np.allclose(scref.data[d][ind], sc3d.data[d][ind])
 
-
+@postclean(working_dir)
 def test_write_read_vectors_3D_defaults():
     dom, topo = init2()
     velo = Field(domain=dom, formula=vec3D, name='velo', is_vector=True)
@@ -184,9 +179,13 @@ def test_write_read_vectors_3D_defaults():
 
     op.finalize()
     filename = ''
-    for v in op.input:
-        filename += v.name
-        filename += '_'
+    names = []
+    for var in op.input:
+        names.append(var.name)
+        names.sort()
+    for nn in names:
+        filename += nn + '_'
+    filename = filename[:-1]
     fullpath = os.path.join(IO.default_path(), filename)
 
     assert os.path.exists(fullpath + '.xmf')
@@ -197,6 +196,8 @@ def test_write_read_vectors_3D_defaults():
     w3d = vorti.discretize(topo)
     ind = topo.mesh.iCompute
 
+    # Copy current values of v3 and w3 into buff1 and buff2, for comparison
+    # after reader.apply, below.
     buff1 = Field(domain=dom, name='buff1', is_vector=True)
     buff2 = Field(domain=dom, name='buff2', is_vector=True)
     b1 = buff1.discretize(topo=topo)
@@ -204,13 +205,16 @@ def test_write_read_vectors_3D_defaults():
     for d in xrange(velo.nb_components):
         b1.data[d][...] = v3d.data[d][...]
         b2.data[d][...] = w3d.data[d][...]
-
+        # reset v3 and w3 to zero.
         v3d.data[d][...] = 0.0
         w3d.data[d][...] = 0.0
 
-    # Read vector fields, using default configuration for output
+    # Read vector fields, using default configuration for input
     # names and location, with a given iteration number.
+    # If everything works fine, reader must read output from
+    # the writer above.
     reader = HDF_Reader(variables={velo: topo, vorti: topo},
+                        io_params=IOParams(filename),
                         restart=simu.currentIteration - 1)
     reader.discretize()
     reader.setup()
@@ -220,11 +224,13 @@ def test_write_read_vectors_3D_defaults():
 
     reader.apply()
     reader.finalize()
+    # Now, v3 and w3 (just read) must be equal to saved values in b1 and b2.
     for d in xrange(v3d.nb_components):
         assert np.allclose(b1.data[d][ind], v3d.data[d][ind])
         assert np.allclose(b2.data[d][ind], w3d.data[d][ind])
 
 
+@postclean(working_dir)
 def test_write_read_vectors_3D():
     dom, topo = init2()
     velo = Field(domain=dom, formula=vec3D, name='velo', is_vector=True)
@@ -233,7 +239,7 @@ def test_write_read_vectors_3D():
     # Write a vector field, using default for output location
     # but with fixed names for datasets
     filename = working_dir + '/testIO_vec'
-    iop = IOParams(filename, fileformat=HDF5)
+    iop = IOParams(filename, fileformat=IO.HDF5)
     op = HDF_Writer(variables={velo: topo, vorti: topo},
                     var_names={velo: 'io_1', vorti: 'io_2'}, io_params=iop)
     simu = Simulation(nbIter=3)
@@ -267,7 +273,7 @@ def test_write_read_vectors_3D():
 
     # Read vector fields, fixed filename, fixed dataset names.
     iop_read = IOParams(working_dir + '/testIO_vec_00001.h5',
-                         fileformat=HDF5)
+                        fileformat=IO.HDF5)
     reader = HDF_Reader(variables={buff1: topo, buff2: topo},
                         io_params=iop_read,
                         var_names={buff1: 'io_2', buff2: 'io_1'})
@@ -282,6 +288,7 @@ def test_write_read_vectors_3D():
         assert np.allclose(b1.data[d][ind], w3d.data[d][ind])
 
 
+@postclean(working_dir)
 def test_write_read_subset_1():
     dom, topo = init2()
     velo = Field(domain=dom, formula=vec3D, name='velo', is_vector=True)
@@ -308,6 +315,7 @@ def test_write_read_subset_1():
     for v in op.input:
         filename += v.name
         filename += '_'
+    filename = filename[:-1]
     fullpath = os.path.join(IO.default_path(), filename)
 
     assert os.path.exists(fullpath + '.xmf')
@@ -321,7 +329,7 @@ def test_write_read_subset_1():
     buff1 = Field(domain=dom, name='buff1', is_vector=True)
 
     # Read vector fields, fixed filename, fixed dataset names.
-    iop = IOParams(filename + '_00000.h5', fileformat=HDF5)
+    iop = IOParams(filename + '_00000.h5', fileformat=IO.HDF5)
     reader = HDF_Reader(variables={buff1: topo},
                         io_params=iop,
                         var_names={buff1: 'io_1'}, subset=mybox)
@@ -335,6 +343,7 @@ def test_write_read_subset_1():
         assert np.allclose(b1.data[d][indsubset], v3d.data[d][indsubset])
 
 
+@postclean(working_dir)
 def test_write_read_subset_2():
     dom, topo = init2()
     velo = Field(domain=dom, formula=vec3D, name='velo', is_vector=True)
@@ -361,6 +370,7 @@ def test_write_read_subset_2():
     for v in op.input:
         filename += v.name
         filename += '_'
+    filename = filename[:-1]
     fullpath = os.path.join(IO.default_path(), filename)
 
     assert os.path.exists(fullpath + '.xmf')
@@ -374,7 +384,7 @@ def test_write_read_subset_2():
     buff1 = Field(domain=dom, name='buff1', is_vector=True)
 
     # Read vector fields, fixed filename, fixed dataset names.
-    iop = IOParams(filename + '_00000.h5', fileformat=HDF5)
+    iop = IOParams(filename + '_00000.h5', fileformat=IO.HDF5)
     reader = HDF_Reader(variables={buff1: topo},
                         io_params=iop,
                         var_names={buff1: 'io_1'}, subset=mybox)
diff --git a/HySoP/hysop/operator/tests/test_multiphase_gradp.py b/HySoP/hysop/operator/tests/test_multiphase_gradp.py
index 43aac523e..477c9e48c 100755
--- a/HySoP/hysop/operator/tests/test_multiphase_gradp.py
+++ b/HySoP/hysop/operator/tests/test_multiphase_gradp.py
@@ -50,9 +50,9 @@ def test_gradp():
     """Testing gradp operator against analytical result"""
     simu = Simulation(tinit=0.0, tend=1.0, timeStep=0.1, iterMax=1)
     box = Box()
-    velo = Field(box, is_vector=True, formula=velo_func)
-    res = Field(box, is_vector=True)
-    true_res = Field(box, is_vector=True, formula=res_func)
+    velo = Field(box, is_vector=True, formula=velo_func, name='v0')
+    res = Field(box, is_vector=True, name='res')
+    true_res = Field(box, is_vector=True, formula=res_func, name='vres')
     d = Discretization([129,129,129])
     op = MultiphaseGradP(velocity=velo, gradp=res, viscosity=VISCOSITY,
                          discretization=d)
diff --git a/HySoP/hysop/operator/tests/test_multiresolution_filter.py b/HySoP/hysop/operator/tests/test_multiresolution_filter.py
index 6d66f2513..91febfd77 100755
--- a/HySoP/hysop/operator/tests/test_multiresolution_filter.py
+++ b/HySoP/hysop/operator/tests/test_multiresolution_filter.py
@@ -37,7 +37,7 @@ def func_periodic_XY(res, x, y, z, t=0):
 
 def filter(d_fine, d_coarse, func, method, atol=1e-8, rtol=1e-5):
     box = Box(length=L, origin=O)
-    f = Field(box, formula=func)
+    f = Field(box, formula=func, name='f0')
     op = MultiresolutionFilter(d_in=d_fine, d_out=d_coarse,
                                variables={f: d_coarse},
                                method=method)
@@ -169,7 +169,7 @@ def func(res, x, y, z, t=0):
 def test_filter_linear():
     """This test compares the GPU linear filter with python implementation"""
     box = Box(length=L, origin=O)
-    f = Field(box, formula=func)
+    f = Field(box, formula=func, name='f0')
     d_fine = Discretization([513, 513, 513])
     d_coarse = Discretization([257, 257, 257], ghosts=[1, 1, 1])
     op = MultiresolutionFilter(d_in=d_fine, d_out=d_coarse,
@@ -195,7 +195,7 @@ def test_filter_linear():
 def test_filter_l2_1():
     """This test compares the GPU linear filter with python implementation"""
     box = Box(length=L, origin=O)
-    f = Field(box, formula=func)
+    f = Field(box, formula=func, name='f0')
     d_fine = Discretization([513, 513, 513])
     d_coarse = Discretization([257, 257, 257], ghosts=[2, 2, 2])
     op = MultiresolutionFilter(d_in=d_fine, d_out=d_coarse,
diff --git a/HySoP/hysop/operator/tests/test_particle_advection.py b/HySoP/hysop/operator/tests/test_particle_advection.py
index 7ca871bea..ad6fc2ef0 100755
--- a/HySoP/hysop/operator/tests/test_particle_advection.py
+++ b/HySoP/hysop/operator/tests/test_particle_advection.py
@@ -1,6 +1,4 @@
-"""
-@file hysop.operator.tests.test_particle_advection
-Testing pure python particle advection with null velocity.
+"""Testing pure python particle advection with null velocity.
 """
 from hysop.domain.box import Box
 from hysop.fields.continuous import Field
@@ -11,7 +9,6 @@ import hysop.tools.numpywrappers as npw
 import numpy as np
 
 
-
 d2d = Discretization([17, 17])
 d3d = Discretization([17, 17, 17])
 
@@ -42,7 +39,9 @@ def setup_list_2D():
     return [scal1, scal2], velo
 
 
-def setup_3D():
+def setup_3d():
+    """Build fields inside a 3d domain
+    """
     box = Box(length=[2., 4., 1.], origin=[-1., -2., 0.])
     scal = Field(domain=box, name='Scalar')
     velo = Field(domain=box, name='Velocity',
@@ -150,7 +149,7 @@ def assertion_list(scal, advec):
 
 
 def test_nullVelocity_2D():
-    """
+    """2D case, advection with null velocity, single resolution.
     """
     scal, velo = setup_2D()
 
@@ -180,7 +179,7 @@ def test_nullVelocity_list_2D():
 def test_nullVelocity_3D():
     """
     """
-    scal, velo = setup_3D()
+    scal, velo = setup_3d()
 
     advec = Advection(velo, scal, discretization=d3d)
     assert assertion(scal, advec)
diff --git a/HySoP/hysop/operator/tests/test_penalization.py b/HySoP/hysop/operator/tests/test_penalization.py
index d26e3b36b..9026e986b 100755
--- a/HySoP/hysop/operator/tests/test_penalization.py
+++ b/HySoP/hysop/operator/tests/test_penalization.py
@@ -3,7 +3,8 @@ from hysop.domain.subsets import HemiSphere, Sphere, Cylinder
 from hysop.domain.porous import Porous
 from hysop.operator.penalization import Penalization, PenalizeVorticity
 from hysop.problem.simulation import Simulation
-from hysop.tools.parameters import Discretization, IOParams
+from hysop.tools.parameters import Discretization
+from hysop.tools.io_utils import IOParams
 from hysop.mpi.topology import Cartesian
 import hysop.tools.numpywrappers as npw
 import numpy as np
diff --git a/HySoP/hysop/operator/tests/test_spectrum.py b/HySoP/hysop/operator/tests/test_spectrum.py
old mode 100644
new mode 100755
index 9de2e6ce7..4acc7264c
--- a/HySoP/hysop/operator/tests/test_spectrum.py
+++ b/HySoP/hysop/operator/tests/test_spectrum.py
@@ -1,10 +1,11 @@
 from hysop.domain.box import Box
-from hysop.constants import np, HDF5
 from hysop.fields.continuous import Field
 from hysop.operator.spectrum import Spectrum
-from hysop.tools.parameters import Discretization, IOParams
+from hysop.tools.parameters import Discretization
 from hysop.problem.simulation import Simulation
 from hysop.operator.hdf_io import HDF_Writer
+import numpy as np
+from hysop.tools.io_utils import IO
 pi = np.pi
 sin = np.sin
 cos = np.cos
diff --git a/HySoP/hysop/problem/tests/test_simulation.py b/HySoP/hysop/problem/tests/test_simulation.py
old mode 100644
new mode 100755
index 65f7ce6c2..8a7ae1264
--- a/HySoP/hysop/problem/tests/test_simulation.py
+++ b/HySoP/hysop/problem/tests/test_simulation.py
@@ -3,14 +3,14 @@
 tests simulation incr and io_utils writer
 """
 from hysop.problem.simulation import Simulation
-from hysop.tools.parameters import IOParams
-from hysop.tools.io_utils import Writer
+from hysop.tools.io_utils import Writer, IOParams, IO
 
 simu = Simulation(tinit=0.0, tend=1.0, nbIter=10)
 
 
 def test_simu_incr():
-    io_params = IOParams(filename='temp_test', frequency=2)
+    io_params = IOParams(filename='temp_test', frequency=2,
+                         fileformat=IO.ASCII)
     wr = Writer(io_params)
     assert wr.do_write(simu.currentIteration)
 
@@ -33,7 +33,8 @@ def test_simu_incr():
 
 
 def test_simu_incr2():
-    io_params = IOParams(filename='temp_test', frequency=3)
+    io_params = IOParams(filename='temp_test', frequency=3,
+                         fileformat=IO.ASCII)
     wr = Writer(io_params)
     assert wr.do_write(simu.currentIteration)
     simu.timeStep = 0.10000000001
diff --git a/HySoP/hysop/tools/io_utils.py b/HySoP/hysop/tools/io_utils.py
old mode 100644
new mode 100755
index 45dc5a4da..f06ce3cd5
--- a/HySoP/hysop/tools/io_utils.py
+++ b/HySoP/hysop/tools/io_utils.py
@@ -1,21 +1,39 @@
 """Tools related to i/o in HySoP.
+
+.. currentmodule hysop.tools.io_utils
+
+* :class:`~IO`
+* :class:`~IOParams`
+* :class:`~Writer`
+* :class:`~XMF`, tools to prepare/write xmf files.
+
 """
 import os
 import scitools.filetable as ft
 import hysop.tools.numpywrappers as npw
 from inspect import getouterframes, currentframe
 import hysop.mpi as mpi
-from re import findall, IGNORECASE
-from hysop.tools.parameters import MPIParams, IOParams
+from re import findall
+from hysop.tools.parameters import MPIParams
+from collections import namedtuple
+from hysop.constants import S_DIR
+import h5py
 
 
 class IO(object):
     """
-    Static class with utilities for file i/o
+    Static class with utilities to set/find the place where i/o files
+    will be read/written.
     """
 
     _default_path = None
 
+    HDF5 = 998
+    """HDF5 format id"""
+
+    ASCII = 997
+    """ascii format id"""
+
     @staticmethod
     def default_path():
         """Get the current default path used for io.
@@ -24,7 +42,7 @@ class IO(object):
         -------
         string
             the default value of the current i/o path. If not set,
-            the basename of the top exec. file is used.
+            the basename of the top exec file is used.
 
         Examples
         --------
@@ -33,46 +51,72 @@ class IO(object):
         results will be saved in ./MyExe/p4.
         The number after 'p' will be the number of mpi processus
         set for the simulation.
-
-
         """
         # Memo for getouterframe usage:
         #    frame, filename, line_number, function_name, lines, index =\
         # inspect.getouterframes(inspect.currentframe())[-1]
         # Warning FP : the behavior of python and ipython is different for
-        # this command.
+        # this command. Moreover, output differs when called from interactive
+        # session.
+
         if IO._default_path is not None:
+            # which must be the case for any interactive session
+            # due to call of set_default_path in hysop.__init__.py.
             return IO._default_path
 
         a = getouterframes(currentframe())
-        ind = -1
+        #ind = -1
         interactive_path = './interactive/p' + str(mpi.main_size)
         interactive_path = os.path.abspath(interactive_path)
         # --- ipython ---
-        from hysop.tools.sys_utils import SysUtils
-        if SysUtils.in_ipython():
-            sublist = [i[1] for i in a]
-            for val in sublist:
-                ll = findall('ipython', val, IGNORECASE)
-                if len(ll) > 0:
-                    ind = sublist.index(val) - 1
-                    break
-            if ind > -1:
-                # -- interactive ipython but call with execfile--
-                if len(findall('io_utils', a[ind][1])) > 0:
-                    return interactive_path
-                a = a[ind]
-            else:
-                # -- interactive ipython without execfile call --
-                return interactive_path
-
-        else:
-            # -- python --
-            a = a[-1]
-            if a[-1] is None or len(findall('py.test', a[1])) > 0:
-                # interactive python
+        #from hysop.tools.sys_utils import SysUtils
+        # if SysUtils.in_ipython():
+        #     # Note FP: because of set_default_path call
+        #     # in __init__.py, this condition must never happen.
+        #     # But we keep the code below, just in case ...
+
+        #     # list of files (fullpath) which contain the callers
+        #     sublist = [i[1] for i in a]
+        #     # look for ipython in callers ...
+        #     # If found, keep index of the file just before
+        #     # first occurence of ipython, i.e. the name
+        #     # of the 'main' file
+        #     for val in sublist:
+        #         ll = findall('ipython', val, IGNORECASE)
+        #         if len(ll) > 0:
+        #             ind = sublist.index(val) - 1
+        #             break
+
+        #     if ind > -1:
+        #         # -- interactive ipython but call with execfile--
+        #         if len(findall('io_utils', a[ind][1])) > 0:
+        #             return interactive_path
+        #         a = a[ind]
+        #     else:
+        #         # -- interactive ipython without execfile call --
+        #         return interactive_path
+
+        # else:
+        # -- python --
+        # if test session, set default path to interactive_path
+        for fname in a:
+            cond1 = len(findall('py.test', fname[1])) > 0
+            cond2 = len(findall('pytest', fname[1])) > 0
+            if cond1 or cond2:
+                IO._default_path = interactive_path
                 return interactive_path
 
+        # else (not in tests)
+        # get name of the file which contains the last caller
+        # and check if not None
+        a = a[-1]
+        if a[-1] is None:
+            # interactive python
+            IO._default_path = interactive_path
+            return interactive_path
+        # Finally (if not interactive, not in tests ...)
+        # use root name of the file which contains
+        # the last caller
         apath = os.path.abspath(os.path.dirname(a[1]))
         sphinxbuild = findall('sphinx-build', a[1])
         if len(sphinxbuild) > 0:
@@ -80,14 +124,14 @@ class IO(object):
         else:
             a = os.path.basename(a[1]).split('.')[-2]
             if a.find('__init__') != -1:
+                IO._default_path = interactive_path
                 return interactive_path
         a = os.path.join(apath, a)
         return os.path.join(a, 'p' + str(mpi.main_size))
 
     @staticmethod
     def check_dir(filename, io_rank=0, comm=None):
-        """
-        Check if the directory of 'filename' exists and creates it if not.
+        """Check if the directory of 'filename' exists and creates it if not.
 
         Parameters
         -----------
@@ -97,8 +141,6 @@ class IO(object):
             processus rank that does the check.
         comm : mpi communicator
             the mpi communicator that does the check.
-
-
         """
         # Create output dir if required
         if comm is None:
@@ -117,19 +159,103 @@ class IO(object):
         pathdir : string
             the new path
 
-
         Notes
         ------
         pN will be add to path name, N being the number of MPI process
         used for the simulation.
 
-
         """
         IO._default_path = pathdir
         IO._default_path = os.path.join(IO._default_path,
                                         'p' + str(mpi.main_size))
         IO.check_dir(IO._default_path)
 
+    @staticmethod
+    def set_datasetname(field_name, topo, direction=None):
+        """Return the dataset name of a given continuous field,
+        saved for a given topology
+        """
+        val = field_name + '_' + str(topo.get_id())
+        if direction is not None:
+            val += S_DIR[direction]
+        return val
+
+    @staticmethod
+    def get_datasetnames(filename):
+        """Return the list of dataset names present
+        in hdf input file
+
+        Parameters
+        ----------
+        filename : string
+            hdf file
+
+        Returns
+        -------
+            a list of strings
+        """
+        hdf_file = h5py.File(filename, 'r')
+        keys = hdf_file.keys()
+        hdf_file.close()
+        return keys
+
+
+class IOParams(namedtuple("IOParams", ['filename', 'filepath',
+                                       'frequency', 'fileformat',
+                                       'io_leader'])):
+    """
+    A struct to handle I/O files parameters
+
+    Parameters
+    -----------
+    filename : string
+        name of the file (absolute or relative path)
+    filepath : string
+        location of the file
+    frequency : int
+        frequency of output or input (e.g. every N times steps)
+    fileformat : int
+        format of the file. See notes for available format. Default=HDF5.
+    io_leader : int
+        rank of the mpi process dealing with the io. Default is 0.
+
+    See examples in hysop.operator.hdf_io
+
+    Notes
+    -----
+    Format parameter must be one of the following :
+      - :class:`~IO.HDF5`
+      - :class:`~IO.ASCII`
+
+    """
+    def __new__(cls, filename, filepath=None, frequency=1,
+                fileformat=None, io_leader=0):
+
+        # Filename is absolute path, filepath arg is ignored.
+        if os.path.isabs(filename):
+            filepath = os.path.dirname(filename)
+
+        else:
+            if filepath is not None:
+                filename = os.path.join(filepath, filename)
+                filepath = os.path.abspath(os.path.dirname(filename))
+            else:
+                filepath = os.path.dirname(filename)
+                if filepath == '':
+                    # Get default output path
+                    filepath = IO.default_path()
+                    filename = os.path.join(filepath, filename)
+                else:
+                    filepath = os.path.abspath(filepath)
+                    filename = os.path.join(filepath,
+                                            os.path.basename(filename))
+        if fileformat is None:
+            fileformat = IO.HDF5
+
+        IO.check_dir(filename)
+        return super(IOParams, cls).__new__(cls, filename, filepath,
+                                            frequency, fileformat, io_leader)
+
 
 class Writer(object):
     """
@@ -138,8 +264,8 @@ class Writer(object):
     Examples
     --------
 
-    >>> from hysop.tools.parameters import IOParams
-    >>> params = IOParams(filename='r.dat')
+    >>> from hysop.tools.io_utils import IOParams, IO
+    >>> params = IOParams(filename='r.dat', fileformat=IO.ASCII)
     >>> wr = Writer(params, buffshape=(1, 2))
     >>> ite = 3 # current iteration number
     >>> if wr.do_write(ite):
@@ -147,10 +273,7 @@ class Writer(object):
     ...    wr.write()
     >>> wr.finalize()
 
-
-    result : buffer is written into r.dat.
-
-
+    result : buffer is written into r.dat
     """
     def __init__(self, io_params, buffshape=None, mpi_params=None,
                  safe_io=True):
@@ -158,7 +281,7 @@ class Writer(object):
 
         Parameters
         ----------
-        io_params : hysop.tools.parameters.IOParams
+        io_params : hysop.tools.io_utils.IOParams
             setup for file ouput (name, location ...)
         buffshape : tuple
             2D numpy.array.shape like tuple, shape of the output/input buffer.
@@ -169,30 +292,27 @@ class Writer(object):
             False --> open at init and close during finalize.
             Cost less but if simu crashes, data are lost.
 
-
         """
         # Absolute path + name for i/o file
         # Note that if filename contains absolute path
         # filepath is ignored
         msg = 'wrong type for io_params arg.'
         assert isinstance(io_params, IOParams), msg
-        self.filename = io_params.filename
-        " I/O file name "
-
-        # How often i/o must occur
-        self._frequency = io_params.frequency
-
-        # Rank of the mpi process that performs i/o (if any)
-        self._io_rank = io_params.io_leader
+        assert io_params.fileformat == IO.ASCII
+        self.io_params = io_params
 
         # A reference communicator, just to identify a
         # process rank for io.
         if mpi_params is None:
             mpi_params = MPIParams()
+        else:
+            msg = 'wrong type for mpi_params arg.'
+            assert isinstance(mpi_params, MPIParams), msg
         self._mpis = mpi_params
 
         # check if output dir exists, create it if not.
-        IO.check_dir(self.filename, self._io_rank, self._mpis.comm)
+        IO.check_dir(self.io_params.filename, self.io_params.io_leader,
+                     self._mpis.comm)
 
         # Shape of the output buffer (must be a 2D numpy array)
         if buffshape is None:
@@ -204,7 +324,7 @@ class Writer(object):
         # The buffer (numpy array) that will be printed to a file
         self.buffer = npw.zeros(self._buffshape)
         " buffer used to save printed data"
-       
+
         # Defines how often
         # output file is written :
         # True --> open/close file everytime
@@ -220,12 +340,11 @@ class Writer(object):
         # Force synchro to be sure that all output dirs
         # have been created.
         self._mpis.comm.barrier()
-        if self._mpis.rank == self._io_rank:
-            self._file = open(self.filename, 'w')
+        if self._mpis.rank == self.io_params.io_leader:
+            self._file = open(self.io_params.filename, 'w')
 
     def do_write(self, ite):
-        """
-        True if output is required
+        """Returns true if output is required
         for iteration ite
 
         Parameters
@@ -233,37 +352,41 @@ class Writer(object):
         ite : int
             current iteration number
 
-
         """
         num = ite + 1  # Number of iterations done
-        return self._mpis.rank == self._io_rank and (num % self._frequency) == 0
+        rk = self._mpis.rank
+        return rk == self.io_params.io_leader and \
+            (num % self.io_params.frequency) == 0
 
     def _fullwrite(self):
-        self._file = open(self.filename, 'a')
+        """open, write and close"""
+        self._file = open(self.io_params.filename, 'a')
         ft.write(self._file, self.buffer)
         self._file.close()
 
     def _partialwrite(self):
+        """just write, no open, nor close"""
         ft.write(self._file, self.buffer)
 
     def finalize(self):
-        if self._mpis.rank == self._io_rank:
+        """close, if required"""
+        if self._mpis.rank == self.io_params.io_leader:
             if not self._file.closed:
                 self._file.close()
 
     def __str__(self):
-        if self._mpis.rank == self._io_rank:
+        if self._mpis.rank == self.io_params.io_leader:
             s = ' === Writer === \n'
-            s += ' - filename = ' + self.filename
+            s += ' - filename = ' + self.io_params.filename
             s += '\n - buffshape = ' + str(self._buffshape)
-            s += '\n - frequ = ' + str(self._frequency)
+            s += '\n - frequ = ' + str(self.io_params.frequency)
             return s
 
 
 class XMF(object):
+    """Static class - Tools to prepare and write xmf file
     """
-    Static class used to define xmf tools.
-    """
+
     @staticmethod
     def _list_format(l):
         """Format a list to the xml output.
@@ -273,7 +396,6 @@ class XMF(object):
         ----------
         l : list to format
 
-
         """
         buff = str(l).replace(',', ' ').replace('[', '')
         return buff.replace(']', '').replace('(', '').replace(')', '')
diff --git a/HySoP/hysop/tools/parameters.py b/HySoP/hysop/tools/parameters.py
old mode 100644
new mode 100755
index 860c639c7..3977a3424
--- a/HySoP/hysop/tools/parameters.py
+++ b/HySoP/hysop/tools/parameters.py
@@ -2,11 +2,15 @@
 
 .. currentmodule hysop.tools
 
+* :class:`~MPIParams`
+* :class:`~Discretization`
+
 """
 
 from collections import namedtuple
 from hysop.mpi.main_var import main_comm, main_rank, MPI
 from hysop.constants import DEFAULT_TASK_ID
+import hysop.tools.numpywrappers as npw
 
 
 class MPIParams(namedtuple('MPIParams', ['comm', 'task_id',
@@ -44,9 +48,6 @@ class MPIParams(namedtuple('MPIParams', ['comm', 'task_id',
                                              rank, on_task)
 
 
-import hysop.tools.numpywrappers as npw
-
-
 class Discretization(namedtuple("Discretization", ['resolution', 'ghosts'])):
     """
     A struct to handle discretization parameters:
@@ -77,65 +78,3 @@ class Discretization(namedtuple("Discretization", ['resolution', 'ghosts'])):
         if result is NotImplemented:
             return result
         return not result
-
-from hysop.constants import HDF5
-import os
-
-
-class IOParams(namedtuple("IOParams", ['filename', 'filepath',
-                                         'frequency', 'fileformat',
-                                         'io_leader'])):
-    """
-    A struct to handle I/O files parameters
-
-    Parameters
-    -----------
-    filename : string
-        name of the file (absolute or relative path)
-    filepath : string
-        location of the file
-    frequency : int
-        frequency of output or input (e.g. every N times steps)
-    fileformat : int
-        format of the file. See notes for available format. Default=HDF5.
-    io_leader : int
-        rank of the mpi process dealing with the io. Default is 0.
-
-
-    See examples in hysop.operator.hdf_io
-
-    Notes
-    -----
-    Format parameter must be one of the following :
-      - hysop.constants.HDF5
-      - hysop.constants.ASCII
-
-    """
-    def __new__(cls, filename, filepath=None, frequency=1,
-                fileformat=None, io_leader=0):
-        from hysop.tools.io_utils import IO
-
-        # Filename is absolute path, filepath arg is ignored.
-        if os.path.isabs(filename):
-            filepath = os.path.dirname(filename)
-
-        else:
-            if filepath is not None:
-                filename = os.path.join(filepath, filename)
-                filepath = os.path.abspath(os.path.dirname(filename))
-            else:
-                filepath = os.path.dirname(filename)
-                if filepath == '':
-                    # Get default output path
-                    filepath = IO.default_path()
-                    filename = os.path.join(filepath, filename)
-                else:
-                    filepath = os.path.abspath(filepath)
-                    filename = os.path.join(filepath,
-                                            os.path.basename(filename))
-        if fileformat is None:
-            fileformat = HDF5
-
-        IO.check_dir(filename)
-        return super(IOParams, cls).__new__(cls, filename, filepath,
-                                            frequency, fileformat, io_leader)
diff --git a/HySoP/hysop/tools/tests/test_parameters.py b/HySoP/hysop/tools/tests/test_parameters.py
index f13a18f3b..e83d4b9d2 100755
--- a/HySoP/hysop/tools/tests/test_parameters.py
+++ b/HySoP/hysop/tools/tests/test_parameters.py
@@ -1,10 +1,8 @@
 # -*- coding: utf-8 -*-
-"""
-Tests for hysop parameters-like variables.
+"""Tests for hysop io
 """
 
-from hysop.tools.parameters import IOParams
-from hysop.tools.io_utils import IO
+from hysop.tools.io_utils import IO, IOParams
 import os
 
 
diff --git a/HySoP/src/fftw/fft3d.f90 b/HySoP/src/fftw/fft3d.f90
index 280b93c02..b9b5d922b 100755
--- a/HySoP/src/fftw/fft3d.f90
+++ b/HySoP/src/fftw/fft3d.f90
@@ -1072,8 +1072,8 @@ contains
   !> Compute spectrum of the given data
   subroutine filter_spectrum_3d(spectrum,wavelengths,length)
 
-    real(mk),dimension(:),intent(out) :: spectrum
-    real(mk),dimension(:),intent(out) :: wavelengths
+    real(mk),dimension(:),intent(inout) :: spectrum
+    real(mk),dimension(:),intent(inout) :: wavelengths
     real(mk),intent(in) :: length
     integer(C_INTPTR_T) :: i,j,k
     real(mk) :: c
-- 
GitLab