From 24f60ca127e80b9ff17fd7cb5cf3320ee993ad09 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Thu, 13 Mar 2014 13:50:25 +0000
Subject: [PATCH] add time-dependent var. and update (some) mpi tests

---
 HySoP/hysop/__init__.py.in                    |  4 +-
 HySoP/hysop/default_methods.py                |  3 +-
 HySoP/hysop/domain/domain.py                  |  5 +-
 HySoP/hysop/f2py/scales2py.f90                |  7 +-
 HySoP/hysop/fields/tests/test_field.py        |  8 ++
 HySoP/hysop/fields/tests/test_variable.py     | 56 +++++++++++++
 .../hysop/{ => fields}/variable_parameter.py  | 59 +++++++++++--
 HySoP/hysop/gpu/tests/test_copy.py            | 10 ++-
 HySoP/hysop/mpi/tests/test_topology.py        |  3 +-
 HySoP/hysop/operator/advection.py             | 18 ++--
 HySoP/hysop/operator/advection_dir.py         | 84 ++++++-------------
 HySoP/hysop/operator/discrete/stretching.py   |  3 +-
 HySoP/hysop/operator/stretching.py            |  4 +-
 HySoP/hysop/operator/tests/test_Stretching.py |  4 +
 .../hysop/operator/tests/test_advec_scales.py | 41 ++++-----
 HySoP/hysop/operator/tests/test_analytic.py   | 15 ++++
 .../operator/tests/test_diff_poisson_3D.py    |  3 +
 .../hysop/operator/tests/test_penalization.py |  8 ++
 HySoP/hysop/operator/tests/test_poisson.py    | 12 ++-
 .../tests/test_velocity_correction.py         |  5 ++
 20 files changed, 237 insertions(+), 115 deletions(-)
 create mode 100644 HySoP/hysop/fields/tests/test_variable.py
 rename HySoP/hysop/{ => fields}/variable_parameter.py (52%)

diff --git a/HySoP/hysop/__init__.py.in b/HySoP/hysop/__init__.py.in
index 7d7e22390..23eec22a2 100755
--- a/HySoP/hysop/__init__.py.in
+++ b/HySoP/hysop/__init__.py.in
@@ -38,8 +38,8 @@ import parmepy.fields.continuous
 Field = parmepy.fields.continuous.Field
 
 ## Variable parameters
-import parmepy.variable_parameter
-VariableParameter = parmepy.variable_parameter.VariableParameter
+import parmepy.fields.variable_parameter
+VariableParameter = parmepy.fields.variable_parameter.VariableParameter
 
 ## Simulation parameters
 import parmepy.problem.simulation
diff --git a/HySoP/hysop/default_methods.py b/HySoP/hysop/default_methods.py
index 8acaa3aed..aa5bbc3e2 100644
--- a/HySoP/hysop/default_methods.py
+++ b/HySoP/hysop/default_methods.py
@@ -9,10 +9,9 @@ from parmepy.numerics.integrators.runge_kutta2 import RK2
 from parmepy.numerics.integrators.runge_kutta3 import RK3
 from parmepy.numerics.interpolation import Linear
 from parmepy.numerics.remeshing import L2_1
+#from parmepy.operator.discrete.stretching import Conservative
 
 
-#class DefaultMethods(object):
-
 ADVECTION = {TimeIntegrator: RK2, Interpolation: Linear,
              Remesh: L2_1, Support: '', Splitting: 'o2', MultiScale: L2_1}
 
diff --git a/HySoP/hysop/domain/domain.py b/HySoP/hysop/domain/domain.py
index 3182ba657..9f950b038 100644
--- a/HySoP/hysop/domain/domain.py
+++ b/HySoP/hysop/domain/domain.py
@@ -65,8 +65,9 @@ class Domain(object):
                 topoCreation = Cartesian.withResolutionFixed
             else:
                 topoCreation = Cartesian.withResolution
-            newTopo = topoCreation(self, topoResolution,
-                                   gridResolution, ghosts=ghosts, comm=comm)
+            newTopo = topoCreation(self, shape=topoResolution,
+                                   globalMeshResolution=gridResolution,
+                                   ghosts=ghosts, comm=comm)
         otherid = self.register(newTopo)
         if otherid < 0:
             otherid = newTopo.getId()
