diff --git a/docs/sphinx/users_guide/fields.rst b/docs/sphinx/users_guide/fields.rst
index ae4729ba07831ee578e0fd7869522f29c599ce7f..f89a5f30d60d4bdc9436383d1ea200b9fd0f49d0 100644
--- a/docs/sphinx/users_guide/fields.rst
+++ b/docs/sphinx/users_guide/fields.rst
@@ -1,207 +1,367 @@
.. _fields:
-Scalar and Vector Fields
-========================
+Continuous and discrete fields
+==============================
-THIS SECTION MUST BE REVIEWED
+.. currentmodule:: hysop.fields
- Tools to build and use 1, 2 or 3D fields (scalar, vector or
- any number of components.)
- Fields are used to described the physical continuous variables involved in
- a problem on a given domain.
+.. contents::
+
+Definition
+----------
+
+In HySoP, a **continuous field** (:class:`~hysop.fields.continuous.Field`) is an abstract object which represents the usual vector field, in mathematical sense, i.e some function which associates a scalar or a vector to each point of the space.
+Such objects are used as input and/or output variables for continuous operators and must be defined with at least:
+
+* a :class:`~hysop.domain.domain.Domain` : the physical domain of definition of the field,
+* a name, compulsory (required for i/o).
+
+Scalar fields are built by default. Use is_vector or nb_component optional arguments to build
+n-components fields (n=domain.dimension for vector fields)::
+
+ from hysop import Field, Box
+
+ # build a domain
+ box = Box()
+ # Scalar field
+ scal = Field(domain=box, name='s1')
+ # Vector field
+ vec = Field(domain=box, name='v1', is_vector=True)
+ # n-component field
+ n = 4
+ vec_n = Field(domain=box, name='vn', nb_components=n)
+ # nb_component = domain.dimension is equivalent to is_vector=True
+
+Obviously one needs to associate one or more topologies (:class:`hysop.mpi.topology.Cartesian`, see :ref:`topologies`) to such a field, to be able to apply operators or numerical methods. Remind that 'topologies' means either the definition of a global space resolution and of the way data are distributed through mpi processes::
+
+ # The continuous field
+ vec = Field(domain=box, name='vn', is_vector=True)
+ # space discretisation and data distribution description
+ d3d = Discretization([65, 65, 65])
+ topo = box.create_topology(d3d)
+ # build a discretisation of vec on d3d and distribute its data according to 'topo'
+ vd = vec.discretize(topo)
+ vd is vec.discreteFields[topo]
+ # --> returns True
+
+Each time 'discretize' function is called, a new entry is added to the continuous field attribute named discreteFields, which maps a topology to a :class:`~hysop.fields.discrete.DiscreteField`. Therefore, a discrete field will handle local arrays of values
+of the field on the local mesh.
+
+To summarize, a **discrete field** (:class:`~hysop.fields.discrete.DiscreteField) is the association of:
+
+* a topology (:class:`hysop.mpi.topology.Cartesian`~),
+* a list of numpy arrays (one per component), values of the field on the local mesh.
+
+For example::
+
+ # The continuous field --> no memory allocation
+ vec = Field(domain=box, name='vn', is_vector=True)
+ # space discretisation and data distribution description
+ d3d = Discretization([65, 65, 65])
+ topo = box.create_topology(d3d)
+ d3d_with_ghosts = Discretization([65, 33, 65], [2, 2, 2])
+ topo_g = box.create_topology(d3d_with_ghosts)
+ # build a discretisation of vec on d3d and distribute its data according to 'topo'
+ vd = vec.discretize(topo)
+ vd is vec.discreteFields[topo]
+ # --> returns True
+ # another data distribution
+ vd2 = vec.discretize(topo_g)
+
+ for i in xrange(vec.nb_components):
+ print vd.data[i].shape
+ print vd2.data[i].shape
- To define a continuous field, one must provide :
- - the domain on which it is defined
- and optionnaly:
- - its number of components
- - a user defined function formula used to initialize the field
- - a name
+will display the local resolution of each component of the discretisations 'topo' and 'topo_g' of vec.
+Their values will depend either on the global resolution, the number of mpi processes on which the program run and
+on the current process rank.
- A continuous field may have several discretisations (i.e. may
- be associated with different hysop topologies).
-
-How to define and discretize a field
--------------------------------------
-.. currentmodule:: hysop.fields
-
-**classes:** :class:`continuous.Field`, :class:`discrete.DiscreteField`
+Notice that the exact behavior of::
-First of all, one must have defined a domain :Code::
-
- from hysop.domain.box import Box
- from hysop.fields.continuous import Field
+ vd = vec.discretize(topo)
- # First, the domain : 3D and 2D boxes
- dom3D = Box()
- dom2D = Box(dimension=2, length=[1.0, 1.0], origin=[0., 0.])
+is
- Definition of a scalar fields on these domains Code::
+* check if a discretisation of vec on topo already exists. If so, returns this discretisation (i.e. vec.discreteFields[topo]),
+* else discretizes vec on topo and returns vec.discreteFields[topo].
- rho = Field(dom3D)
- rho2 = Field(dom3D, name="second_field")
- scal_2D = Field(dom2D)
- ...
-Vector fields :Code::
- velocity = Field(dom3D, is_vector=True)
- # ---> velocity is composed of 3 components.
- velo2D = Field(dom2D, is_vector=True)
- #---> velocity is composed of 2 components.
+**Remark:**
-Other fields :Code::
+Usually, in a 'complete' problem, the user does not have
+to bother with the topologies' stuff and the fields discretisations.
+All this is hidden from user : one just has to define
+some operators and the fields they work with, how those fields are initialized (see below)
+, bricks of the global problem, associated with
+some required global resolutions and some solving methods.
+Then the setup process of the problem/operators will automatically
+performs the topologies creation and the discretisation of the fields.
+See operators/problem documentation for more details or
+the examples in HySoP/Examples.
- Jacobian = Field(dom3D, nb_components = 9, name='somename')
- # ---> Creates a 9 components field defined on a 3D domain.
+
+How to initialize a field?
+--------------------------
-Now, let's turn to fields discretisation. We suppose that
-two hysop topologies have been defined (see hysop.mpi.topology)
-on domain dom3D:
-- topo1 : with a local (for each mpi process) resolution "r1"
-- topo2 : r2 Code::
+Either fields values are set to zero (default behavior) or computed with a user-defined function.
+To do so, continuous field has an 'initialize' method. It's behavior depends on the way the field
+has been defined and on the arguments of 'initialize'.
- rho.discretize(topo1)
- # --> allocate memory for numpy array of shape r1
- velocity.discretize(topo2)
- #--> allocate memory for 3 numpy arrays of shape r2
- velocity.discretize(topo1)
- #--> allocate memory for 3 numpy arrays of shape r1
-
-To handle the underlying numpy arrays do :Code::
- velo_d = velocity.discreteFields[topo1]
- #--> the discrete field for velocity on resolution r1
- print velo2_d.data[0] first component of velocity on topo1.
-
-How to initialize a field
--------------------------
+*Notice that you can find examples of all the different ways to initialize fields in fields/tests/test_field.py.*
-The simplest way :Code::
+If no function has been associated with the field
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
- #initialize on topoName
- velocity->initialize(topo=topoName)
- #initialize on all discretizations declared for the field:
- velocit->initialize()
-
-All components of the field on topo1 will be set to 0.0.
+::
+
+ # case 1: set all components values of all discrete fields to zero.
+ vec.initialize()
+
+ # case 2: set all components values of vec.discreteFields[some_topo] to zero
+ vec.initialize(time=0., topo=some_topo)
-It's also possible to set an analytical formula (depending on
-domain coordinates) to initialize a field?
+Remarks:
-There are two possibilities to define a user function:
-- either a function with numpy arrays arguments, see func1 below,
-- or a function that will be vectorized during initialize call, see func2.
+* time argument in case 2 is not used but must be provided.
+* in case 2, a call to initialize will proceed with the discretisation of vec on some_topo
+ if it does not already exist.
-Notice that :
-- func_1 required arguments : (res, x, y, ...) with res a list of numpy
-arrays which corresponds to the return values of the function and
-x, y ... (depending on the domain dimension) the grid coordinates.
-- func_2 required arguments: (x,y ...) only but will probably be less
-efficient and may result in memory reallocation.
-- optional extra args can be added in the function (see example below)
+Set a function to compute field values
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-First example : a function of numpy arrays (func_1 case):
+A user-defined function can be connected to a field::
-Code::
+ # during definition
+ vec = Field(domain=box, name='vn', is_vector=True, formula=myfunction)
+ # or later
+ vec2 = Field(domain=box, name='v2', is_vector=True)
+ vec2.set_formula(myfunction, vectorize_formula=True)
+
+myfunction must be a function of time and space, returning n numpy arrays, n being the
+number of components of the field. The exact possible signatures of this function will
+be shown later.
- def func_1_3D(res, x, y, z):
- res[0][...] = np.sin(x)
- res[1][...] = np.sin(y)
- res[2][...] = np.sin(z)
- return res
- def func_1_2D(res, x, y):
- res[0][...] = np.sin(x)
- res[1][...] = np.sin(y)
+Therefore, to initialize vec, run::
+
+ # case 1: use myfunction to compute all components values of all discrete fields
+ # at time t=0.
+ vec.initialize()
+ # case 2: use myfunction to compute all components values of all discrete fields
+ # at time t
+ t = 2.4
+ vec.initialize(time=t)
+ # case 3: use myfunction to compute all components values of vec.discreteFields[some_topo]
+ # at time t
+ vec.initialize(time=t, topo=some_topo)
+
+
+User-defined function must be::
+
+ # case 1: 'vectorize' case
+ def myfunction(x, y, z, t):
+ f_x = ...
+ f_y = ...
+ f_z = ...
+ return f_x, f_y, f_z
+
+
+ # case 2:
+
+ def myfunction(res, x, y, z, t):
+ res[0] = ...
+ res[1] = ...
+ res[2] = ...
return res
- velocity = Field(dom3D, is_vector=True, formula=func_1_3D)
- velo2D = Field(dom2D, is_vector=True, formula=func_1_2D)
- velocity.discretize(topo)
- velo2D.discretize(topo2D)
- velocity.initialize()
- velo2D.initialize()
-Obviously, the number of components return by the user function must
-corresponds to the number of components of the field.
-Second example : a function of tha will be vectorized (func_2 case).
-The code is almost the same:Code::
+Remarks:
- def func_2_3D(x, y, z):
- rx = np.sin(x)
- ry = np.sin(y)
- rz = np.sin(z)
- return rx, ry, rz
+* in 2d, z argument is removed, in 1d, only x remains (with t).
+* Res must be a list of numpy arrays of length equal to the number of
+ component of the field to be initialized.
+* The number of returned values in case 1 must be equal to the number of
+ component of the field to be initialized.
+* Case1 function takes scalar arguments as inputs and so will be vectorized
+ during initialization (about vectorization in numpy, see http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.vectorize.html), while case2 function works on numpy arrays. Case 2 must be prefered when possible for performance reasons.
+* When vectorization is needed, vectorize_formula=True must be set during definition of the field.
- velocity = Field(dom3D, is_vector=True, formula=func_2_3D)
- velocity.discretize(topo)
- velocity.initialize(vectorize_formula=True)
+Here is an example where a vector field 'velo' on a 2d domain is initialized with a vectorized function
+such that
- #--> we ask for a vectorization of func_2 so that
- # it can use arrays of coordinates as input and return
- # numpy arrays to compute velocity discrete field.
-
- Third example : add extra parameters.\n
- Let's assume that we need extra arguments to compute the field:Code::
+.. math::
+
+ velo_x(x,y,t) &=& cos(x)\\
+ velo_y(x,y,t) &=& cos(y) + tx
+
+and a vector field :math:`\omega` on a 3d domain with a 'numpy-ready' function. Notice that this last case is what we call 'Taylor-Green' in many of our examples.
+
+.. math::
+
+ \omega_x(x, y, z, t) = - cos(x)sin(y)sin(z)\\
+ \omega_y(x, y, z, t) = - sin(x)cos(y)sin(z)\\
+ \omega_z(x, y, z, t) = 2sin(x)sin(y)cos(z)
+
+The corresponding implementation will be::
+
+ # user defined functions
+ def func2d(x, y, t):
+ vx = cos(x)
+ vy = cos(y) + x*t
+ return vx, vy
+
+ def w_TG(res, x, y, z, t):
+ 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
+
+ # build a 2d domain
+ box2d = Box(length=[1., 1.])
+ # vector field
+ v2 = Field(domain=box2d, name='v2', formula=func2d, is_vector=True,
+ vectorize_formula=True)
+
+ # space discretisation and data distribution description
+ d2d = Discretization([65, 65])
+ topo2d = box2d.create_topology(d2d)
+ v2d = v2.discretize(topo2d)
+ t=3.2
+ v2.initialize(time=t, topo=topo2d)
+ # --> call func2d for each grid point and time=0
+ # to compute v2d.data[0] and v2d.data[1]
+
+ # 3d
+ box3d = Box()
+ v3 = Field(domain=box3d, name='v3', formula=w_TG)
+ # space discretisation and data distribution description
+ d3d = Discretization([65, 65, 129])
+ d3d_bis = Discretization([23, 65, 19], [2, 2, 2])
+ topo3d = box3d.create_topology(d3d)
+ topo3d_bis = box3d.create_topology(d3d_bis)
+ v3d = v3.discretize(topo3d)
+ v3d_bis = v3.discretize(topo3d_bis)
+ t=3.2
+ v3.initialize(time=t)
+ # --> call w_TG for each grid point and time=t
+ # to compute v3d.data[0] and v3d.data[1], v3d.data[2]
+ # AND v3d_bis components.
+
+
+Notice that the local mesh resolution or the data distribution are hidden from the user-defined function.
+
+
+Add extra-parameters for field'initialization
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+In addition to the standard, required, arguments of user-defined function specified above, one can
+add some parameters in the function definitions. Consider for example the definition of velo above
+and suppose that we need to introduce a dependency with the lengthes of the domain, Lx and Ly, such that:
+
+.. math::
+
+ velo_x(x,y,t) &=& L_x cos(x)\\
+ velo_y(x,y,t) &=& L_y cos(y) + tx
+
+This is possible using::
+
+ def func2d(x, y, t, lx, ly):
+ vx = lx * cos(x)
+ vy = ly * cos(y) + x*t
+ return vx, vy
+
+ Lx = 1.
+ Ly = 3.
+ extras = (Lx, Ly,) # Must be a tuple!
+ vec = Field(domain=box, is_vector=True, formula=func2d,
+ name='v2', vectorize_formula=True)
+ vec.set_formula_parameters(extras)
+ # ...
+ vec.initialize() # --> use lx=1, ly=3
+ vec.set_formula_parameters((3., 4.))
+ vec.initialize() # --> use lx=3., ly=4.
- def func_1_3D(res, x, y, z, alpha, t):
- res[0][...] = np.sin(alpha * x + t)
- res[1][...] = np.sin(alpha * y)
- res[2][...] = np.sin(z * t)
- return res
- velocity = Field(dom3D, is_vector=True, formula=func_2_3D)
-
- velocity.discretize(topo)
- #We set the values of extra args :
- alpha = ...
- t = ...
- velocity.set_formula_parameters(alpha, t)
- velocity.initialize()
-
-Remarks :
-- extra param in set_formula_parameters must be in the same order as in the func_1_3D arg list.
-- set_formula_parameters must be called for any new value of alpha/t
-- if the field has to be initialized at each step, or at least several times,
- you'd rather use hysop.operator.analytic.Analytic.
-- if you only need to compute the value of the field at a specific point, you can just call:Code::
-
- res = velocity.value(res, 1.,2.,0., alpha, t)
+Other ways to initialize a field
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-- it's possible to reset a formula with :Code::
+* Copy an existing field::
- velocity.set_formula(newFormula)
+ vec = Field(domain=dom, is_vector=True, name='v1')
+ vec2 = Field(domain=dom, is_vector=True, name='v2')
+ vec.discretize(topo)
+ vec.initialize()
+
+ vec2.copy(vec, topo)
+ # --> discretize vec2 on topo and copy vec values for topo
+
+* Random initialization::
+
+ vec = Field(domain=dom, is_vector=True, name='v1')
+ vec.randomize(topo)
+ # --> discretize vec on topo and init discrete field with
+ # random values.
+
Fields' data access
---------------------
+-------------------
-Some examples on how to get access or set values of the fields:Code::
+Some examples on how to get access or set values of the fields::
- Get a discretisation for a given continuous field:
- v3 = velocity.discreteFields[topo1]
- #Note that in this case v3 == v1
+ # Get a discretisation for a given continuous field:
+ velo = Field(domain=box, is_vector=True, name='v1')
+ # discretize and get discrete field:
+ vd = velo.discretize(some_topo)
# Get all the discretisations :
- # for topo in velocity.discreteFields.keys()
- print velocity.discreteFields[topo]
-
- # Get second component value of v2 at point of index 8,8,8
- # in the local mesh of topo2:
- print v2[1][8,8,8]
- # use vorticity formula to initialize v1 values :
- v1.initialize(formula=compute_ref)
- # norm 2 of v1
- nx, ny, nz = v1.norm()
-
-Remark : usually, in a 'complete' problem, the user does not have
-to bother with the topologies' stuff and the fields discretisations.
-All this is hidden from user : one just has to define
-some operators, bricks of the global problem, associated with
-some required global resolutions and some solving methods.
-Then the setup process of the problem/operators will automatically
-performs the topologies creation and the discretisation of the fields.
-See operators/problem documentation for more details or
-the examples in HySoP/Examples.
+ for topo in velo.discreteFields.keys():
+ current_field = velo.discreteFields[topo]
+ # Get the value of the second component of velo at point of index 8,8,8
+ # in some_topo discretization
+ print vd[1][8,8,8]
+
+
+Writing/reading fields to/from file
+-----------------------------------
+
+Fields can be written to or read from hdf files::
+
+ vec = Field(domain=box, is_vector=True, name='v1')
+ vec.discretize(topo)
+ simu = Simulation(...)
+ #...
+ # Write to file
+ vec.hdf_dump(topo, simu)
+
+ # Read from file
+ vec.hdf_load(topo)
+
+
+Other useful functions
+----------------------
+
+
+* :math:`norm(v)_d = (\sum_{i,j,...} (|v_d(i,j,...)|^2) )^{0.5}`
+
+ Sum over all grid points of the values of component d of the discrete field. Computed with::
+
+ v1 = Field(...)
+ v1.discretize(topo)
+ res = v1.norm(topo)
+ # res an array of size v1.nb_components
+
+* :math:`normh(v)_d = (h_d\sum_{i,j,...} (|v_d(i,j,...)|^2) )^{0.5}`
+
+ h_d being the space discretisation step in direction d::
+
+ v1 = Field(...)
+ v1.discretize(topo)
+ res = v1.normh(topo)
+ # res an array of size v1.nb_components
+
+
diff --git a/docs/sphinx/users_guide/index.rst b/docs/sphinx/users_guide/index.rst
index 8d5c393d31f07dd5853f9b6ae371e9e4fca1821a..bc4a54027e55a7dd707211f1a4a6af4e4a5613f3 100644
--- a/docs/sphinx/users_guide/index.rst
+++ b/docs/sphinx/users_guide/index.rst
@@ -7,15 +7,13 @@ HySoP User Guide
.. toctree::
- :maxdepth: 2
+ :maxdepth: 8
- introduction
- domains
+ space
fields
operators
numerics
- redistribute
- memo_sphinx
+ tools
biblio
diff --git a/docs/sphinx/users_guide/io_utils.rst b/docs/sphinx/users_guide/io_utils.rst
new file mode 100644
index 0000000000000000000000000000000000000000..e6f0a4f9f1bc76c6ba09e520ecf3548d49f2e5a9
--- /dev/null
+++ b/docs/sphinx/users_guide/io_utils.rst
@@ -0,0 +1,154 @@
+.. _io_utils:
+
+
+I/O management: tools and conventions in HySoP
+==============================================
+
+Different tools are available in HySoP to deal with file reading or writing, for:
+
+* scalar or vector fields, in hdf files,
+* basic text (ascii) files, for example to write the time evolution of a specific variable.
+
+
+.. _file_descriptor:
+
+Files descriptors: :class:`~hysop.tools.io_utils.IOParams`
+----------------------------------------------------------
+
+For I/O operations, especially in a parallel context, some specfici parameters must be properly defined.
+In HySoP, :class:`~hysop.tools.io_utils.IOParams` object is devoted to this task.
+
+Whenever you need to access to a file, you must use IOParams. Most of the time, you can rely on the default behavior and just specify a file name, like this::
+
+ from hysop import IOParams, IO
+ iop = IOParams('myfile.h5')
+ # use iop as input argument for i/o operators
+ # ...
+ # or to explicitely set the file type:
+ iop2 = IOParams('f2.h5', fileformat=IO.HDF5)
+ iop3 = IOParams('f3.txt', fileformat=IO.ASCII)
+
+In any case, a default path will be set, the file created if needed, an mpi communicator associated with it and so on.
+
+
+Let us now get more into details about IOParams. Such an object defines uniquely:
+
+* **IOParams.filename**, name of the output/input file (relative or absolute)
+* **IOParams.filepath**: full path to filename. Ignored if filename is absolute.
+* **IOParams.frequency**: frequency of access to this file (e.g. if equal to N, print or read every N time steps),
+ default = 1
+* **IOParams.fileformat**: the chosen format (ASCII, HDF ...), default=IO.HDF5
+* **IOParams.io_leader**, the rank (in communicator of the object using IOParams) of the mpi process that will lead the read/write process, default=0.
+
+Filename is the only required parameter, other ones have default values as shown above.
+
+If filename is 'some_relative_path/name', then filepath=$PWD/some_relative_path/name.
+
+If filename is 'just_a_name', then filepath=default_path/just_a_name
+where default_path depends on the way you run python (interactivly, ipython ...)
+and is equal to :class:`~hysop.tools.io_utils.IO.default_path()` (see table below). To check your current value, try::
+
+ from hysop.tools.io_utils import IO
+ print IO.default_path()
+
+
+=================== =======================
+Type of execution Default path values
+=================== =======================
+interactive $PWD/interactive/p1
+command line $PWD/pythonfile_name/pN
+=================== =======================
+
+N stands for the number of mpi processes used, and $PWD for current directory.
+For example, if you run in /tmp/hysop::
+
+ mpirun -np 4 myscript.py
+
+default path will be /tmp/hysop/myscript/p4
+
+
+Anyway, since the behavior of default_path strongly depends on your python version, the way you run your scropt and others things, it may be a good idea to explicitely set default path at the beginning of your
+python main program, using::
+
+ from hysop import IO
+ IO.set_default_path('somewhere')
+
+Notice that even in this case, pN will be appended to 'somewhere'.
+
+Examples::
+
+ from hysop import IOParams
+ # minimal config
+ io1 = IOParams('toto.h5')
+ print io1
+
+ io2 = IOParams('/tmp/myexamples/tutu.h5')
+ io3 = IOParams('tutu.h5', filepath='/tmp/myexamples/', frequency=4)
+
+
+As usual, check tools/tests/test_io.py for an overview of all behaviors.
+
+
+About HDF5 files in HySoP
+-------------------------
+
+
+HDF5 is, among other things, a portable file format well fitted and efficient to deal with parallel I/O for high volume of data. In HySoP, the python interface h5py is used.
+
+References :
+
+* https://www.hdfgroup.org/HDF5/
+* http://www.h5py.org
+
+
+Here is a short reminder of the basic things you need to know about hdf5 in HySoP.
+
+
+A **dataset** is a multidimensionnal array of data, in our case the values of a scalar field (or of a component of a vector field) on a grid. Each dataset may have some attributes (name ...)
+
+An **hdf file** is a container of datasets, with all the associated meta-data. It is thus important to notice that a single hdf file contains all the informations concerning the mesh on which datasets are defined, the name of the fields, \ldots and is in some way self-sufficient.
+
+
+File name convention
+^^^^^^^^^^^^^^^^^^^^
+
+In HySoP, hdf files are always named 'somename.h5'. Each file can contain several datasets but corresponds to
+one and only one instant of simulation.
+
+When fields are written at several time instants (see for example :class:`~hysop.operator.hdf_io.HDF_Writer`)
+during simulation, hdf files are named 'somename_XXXXX.h5' and a 'somename.xmf' file is created.
+
+XXXXX is the iteration number which corresponds to the time when the file has been written. The xmf file gathered all h5 files, i.e the complete collection of field values at all time steps and can be post-processed with Paraview or any similar tool.
+
+
+h5 file can be 'visualized' with h5dump or HDFView or uploaded into python (see h5py doc).
+
+Saving fields to HDF files
+==========================
+
+There are two ways to write the components of a field into an hdf file. Either with the field method :class:`~hysop.fields.continuous.Filed.hdf_dump` or
+with the operator :class:`~hysop.operator.hdf_io.HDF_Writer`.
+
+Let us start with an example::
+
+ from hysop import *
+ dom = Box()
+ d3d = Discretization([65, 65, 65])
+ f = Field(dom, name='f1', is_vector=True)
+ topo = dom.create_topology(d3d)
+ simu = Simulation(nb_iter=4)
+ simu.initialize()
+ f.hdf_dump(topo, simu)
+
+Executing this will create two files named f1.xmf and f1_00000.h5
+
+The output is done for:
+
+* a given space discretisation of the field (i.e. a :class:`~hysop.mpi.topology.Cartesian`), topo parameter in the example,
+* the current values of the field, i.e. for the current simulation time, from simu parameter in the example
+* some IOParams (see :ref:`file_descriptors`), not specified in the example, default values are used (see below)
+
+
+
+
+
diff --git a/docs/sphinx/users_guide/mpi_utils.rst b/docs/sphinx/users_guide/mpi_utils.rst
new file mode 100644
index 0000000000000000000000000000000000000000..1386f02215108417f4f83b08613cc3c911c074c7
--- /dev/null
+++ b/docs/sphinx/users_guide/mpi_utils.rst
@@ -0,0 +1,116 @@
+.. _mpi_utils:
+
+MPI interface and tools in HySoP
+================================
+
+HySoP is able to run on multi-core/multi-nodes architectures, with distributed memory systems and is based on Message Passing Interface (MPI) standard.
+For a standard usage of the software, you do not need a deep understanding of MPI standard, and even no understanding at all.
+
+Anyway, this part is a reminder of a few things to know for a more advanced usage of HySoP parallel capabilities.
+You can skip most of it if you are at ease with mpi.
+
+Glossary
+--------
+
+HySoP is able to run on hybrid systems, multi-nodes, multi-cores, multi-gpu.
+Let us then start with some clarifications about the vocabulary used in this documentation, concerning those hybrid architectures.
+
+**cluster** : a set of connected computers, all the available ressources for your simulation.
+
+**core** : a processing unit with its own memory (cache)
+
+**socket** or processor: a single computing component with two or more cores
+
+**node** : a single computer inside a cluster.
+
+An example of a cluster is shown below.
+
+.. image:: /figures/cluster_general.*
+
+What is important to notice is that there are different memory areas (cache, RAM ...), shared or not between cores, sockets ...
+
+
+MPI implementation
+------------------
+
+HySoP is based on mpi4py which is the Python package implementing the MPI standard.
+
+References:
+
+* MPI standard : http://www.mpi-forum.org/docs/docs.html
+* About mpi4py : http://pythonhosted.org/mpi4py/
+
+Once again, a glossary is required to clearly defined standard concepts used in this documentation, concerning MPI.
+
+**MPI process**: an instance of a running program. When executing a python programme as::
+
+ mpirun -np 4 python yourscript.py
+
+it means you launch 4 autonomous 'processes', each of them executing its own code from the instructions in yourscript.py. Each process executes in its own address space, its own memory, but has the ability to communicate with the others via mpi communication routines calls. Communicate means more or less to read or write some data from/to another process memory. The number of mpi processes is not supposed to correspond to the number of real physical cores available on the cluster but it's certainly not a good idea to use a number of mpi processes greater than the number of available physical cores. The mapping between processes and cores (that may be anywhere on the cluster) is done by default by MPI but may be driven by user.
+
+**Communicator**: a group (i.e. an ordered set) of MPI processes, inside a 'context' of communication.
+The context defines the properties of the communicator and how communications occur between the processes of the communicator. Mind that a message sent in one context cannot be received in another context.
+
+**Rank** : integer id of a process, relative to a given communicator. Mind that the same process may have different
+ranks in two different communicators. Ranks always start at 0.
+
+**Task** : a 'job' to be done by a set of processes. For example, consider the case where you need to solve
+an advection problem, i.e. find the values of a scalar field in a given velocity field and another problem (Navier-Stokes not to mention it) used to find the velocity field. Then, two **tasks** can be defined: one for the advection problem, attributed to one mpi process which may drive a GPU, and another task, attributed to all other processes running in parallel to solve Navier-Stokes.
+
+
+hysop and mpi
+-------------
+
+When running a simulation with hysop, a default communicator is used, including all the processes attributed at runtime
+(by mpirun command), namely main_comm. The rank of each process in main_comm is main_rank::
+
+ from hysop.mpi import main_comm, main_rank, main_size
+ print 'I am process number ', main_rank, ' among ', main_size
+
+Running a script containing the previous lines with::
+
+ mpirun -np 4 python myscript.py
+
+ended in::
+
+ I am process number 1 among 4
+ I am process number 2 among 4
+ I am process number 3 among 4
+ I am process number 0 among 4
+
+
+The mpi4py implementation is hidden into hysop.mpi package. To call any mpi4py function, try::
+
+ from hysop.mpi import MPI as MPI
+ # example : get main communicator size:
+ MPI.COMM_WORLD.Get_size()
+
+
+mpi parameters
+--------------
+
+To pass mpi-related arguments in many functions and classes in HySoP, a specific object is used, namely :class:`~hysop.tools.parameters.MPIParams`.
+
+with the following attributes:
+
+* comm : a communicator, default=main_comm
+* task_id : ... id of a task, default=DEFAULT_TASK_ID
+* rank : current process rank in comm, default=main_rank
+* on_task: boolean, true if the task_id of the current process is task_id, default=True
+
+
+Example::
+
+ from hysop import MPIParams
+ # default
+ mp = MPIParams()
+
+ # new sub-communicator:
+ comm = ... # create a subcomm
+ mp = MPIParams(comm)
+
+
+Most of the time, you will not have to set this parameter and can rely on default behavior.
+It's only useful when several tasks are required.
+
+Communicator splitting, task creation and so on are domain'jobs and the process is detailed in :ref:`domains`.
diff --git a/docs/sphinx/users_guide/tools.rst b/docs/sphinx/users_guide/tools.rst
new file mode 100644
index 0000000000000000000000000000000000000000..6e2187f802b3954602b69121af4b1726fe7e23e6
--- /dev/null
+++ b/docs/sphinx/users_guide/tools.rst
@@ -0,0 +1,8 @@
+.. _tools::
+
+
+.. toctree::
+ :maxdepth: 5
+
+ io_utils
+ mpi_utils
diff --git a/docs/sphinx/users_guide/work_arrays.rst b/docs/sphinx/users_guide/work_arrays.rst
new file mode 100644
index 0000000000000000000000000000000000000000..48da58af85310d1eba330a33a876a1defcb9cdb2
--- /dev/null
+++ b/docs/sphinx/users_guide/work_arrays.rst
@@ -0,0 +1,51 @@
+.. _work_spaces:
+
+
+Internal work arrays management
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. currentmodule:: hysop.operator.computational
+
+Internal temporary buffers may be required by computational operators (buffers for Runge-Kutta scheme for example). By default, each operator handles its own work arrays and manage their memory.
+
+But, to improve performance and limit the required memory amount, user can provide some external
+storage for those arrays, for example to share them between several operators.
+
+Each operator has a :py:meth:`~Computational.get_work_properties` function, returning information on the required memory for internal stuff.
+This function returns a dictionnary with the following keys:
+
+ * 'rwork' : parameters for internal real arrays
+ * 'iwork' : parameters for internal integer arrays
+
+To each key corresponds a value which may be either 'None' if no internal memory is required
+or a list of tuples, each tuple being the shape of the required array.
+
+For example:
+
+ >>> op = ...
+ >>> op.discretize()
+ >>> work_params = op.get_work_properties()
+ >>> print work_params
+ {'rwork': [(12, 12), (45, 12, 33)], 'iwork': None}
+
+
+means that op requires two real arrays, of shape (12,12) and (45, 12, 33) while no
+integer arrays is needed.
+
+Note that this function must be called after discretization, else the operator may lack some
+information about local resolution and data distribution.
+
+Then, if needed, one may set internal arrays for the operator with setup optional parameters, as in the example below:
+
+ >>> myrwork = []
+ >> for shape in work.params['rwork']
+ ... myrwork.append(npw.zeros(shape))
+
+ >>> op.setup(rwork=myrwork)
+ >>> # And this work may be used by another operator
+ >>> op2.setup(rwork=myrwork
+
+
+Warning: obviously, data in work arrays may ne modified at any time by the operator.
+
+