From 5c7f5d031f9073b24d3587670d2bd5d7d43340fd Mon Sep 17 00:00:00 2001
From: Eric Dalissier <Eric.Dalissier@imag.fr>
Date: Tue, 4 Jun 2013 12:28:08 +0000
Subject: [PATCH] synchronization of field of array in stretching and
 differential operator

test stretching and test synchronize ok
---
 HySoP/hysop/numerics/differentialOperator.py  | 29 +++---
 HySoP/hysop/numerics/fct2op.py                | 99 +++++++++++--------
 .../numerics/integrators/runge_kutta3.py      | 11 ++-
 HySoP/hysop/operator/discrete/stretching.py   | 23 +++--
 .../operator/discrete/synchronizeGhosts.py    | 20 ++--
 HySoP/hysop/operator/synchronizeGhosts.py     | 39 +++++---
 HySoP/hysop/operator/tests/test_Stretching.py | 39 +++++---
 .../operator/tests/test_synchronizeGhosts.py  |  6 +-
 HySoP/hysop/problem/problem.py                |  2 +-
 9 files changed, 160 insertions(+), 108 deletions(-)

diff --git a/HySoP/hysop/numerics/differentialOperator.py b/HySoP/hysop/numerics/differentialOperator.py
index 9c8fca832..6dce5f092 100755
--- a/HySoP/hysop/numerics/differentialOperator.py
+++ b/HySoP/hysop/numerics/differentialOperator.py
@@ -5,7 +5,7 @@ Discrete stretching representation
 """
 from parmepy.constants import np, PARMES_REAL, ORDER
 from parmepy.numerics.method import NumMethod
-from parmepy.tools.cpu_data_transfer import Synchronize
+#from parmepy.tools.cpu_data_transfer import Synchronize
 
 
 class DifferentialOperator(NumMethod):
@@ -37,7 +37,7 @@ class DifferentialOperator(NumMethod):
 #    def setUp(self):
 #        pass
 
-    def __call__(self):
+    def __call__(self, synchroOp=None):
         """
         Apply operator.
         """
@@ -50,7 +50,6 @@ class DifferentialOperator(NumMethod):
         ind0 = np.arange(ind0a, ind0b)
         ind1 = np.arange(ind1a, ind1b)
         ind2 = np.arange(ind2a, ind2b)
-
         if self.choice.find('divWU') >= 0:
 
             # Ghosts synchronization
@@ -75,6 +74,7 @@ class DifferentialOperator(NumMethod):
 #            OpSynchronize = Synchronize(self.topology)
 #            OpSynchronize.apply(self.field2, self.field1)
 
+            synchroOp.apply()
             if self.method.find('FD_order2') >= 0:
                 raise ValueError("2nd order scheme Not yet implemented")
 ##           X components of temp and result
@@ -271,9 +271,9 @@ class DifferentialOperator(NumMethod):
         elif self.choice.find('gradV') >= 0:
 
             # Ghosts synchronization
-            OpSynchronize = Synchronize(self.topology)
-            OpSynchronize.apply(self.field2)
-
+            synchroOp.apply()
+#            OpSynchronize = Synchronize(self.topology)
+#            OpSynchronize.apply(self.field2)
             if self.method.find('FD_order2') >= 0:
                 raise ValueError("2nd order scheme Not yet implemented")
 
@@ -428,8 +428,9 @@ class DifferentialOperator(NumMethod):
         elif (self.choice.find('gradS') >= 0):
 
             # Ghosts synchronization
-            OpSynchronize = SynchronizeS(self.topology)
-            OpSynchronize.apply(self.field2)
+            #OpSynchronize = SynchronizeS(self.topology)
+            #OpSynchronize.apply(self.field2)
+            synchroOp.apply()
 
             if self.method.find('FD_order2') >= 0:
                 raise ValueError("2nd order scheme Not yet implemented")
@@ -467,9 +468,10 @@ class DifferentialOperator(NumMethod):
         elif self.choice.find('laplacianV') >= 0:
 
             # Ghosts synchronization
-            OpSynchronize = Synchronize(self.topology)
-            OpSynchronize.apply(self.field2)
-
+#            OpSynchronize = Synchronize(self.topology)
+#            OpSynchronize.apply(self.field2)
+            
+            synchroOp.apply()
             if self.method.find('FD_order2') >= 0:
                 raise ValueError("2nd order scheme Not yet implemented")
 
@@ -616,9 +618,10 @@ class DifferentialOperator(NumMethod):
         elif self.choice.find('curl') >= 0:
 
             # Ghosts synchronization
-            OpSynchronize = Synchronize(self.topology)
-            OpSynchronize.apply(self.field1)
+#            OpSynchronize = Synchronize(self.topology)
+#            OpSynchronize.apply(self.field1)
 
+            synchroOp.apply()
             if self.method.find('FD_order2') >= 0:
                 raise ValueError("2nd order scheme Not yet implemented")
 
diff --git a/HySoP/hysop/numerics/fct2op.py b/HySoP/hysop/numerics/fct2op.py
index 43b783d77..9cbf48054 100644
--- a/HySoP/hysop/numerics/fct2op.py
+++ b/HySoP/hysop/numerics/fct2op.py
@@ -11,9 +11,9 @@ import numpy as np
 
 class Fct2Op(object):
     """