diff --git a/HySoP/hysop/f2py/scales2py.f90 b/HySoP/hysop/f2py/scales2py.f90
index 264aeeff0..ab5b440ca 100755
--- a/HySoP/hysop/f2py/scales2py.f90
+++ b/HySoP/hysop/f2py/scales2py.f90
@@ -20,7 +20,7 @@ contains
   !! @param[in] mpi communicator from python
   !! @param[out] datashape local dimension of the input/output field
   !! @param[out] offset absolute index of the first component of the local field
-  subroutine init_advection_solver(ncells,lengths,topodims,main_comm,datashape,offset,dim,order,stab_coeff,dim_split)
+  subroutine init_advection_solver(ncells,lengths,topodims,main_comm,datashape,offset,dim,order,dim_split)
     integer, intent(in) :: dim
     integer, dimension(dim),intent(in) :: ncells
     real(pk),dimension(dim), intent(in) :: lengths
@@ -29,7 +29,8 @@ contains
     integer(ik), dimension(dim), intent(out) :: datashape
     integer(ik), dimension(dim), intent(out) :: offset
     character(len=*), optional, intent(in)  ::  order, dim_split
-    real(pk), optional, intent(out) ::  stab_coeff
+    !! real(pk), optional, intent(out) ::  stab_coeff
+    real(pk) ::  stab_coeff
     !f2py integer optional , depend(ncells) :: dim=len(ncells)
     !f2py intent(hide) dim
     !f2py character(*) optional, intent(in) :: order = 'p_O2'
@@ -54,7 +55,7 @@ contains
     call advec_init(order,stab_coeff,dim_split=dim_split)
 
     ! get the local resolution (saved in scales global variable "N_proc")
-    datashape = N_proc + 1
+    datashape = N_proc
     ! get offset (i.e. global index of the lowest point of the current subdomain == scales global var "coord")
     offset = coord
 
diff --git a/HySoP/hysop/fields/tests/test_field.py b/HySoP/hysop/fields/tests/test_field.py
index 020e8df0f..9af2b3f25 100644
--- a/HySoP/hysop/fields/tests/test_field.py
+++ b/HySoP/hysop/fields/tests/test_field.py
@@ -83,3 +83,11 @@ def test_integrate_onSurf():
     res = myf.integrateOnSurface(surf, topo)
     sref = dom.length[1] * dom.length[2]
     assert(abs(res - sref) < 1e-6)
