diff --git a/HySoP/hysop/domain/domain.py b/HySoP/hysop/domain/domain.py
index baa7b714aa6929b2981b04c604ee4c5ac2c0caf9..505dfdd39ad1c4ff62feba3ceefa2720d160f833 100644
--- a/HySoP/hysop/domain/domain.py
+++ b/HySoP/hysop/domain/domain.py
@@ -34,7 +34,7 @@ class Domain(object):
     def getOrCreateTopology(self, topoDim, gridResolution,
                             topoResolution=None, ghosts=None,
                             precomputed=False, localres=None,
-                            offset=None):
+                            offset=None, fixedResolution=False):
         """
         This routines checks if a topology is present in the list
         of topologies of the domain.
@@ -58,8 +58,12 @@ class Domain(object):
                                                           ghosts=ghosts)
 
         else:
-            newTopo = Cartesian.withResolution(self, topoResolution,
-                                               gridResolution, ghosts=ghosts)
+            if fixedResolution:
+                topoCreation = Cartesian.withResolutionFixed
+            else:
+                topoCreation = Cartesian.withResolution
+            newTopo = topoCreation(self, topoResolution,
+                                   gridResolution, ghosts=ghosts)
         otherid = self.register(newTopo)
         if otherid < 0:
             otherid = newTopo.getId()
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
index e3913378a2b5115aec7846f1ec005e3b41cfa8e2..c980f51bfc54268b9dc9c9eff71f837b877d9f74 100644
--- a/HySoP/hysop/fields/continuous.py
+++ b/HySoP/hysop/fields/continuous.py
@@ -10,7 +10,8 @@ from parmepy.mpi import main_rank
 
 
 class Field(object):
-    """ Continuous field defined on a physical domain.
+    """
+    Continuous field defined on a physical domain.
 
     This object handles a dictionnary of discrete fields
     (from 0 to any number).
@@ -49,6 +50,7 @@ class Field(object):
         (default = 0 all over the domain)
         @param name : name of the variable.
         @param isVector : true if the field is a vector.
+        @param nbComponents : Components number (1 for scalar fields).
         """
         ## Physical domain of definition of the variable.
         self.domain = domain
diff --git a/HySoP/hysop/gpu/cl_src/advection/basic.cl b/HySoP/hysop/gpu/cl_src/advection/basic.cl
index 2bc92e4f061ddde164d56f8ec3d522bea55a6bdc..2a0eecef52ea4c3281e9128adfc3438b9040bc55 100644
--- a/HySoP/hysop/gpu/cl_src/advection/basic.cl
+++ b/HySoP/hysop/gpu/cl_src/advection/basic.cl
@@ -1,5 +1,5 @@
 /**
- * @file basic.cl
+ * @file advection/basic.cl
  * Advection function, vectorized version, no use of builtins functions.
  */
 
diff --git a/HySoP/hysop/gpu/cl_src/advection/basic_noVec.cl b/HySoP/hysop/gpu/cl_src/advection/basic_noVec.cl
index 2be4a657f41019424b2bf158765337defb1812fb..e062208bd7f009a70ecdb7b4c72682701907f1dc 100644
--- a/HySoP/hysop/gpu/cl_src/advection/basic_noVec.cl
+++ b/HySoP/hysop/gpu/cl_src/advection/basic_noVec.cl
@@ -1,5 +1,5 @@
 /**
- * @file basic_noVec.cl
+ * @file advection/basic_noVec.cl
  * Advection function, basic version
  */
 