-    Fct2Op 
-    
-    make the link between differentialOperator and the integration time scheme (Euler,RK,...)
+    Fct2Op
+    make the link between differentialOperator and
+    the integration time scheme (Euler,RK,...)
     """
 
     def __init__(self, field1, field2, method='', choice='', topology=None):
@@ -21,8 +21,9 @@ class Fct2Op(object):
         Constructor.
         store field2 to work with it in the apply function
 
-        @param field1 : vorticity field 
-        @param field2 : velocity field or grad(U) depending on 'choice' parameter
+        @param field1 : vorticity field
+        @param field2 : velocity field or grad(U)
+            depending on 'choice' parameter
         @param choice : can be 'gradV' for GradU.w or 'divWU' for div(wu)
         """
         self.field1 = field1
@@ -31,13 +32,13 @@ class Fct2Op(object):
         self.choice = choice
         self.topology = topology
 
-    def apply(self, t , y):
+    def apply(self, t, y):
         """
         Apply operator.
-        
         Fct2Op.apply is used like a def function
         """
-        if self.choice == 'gradV' : # field2 = gradU
+#        self.synchro.apply()
+        if self.choice == 'gradV':  # field2 = gradU
             temp0 = y[0][...] * 0.
             temp1 = y[0][...] * 0.
             temp2 = y[0][...] * 0.
@@ -48,48 +49,60 @@ class Fct2Op(object):
             ind1b = self.topology.mesh.local_end[1] + 1
             ind2a = self.topology.mesh.local_start[2]
             ind2b = self.topology.mesh.local_end[2] + 1
-            ind0 = np.arange(ind0a,ind0b)
-            ind1 = np.arange(ind1a,ind1b)
-            ind2 = np.arange(ind2a,ind2b)
+            ind0 = np.arange(ind0a, ind0b)
+            ind1 = np.arange(ind1a, ind1b)
+            ind2 = np.arange(ind2a, ind2b)
 
             ## Perform Matrix/Scalar product on local mesh