+
+# This may be useful to run mpi tests
+if __name__ == "__main__":
+    test_continuous()
+    test_analytical()
+    test_analytical_reset()
+    test_integrate_onSurf()
+    test_discretization()
diff --git a/HySoP/hysop/fields/tests/test_variable.py b/HySoP/hysop/fields/tests/test_variable.py
new file mode 100644
index 000000000..d3fac2c7b
--- /dev/null
+++ b/HySoP/hysop/fields/tests/test_variable.py
@@ -0,0 +1,56 @@
+"""
+Testing parmepy.field.variable_parameter.Variable_parameter
+"""
+from parmepy.fields.variable_parameter import VariableParameter
+from parmepy.problem.simulation import Simulation
+from math import sin, cos
+import numpy as np
+
+
+def test_constantVar():
+    uinf = 1.
+    var = VariableParameter(data=uinf, name='uinf')
+    assert var['uinf'] == 1.
+    var['uinf'] *= 3
+    assert var['uinf'] == 3.
+
+    v2 = VariableParameter({'uinf': uinf})
+    assert v2['uinf'] == 1.
+
+
+def func(simu):
+    time = simu.tk
+    return np.asarray((sin(time), cos(time)))
+
+
+def test_timeVar():
+    var = VariableParameter(formula=func)
+    simu = Simulation(tinit=0., tend=1., timeStep=0.1)
+    var.update(simu)
+    assert np.allclose(var['func'], [0., 1.])
+    assert var.name == 'func'
+    simu.advance()
+    var.update(simu)
+    assert np.allclose(var['func'], [sin(0.1), cos(0.1)])
+    var = VariableParameter(formula=func, name='nn')
+    assert var.name == 'nn'
+    assert var.formula is not None
+
+
+def test_timeVar2():
+    var = VariableParameter(formula=func, name='toto', data={'alpha': 1.})
+    simu = Simulation(tinit=0., tend=1., timeStep=0.1)
+    assert var.name == 'alpha'
+    assert np.allclose(var['alpha'], [1.])
+    var.update(simu)
+    assert np.allclose(var['alpha'], [sin(0.), cos(0.)])
+    simu.advance()
+    var.update(simu)
+    assert np.allclose(var['alpha'], [sin(0.1), cos(0.1)])
+
+
+# This may be useful to run mpi tests
+if __name__ == "__main__":
+    test_constantVar()
+    test_timeVar()
+    test_timeVar2()
diff --git a/HySoP/hysop/variable_parameter.py b/HySoP/hysop/fields/variable_parameter.py
similarity index 52%
rename from HySoP/hysop/variable_parameter.py
rename to HySoP/hysop/fields/variable_parameter.py
index 3cd0bea10..a43388cd6 100644
--- a/HySoP/hysop/variable_parameter.py
+++ b/HySoP/hysop/fields/variable_parameter.py
@@ -7,7 +7,7 @@
 #
 #
 """
-@file variables/variables.py
+@file variable_parameter.py
 
 Tools to define parameters that do not depend on space but may
 be used as input-output parameter during simulation.
@@ -26,6 +26,32 @@ op.apply(simu) ===> change simu.time_step
 
 \endcode
 
+A python function may be set to compute the parameter.
+Example : if you want to define a time-dependent parameter
+alpha, such that alpha = [sin(t), cos(t)], do :
+
+\code
+def myfunc(simu):
+    return sin(simu.time), cos(simu.time)
+
+
+alpha = VariableParameter(formula=myfunc, name='alpha')
+
+simu = Simulation(tstart = 0., tend=1., timeStep=0.1)
+
+
+print alpha
+# --> return {'alpha': None}
+
+alpha.update(simu)
+print alpha
+# --> return {'alpha': (0., 1.)}
+
+simu.advance()
+alpha.update(simu)
+print alpha
+# --> return {'alpha': (0.19866933079506122, 0.9800665778412416)}
+
 """
 
 
@@ -36,7 +62,7 @@ class VariableParameter(object):
     VariableParameter has a member data which is mutable.
     """
 
-    def __init__(self, data=None, name=None):
+    def __init__(self, data=None, name=None, formula=None):
         """
         Creates a dictionnary with data
         @param data: the data used as parameter
@@ -47,14 +73,35 @@ class VariableParameter(object):
         """
         ## name of the variable
         self.name = name
+        if self.name is None and formula is not None:
+            self.name = formula.__name__
+        ## Formula used to compute data (a python function)
+        self.formula = formula
+        if self.formula is None:
+            self.formula = self._constant
         ## data values
         if isinstance(data, dict):
+            msg = 'A dictionnary is used to initialize the variable.'
+            msg += 'Name parameter will be ignored.'
+            print msg
             self.data = data
+            self.name = data.keys()[0]
         else:
-            if name is None:
+            if self.name is None:
                 msg = "Name arg is required when data is not a dict."
                 raise AttributeError(msg)
-            self.data = {name: data}
+            self.data = {self.name: data}
+
+    def _constant(self, simu=None):
+        pass
+
+    def update(self, simu=None):
+        """
+        Apply formula to compute data for
+        a given simulation (current time ...)
+        @param simu : a parmepy.problem.simulation.Simulation
+        """
+        self.data[self.name] = self.formula(simu)
 
     def __getitem__(self, key):
         """ Access to the content of the data[key]
@@ -70,6 +117,4 @@ class VariableParameter(object):
         self.data[key] = value
 
     def __str__(self):
-        s = "Variable parameter, which contains: "
-        s += str(self.data)
-        return s
+        return str(self.data)
diff --git a/HySoP/hysop/gpu/tests/test_copy.py b/HySoP/hysop/gpu/tests/test_copy.py
index 0dd2ea6f1..087cf9155 100644
--- a/HySoP/hysop/gpu/tests/test_copy.py
+++ b/HySoP/hysop/gpu/tests/test_copy.py
@@ -7,10 +7,12 @@ from parmepy.constants import ORDER, np, PARMES_REAL
 from parmepy.gpu.tools import get_opencl_environment
 from parmepy.gpu.gpu_kernel import KernelLauncher
 