diff --git a/HySoP/hysop/gpu/cl_src/advection/basic_rk4_noVec.cl b/HySoP/hysop/gpu/cl_src/advection/basic_rk4_noVec.cl
index b571c4bc4702a3bd3669cd0dcb5b46bfe1e2716b..0e2317f43b73fa845462230f7cd87536cda0cc81 100644
--- a/HySoP/hysop/gpu/cl_src/advection/basic_rk4_noVec.cl
+++ b/HySoP/hysop/gpu/cl_src/advection/basic_rk4_noVec.cl
@@ -14,7 +14,7 @@ float advection(uint i, float dt, float dx, float invdx, __local float* gvelo_lo
  * @param dt Time step.
  * @param dx Space step.
  * @param invdx 1/dx.
- * @param gvelo Global velocity field.
+ * @param gvelo_loc Local velocity field.
  * @return Particle position
  *
  * @remark NB_I, NB_II, NB_III : points number in directions from 1st varying index to last.
diff --git a/HySoP/hysop/gpu/cl_src/advection/builtin_noVec.cl b/HySoP/hysop/gpu/cl_src/advection/builtin_noVec.cl
index a6ba8b975cc467646bb0afdeb113eaa3b08fed2a..6376911f4d7caf95961067468d6209d47b1ada8c 100644
--- a/HySoP/hysop/gpu/cl_src/advection/builtin_noVec.cl
+++ b/HySoP/hysop/gpu/cl_src/advection/builtin_noVec.cl
@@ -14,7 +14,7 @@ float advection(uint i, float dt, float dx, float invdx, __local float* gvelo_lo
  * @param dt Time step.
  * @param dx Space step.
  * @param invdx 1/dx.
- * @param gvelo Global velocity field.
+ * @param gvelo_loc Local velocity field.
  * @return Particle position
  *
  * @remark NB_I, NB_II, NB_III : points number in directions from 1st varying index to last.
diff --git a/HySoP/hysop/gpu/cl_src/advection/builtin_rk4_noVec.cl b/HySoP/hysop/gpu/cl_src/advection/builtin_rk4_noVec.cl
index da79c8950fae3ce5678949ac7a080cec1030df89..f039542d6b5d285bfe55f6b7eb0bf6e652998c45 100644
--- a/HySoP/hysop/gpu/cl_src/advection/builtin_rk4_noVec.cl
+++ b/HySoP/hysop/gpu/cl_src/advection/builtin_rk4_noVec.cl
@@ -14,7 +14,7 @@ float advection(uint i, float dt, float dx, float invdx, __local float* gvelo_lo
  * @param dt Time step.
  * @param dx Space step.
  * @param invdx 1/dx.
- * @param gvelo Global velocity field.
+ * @param gvelo_loc Local velocity field.
  * @return Particle position
  *
  * @remark NB_I, NB_II, NB_III : points number in directions from 1st varying index to last.
diff --git a/HySoP/hysop/gpu/cl_src/remeshing/basic.cl b/HySoP/hysop/gpu/cl_src/remeshing/basic.cl
index 73596e388c365dbfb3b9e7ea3b6dfb6bdc5169a2..0408a80f7f1b7196f6ce752157936f05d8e6080b 100644
--- a/HySoP/hysop/gpu/cl_src/remeshing/basic.cl
+++ b/HySoP/hysop/gpu/cl_src/remeshing/basic.cl
@@ -1,5 +1,5 @@
 /**
- * @file basic.cl
+ * @file remeshing/basic.cl
  * Remeshing function, vectorized version.
  */
 
diff --git a/HySoP/hysop/gpu/cl_src/remeshing/basic_noVec.cl b/HySoP/hysop/gpu/cl_src/remeshing/basic_noVec.cl
index a9d67638dc6b2152773b4f86edab625b030447e1..29d0c7ee456f61fd8f9ce6653119a478d2e103e2 100644
--- a/HySoP/hysop/gpu/cl_src/remeshing/basic_noVec.cl
+++ b/HySoP/hysop/gpu/cl_src/remeshing/basic_noVec.cl
@@ -1,5 +1,5 @@
 /**
- * @file basic_noVec.cl
+ * @file remeshing/basic_noVec.cl
  * Remeshing function, vectorized version.
  */
 
diff --git a/HySoP/hysop/gpu/cl_src/remeshing/weights_noVec_builtin.cl b/HySoP/hysop/gpu/cl_src/remeshing/weights_noVec_builtin.cl
index 0a703d65a13b3f89d7fb00a2c5efa9f815bb1ce2..9314559ea34fd50d91f3dfe83d85e815e18b4e73 100644
--- a/HySoP/hysop/gpu/cl_src/remeshing/weights_noVec_builtin.cl
+++ b/HySoP/hysop/gpu/cl_src/remeshing/weights_noVec_builtin.cl
@@ -1,5 +1,5 @@
 /**
- * @file weights_npVec_builtin.cl
+ * @file weights_noVec_builtin.cl
  * Remeshing formulas, vectorized version, use of builtin OpenCL fma.
  * Polynomials under Horner form.
  */
diff --git a/HySoP/hysop/operator/advection.py b/HySoP/hysop/operator/advection.py
index c5a3e3dfec040f79314252087a9d5e3c3b5721a3..a1de2517fc697c53ddda3a88e23623376f45e927 100644
--- a/HySoP/hysop/operator/advection.py
+++ b/HySoP/hysop/operator/advection.py
@@ -9,7 +9,7 @@ from parmepy.tools.timers import Timer
 from parmepy.fields.continuous import Field
 from parmepy.mpi import main_size
 from parmepy.operator.redistribute import Redistribute
-
+from parmepy.operator.advectionDir import AdvectionDir
 
 class Advection(Operator):
     """
@@ -61,16 +61,16 @@ class Advection(Operator):
         v = [velocity]
         if isinstance(advectedFields, list):
             [v.append(f) for f in advectedFields]
-            self.output = advectedFields
         else:
             v.append(advectedFields)
-            self.output = [advectedFields]
         Operator.__init__(self, v, method, topo=topo,
                           ghosts=ghosts, name="Advection")
         ## Transport velocity
         self.velocity = velocity
         ## Transported fields
-        self.advectedFields = advectedFields
+        self.advectedFields = [v for v in self.variables
+                               if not v is self.velocity]
+        self.output = self.advectedFields
         ## Resolutions of the fields. Dictionnary with
         ## fields as keys.
         ## Example :
@@ -78,29 +78,12 @@ class Advection(Operator):
         self.resolutions = resolutions
         ## Extra parameters (depend on the method)
         self.config = other_config
-        self.input = [self.velocity]
+        self.input = [v for v in self.variables]
         self._isSCALES = False
         self._old_dir = 0
-        dim = self.velocity.dimension
         self.splitting = []
-        if splittingConfig == 'o2_FullHalf':
-            ## Half timestep in all directions
-            [self.splitting.append((i, 0.5)) for i in xrange(dim)]
-            [self.splitting.append((dim - 1 - i, 0.5)) for i in xrange(dim)]
-        elif splittingConfig == 'o1':
-            [self.splitting.append((i, 1.)) for i in xrange(dim)]
-        elif splittingConfig == 'x_only':
-            self.splitting.append((0, 1.))
-        elif splittingConfig == 'y_only':
-            self.splitting.append((1, 1.))
-        elif splittingConfig == 'z_only':
-            self.splitting.append((2, 1.))
-        else:
-            ## Half timestep in all directions but last
-            [self.splitting.append((i, 0.5)) for i in xrange(dim - 1)]
-            self.splitting.append((dim - 1, 1.))
-            [self.splitting.append((dim - 2 - i, 0.5))
-             for i in xrange(dim - 1)]
+        self.splittingConfig = splittingConfig
+        self.advecDir = [None] * self.velocity.dimension
 
     @debug
     def setUp(self):
@@ -131,10 +114,9 @@ class Advection(Operator):
           - 'm8prime' : M8prime formula
           - 'l6star' : Lambda6star formula
         """
-
         # --- Advection solver from SCALES ---
         if self.method.find('scales') >= 0:
-            if not self.advectedFields.dimension == 3:
+            if not self.advectedFields[0].dimension == 3:
                     raise ValueError("Scales Advection not implemented in 2D.")
             if self.ghosts is not None:
                 if (self.ghosts > 0).any():
@@ -154,7 +136,7 @@ class Advection(Operator):
 
             # - Create the topologies (get param from scales) -
             # Scales nbcells equals resolutions - 1
-            nbcells = np.asarray(self.resolutions[self.advectedFields],
+            nbcells = np.asarray(self.resolutions[self.advectedFields[0]],
                                  dtype=PARMES_INDEX) - 1
             topodims = [1, 1, main_size]
             scalesres, scalesoffset, stab_coeff = \
@@ -182,14 +164,13 @@ class Advection(Operator):
                     # ... and the corresponding discrete field
                     self.discreteFields[v] = v.discretize(topo)
 
-            advectedFields = [self.discreteFields[f]
-                              for f in self.discreteFields.keys()
-                              if f is not self.velocity]
+            advectedDiscreteFields = [self.discreteFields[f]
+                                      for f in self.advectedFields]
 
             # - Create the discreteOperator from the list of discrete fields -
             self.discreteOperator = ScalesAdvection(
                 self.discreteFields[self.velocity],
-                advectedFields, method=self.method,
+                advectedDiscreteFields, method=self.method,
                 **self.config)
 
             # -- Final set up --
@@ -198,115 +179,52 @@ class Advection(Operator):
 
         # --- GPU or pure-python advection ---
         else:
-            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:
-                        topo = self.domain.getOrCreateTopology(
-                            self.domain.dimension, self.resolutions[v])
-                        self.discreteFields[v] = v.discretize(topo)
-
-                advectedFields = [self.discreteFields[f]
-                                  for f in self.discreteFields.keys()
-                                  if f is not self.velocity]
-                particles_advectedFields = [
-                    Field(adF.topology.domain, "Particle_AdvectedFields",
-                          isVector=adF.isVector).discretize(adF.topology)
-                    for adF in advectedFields]
-
-                if self.method.find('gpu_2k') >= 0:
-                    from parmepy.gpu.gpu_particle_advection_2k \
-                        import GPUParticleAdvection2k as ParticleAdvection
-                    particles_positions = \
-                        Field(advectedFields[0].topology.domain,
-                              "Particle_Position",  isVector=False
-                              ).discretize(advectedFields[0].topology)
-                elif self.method.find('gpu_1k') >= 0:
-                    from parmepy.gpu.gpu_particle_advection_1k \
-                        import GPUParticleAdvection1k as ParticleAdvection
-                    particles_positions = None
-                else:
-                    from parmepy.operator.discrete.particle_advection \
-                        import ParticleAdvection
-                    particles_positions = \
-                        Field(advectedFields[0].topology.domain,
-                              "Particle_Position",  isVector=False
-                              ).discretize(advectedFields[0].topology)
-
-                self.discreteOperator = [
-                    ParticleAdvection(
-                        self.discreteFields[self.velocity],
-                        advectedFields, d, method=self.method,
-                        part_position=particles_positions,
-                        part_advectedFields=particles_advectedFields,
-                        **self.config)
-                    for d in xrange(self.velocity.dimension)]
-                # -- Final set up --
-                for dop in self.discreteOperator:
-                    dop.setUp()
-                self._isUpToDate = True
+            dim = self.velocity.dimension
+            # Splitting configuration
+            if self.splittingConfig == 'o2_FullHalf':
+                ## Half timestep in all directions
+                [self.splitting.append((i, 0.5)) for i in xrange(dim)]
+                [self.splitting.append((dim - 1 - i, 0.5)) for i in xrange(dim)]
+            elif self.splittingConfig == 'o1':
+                [self.splitting.append((i, 1.)) for i in xrange(dim)]
+            elif self.splittingConfig == 'x_only':
+                self.splitting.append((0, 1.))
+            elif self.splittingConfig == 'y_only':
+                self.splitting.append((1, 1.))
+            elif self.splittingConfig == 'z_only':
+                self.splitting.append((2, 1.))
             else:
-                if self._predefinedTopo is not None:
-                    raise ValueError("User-defined topology is not\
-                    (yet) allowed for advection.")
-                self.discreteOperator = self.domain.dimension * [None]
-                for d in xrange(self.domain.dimension):
-                    topodims = np.ones((self.domain.dimension))
-                    if d == self.domain.dimension - 1:
-                        topodims[0] = main_size
-                    else:
-                        topodims[-1] = main_size
-                    for v in self.variables:
-                        topo = self.domain.getOrCreateTopology(
-                            self.domain.dimension, self.resolutions[v],
-                            topoResolution=topodims)
-                        self.discreteFields[v] = v.discretize(topo)
-
-                    advectedFields = [self.discreteFields[f]
-                                      for f in self.discreteFields.keys()
-                                      if f is not self.velocity]
-                    particles_advectedFields = [
-                        Field(adF.topology.domain, "Particle_AdvectedFields",
-                              isVector=adF.isVector).discretize(adF.topology)
-                        for adF in advectedFields]
-
-                    if self.method.find('gpu_2k') >= 0:
-                        from parmepy.gpu.gpu_particle_advection_2k \
-                            import GPUParticleAdvection2k as ParticleAdvection
-                        particles_positions = \
-                            Field(advectedFields[0].topology.domain,
-                                  "Particle_Position",  isVector=False
-                                  ).discretize(advectedFields[0].topology)
-                    elif self.method.find('gpu_1k') >= 0:
-                        from parmepy.gpu.gpu_particle_advection_1k \
-                            import GPUParticleAdvection1k as ParticleAdvection
-                        particles_positions = None
-                    else:
-                        print "Using default advection operator"
-                        from parmepy.operator.discrete.particle_advection \
-                            import ParticleAdvection
-                        particles_positions = \
-                            Field(advectedFields[0].topology.domain,
-                                  "Particle_Position",  isVector=False
-                                  ).discretize(advectedFields[0].topology)
+                ## Half timestep in all directions but last
+                [self.splitting.append((i, 0.5)) for i in xrange(dim - 1)]
+                self.splitting.append((dim - 1, 1.))
+                [self.splitting.append((dim - 2 - i, 0.5))
+                 for i in xrange(dim - 1)]
 
-                    self.discreteOperator[d] = ParticleAdvection(
-                        self.discreteFields[self.velocity],
-                        advectedFields, d, method=self.method,
-                        part_position=particles_positions,
-                        part_advectedFields=particles_advectedFields,
-                        **self.config)
+            particles_advectedFields = [
+                Field(adF.domain, "Particle_AdvectedFields",
+                      isVector=adF.isVector) for adF in self.advectedFields]
+            if self.method.find('gpu_1k') >= 0:
+                particles_positions = None
+            else:
+                particles_positions = \
+                    Field(self.advectedFields[0].domain,
+                          "Particle_Position",  isVector=False
+                          )
 
-                # -- Final set up --
-                for dop in self.discreteOperator:
-                    dop.setUp()
+            # Directional continuous Advection operators
+            for i in xrange(dim):
+                self.advecDir[i] = AdvectionDir(
+                    self.velocity, self.advectedFields, i,
+                    particles_advectedFields, particles_positions,
+                    self.resolutions,
+                    method=self.method, topo=self._predefinedTopo,
+                    ghosts=self.ghosts, **self.config)
+            for i in xrange(dim):
+                self.advecDir[i].setUp()
 
-                # # Build bridges
-                self.bridges = self.domain.dimension * [None]
+            # Build bridges
+            self.bridges = self.domain.dimension * [None]
+            if main_size > 1:
                 for dfrom in xrange(self.domain.dimension):
                     self.bridges[dfrom] = self.domain.dimension * [None]
                     for dto in xrange(self.domain.dimension):
@@ -314,12 +232,11 @@ class Advection(Operator):
                             self.bridges[dfrom][dto] = None
                         else:
                             self.bridges[dfrom][dto] = Redistribute(
-                                self.variables,
-                                self.discreteOperator[dfrom],
-                                self.discreteOperator[dto],
-                                method='FromDiscreteOperators')
+                                self.input,
+                                self.advecDir[dfrom],
+                                self.advecDir[dto])
                             self.bridges[dfrom][dto].setUp()
-                self._isUpToDate = True
+
 
     @debug
     def apply(self, simulation=None):
@@ -336,7 +253,7 @@ class Advection(Operator):
         else:
             if main_size == 1:  # Sequential
                 for split_id, split in enumerate(self.splitting):
-                    self.discreteOperator[split[0]].apply(
+                    self.advecDir[split[0]].apply(
                         simulation, split[1], split_id, self._old_dir)
                     self._old_dir = split[0]
             else:  # Parallel, must use bridges to exchange data
@@ -344,7 +261,7 @@ class Advection(Operator):
                     if not self._old_dir == split[0]:
                         self.bridges[self._old_dir][split[0]].apply()
                         self.bridges[self._old_dir][split[0]].wait()
-                    self.discreteOperator[split[0]].apply(
+                    self.advecDir[split[0]].apply(
                         simulation, split[1], split_id, self._old_dir)
                     self._old_dir = split[0]
 
@@ -354,7 +271,7 @@ class Advection(Operator):
         Memory cleaning.
         """
         self.timer = Timer(self)
-        for dop in self.discreteOperator:
+        for dop in self.advecDir:
             dop.finalize()
             self.timer = self.timer + dop.timer
         print self.timer
diff --git a/HySoP/hysop/operator/advection_dir.py b/HySoP/hysop/operator/advection_dir.py
new file mode 100644
index 0000000000000000000000000000000000000000..b82f6a2126d915b1bb3229e6016d737571173158
--- /dev/null
+++ b/HySoP/hysop/operator/advection_dir.py
@@ -0,0 +1,144 @@
+"""
+@file advectionDir.py
+
+Advection of a field in a single direction.
+"""
+from parmepy.constants import debug, np, PARMES_INDEX
+from parmepy.operator.continuous import Operator
+from parmepy.mpi import main_size
+from parmepy.fields.continuous import Field
+from parmepy.tools.timers import Timer
+
+class AdvectionDir(Operator):
+    """
+    Advection of a field,
+    \f{eqnarray*}
+    X = Op(X,velocity)
+    \f} for :
+    \f{eqnarray*}
+    \frac{\partial{X}}{\partial{t}} + velocity.\nabla X = 0
+    \f}
+    Note : we assume incompressible flow.
+
+    Computations are performed in a given direction.
+    """
+
+    @debug
+    def __init__(self, velocity, advectedFields, d, particle_fields,
+                 particle_positions, resolutions,
+                 method='', topo=None, ghosts=None, **other_config):
+        v = [velocity]
+        [v.append(f) for f in advectedFields]
+        self.output = advectedFields
+        Operator.__init__(self, v, method, topo=topo,
+                          ghosts=ghosts, name="AdvectionDir_" + str(d))
+        ## Transport velocity
+        self.velocity = velocity
+        ## direction to advect
+        self.dir = d
+        ## Resolutions of the fields. Dictionnary with
+        ## fields as keys.
+        ## Example :
+        ## resolutions[velocity] may return [32,32,32].
+        self.resolutions = resolutions
+        ## Extra parameters (depend on the method)
+        self.config = other_config
+        self.input = [self.velocity]
+        self.particle_fields = particle_fields
+        self.particle_positions = particle_positions
+
+    @debug
+    def setUp(self):
+        """
+        Discretisation according to the chosen method.
+        Available methods :
+        - 'gpu' : OpenCL kernels (2d and 3d, single field, scalar or vector)
+          - '1k' : Single OpenCL kernel for advection and remeshing
+          - '2k' : Separate kernels
+          - 'm4prime' : M4prime formula
+          - 'm6prime' : M6prime formula
+          - 'm8prime' : M8prime formula
+          - 'l6star' : Lambda6star formula
+        - other : Pure python (2d and 3d, list of vector and/or scalar)
+          - 'rk2' : Runge Kutta 2nd order advection
+          - 'm4prime' : M4prime formula
+          - 'm6prime' : M6prime formula
+          - 'm8prime' : M8prime formula
+          - 'l6star' : Lambda6star formula
+        """
+        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:
+                    topo = self.domain.getOrCreateTopology(
+                        self.domain.dimension, self.resolutions[v])
+                    self.discreteFields[v] = v.discretize(topo)
+        else:
+            if self._predefinedTopo is not None:
+                raise ValueError("User-defined topology is not\
+                (yet) allowed for advection.")
+
+            topodims = np.ones((self.domain.dimension))
+            if self.dir == self.domain.dimension - 1:
+                topodims[0] = main_size
+            else:
+                topodims[-1] = main_size
+
+            for v in self.variables:
+                topo = self.domain.getOrCreateTopology(
+                    self.domain.dimension, self.resolutions[v],
+                    topoResolution=topodims, fixedResolution=True)
+                self.discreteFields[v] = v.discretize(topo)
+
+        advectedDiscreteFields = [self.discreteFields[v]
+                                  for v in self.variables
+                                  if not v is self.velocity]
+
+        particles_advectedDiscreteFields = [
+            v.discretize(advectedDiscreteFields[0].topology)
+            for v in self.particle_fields]
+        particles_positionsDiscreteField = None
+        if self.particle_positions is not None:
+            particles_positionsDiscreteField = \
+                self.particle_positions.discretize(
+                advectedDiscreteFields[0].topology)
+
+        if self.method.find('gpu_2k') >= 0:
+            from parmepy.gpu.gpu_particle_advection_2k \
+                import GPUParticleAdvection2k as ParticleAdvection
+        elif self.method.find('gpu_1k') >= 0:
+            from parmepy.gpu.gpu_particle_advection_1k \
+                import GPUParticleAdvection1k as ParticleAdvection
+        else:
+            from parmepy.operator.discrete.particle_advection \
+                import ParticleAdvection
+
+        self.discreteOperator = ParticleAdvection(
+            self.discreteFields[self.velocity],
+            advectedDiscreteFields, self.dir, method=self.method,
+            part_position=particles_positionsDiscreteField,
+            part_advectedFields=particles_advectedDiscreteFields,
+            **self.config)
+
+        # -- Final set up --
+        self.discreteOperator.setUp()
+        self._isUpToDate = True
+
+    @debug
+    def apply(self, simulation, dtCoeff, split_id, old_dir=None):
+        self.discreteOperator.apply(
+            simulation, dtCoeff, split_id, old_dir)
+
+    @debug
+    def finalize(self):
+        """
+        Memory cleaning.
+        """
+        self.timer = Timer(self)
+        self.discreteOperator.finalize()
+        self.timer = self.timer + self.discreteOperator.timer
+        print self.timer
diff --git a/HySoP/hysop/operator/monitors/printer.py b/HySoP/hysop/operator/monitors/printer.py
index e9c121495ba240ffde58c9258c0c54317c335c51..c23a7f3be57207b9123b9b36fff2ae7e458c258e 100644
--- a/HySoP/hysop/operator/monitors/printer.py
+++ b/HySoP/hysop/operator/monitors/printer.py
@@ -153,7 +153,7 @@ class Printer(Monitoring):
                                             df[1][i, j]))
                                 else:
                                     f.write("{0:8.12} ".format(
-                                            df[i, j]))
+                                            df[0][i, j]))
                             f.write("\n")
             elif self.variables[0].domain.dimension == 3:
                 df = self.variables[0].discreteFields.values()[0]
@@ -178,7 +178,7 @@ class Printer(Monitoring):
                                             df[2][i, j, k]))
                                 else:
                                     f.write("{0:8.12} ".format(
-                                            df[i, j, k]))
+                                            df[0][i, j, k]))
                             f.write("\n")
             else:
                 df = self.variables[0].discreteFields.values()[0]
diff --git a/HySoP/hysop/operator/tests/test_advec_scales.py b/HySoP/hysop/operator/tests/test_advec_scales.py
index e2a58695815dfd4079ed48955667d82020871311..72dc29bd49dfa1b3babce351d59e87aab507a86a 100755
--- a/HySoP/hysop/operator/tests/test_advec_scales.py
+++ b/HySoP/hysop/operator/tests/test_advec_scales.py
@@ -9,7 +9,7 @@ from parmepy.operator.advection import Advection
 from parmepy.problem.simulation import Simulation
 
 
-def test_nullVelocity_m4():
+def _nullVelocity_m4():
     """Basic test with random velocity. Using M4prime"""
     box = Box(3, length=[1., 1., 1.], origin=[0., 0., 0.])
     scal = Field(domain=box, name='Scalar')
@@ -44,7 +44,7 @@ def test_nullVelocity_m4():
     assert np.allclose(scal_ref_d.data[0], scal_d.data[0])
 
 
-def test_nullVelocity_vec_m4():
+def _nullVelocity_vec_m4():
     """Basic test with random velocity and vector field. Using M4prime"""
     box = Box(3, length=[1., 1., 1.], origin=[0., 0., 0.])
     scal = Field(domain=box, name='Scalar', isVector=True)
diff --git a/HySoP/hysop/problem/tests/test_transport.py b/HySoP/hysop/problem/tests/test_transport.py
index 5b5eb07dcba5561640cfa9fdba78b7b6de77ab7e..3c390693f51f82604085f29af3eb10458c39c9b7 100644
--- a/HySoP/hysop/problem/tests/test_transport.py
+++ b/HySoP/hysop/problem/tests/test_transport.py
@@ -46,7 +46,7 @@ def rotating_velocity_2D(x, y):
     return -c * y, c * x
 
 
-def assert_nullVelocity(dim, boxLength, boxMin, nbElem, finalTime, timeStep,
+def assertion(dim, boxLength, boxMin, nbElem, finalTime, timeStep,
                         s, v, rtol=1e-05, atol=1e-08):
     box = Box(dim, length=boxLength, origin=boxMin)
     scal = Field(domain=box, formula=s, name='Scalar')
@@ -55,7 +55,7 @@ def assert_nullVelocity(dim, boxLength, boxMin, nbElem, finalTime, timeStep,
                       resolutions={velo: nbElem,
                                    scal: nbElem},
                       )
-    simu = Simulation(tend=finalTime, iterMax=1)
+    simu = Simulation(tinit=0.0, tend=finalTime, timeStep=timeStep, iterMax=1)
     pb = TransportProblem([advec], simu)
     pb.setUp()
     initial_scalar = scal.discreteFields.values()[0].data[0].copy()
@@ -72,7 +72,7 @@ def test_nullVelocity_2D():
     nbElem = [nb, nb]
     timeStep = 0.01
     finalTime = timeStep
-    assert assert_nullVelocity(dim, boxLength, boxMin,
+    assert assertion(dim, boxLength, boxMin,
                                nbElem, finalTime, timeStep,
                                lambda x, y: np.random.random(),
                                lambda x, y: (0., 0.))
@@ -86,7 +86,7 @@ def test_nullVelocity_3D():
     nbElem = [nb, nb, nb]
     timeStep = 0.01
     finalTime = timeStep
-    assert assert_nullVelocity(dim, boxLength, boxMin,
+    assert assertion(dim, boxLength, boxMin,
                                nbElem, finalTime, timeStep,
                                lambda x, y, z: np.random.random(),
                                lambda x, y, z: (0., 0., 0.))
@@ -100,7 +100,7 @@ def test_gaussian_2D():
     nbElem = [nb, nb]
     timeStep = 0.01
     finalTime = timeStep
-    assert assert_nullVelocity(dim, boxLength, boxMin,
+    assert assertion(dim, boxLength, boxMin,
                                nbElem, finalTime, timeStep,
                                gaussian_scalar_2D, rotating_velocity_2D,
                                rtol=1e-04, atol=1e-05)
@@ -114,7 +114,7 @@ def test_cosinus_translation_2D():
     nbElem = [nb, nb]
     timeStep = 1.
     finalTime = 1.
-    assert assert_nullVelocity(dim, boxLength, boxMin,
+    assert assertion(dim, boxLength, boxMin,
                                nbElem, finalTime, timeStep,
                                cosinus_product_2D,
                                lambda x, y: (1., 2.))
@@ -128,7 +128,7 @@ def test_cosinus_translation_3D():
     nbElem = [nb, nb, nb]
     timeStep = 1.
     finalTime = timeStep
-    assert assert_nullVelocity(dim, boxLength, boxMin,
+    assert assertion(dim, boxLength, boxMin,
                                nbElem, finalTime, timeStep,
                                cosinus_product_3D,
                                lambda x, y, z: (1., 2., 0.5))