-            temp0[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]*\
-                                                            y[0][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] + \
-                                                        self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]* \
-                                                            y[1][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] + \
-                                                        self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]* \
-                                                            y[2][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
-            temp1[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = self.field2[3][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]* \
-                                                            y[0][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] + \
-                                                        self.field2[4][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]* \
-                                                            y[1][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] + \
-                                                        self.field2[5][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]* \
-                                                            y[2][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
-            temp2[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = self.field2[6][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]* \
-                                                            y[0][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] + \
-                                                        self.field2[7][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]* \
-                                                            y[1][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] + \
-                                                        self.field2[8][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]* \
-                                                            y[2][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
-#            result = np.concatenate((np.array([temp0]),np.array([temp1]),np.array([temp2])))
+            temp0[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] =\
+                self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
+                y[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
+                self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
+                y[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
+                self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
+                y[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b]
+            temp1[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] =\
+                self.field2[3][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
+                y[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
+                self.field2[4][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
+                y[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
+                self.field2[5][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
+                y[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b]
+            temp2[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] =\
+                self.field2[6][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
+                y[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
+                self.field2[7][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
+                y[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
+                self.field2[8][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
+                y[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b]
+#            result = np.concatenate((
+#                np.array([temp0]),np.array([temp1]),np.array([temp2])))
             return temp0, temp1, temp2
 
-        if self.choice == 'divWU' : # field2 = velocity field
-            result = DifferentialOperator(y, self.field2, self.method, self.choice).__call__()
+        if self.choice == 'divWU':  # field2 = velocity field
+            result = DifferentialOperator(
+                y, self.field2, self.method, self.choice).__call__()
             return result
 
-        if self.choice == 'hessianU' : # field2 = velocity field
-            gradU, maxgersh = DifferentialOperator(self.field1, self.field2, self.method, choice='gradV').__call__()
-            GradUline1 = np.concatenate(([gradU[0]],[gradU[1]],[gradU[2]]))
-            GradUline2 = np.concatenate(([gradU[3]],[gradU[4]],[gradU[5]]))
-            GradUline3 = np.concatenate(([gradU[6]],[gradU[7]],[gradU[8]]))
-            result1, maxgersh = DifferentialOperator(self.field1, GradUline1, self.method, choice='gradV').__call__()
-            result2, maxgersh = DifferentialOperator(self.field1, GradUline2, self.method, choice='gradV').__call__()
-            result3, maxgersh = DifferentialOperator(self.field1, GradUline3, self.method, choice='gradV').__call__()
+        if self.choice == 'hessianU':  # field2 = velocity field
+            gradU, maxgersh = DifferentialOperator(
+                self.field1, self.field2,
+                self.method, choice='gradV').__call__()
+            GradUline1 = np.concatenate(([gradU[0]], [gradU[1]], [gradU[2]]))
+            GradUline2 = np.concatenate(([gradU[3]], [gradU[4]], [gradU[5]]))
+            GradUline3 = np.concatenate(([gradU[6]], [gradU[7]], [gradU[8]]))
+            result1, maxgersh = DifferentialOperator(
+                self.field1, GradUline1,
+                self.method, choice='gradV').__call__()
+            result2, maxgersh = DifferentialOperator(
+                self.field1, GradUline2,
+                self.method, choice='gradV').__call__()
+            result3, maxgersh = DifferentialOperator(
+                self.field1, GradUline3,
+                self.method, choice='gradV').__call__()
             return result1, result2, result3
-        else :
-            print 'choice', choice , 'in Fct2Op is not defined'
-
+        else:
+            print 'choice', choice, 'in Fct2Op is not defined'
 
     def __str__(self):
         s = "fct2Op(DiscreteOperator). " + DiscreteOperator.__str__(self)
diff --git a/HySoP/hysop/numerics/integrators/runge_kutta3.py b/HySoP/hysop/numerics/integrators/runge_kutta3.py
index c5de5014f..d417ab7c2 100755
--- a/HySoP/hysop/numerics/integrators/runge_kutta3.py
+++ b/HySoP/hysop/numerics/integrators/runge_kutta3.py
@@ -41,8 +41,9 @@ class RK3(ODESolver):
 #        # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
 #        K2 = dt * np.asarray(f(t + dt / 3., y + K1[:, :, :, :] / 3.),
 #                             dtype=PARMES_REAL, order=ORDER)
-#        K3 = dt * np.asarray(f(t + dt * 2. / 3., y + K2[:, :, :, :] * 2. / 3.),
-#                             dtype=PARMES_REAL, order=ORDER)
+#        K3 = dt * np.asarray(
+#            f(t + dt * 2. / 3., y + K2[:, :, :, :] * 2. / 3.),
+#            dtype=PARMES_REAL, order=ORDER)
 #        # result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
 #        result = y + 1. / 4. * (K1 + 3. * K3)
 #        return result
@@ -51,15 +52,15 @@ class RK3(ODESolver):
         K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
 
         RHS2 = f(t + dt / 3., y + K1 / 3.)
-        K2 =  np.array([dt * RHS2[0], dt * RHS2[1], dt * RHS2[2]])
+        K2 = np.array([dt * RHS2[0], dt * RHS2[1], dt * RHS2[2]])
 
         RHS3 = f(t + dt * 2. / 3., y + K2 * 2. / 3.)
         K3 = np.array([dt * RHS3[0], dt * RHS3[1], dt * RHS3[2]])
 
         # result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
         return y[0] + 1. / 4. * (K1[0] + 3. * K3[0]), \
-               y[1] + 1. / 4. * (K1[1] + 3. * K3[1]), \
-               y[2] + 1. / 4. * (K1[2] + 3. * K3[2])
+            y[1] + 1. / 4. * (K1[1] + 3. * K3[1]), \
+            y[2] + 1. / 4. * (K1[2] + 3. * K3[2])
 
     def __str__(self):
         """ToString method"""
diff --git a/HySoP/hysop/operator/discrete/stretching.py b/HySoP/hysop/operator/discrete/stretching.py
index 3d55cbc29..46ab540e3 100755
--- a/HySoP/hysop/operator/discrete/stretching.py
+++ b/HySoP/hysop/operator/discrete/stretching.py
@@ -14,6 +14,7 @@ from parmepy.numerics.integrators.runge_kutta2 import RK2
 from parmepy.numerics.integrators.runge_kutta3 import RK3
 from parmepy.numerics.integrators.runge_kutta4 import RK4
 from parmepy.tools.timers import timed_function
+from parmepy.operator.discrete.synchronizeGhosts import SynchronizeGhosts_d
 
 
 class Stretching_d(DiscreteOperator):
@@ -103,8 +104,12 @@ class Stretching_d(DiscreteOperator):
                     self.vorticity, self.velocity,
                     method='FD_order4',
                     choice='gradV')
-
-                gradResult, maxArray, divergence = self.stretchOp()
+                #une synchro
+                self.synchro = SynchronizeGhosts_d(
+                    self.vorticity.data,
+                    self.vorticity.topology, 'greaterSend')
+                self.synchro.setUp()
+                gradResult, maxArray, divergence = self.stretchOp(self.synchro)
 
                 print "maxadv, maxgersh", maxArray[1], maxArray[0]
                 print "dtAdvec", 0.5 / maxArray[1]
@@ -115,11 +120,12 @@ class Stretching_d(DiscreteOperator):
                 self.dtstretch = dt / float(self.ndt)
 
                 ## fct2Op
-                self.function = Fct2Op(self.vorticity.data,
+                self.function = Fct2Op(self.vorticity,
                                        gradResult,
                                        method='FD_order4',
                                        choice='gradV',
                                        topology=self.vorticity.topology)
+                self.synchro.apply()
 
             ## Time integration
             t = simulation.time
@@ -128,11 +134,12 @@ class Stretching_d(DiscreteOperator):
 #                self.vorticity.data[...] = methodInt(self.function.apply,
 #                                                     t, self.dtstretch,
 #                                                     self.vorticity.data)
-                self.vorticity.data[0][...], \
-                self.vorticity.data[1][...], \
-                self.vorticity.data[2][...] = methodInt(self.function.apply,
-                                                        t, self.dtstretch,
-                                                        self.vorticity.data)
+                self.vorticity.data[0][...], self.vorticity.data[1][...], \
+                    self.vorticity.data[2][...] = methodInt(
+                        self.function.apply,
+                        t, self.dtstretch,
+                        self.vorticity.data)
+                self.synchro.apply()
 #            print 'norm omega', LA.norm(self.vorticity.data)
 
     def stabilityTest(self, dt, maxi, num_method, cststretch):
diff --git a/HySoP/hysop/operator/discrete/synchronizeGhosts.py b/HySoP/hysop/operator/discrete/synchronizeGhosts.py
index e8c71b03e..9b3e36dcc 100644
--- a/HySoP/hysop/operator/discrete/synchronizeGhosts.py
+++ b/HySoP/hysop/operator/discrete/synchronizeGhosts.py
@@ -25,8 +25,8 @@ class SynchronizeGhosts_d(DiscreteOperator):
         @param topology : the topology common to all fields.
         @param transferMethod : which type of exchange is used.
         """
-        #DiscreteOperator.__init__(self, fieldslist, method=transferMethod,
-        #                          name="Ghosts Synchronization")
+        DiscreteOperator.__init__(self, fieldslist, method=transferMethod,
+                                  name="Ghosts Synchronization")
         self.input = fieldslist
         self.output = fieldslist
         self.compute_time = 0.
@@ -49,7 +49,8 @@ class SynchronizeGhosts_d(DiscreteOperator):
                 self.topology.dims,
                 periods=self.topology.periods)
         ghosts = self.topology.ghosts
-        resolution = self.topology.localGridResolution
+        resolution = self.topology.mesh.resolution
+        self.fieldslist = np.asarray([self.fieldslist])
         if (self.transferMethod == 'SendRecv'):
             # Allocation of space for reciev data
             self.dataRecvX = np.zeros((
@@ -117,14 +118,14 @@ class SynchronizeGhosts_d(DiscreteOperator):
     @debug
     @timed_function
     def apply(self, simulation=None):
-        self.resolution = self.topology.localGridResolution
-        resolution = self.topology.localGridResolution
+        self.resolution = self.topology.mesh.resolution
+        resolution = self.topology.mesh.resolution
         ghosts = self.topology.ghosts
         dim = self.topology.dim
         # Do the send/recv as defined in setup.
         #### if args is None: \todo : A verifier.
-        args = self.fieldslist
-        self.fieldslist = np.asarray([self.fieldslist])
+#        args = self.fieldslist
+#        self.fieldslist = np.asarray([self.fieldslist])
         if (self.transferMethod == 'SendRecv'):
             for f in fieldslist:
                 # DOWN
@@ -233,6 +234,7 @@ class SynchronizeGhosts_d(DiscreteOperator):
                         # WEST
                         if (self.tabSort[i, 1] == 0):
                             for j in xrange(dim):
+                                #print '\nTEST', resolution,  f[j][...].shape, '\t' , ghosts
                                 f[j][resolution[0] - ghosts[0]: resolution[0],
                                      :, :] =\
                                     f[j][ghosts[0]:2 * ghosts[0], :, :]
@@ -333,6 +335,10 @@ class SynchronizeGhosts_d(DiscreteOperator):
                     tabfinal[deb, :] = tab[deb, :]
                 else:
                     tabNew = tab[deb:i, :]
+#                    for l in xrange(int(tabNew[:,1].shape[0] / 2)):
+#                        tmp = tabNew[2 * l, 1]
+#                        tabNew[2 * l, 1] = tabNew[2 * l + 1,1]
+#                        tabNew[2 * l + 1, 1] = tmp
                     tabNew = tabNew[np.invert(tabNew[:, 1].argsort())]
                     tabfinal[deb:i, :] = tabNew
                 deb = i
diff --git a/HySoP/hysop/operator/synchronizeGhosts.py b/HySoP/hysop/operator/synchronizeGhosts.py
index 9ecf5140e..652a8e99a 100644
--- a/HySoP/hysop/operator/synchronizeGhosts.py
+++ b/HySoP/hysop/operator/synchronizeGhosts.py
@@ -6,7 +6,8 @@ Update ghost points for some fields defined on a specific topology.
 from parmepy.operator.continuous import Operator
 from parmepy.operator.discrete.synchronizeGhosts import SynchronizeGhosts_d
 from parmepy.constants import debug
-
+from parmepy.fields.vector import VectorField
+from parmepy.fields.scalar import ScalarField
 
 class SynchronizeGhosts(Operator):
     """
@@ -15,7 +16,7 @@ class SynchronizeGhosts(Operator):
 
     @debug
     def __init__(self, fieldslist, resolution=None,
-                 method='', **other_config):
+                 method='', topology=None, **other_config):
         """
         Defines a way to send/recv values at ghosts points for a list
         of fields, discretized on a given topology.
@@ -31,10 +32,11 @@ class SynchronizeGhosts(Operator):
         self.fieldslist = fieldslist
         self.resolution = resolution
         self.method = method
+        self.topology = topology
         self.config = other_config
 
     @debug
-    def setUp(self, ghosts):
+    def setUp(self, ghosts=None):
         """
         SynchroGhost operator discretization method.
         Create an discrete SynchroGhost operator from given specifications.
@@ -43,15 +45,28 @@ class SynchronizeGhosts(Operator):
         Operator.setUp(self)
         for v in self.variables:
         # the topology for v ...
-            topo = self.domain.getOrCreateTopology(
-                self.domain.dimension,
-                self.resolution[v],
-                ghosts=ghosts)
-            topo.ghosts = ghosts
-        # ... and the corresponding discrete field
-            self.discreteFields[v] = v.discretize(topo)
-            self.discreteOperator = \
-                SynchronizeGhosts_d(self.discreteFields[v], topo, self.method)
+            if self.topology is not None:
+                topo = self.topology
+            else:
+                if ghosts is None:
+                    topo = self.domain.getOrCreateTopology(
+                        self.domain.dimension,
+                        self.resolution[v])
+                else:
+                    topo = self.domain.getOrCreateTopology(
+                        self.domain.dimension,
+                        self.resolution[v],
+                        ghosts=ghosts)
+                topo.ghosts = ghosts
+            # ... and the corresponding discrete field
+                if v.__class__ is VectorField or v.__class__ is ScalarField:
+                    self.discreteOperator = \
+                        SynchronizeGhosts_d(v, topo, self.method)
+                else:
+                    self.discreteFields[v] = v.discretize(topo)
+                    self.discreteOperator = \
+                        SynchronizeGhosts_d(self.discreteFields[v],
+                        topo, self.method)
 
         self.discreteOperator.setUp()
         self._isUpToDate = True
diff --git a/HySoP/hysop/operator/tests/test_Stretching.py b/HySoP/hysop/operator/tests/test_Stretching.py
index ab92bce44..582613a36 100755
--- a/HySoP/hysop/operator/tests/test_Stretching.py
+++ b/HySoP/hysop/operator/tests/test_Stretching.py
@@ -3,9 +3,11 @@ import time
 
 import parmepy as pp
 import numpy as np
-from parmepy.fields.analytical import AnalyticalField
+from parmepy.fields.continuous import Field
 from parmepy.operator.stretching import Stretching
 from parmepy.problem.navier_stokes import NSProblem
+from parmepy.problem.simulation import Simulation
+from parmepy.mpi.main_var import main_comm, main_rank
 
 
 def computeVel(x, y, z):
@@ -65,24 +67,29 @@ def test_Stretching():
     box = pp.Box(dim, length=boxLength, origin=boxMin)
 
     ## Fields
-    velo = AnalyticalField(domain=box, formula=computeVel,
-                           name='Velocity', vector=True)
-    vorti = AnalyticalField(domain=box, formula=computeVort,
-                            name='Vorticity', vector=True)
-
+    velo = Field(
+        domain=box, formula=computeVel,
+        name='Velocity', isVector=True)
+    vorti = Field(
+        domain=box, formula=computeVort,
+        name='Vorticity', isVector=True)
+
+    velo.initialize()
+    vorti.initialize()
     ## Operators
     stretch = Stretching(velo, vorti,
                          resolutions={velo: nbElem,
                                       vorti: nbElem},
                          method='FD_order4 RK3',
-                         propertyOp='divConservation'
-                        )
-
+                         propertyOp='divConservation')
+    stretch.setUp()
+    simulation = Simulation(0, timeStep, 20)
     ##Problem
-    pb = NSProblem([stretch])
+    pb = NSProblem([stretch], simulation)
 
     ## Setting solver to Problem
-    pb.setUp(finalTime, timeStep)
+#    pb.setUp(simulation)
+    pb.setUp()
 
     t1 = time.time()
 
@@ -90,11 +97,11 @@ def test_Stretching():
     pb.solve()
 
     tf = time.time()
-
-    print "\n"
-    print "Total time : ", tf - t0, "sec (CPU)"
-    print "Init time : ", t1 - t0, "sec (CPU)"
-    print "Solving time : ", tf - t1, "sec (CPU)"
+    if main_rank == 0:
+        print "\n"
+        print "Total time : ", tf - t0, "sec (CPU)"
+        print "Init time : ", t1 - t0, "sec (CPU)"
+        print "Solving time : ", tf - t1, "sec (CPU)"
 
 
 if __name__ == "__main__":
diff --git a/HySoP/hysop/operator/tests/test_synchronizeGhosts.py b/HySoP/hysop/operator/tests/test_synchronizeGhosts.py
index ddde726f4..66f79bf85 100644
--- a/HySoP/hysop/operator/tests/test_synchronizeGhosts.py
+++ b/HySoP/hysop/operator/tests/test_synchronizeGhosts.py
@@ -9,7 +9,7 @@ from parmepy.mpi.topology import Cartesian
 #from parmepy.fields.analytical import AnalyticalField
 from parmepy.mpi.main_var import main_comm, main_size, main_rank, MPI
 from parmepy.operator.analytic import Analytic
-
+from parmepy.problem.simulation import Simulation
 
 def computeVel(x, y, z):
     vx = x
@@ -31,18 +31,18 @@ def test_Synchro():
     field = Field(domain=dom, formula=computeVel, name='velocity', isVector=True)
 #    velo = AnalyticalField(domain=dom, formula=computeVel,
 #                           name='Velocity', isVector=True)
-
 #    Synchrotest = SynchronizeGhosts(velo, resolution={velo: nbElem},
     Synchrotest = SynchronizeGhosts(field, resolution={field: nbElem},
                                     method='greaterSend')
 
+    simu = Simulation(tinit=0., tend=1., timeStep=0.5)
     if main_rank ==0 :
         print "Setup Synchro"
     Synchrotest.setUp(nbGhosts)
     field.initialize()
     if main_rank ==0 :
         print "Apply Synchro"
-    Synchrotest.apply(t=None, dt=None, ite=None)
+    Synchrotest.apply(simu)
 
 if __name__ == "__main__":
 #    print "Rank", main_rank
diff --git a/HySoP/hysop/problem/problem.py b/HySoP/hysop/problem/problem.py
index 0c1a9e9f0..9d8df27c2 100644
--- a/HySoP/hysop/problem/problem.py
+++ b/HySoP/hysop/problem/problem.py
@@ -93,7 +93,7 @@ class Problem(object):
                 if not isinstance(op, Redistribute):
                     op.setUp()
 
-        if  __VERBOSE__:
+        if __VERBOSE__:
             print "==== Variables initialization ===="
 
         # Build variables list to initialize
-- 
GitLab