+DEVICE_NUMBER = 1
+
 
 def test_copy2D():
     resolution = (256, 256)
-    cl_env = get_opencl_environment(0, 0, 'gpu', PARMES_REAL)
+    cl_env = get_opencl_environment(0, DEVICE_NUMBER, 'gpu', PARMES_REAL)
     vec = 2
     src_copy = 'kernels/copy.cl'
     build_options = ""
@@ -48,7 +50,7 @@ def test_copy2D():
 def test_copy2D_rect():
     resolution = (256, 512)
     resolutionT = (512, 256)
-    cl_env = get_opencl_environment(0, 0, 'gpu', PARMES_REAL)
+    cl_env = get_opencl_environment(0, DEVICE_NUMBER, 'gpu', PARMES_REAL)
     vec = 2
     src_copy = 'kernels/copy.cl'
     build_options = ""
@@ -109,7 +111,7 @@ def test_copy2D_rect():
 
 def test_copy3D():
     resolution = (64, 64, 64)
-    cl_env = get_opencl_environment(0, 0, 'gpu', PARMES_REAL)
+    cl_env = get_opencl_environment(0, DEVICE_NUMBER, 'gpu', PARMES_REAL)
     vec = 4
     src_copy = 'kernels/copy.cl'
     build_options = ""
@@ -151,7 +153,7 @@ def test_copy3D_rect():
     resolution_x = (16, 32, 64)
     resolution_y = (32, 16, 64)
     resolution_z = (64, 16, 32)
-    cl_env = get_opencl_environment(0, 0, 'gpu', PARMES_REAL)
+    cl_env = get_opencl_environment(0, DEVICE_NUMBER, 'gpu', PARMES_REAL)
     vec = 4
     src_copy = 'kernels/copy.cl'
 
diff --git a/HySoP/hysop/mpi/tests/test_topology.py b/HySoP/hysop/mpi/tests/test_topology.py
index 4879ae9ee..00ef3e6b6 100644
--- a/HySoP/hysop/mpi/tests/test_topology.py
+++ b/HySoP/hysop/mpi/tests/test_topology.py
@@ -58,14 +58,13 @@ def test_operator_equal():
     assert not topo == topo2
     topo2 = Cartesian(dom, 2, [11, 33, 17])
     assert not topo == topo2
-    
+
     topo2 = Cartesian(dom, 2, shape=topoDims[-1::-1],
                       globalMeshResolution=resolTopo)
     if topo.rank == 0:
         print topoDims.reverse()
         print topo
         print topo2
-        
     assert topo == topo2
 
 
diff --git a/HySoP/hysop/operator/advection.py b/HySoP/hysop/operator/advection.py
index 21e4cdae2..5037a1fa6 100644
--- a/HySoP/hysop/operator/advection.py
+++ b/HySoP/hysop/operator/advection.py
@@ -83,10 +83,7 @@ class Advection(Operator):
         ## Transported fields
         self.output = self.advectedFields
         self.input = [var for var in self.variables]
-        ## Resolutions of the fields. Dictionnary with
-        ## fields as keys.
-        ## Example :
-        ## resolutions[velocity] may return [32,32,32].
+        ## Global resolutions of the fields.
         self.resolutions = resolutions
         self._isMultiScale = False
         v_resol = self.resolutions[self.velocity]
@@ -203,7 +200,7 @@ class Advection(Operator):
         # --- Advection solver from SCALES ---
         if self._isSCALES:
             if not self._dim == 3:
-                    raise ValueError("Scales Advection not implemented in 2D.")
+                raise ValueError("Scales Advection not implemented in 2D.")
             if self.ghosts is not None:
                 if (self.ghosts > 0).any():
                     raise ValueError("Scales Advection not\
@@ -236,12 +233,12 @@ class Advection(Operator):
                 comm = self._predefinedTopo.comm
             elif self._comm is not None:
                 main_size = self._comm.Get_size()
-                comm =self._comm
+                comm = self._comm
             else:
                 from parmepy.mpi.main_var import main_size
                 from parmepy.mpi.main_var import main_comm as comm
             topodims = [1, 1, main_size]
-            scalesres, scalesoffset, stab_coeff = \
+            scalesres, scalesoffset = \
                 scales.init_advection_solver(nbcells,
                                              self.domain.length,
                                              topodims, comm.py2f(),
@@ -255,15 +252,16 @@ class Advection(Operator):
                     not scales compliant.'
                 for v in self.variables:
                     self.discreteFields[v] = v.discretize(topo)
-            else:  # Note Franck : is it reasonable to allow different
-                # topologies in the same
-                # advection operator???
+            else:
                 for v in self.variables:
                     # the topology for v ...
                     topo = \
                         self.domain.getOrCreateTopology(self._dim,
                                                         self.resolutions[v],
                                                         topodims,
+                                                        precomputed=True,
+                                                        offset=scalesoffset,
+                                                        localres=scalesres,
                                                         ghosts=self.ghosts,
                                                         comm=self._comm)
                     # ... and the corresponding discrete field
diff --git a/HySoP/hysop/operator/advection_dir.py b/HySoP/hysop/operator/advection_dir.py
index 97b76f9c7..919fe9b3d 100644
--- a/HySoP/hysop/operator/advection_dir.py
+++ b/HySoP/hysop/operator/advection_dir.py
@@ -42,6 +42,10 @@ class AdvectionDir(Operator):
         Operator.__init__(self, v, method, topo=topo,
                           ghosts=ghosts, task_id=task_id, comm=comm,
                           name_suffix=name_suffix + S_DIR[d])
+        if self._predefinedTopo is not None:
+            raise ValueError("User defined topology is not\
+               allowed for advecDir.")
+
         self.output = advectedFields
         self.input = [var for var in self.variables]
         ## Transport velocity
@@ -54,8 +58,7 @@ class AdvectionDir(Operator):
         ## resolutions[velocity] may return [32,32,32].
         self.resolutions = resolutions
         self._isMultiScale = isMultiScale
-        self._v_ghosts = np.array([0, ] * self.domain.dimension,
-                                  dtype=PARMES_INDEX)
+        self._v_ghosts = None
         if self._isMultiScale:
             if self.method[Support].find('gpu') < 0:
                 raise ValueError("Multiscale advection is not supported in "
@@ -86,68 +89,31 @@ class AdvectionDir(Operator):
         Discretisation according to the chosen method.
         Available methods : See Advection.setUp
         """
-        if self._predefinedTopo is not None:
-            main_size = self._predefinedTopo.size
-        elif self._comm is not None:
+        if self._comm is not None:
             main_size = self._comm.Get_size()
         else:
             from parmepy.mpi import main_size
-        if main_size == 1:
-            # - topologies creation and variables discretization -
-            if self._predefinedTopo is not None:
-                topo = self._predefinedTopo
-                for v in self.variables:
-                    self.discreteFields[v] = v.discretize(topo)
-            else:
-                for v in self.variables:
-                    if v == self.velocity:
-                        topo = self.domain.getOrCreateTopology(
-                            self.domain.dimension, self.resolutions[v],
-                            ghosts=self._v_ghosts, comm=self._comm)
-                    else:
-                        topo = self.domain.getOrCreateTopology(
-                            self.domain.dimension, self.resolutions[v],
-                            comm=self._comm)
-                    self.discreteFields[v] = v.discretize(topo)
+
+        topodims = np.ones((self.domain.dimension))
+        # MPI topology depends on direction
+        if self.dir == self.domain.dimension - 1:
+            # Cut in first dir for last dir computations
+            topodims[0] = main_size
         else:
-            topodims = np.ones((self.domain.dimension))
-            if self._predefinedTopo is not None:
-                if self._predefinedTopo.dim != 1:
-                    raise ValueError("User-defined topology is not\
-                (yet) allowed for advection.")
-                else:
-                    comm = self._predefinedTopo._comm_origin
-                    shape = self._predefinedTopo.shape
-                    size = self._predefinedTopo.size
-                    if self.dir == self.domain.dimension - 1 and \
-                        shape[-1] > 1:
-                        topodims[0] = size
-                    elif self.dir != self.domain.dimension - 1 and \
-                        shape[self.dir] > 1:
-                        topodims[-1] = size
-                    else:
-                        topodims[...] = shape
+            # Cut in last dir
+            topodims[-1] = main_size
+        for v in self.variables:
+            if v == self.velocity:
+                topo = self.domain.getOrCreateTopology(
+                    1, self.resolutions[v], topoResolution=topodims,
+                    fixedResolution=True, ghosts=self._v_ghosts,
+                    comm=self._comm)
             else:
-                # MPI topology depends on direction
-                if self.dir == self.domain.dimension - 1:
-                    # Cut in first dir for last dir computations
-                    topodims[0] = main_size
-                else:
-                    # Cut in last dir
-                    topodims[-1] = main_size
-
-            for v in self.variables:
-                if v == self.velocity:
-                    topo = self.domain.getOrCreateTopology(
-                        self.domain.dimension, self.resolutions[v],
-                        topoResolution=topodims, fixedResolution=True,
-                        ghosts=self._v_ghosts, comm=self._comm)
-                else:
-                    topo = self.domain.getOrCreateTopology(
-                        self.domain.dimension, self.resolutions[v],
-                        topoResolution=topodims, fixedResolution=True,
-                        comm=self._comm)
-                self.discreteFields[v] = v.discretize(topo)
+                topo = self.domain.getOrCreateTopology(
+                    1, self.resolutions[v], topoResolution=topodims,
+                    fixedResolution=True, ghosts=self.ghosts,
+                    comm=self._comm)
+            self.discreteFields[v] = v.discretize(topo)
 
     @staticmethod
     def getWorkLengths(method=None, domain_dim=None):
diff --git a/HySoP/hysop/operator/discrete/stretching.py b/HySoP/hysop/operator/discrete/stretching.py
index 85c213a3b..146c5418b 100755
--- a/HySoP/hysop/operator/discrete/stretching.py
+++ b/HySoP/hysop/operator/discrete/stretching.py
@@ -17,7 +17,6 @@ import parmepy.tools.numpywrappers as npw
 from parmepy.numerics.updateGhosts import UpdateGhosts
 from parmepy.mpi import MPI
 from abc import ABCMeta, abstractmethod
-import parmepy.default_methods as default
 import math
 ceil = math.ceil
 
@@ -45,6 +44,7 @@ class Stretching(DiscreteOperator):
         ## vorticity discrete field
         self.vorticity = vorticity
         if method is None:
+            import parmepy.default_methods as default
             method = default.STRETCHING
         DiscreteOperator.__init__(self, [self.velocity, self.vorticity],
                                   method=method)
@@ -196,6 +196,7 @@ class GradUW(Stretching):
         @param method : the dict of parameters for the operator.
         """
         if method is None:
+            import parmepy.default_methods as default
             method = default.STRETCHING
         assert TimeIntegrator in method,\
             'A time integrator is required for the stretching.'
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index ce310da71..9612dc460 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -11,7 +11,6 @@ from parmepy.methods_keys import TimeIntegrator, Formulation, \
 from parmepy.numerics.finite_differences import FD_C_4
 from parmepy.operator.continuous import Operator
 from parmepy.operator.discrete.stretching import Conservative, GradUW
-import parmepy.default_methods as default
 import numpy as np
 
 
@@ -48,7 +47,9 @@ class Stretching(Operator):
         ## Grid resolution for each variable (dictionnary)
         self.resolutions = resolutions
         ## Numerical methods for time and space discretization
+
         if method is None:
+            import parmepy.default_methods as default
             method = default.STRETCHING
         self.method = method
         assert Formulation in self.method.keys()
@@ -104,6 +105,7 @@ class Stretching(Operator):
                       SpaceDiscretisation: FD_C_4}
         """
         if method is None:
+            import parmepy.default_methods as default
             method = default.STRETCHING
         assert Formulation in method,\
             'A formulation is required for the stretching.'
diff --git a/HySoP/hysop/operator/tests/test_Stretching.py b/HySoP/hysop/operator/tests/test_Stretching.py
index 2d5f40e48..1974dc733 100755
--- a/HySoP/hysop/operator/tests/test_Stretching.py
+++ b/HySoP/hysop/operator/tests/test_Stretching.py
@@ -86,3 +86,7 @@ def test_Stretching():
     stretch.setUp()
     simulation = Simulation(tinit=0, tend=20, timeStep=timeStep)
     stretch.apply(simulation)
+
+
+if __name__ == "__main__":
+    test_Stretching()
diff --git a/HySoP/hysop/operator/tests/test_advec_scales.py b/HySoP/hysop/operator/tests/test_advec_scales.py
index c2df44665..a3fa3b831 100755
--- a/HySoP/hysop/operator/tests/test_advec_scales.py
+++ b/HySoP/hysop/operator/tests/test_advec_scales.py
@@ -39,7 +39,6 @@ def test_nullVelocity_m4():
     advec_py.discretize()
     advec.setUp()
     advec_py.setUp()
-
     scal_d = scal.discreteFields.values()[0]
     scal_ref_d = scal_ref.discreteFields.values()[0]
 
@@ -47,10 +46,10 @@ def test_nullVelocity_m4():
         np.random.random(scal_d.data[0].shape),
         dtype=PARMES_REAL, order=ORDER)
     scal_ref_d.data[0][...] = scal_d.data[0][...]
-
+    assert (velo.norm() == 0).all()
     advec.apply(Simulation(tinit=0., tend=0.1, nbIter=1))
     advec_py.apply(Simulation(tinit=0., tend=0.1, nbIter=1))
-
+    print np.max(np.abs(scal_ref_d.data[0] - scal_d.data[0]))
     assert np.allclose(scal_ref_d.data[0], scal_d.data[0])
 
 
@@ -233,29 +232,19 @@ def test_nullVelocity_vec_m8():
     advec_py.discretize()
     advec.setUp()
     advec_py.setUp()
-
+    assert (velo.norm() == 0).all()
     scal_d = scal.discreteFields.values()[0]
     scal_ref_d = scal_ref.discreteFields.values()[0]
-
-    scal_d.data[0][...] = np.asarray(
-        np.random.random(scal_d.data[0].shape),
-        dtype=PARMES_REAL, order=ORDER)
-    scal_d.data[1][...] = np.asarray(
-        np.random.random(scal_d.data[1].shape),
-        dtype=PARMES_REAL, order=ORDER)
-    scal_d.data[2][...] = np.asarray(
-        np.random.random(scal_d.data[2].shape),
-        dtype=PARMES_REAL, order=ORDER)
-    scal_ref_d.data[0][...] = scal_d.data[0][...]
-    scal_ref_d.data[1][...] = scal_d.data[1][...]
-    scal_ref_d.data[2][...] = scal_d.data[2][...]
+    for i in xrange(box.dimension):
+        scal_d.data[i][...] = \
+            npw.realarray(np.random.random(scal_d.data[i].shape))
+        scal_ref_d.data[i][...] = scal_d.data[i][...]
 
     advec.apply(Simulation(tinit=0., tend=0.075, nbIter=1))
     advec_py.apply(Simulation(tinit=0., tend=0.075, nbIter=1))
 
-    assert np.allclose(scal_ref_d.data[0], scal_d.data[0])
-    assert np.allclose(scal_ref_d.data[1], scal_d.data[1])
-    assert np.allclose(scal_ref_d.data[2], scal_d.data[2])
+    for i in xrange(box.dimension):
+        assert np.allclose(scal_ref_d.data[i], scal_d.data[i])
 
 
 def _randomVelocity_m4():
@@ -580,3 +569,15 @@ def test_randomVelocity_vec_m8():
     assert np.allclose(scal_ref_d.data[0], scal_d.data[0], atol=1e-07)
     assert np.allclose(scal_ref_d.data[1], scal_d.data[1], atol=1e-07)
     assert np.allclose(scal_ref_d.data[2], scal_d.data[2], atol=1e-07)
+
+if __name__ == "__main__":
+    test_nullVelocity_m4()
+    test_nullVelocity_vec_m4()
+    test_nullVelocity_m6()
+    test_nullVelocity_vec_m6()
+    test_nullVelocity_m8()
+    test_nullVelocity_vec_m8()
+    test_randomVelocity_m6()
+    test_randomVelocity_vec_m6()
+    test_randomVelocity_m8()
+    test_randomVelocity_vec_m8()
diff --git a/HySoP/hysop/operator/tests/test_analytic.py b/HySoP/hysop/operator/tests/test_analytic.py
index 5c5f159a5..d7d3cae67 100644
--- a/HySoP/hysop/operator/tests/test_analytic.py
+++ b/HySoP/hysop/operator/tests/test_analytic.py
@@ -457,3 +457,18 @@ def test_analytical_op_6():
     for i in xrange(caf.nbComponents):
         assert np.allclose(cafd[i], refd.data[i])
         assert id(cafd.data[i]) == ids[i]
+
+if __name__ == "__main__":
+    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_analytical_op_1()
+    test_analytical_op_3()
+    test_analytical_op_4()
+    test_analytical_op_5()
+    test_analytical_op_6()
diff --git a/HySoP/hysop/operator/tests/test_diff_poisson_3D.py b/HySoP/hysop/operator/tests/test_diff_poisson_3D.py
index 978c6477f..3995bf406 100755
--- a/HySoP/hysop/operator/tests/test_diff_poisson_3D.py
+++ b/HySoP/hysop/operator/tests/test_diff_poisson_3D.py
@@ -72,3 +72,6 @@ def test_Diff_Poisson():
                       iterMax=1000000)
     diffusion.apply(simu)
     poisson.apply(simu)
+
+if __name__ == "__main__":
+    test_Diff_Poisson()
diff --git a/HySoP/hysop/operator/tests/test_penalization.py b/HySoP/hysop/operator/tests/test_penalization.py
index 2b3d68c5c..f43a2f00b 100644
--- a/HySoP/hysop/operator/tests/test_penalization.py
+++ b/HySoP/hysop/operator/tests/test_penalization.py
@@ -177,3 +177,11 @@ def testPenalVec3D():
     veloRef.discretize(topo)
     veloRef.load('ref_velo3D_PenalHspherePlane', fieldname='VeloRef')
     assert (veloRef.norm() == velo.norm()).all()
+
+if __name__ == "__main__":
+    testPenalScal2D()
+    ## testPenalScal3D()
+    ## testPenalScal2D_2()
+    ## testPenalScal3D_2()
+    ## testPenalVec2D()
+    ## testPenalVec3D()
diff --git a/HySoP/hysop/operator/tests/test_poisson.py b/HySoP/hysop/operator/tests/test_poisson.py
index 6880b2a48..2cf8231b6 100755
--- a/HySoP/hysop/operator/tests/test_poisson.py
+++ b/HySoP/hysop/operator/tests/test_poisson.py
@@ -76,7 +76,6 @@ def test_Poisson3D():
     ref = pp.Field(domain=dom, name='reference', isVector=True)
     refOp = Analytic(ref, formula=computeRef,
                      resolutions={ref: resol})
-
     simu = Simulation(nbIter=10)
     refOp.discretize()
     refOp.setUp()
@@ -84,12 +83,17 @@ def test_Poisson3D():
     poisson.setUp()
     refOp.apply(simu)
     topo = poisson.discreteFields[vorticity].topology
+    ## from parmepy.operator.redistribute import Redistribute
+    ## distr = Redistribute([ref], refOp, poisson)
+    ## distr.discretize()
+    ## distr.setUp()
+    ## distr.apply(simu)
     vorticity.setTopoInit(topo)
     vorticity.initialize()
     poisson.apply(simu)
     assert np.allclose(ref.norm(), velocity.norm())
     poisson.finalize()
-    refD = ref.discreteFields.values()[0]
+    refD = ref.discretization(topo)
     vd = velocity.discreteFields.values()[0]
 
     assert np.allclose(ref.norm(), velocity.norm())
@@ -129,3 +133,7 @@ def test_Poisson3D():
 ## #    assert np.allclose(ref.norm(), velocity.norm())
 ## #    for i in range(dom.dimension):
 ## #        assert np.allclose(vd[i], refD[i])
+
+# This may be useful to run mpi tests
+if __name__ == "__main__":
+    test_Poisson3D()
diff --git a/HySoP/hysop/operator/tests/test_velocity_correction.py b/HySoP/hysop/operator/tests/test_velocity_correction.py
index 06c2ff23d..92ac828b3 100755
--- a/HySoP/hysop/operator/tests/test_velocity_correction.py
+++ b/HySoP/hysop/operator/tests/test_velocity_correction.py
@@ -146,3 +146,8 @@ def test_velocity_correction_2D():
     for i in xrange(1, dim):
         flowrate2 = velo.integrateOnSurface(sref, topo, component=i)
         assert flowrate2 < tol
+
+# This may be useful to run mpi tests
+if __name__ == "__main__":
+    test_velocity_correction_2D()
+    test_velocity_correction_3D()
-- 
GitLab