Update stretching (two formulations) and integrators

HySoP/CMakeLists.txt

@@ -180,7 +180,7 @@ else()
 endif()
 
 # Target to remove parmepy from install-path.
-add_custom_target(python-cleaninstall COMMAND rm -rf ${${PROJECT_NAME}_INSTALL_DIR}/${PROJECT_NAME}*
+add_custom_target(python-cleaninstall COMMAND rm -rf ${CMAKE_INSTALL_PREFIX}/${PROJECT_NAME}*
   COMMENT "remove parmepy package and its dependencies")
 
 # Target to clean sources (remove .pyc files) and build dir.
@@ -246,9 +246,17 @@ if(VERBOSE_MODE)
   message(STATUS " Project uses Scales : ${WITH_SCALES}")
   message(STATUS " Project uses FFTW : ${WITH_FFTW}")
   message(STATUS " Project uses GPU : ${WITH_GPU}")
-  message(STATUS " Project will be built in ${${PROJECT_NAME}_PYTHON_BUILD_DIR}.")
+  message(STATUS " Project will be built in build/${${PROJECT_NAME}_PYTHON_BUILD_DIR}.")
   message(STATUS " Python packages will be installed in ${CMAKE_INSTALL_PREFIX}")
   message(STATUS " ${PROJECT_NAME} debug mode : ${DEBUG}")
   message(STATUS "====================== ======= ======================")
+  message(STATUS " ")
+  message(STATUS "Try :")
+  message(STATUS " 'make -jN' to build the project, N being the number of available processes.")
+  message(STATUS " 'make python-install' to install python modules and their dependencies. ")
+  message(STATUS " 'make doc' to generate doxygen documentation for parmepy.")
+  message(STATUS " 'make test' to run some test (after the build!).")
+  message(STATUS " 'make clean' to clean build directory.")
+  message(STATUS " 'make python-cleaninstall' to clean install directory.")
 endif()
 

HySoP/hysop/constants.py

@@ -31,6 +31,10 @@ ZDIR = 2
 PERIODIC = 0
 ## Directions string
 S_DIR = ["_X", "_Y", "_Z"]
+## Stretching formulation (div(w:u))
+CONSERVATIVE = 1
+## Stretching formulation ([grad(u)][w])
+GRADUW = 0
 
 
 #define debug decorator:

HySoP/hysop/numerics/differential_operations.py

-# -*- codingind utf-8 -*-
+# -*- coding: utf-8 -*-
 """
 @file differential_operations.py
 Library of functions used to perform classical vector calculus
@@ -7,6 +7,7 @@ Library of functions used to perform classical vector calculus
 from parmepy.constants import debug, XDIR, YDIR, ZDIR
 from abc import ABCMeta, abstractmethod
 from finite_differences import FD_C_4
+import numpy as np
 
 
 class DifferentialOperation(object):
@@ -18,15 +19,17 @@ class DifferentialOperation(object):
 
     @debug
     @abstractmethod
-    def __init__(self):
+    def __init__(self, topo, method=None):
         pass
-
+    
 
 class Curl(DifferentialOperation):
     """
     Computes \f$ \nabla \ times (f) \f$, f being a vector field.
     """
     def __init__(self, topo, method=None):
+
+        DifferentialOperation.__init__(self, topo)
         if method is not None and method.find('FD_order2') >= 0:
                 raise ValueError("2nd order scheme Not yet implemented")
         else:
@@ -36,19 +39,21 @@ class Curl(DifferentialOperation):
             self.fd_scheme = FD_C_4(topo.mesh.space_step)
             self.indices = topo.mesh.iCompute
             self.fd_scheme.computeIndices(self.indices)
+        self.nbComponents = topo.domain.dimension
 
     def __call__(self, variable, data=None):
         return self.fcall(variable, data)
 
     def FDCentral4(self, variable, data=None):
         """
+        @param variable : a list of numpy arrays
         """
         # Get slices for subarrays ...
        # init (if required) return data
         if data is None:
             data = []
-            for i in xrange(variable.dimension):
-                data.append(variable.data[i].copy())
+            for i in xrange(self.nbComponents):
+                data.append(variable[i].copy())
                 #data[i][...] = 0.0
 
         #data[0][...] = 0.0
@@ -100,6 +105,7 @@ class DivV(DifferentialOperation):
         self.fd_scheme.computeIndices(self.indices)
         assert (topo.ghosts >= 2).all(),\
             'you need a ghost layer for FD4 scheme.'
+        self.nbComponents = topo.domain.dimension
 
     def __call__(self, var1, scal, data=None, work=None):
         return self.fcall(var1, scal, data)
@@ -109,9 +115,9 @@ class DivV(DifferentialOperation):
         Compute central finite difference scheme, order 4.
         """
         if data is None:
-            data = var1.data[0].copy()
+            data = np.zeros_like(var1[0])
         if work is None:
-            work = var1.data[0].copy()
+            work = np.zeros_like(var1[0])
 
         #data[...] = 0.0
 
@@ -129,9 +135,9 @@ class DivV(DifferentialOperation):
         """
         """
         if data is None:
-            data = var1.data[0].copy()
-
-        data[...] = 0.0
+            data = np.zeros_like(var1[0])
+        else:
+            data[...] = 0.0
 
         #  d/dx (scal * v1x)
         self.fd_scheme.accumulate(var1[XDIR] * scal, XDIR, data)
@@ -154,6 +160,7 @@ class DivT(DifferentialOperation):
             self.fcall = self.FDCentral4
             self.div = DivV(topo, method='FD_C4')
             self.divZ = DivV(topo, method='FD_C4_nowork')
+        self.nbComponents = topo.domain.dimension
 
     def __call__(self, var1, var2, data=None):
         return self.fcall(var1, var2, data)
@@ -164,8 +171,9 @@ class DivT(DifferentialOperation):
         # init (if required) return data
         if data is None:
             data = []
-            for i in xrange(var1.dimension):
-                data.append(var1.data[i].copy())
+            for i in xrange(self.nbComponents):
+                data.append(np.zeros_like(var1[i]))
+
         # we use data[Y] and data[Z] as work ...
         data[XDIR] = self.div(var1, var2[XDIR],
                               data=data[XDIR], work=data[YDIR])
@@ -174,3 +182,124 @@ class DivT(DifferentialOperation):
         # No work vector for zdir
         data[ZDIR] = self.divZ(var1, var2[ZDIR], data=data[ZDIR])
         return data
+
+
+class GradS(DifferentialOperation):
+    """
+    Computes \f$ \nabla(\rho) \f$, \f$ \rho\f$ a scalar field
+    """
+    def __init__(self, topo, method=None):
+        if method is not None and method.find('FD_order2') >= 0:
+                raise ValueError("2nd order scheme Not yet implemented")
+        elif method.find('FD_C4') >= 0:
+            self.fcall = self.FDCentral4
+            self.fd_scheme = FD_C_4((topo.mesh.space_step))
+        else:
+            raise ValueError("FD scheme Not yet implemented")
+
+        self.indices = topo.mesh.iCompute
+        self.fd_scheme.computeIndices(self.indices)
+        assert (topo.ghosts >= 2).all(),\
+            'you need a ghost layer for FD4 scheme.'
+        ## number of components for vector fields
+        self.nbComponents = topo.domain.dimension
+
+    def __call__(self, scal, data=None, work=None):
+        return self.fcall(scal, data)
+
+    def FDCentral4(self, scal, data=None):
+        """
+        Compute central finite difference scheme, order 4.
+        @param scal a numpy array
+        """
+        if data is None:
+            data = []
+            for i in xrange(self.nbComponents):
+                data.append(np.zeros_like(scal))
+        #data[...] = 0.0
+
+        #  d/dx (scal), saved in data
+        self.fd_scheme.inplace(scal, XDIR, data[XDIR])
+        #  d/dy (scal), saved in data
+        self.fd_scheme.inplace(scal, YDIR, data[YDIR])
+        #  d/dz(scal), saved in data
+        self.fd_scheme.inplace(scal, ZDIR, data[ZDIR])
+        return data
+
+
+class GradVxW(DifferentialOperation):
+    """
+    Computes \f$ [\nabla(V)][W] \f$, with
+    \f$V\f$ and \f$W\f$ some vector fields.
+    """
+    def __init__(self, topo, method=None):
+        if method is not None and method.find('FD_order2') >= 0:
+                raise ValueError("2nd order scheme Not yet implemented")
+        else:
+            self.fcall = self.FDCentral4
+            self.grad = GradS(topo, method='FD_C4')
+
+            self.divZ = DivV(topo, method='FD_C4_nowork')
+        self.nbComponents = topo.domain.dimension
+
+        if self.nbComponents == 1:
+            self.worklength = 1
+        elif self.nbComponents == 2:
+            self.worklength = 3
+        else:
+            self.worklength = 5
+        self.prod = BlockProd()
+
+    def __call__(self, var1, var2, diagnostics=None, data=None):
+        return self.fcall(var1, var2, diagnostics, data)
+
+    def FDCentral4(self, var1, var2, diagnostics=None, data=None):
+        """
+        """
+        # init (if required) return data
+        if data is None:
+            data = []
+            for i in xrange(self.worklength):
+                data.append(np.zeros_like(var1[0]))
+        else:
+            assert len(data) == self.worklength
+        if diagnostics is not None:
+            diagnostics[:] = 0.0
+
+        for direc in xrange(self.nbComponents):
+            data[direc:self.nbComponents + direc] = \
+                self.grad(var1[direc],
+                          data=data[direc:self.nbComponents + direc])
+
+            if diagnostics is not None:
+                diagnostics[0] = max(diagnostics[0],
+                                     sum([abs(data[i])
+                                          for i in xrange(direc,
+                                                          self.nbComponents
+                                                          + direc)]))
+
+                # compute grad(Vx).var2
+                data[direc] = self.prod(data[direc:self.nbComponents + direc],
+                                        var2, data[direc])
+        return data
+
+
+class BlockProd(object):
+    """
+    Block product of two lists of numpy arrays
+    """
+    def __init__(self):
+        pass
+
+    def __call__(self, list1, list2, data=None):
+        assert len(list1) == len(list2)
+
+        if data is not None:
+            assert data.shape == list1[0].shape
+            data[...] = 0.0
+        else:
+            data = np.zeros_like(list1[0])
+
+        for i in xrange(len(list1)):
+            data += list1[i] * list2[i]
+        return data

HySoP/hysop/numerics/finite_differences.py

@@ -29,10 +29,11 @@ class FD_C_4(FiniteDifference):
     def __init__(self, step):
         """
         Constructor
-        @param step : grid step size
+        @param step : a numpy array with grid step size in each dir.
         """
         #  dim of the field on which scheme will be applied
         # (i.e dim of the domain)
+        step = np.asarray(step)
         self._dim = step.size
         # list of indices for index - 1
         self._m1 = []
@@ -73,7 +74,7 @@ class FD_C_4(FiniteDifference):
 
     def inplace(self, tab, cdir, result):
         """
-        Compute \f$ \frac{\partial d tab()}{\partial cdir} \f$ using fd scheme,
+        Compute \f$ \frac{\partial tab()}{\partial cdir} \f$ using fd scheme,
         The result is saved in input/ouput parameter result (in place!)
         @param tab : a numpy array
         @param cdir : direction of differentiation
@@ -88,7 +89,7 @@ class FD_C_4(FiniteDifference):
 
     def accumulate(self, tab, cdir, result):
         """
-        Compute \f$ \frac{\partial d tab()}{\partial cdir} \f$ using fd scheme,
+        Compute \f$ \frac{\partial tab()}{\partial cdir} \f$ using fd scheme,
         The result is ADDED to input/ouput parameter result.
         @param tab : a numpy array
         @param cdir : direction of differentiation

HySoP/hysop/numerics/integrators/euler.py

@@ -4,7 +4,6 @@
 
 Forward Euler method.
 """
-from parmepy.constants import np, PARMES_REAL, ORDER
 from parmepy.numerics.integrators.odesolver import ODESolver
 
 
@@ -17,35 +16,18 @@ class Euler(ODESolver):
     y(t_{n+1}) = y(t_n) + dt*f(y(t_n))
 
     """
-    def __init__(self, f=None, dim=1):
+    def __init__(self, dim, f=None):
         """
-        Constructor.
-
         @param f function f(t,y) : Right hand side of the equation to solve.
         @param dim : dimensions
-        @param conditions : function to apply boundary conditions.
         """
-        ODESolver.__init__(self, f)
+        ODESolver.__init__(self, dim, f)
         self.name = 'euler'
 
-    def __call__(self, f, t, dt, y):
-        """
-        Integration step for Euler Method.
-         up = f[t, u(t,x,y,z)].
-
-        @param u : function of position at time t.
-        @param f : function of function of position at time t.
-        @param t : current time.
-        @param dt : time step.
-#        """
-#        result = y + dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-#        return result
-
-        RHS = f(t, y)
-
-        return y[0] + dt * RHS[0], \
-               y[1] + dt * RHS[1], \
-               y[2] + dt * RHS[2]
+    def __call__(self, y, dt, t):
+
+        result = self.f(t, y)
+        return [y[i] + dt * result[i] for i in xrange(self._dim)]
 
     def __str__(self):
         """ToString method"""

HySoP/hysop/numerics/integrators/odesolver.py

@@ -12,46 +12,40 @@ class ODESolver(NumMethod):
     """
     Abstract description for ODE solvers.
     Solve the system :
-    y'(t) = f(t,y)
+    \f[ \frac{dy(t)}{dt} = f(t,y) \f]
 
-    By default, f = 0.
     """
 
     __metaclass__ = ABCMeta
 
     @abstractmethod
-    def __init__(self, f=None):
+    def __init__(self, dim, f=None):
         """
-        Constructor.
-
-        @param f function f(t,y) : Right hand side of the equation to solve.
+        @param dim, vector field (f) dimension (number of component)
+        @param f : function f(t,y), right hand side of the equation to solve.
+        WARNING : y arg must be a list of numpy arrays and f must return
+        a list of numpy arrays.
         """
         super(ODESolver, self).__init__()
         ## RHS.
         if f is not None:
             self.f = f
         else:
-            #[0. for x in y]
-            self.f = lambda t, y: np.vectorize(f)(t, y)
+            self.f = lambda t, y: [np.zeros_like(y[i]) for i in xrange(dim)]
+
+        self._dim = dim
 
     @abstractmethod
-    def __call__(self, y, t, dt):
+    def __call__(self, y, dt, t):
         """
-        Abstract method, apply operaton on a variable.
-        Must be implemented by sub-class.
-
         @param y : position at time t.
         @param t : current time.
         @param dt : time step.
         @return y at t+dt.
         """
-        pass
 
 if __name__ == "__main__":
     print __doc__
     print "- Provided class : ODESolver (abstract)"
     print ODESolver.__doc__
 
-
-
-

HySoP/hysop/numerics/integrators/runge_kutta2.py

@@ -3,7 +3,6 @@
 
 RK2 method interface.
 """
-from parmepy.constants import np
 from parmepy.numerics.integrators.odesolver import ODESolver
 
 
@@ -15,98 +14,79 @@ class RK2(ODESolver):
     y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
 
     """
-    def __init__(self, f=None, dim=3):
+    def __init__(self, dim, f=None):
         """
-        Constructor.
-
+        @param dim, vector field (f) dimension (number of component)
         @param f function f(t,y) : Right hand side of the equation to solve.
-        @param dim : dimensions
         """
-        ODESolver.__init__(self, f)
+        ODESolver.__init__(self, dim, f)
         self.name = 'RK2'
 
-    def __call__(self, f, t, dt, y):
-
-        """
-        Integration step for RK2 Method.
-         up = f[t, Y(t,space)].
-
-        @param f : right-hand side.
-        @param t : current time.
-        @param dt : time step.
-        @param y : function Y(t,space)
-        """
-
-#        # K1 = dt *f[t,y(t,x,y,z)]
-#        K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-#        # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
-
-#        result = y + dt * np.asarray(f(t + dt / 2., y + K1 / 2.),
-#                                     dtype=PARMES_REAL, order=ORDER)
-#        return result
-
-        RHS1 = f(t, y)
-        K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
-
-        RHS2 = f(t + dt / 2., y + K1 / 2.)
-
-        return y[0] + dt * RHS2[0], \
-               y[1] + dt * RHS2[1], \
-               y[2] + dt * RHS2[2]
-
-
-class RK2_scalar(ODESolver):
-    """
-    ODESolver implementation for solving an equation system with RK2 method
-    for scalar quantities.
-    y'(t) = f(t,y)
-
-    y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
-
-    """
-    def __init__(self, f=None):
-        """
-        Constructor.
-
-        @param f function f(t,y) : Right hand side of the equation to solve.
-        @param dim : dimensions
-        """
-        ODESolver.__init__(self, f)
-        self.name = 'RK2'
-
-    def __call__(self, t, dt, y, res=None,
-                 work=None,  y_temp=None,
-                 yprime=None, **f_kwargs):
-
-        """
-        Integration step for RK2 Method.
-         up = f[t, Y(t,space)].
-
-        @param t : current time.
-        @param dt : time step.
-        @param y : function Y(t,space).
-        @param res : result (may be allocated when not None).
-        @param work : allocated working array.
-        @param yprime : y'(t) already computed (avoid 1st call to f)
-        @param f_kwargs : Parameters for right-hand side computing.
-        """
-        if res is None:
-            res = np.empty_like(y)
-        if work is None:
-            work = np.empty_like(y)
-        if yprime is None:
-            yprime = np.empty_like(y)
-            yprime[...] = self.f(t, y, res=yprime, **f_kwargs)
-
-        y_temp[...] = y + 0.5 * dt * yprime
-        work[...] = self.f(t + 0.5 * dt, y_temp, res=work, **f_kwargs)
-
-        res[...] = y
-        res[...] += dt * work
-
-    def __str__(self):
-        """ToString method"""
-        return "RK2 Method"
+    def __call__(self, y, dt, t):
+        # k1 = f(t,y)
+        # k2 = f(t + dt/2, y + dt/2 * k1)
+        # result = y + dt * k2
+
+        assert len(y) == self._dim
+        # k1 ... in k2
+        k2 = self.f(t, y)
+        # k2
+        k2 = [y[i] + dt / 2 * k2[i] for i in xrange(self._dim)]
+        k2 = self.f(t + dt / 2., k2)
+
+        return [y[i] + dt * k2[i] for i in xrange(self._dim)]
+
+
+## class RK2_scalar(ODESolver):
+##     """
+##     ODESolver implementation for solving an equation system with RK2 method
+##     for scalar quantities.
+##     y'(t) = f(t,y)
+
+##     y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
+
+##     """
+##     def __init__(self, f=None):
+##         """
+##         Constructor.
+
+##         @param f function f(t,y) : Right hand side of the equation to solve.
+##         @param dim : dimensions
+##         """
+##         ODESolver.__init__(self, f)
+##         self.name = 'RK2'
+
+##     def __call__(self, f, t, dt, y, res=None, work=None, yprime=None,
+##                  **f_kwargs):
+
+##         """
+##         Integration step for RK2 Method.
+##          up = f[t, Y(t,space)].
+
+##         @param f : right-hand side.
+##         @param t : current time.
+##         @param dt : time step.
+##         @param y : function Y(t,space).
+##         @param res : result (may be allocated when not None).
+##         @param work : allocated working array.
+##         @param yprime : y'(t) already computed (avoid 1st call to f)
+##         @param f_kwargs : Parameters for right-hand side computing.
+##         """
+##         if res is None:
+##             res = np.empty_like(y)
+##         if work is None:
+##             work = np.empty_like(y)
+##         if yprime is None:
+##             work[...] = f(t, y, res=work, **f_kwargs)
+##         else:
+##             work[...] = yprime
+##         work[...] *= 0.5 * dt
+##         res[...] = y
+##         res[...] += dt * f(t + 0.5 * dt, y + work, **f_kwargs)
+
+##     def __str__(self):
+##         """ToString method"""
+##         return "RK2 Method"
 
 
 if __name__ == "__main__":

HySoP/hysop/numerics/integrators/runge_kutta3.py

@@ -3,64 +3,41 @@
 
 RK3 method interface.
 """
-from parmepy.constants import np, PARMES_REAL, ORDER
 from parmepy.numerics.integrators.odesolver import ODESolver
 
 
 class RK3(ODESolver):
     """
     ODESolver implementation for solving an equation system with RK3 method.
-    y'(t) = f(t,y)
-
-    y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
-
     """
-    def __init__(self, f=None, dim=3):
+    def __init__(self, dim, f=None):
         """
-        Constructor.
-
         @param f function f(t,y) : Right hand side of the equation to solve.
         @param dim : dimensions
-        @param conditions : function to apply boundary conditions.
         """
-        ODESolver.__init__(self, f)
+        ODESolver.__init__(self, dim, f)
         self.name = 'RK3'
 
-    def __call__(self, f, t, dt, y):
-        """
-        Integration step for RK2 Method.
-         up = f[t, Y(t,space)].
-
-        @param f : right-hand side.
-        @param t : current time.
-        @param dt : time step.
-        @param y : function Y(t,space)
-        """
-#        # K1 = dt *f[t,y(t,x,y,z)]
-#        K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-#        # 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)
-#        # result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
-#        result = y + 1. / 4. * (K1 + 3. * K3)
-#        return result
+    def __call__(self, y, dt, t):
 
-        RHS1 = f(t, y)
-        K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
+        # k1 = f(t,y)
+        # k2 = f(t+dt/3, y + dt/3 *k1))
+        # k3 = f(t + 2/3 *dt , y + 2/3 dt * k2))
+        # result = y + 0.25 * dt * (k1 + 3 * k3)
 
-        RHS2 = f(t + dt / 3., y + K1 / 3.)
-        K2 = np.array([dt * RHS2[0], dt * RHS2[1], dt * RHS2[2]])
+        assert len(y) == self._dim
 
-        RHS3 = f(t + dt * 2. / 3., y + K2 * 2. / 3.)
-        K3 = np.array([dt * RHS3[0], dt * RHS3[1], dt * RHS3[2]])
+        k1 = self.f(t, y)
 
+        # k2 ... in k3
+        k3 = [y[i] + dt / 3 * k1[i] for i in xrange(self._dim)]
+        k3 = self.f(t + dt / 3., k3)
+        # k3 ...
+        k3 = [y[i] + 2 * dt / 3 * k3[i] for i in xrange(self._dim)]
+        k3 = self.f(t + dt * 2. / 3., k3)
         # 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])
+        return [y[i] + 0.25 * dt * (k1[i] + 3 * k3[i])
+                for i in xrange(self._dim)]
 
     def __str__(self):
         """ToString method"""

HySoP/hysop/numerics/integrators/runge_kutta4.py

@@ -4,75 +4,59 @@
 
 RK4 method interface.
 """
-from parmepy.constants import np, PARMES_REAL, ORDER
 from parmepy.numerics.integrators.odesolver import ODESolver
 
 
 class RK4(ODESolver):
     """
     ODESolver implementation for solving an equation system with RK4 method.
-    y'(t) = f(t,y)
-
-    y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
-
     """
-    def __init__(self, f=None, dim=3):
+    def __init__(self, dim, f=None):
         """
-        Constructor.
-
         @param f function f(t,y) : Right hand side of the equation to solve.
         @param dim : dimensions
-        @param conditions : function to apply boundary conditions.
         """
-        ODESolver.__init__(self, f)
+        ODESolver.__init__(self, dim, f)
         self.name = 'RK4'
 
-    def __call__(self, f, t, dt, y):
-        """
-        Integration step for RK4 Method.
-         up = f[t, Y(t,space)].
-
-        @param f : right-hand side.
-        @param t : current time.
-        @param dt : time step.
-        @param y : function Y(t,space)
-        """
-
-#        # K1 = dt *f[t, y(t, x, y, z)]
-#        K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-#        K2 = dt * np.asarray(f(t + dt / 2., y + K1[:, :, :, :] / 2.),
-#                             dtype=PARMES_REAL, order=ORDER)
-#        K3 = dt * np.asarray(f(t + dt / 2., y + K2[:, :, :, :] / 2.),
-#                             dtype=PARMES_REAL, order=ORDER)
-#        K4 = dt * np.asarray(f(t + dt, y + K3[:, :, :, :]),
-#                             dtype=PARMES_REAL, order=ORDER)
-#        # result = y(t, x, y, z)]  +  dt *f[t + dt/2 ,  y(t, x, y, z) + k1/2.]
-#        result = y + 1. / 6. * (K1[:, :, :, :] + 2. * K2[:, :, :, :] +
-#                                2. * K3[:, :, :, :] + K4[:, :, :, :])
-#        return result
-
-        RHS1 = f(t, y)
-        K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
-
-        RHS2 = f(t + dt / 2., y + K1 / 1.)
-        K2 =  np.array([dt * RHS2[0], dt * RHS2[1], dt * RHS2[2]])
-
-        RHS3 = f(t + dt * 2., y + K2 * 2.)
-        K3 = np.array([dt * RHS3[0], dt * RHS3[1], dt * RHS3[2]])
-
-        RHS4 = f(t + dt, y + K3)
-        K4 = np.array([dt * RHS4[0], dt * RHS4[1], dt * RHS4[2]])
-
-        # result = y(t, x, y, z)]  +  dt *f[t + dt/2 ,  y(t, x, y, z) + k1/2.]
-        return y[0] + 1. / 6. * (K1[0] + 2. * K2[0] + 2. * K3[0] + K4[0]), \
-               y[1] + 1. / 6. * (K1[1] + 2. * K2[1] + 2. * K3[1] + K4[1]), \
-               y[2] + 1. / 6. * (K1[2] + 2. * K2[2] + 2. * K3[2] + K4[2])
+    def __call__(self, y, dt, t):
+
+        # k1 = f(t,y)
+        # k2 = f(t + dt/2, y + dt/2 * k1)
+        # k3 = f(t + dt/2, y + dt/2 * k2)
+        # k4 = f(t + dt, y + dt * k3)
+        # result = y + dt/6( k1 + 2 * k2 + 2 * k3 + k4)
+
+        # k1 in result
+        result = self.f(t, y)
+        # k2 in work
+        work = [y[i] + dt / 2 * result[i] for i in xrange(self._dim)]
+        work = self.f(t + dt / 2, work)
+
+        # k1 + 2 * k2
+        for i in xrange(self._dim):
+            result[i] += 2 * work[i]
+        # k3 in work
+        work = [y[i] + dt / 2 * work[i] for i in xrange(self._dim)]
+        work = self.f(t + dt / 2, work)
+        # k1 + 2 * k2 + 2 * k3
+        for i in xrange(self._dim):
+            result[i] += 2 * work[i]
+        # k4 in work
+        work = [y[i] + dt * work[i] for i in xrange(self._dim)]
+        work = self.f(t + dt, work)
+        # k1 + 2 * k2 + 2 * k3 + k4
+        for i in xrange(self._dim):
+            result[i] += work[i]
+            result[i] *= dt / 6
+            result[i] += y[i]
+
+        return result
 
     def __str__(self):
         """ToString method"""
         return "RK4 Method"
 
-
 class RK4_scalar(ODESolver):
     """
     ODESolver implementation for solving an equation system with RK4 method
@@ -132,7 +116,6 @@ class RK4_scalar(ODESolver):
         """ToString method"""
         return "RK4 Method"
 
-
 if __name__ == "__main__":
     print __doc__
     print "- Provided class : RK4"

HySoP/hysop/operator/continuous.py

@@ -75,7 +75,6 @@ class Operator(object):
         ## Number of points in the ghost layer
         if ghosts is not None:
             self.ghosts = np.asarray(ghosts)
-            print "ooioii", self.ghosts, self.ghosts.size, self.domain.dimension
             assert self.ghosts.size == self.domain.dimension
         else:
             self.ghosts = None

HySoP/hysop/operator/curl.py

@@ -6,6 +6,7 @@ Curl of a field.
 from parmepy.constants import debug, np
 from parmepy.operator.continuous import Operator
 from parmepy.operator.discrete.curl import CurlFD
+from parmepy.operator.updateGhosts import UpdateGhosts
 
 
 class Curl(Operator):
@@ -21,7 +22,7 @@ class Curl(Operator):
         @param invar : the input vector field
         @param outvar : curl of input vector field
         @param resolution : resolution (global) of the fields
-        @param method : Method used (default = finite difference, 4th order)
+        @param method : method used (default = finite difference, 4th order)
         @param topo : a predefined topology to discretize invar/outvar
         """
 
@@ -35,15 +36,15 @@ class Curl(Operator):
         self.invar = invar
         self.outvar = outvar
         if self.method is None:
-            self.method = 'FD_order4'
+            self.method = 'FD_C4'
+        ## An operator to update ghost points of input var.
+        self.ghostsOperator = None
 
     @debug
     def setUp(self):
-
-        print "set up for curl ..."
         if self.method.find('FD_order2') >= 0:
             nbGhosts = 1
-        elif self.method.find('FD_order4') >= 0:
+        elif self.method.find('FD_C4') >= 0:
             nbGhosts = 2
         else:
             nbGhosts = 0
@@ -75,8 +76,18 @@ class Curl(Operator):
                                        self.discreteFields[self.outvar],
                                        method=self.method)
         self.discreteOperator.setUp()
+        # ghost op
+        self.ghostsOperator = UpdateGhosts(self.input, topo=topo)
+
         self._isUpToDate = True
 
+    def apply(self, simulation=None):
+        # update ghost points, if required
+        if self.ghostsOperator is not None:
+            self.ghostsOperator.apply()
+        # computation ...
+        self.discreteOperator.apply(simulation)
+
 if __name__ == "__main__":
     print __doc__
     print "- Provided class : Curl"

HySoP/hysop/operator/discrete/stretching.py

@@ -2,50 +2,64 @@
 """
 @file operator/discrete/stretching.py
 
-Discrete stretching representation
+Discretisation of the stretching operator for two different formulations.
 """
 
 from parmepy.constants import debug
 from discrete import DiscreteOperator
-from parmepy.numerics.differentialOperator import DifferentialOperator
-from parmepy.numerics.fct2op import Fct2Op
 from parmepy.numerics.integrators.euler import Euler
 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
+from parmepy.numerics.differential_operations import DivT, GradVxW
+import numpy as np
+import math
+ceil = math.ceil
 
 
-class Stretching_d(DiscreteOperator):
+class ConservativeStretching(DiscreteOperator):
     """
-    Stretching operator representation.
+    Discretisation of the following problem :
+    \f{eqnarray*} \frac{\partial\omega}{\partial t} = \nabla.(\omega:v) \f}
     """
 
     @debug
-    def __init__(self, velocity, vorticity,
-                 method='FD_order4 RK3', propertyOp='divConservation'):
+    def __init__(self, velocity, vorticity, method=None):
         """
-        Constructor.
-        Create a Stretching operator on a given continuous domain.
-        Work on a given velocity and vorticity to compute the stretching term.
-
-        @param stretch : Stretching operator.
-        @param method : Method to use
+        @param velocity : discrete field
+        @param vorticity : discrete field
+        @param method : numerical method for space/time discretizations
         """
-        DiscreteOperator.__init__(self, [velocity, vorticity], method,
-                                  name="Stretching")
-        ## Velocity.
+        ## velocity discrete field
         self.velocity = velocity
-        ## Vorticity.
+        ## vorticity discrete field
         self.vorticity = vorticity
-        #        ## input fields
-        #        self.input = [self.velocity, self.vorticity]
-        self.propertyOp = propertyOp
+
+        DiscreteOperator.__init__(self, [self.velocity, self.vorticity],
+                                  method,
+                                  name="ConservativeStretching")
+
         self.input = [self.velocity, self.vorticity]
         self.output = [self.vorticity]
 
+        # \todo multiresolution case
+        assert self.velocity.topology == self.vorticity.topology,\
+            'Multiresolution case not yet implemented.'
+
+        # We must have a topo for this operator ...
+        # A function to compute the gradient of a vector field
+        self.divFunc = DivT(self.velocity.topology,
+                            method=method)
+
+        # Time integration scheme.
+        def rhs(t, y):
+            return self.divFunc(self.velocity.data, y)
+
+        # TOdo : check fd call with numpy and/or parmes fields.
+        # stability constant
         self.cststretch = 0.
+        # set time integration method
         if self.method.find('Euler') >= 0:
             self.num_method = Euler
             self.cststretch = 2.0
@@ -62,99 +76,124 @@ class Stretching_d(DiscreteOperator):
             self.num_method = RK3
             self.cststretch = 2.5127
 
+        self.timeIntegrator = self.num_method(self.vorticity.nbComponents, rhs)
+
     @timed_function
     def apply(self, simulation):
+        if simulation is None:
+            raise ValueError("Missing simulation value for computation.")
+
+        # time step
+        dt = simulation.timeStep
+        # current time
+        t = simulation.time
+        # Call time integrator
+        self.vorticity.data[0][...], self.vorticity.data[1][...], \
+            self.vorticity.data[2][...] = \
+            self.timeIntegrator(self.vorticity.data, dt, t)
+
+
+class GradStretching(DiscreteOperator):
+    """
+    Discretisation of the following problem :
+    \f{eqnarray*} \frac{\partial\omega}{\partial t}=[\nabla(v)][\omega]\f}
+    """
+
+    @debug
+    def __init__(self, velocity, vorticity, method=None):
         """
-        Apply Stretching operator.
-
-        Stretching algorithm:
-        @li 1. discretization of the divergence operator (div wu)
-        or w grad(u) : \n
-        @li 2. Time integration. Performs a method choose to resolove
-        the equation dw/dt = div(wu) or dw/dt = w grad(u).
-         -define the function of ODEsolver as the result of the first step
-         -integrate with the chosen method (RK4. RK2. euler)
+        @param velocity : discrete field
+        @param vorticity : discrete field
+        @param method : numerical method for space/time discretizations
         """
+        ## velocity discrete field
+        self.velocity = velocity
+        ## vorticity discrete field
+        self.vorticity = vorticity
+
+        DiscreteOperator.__init__(self, [self.velocity, self.vorticity],
+                                  method,
+                                  name="ConservativeStretching")
+
+        self.input = [self.velocity, self.vorticity]
+        self.output = [self.vorticity]
+
+        # \todo multiresolution case
+        assert self.velocity.topology == self.vorticity.topology,\
+            'Multiresolution case not yet implemented.'
+
+        # We must have a topo for this operator ...
+        # A function to compute the gradient of a vector field
+        self.graduv = GradVxW(self.velocity.topology, method=method)
+
+        # Time integration scheme.
+        def rhs(t, y):
+            return self.graduv(self.velocity.data, y)
+
+        # TOdo : check fd call with numpy and/or parmes fields.
+        # stability constant
+        self.cststretch = 0.
+        # set time integration method
+        if self.method.find('Euler') >= 0:
+            self.num_method = Euler
+            self.cststretch = 2.0
+        elif self.method.find('RK2') >= 0:
+            self.num_method = RK2
+            self.cststretch = 2.0
+        elif self.method.find('RK3') >= 0:
+            self.num_method = RK3
+            self.cststretch = 2.5127
+        elif self.method.find('RK4') >= 0:
+            self.num_method = RK4
+            self.cststretch = 2.7853
+        else:
+            self.num_method = RK3
+            self.cststretch = 2.5127
+
+        ## time integrator
+        self.timeIntegrator = self.num_method(self.vorticity.nbComponents, rhs)
+        ## a vector of float with diagnostics computed
+        ## from GradVxW (max div ...)
+        self.diagnostics = np.ones((5))
+
+    @timed_function
+    def apply(self, simulation):
         if simulation is None:
             raise ValueError("Missing simulation value for computation.")
 
+        # time step
         dt = simulation.timeStep
-        if self.num_method is not None:
-            condA = self.propertyOp.find('divConservation') < 0
-            if condA and self.propertyOp.find('massConservation') < 0:
-                print 'No or wrong stretching operator property specified, \
-                use default divConservation'
-                self.propertyOp = 'divConservation'
-
-            if self.propertyOp.find('massConservation') >= 0:
-                ## Determination of Stretching time step (stability condition)
-                self.ndt = 1
-                self.dtstretch = dt
-
-                ## fct2Op
-                self.function = Fct2Op(self.vorticity,
-                                       self.velocity,
-                                       choice='divWU',
-                                       topology=self.velocity.topology)
-
-            else:  # self.propertyOp == 'divConservation'
-
-                ## Space discretization of grad(u)
-                self.stretchOp = DifferentialOperator(
-                    self.vorticity, self.velocity,
-                    method='FD_order4',
-                    choice='gradV')
-                #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]
-
-                ## Determination of Stretching time step (stability condition)
-                self.ndt = self.stabilityTest(dt, maxArray[0],
-                                              self.num_method, self.cststretch)
-                self.dtstretch = dt / float(self.ndt)
-
-                ## fct2Op
-                self.function = Fct2Op(self.vorticity,
-                                       gradResult,
-                                       method='FD_order4',
-                                       choice='gradV',
-                                       topology=self.vorticity.topology)
-                self.synchro.apply()
-
-            ## Time integration
-            t = simulation.time
-            for i in range(self.ndt):
-                methodInt = self.num_method(self.function)
-#                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.synchro.apply()
-#            print 'norm omega', LA.norm(self.vorticity.data)
-
-    def stabilityTest(self, dt, maxi, num_method, cststretch):
-        dtstretch = dt
-        dtstretch = min(dtstretch, cststretch / maxi)
-        if dt == dtstretch:
-            print "dt, ndt , dtstretch, upperbound", dt, int(dt / dtstretch),\
-                  dt / float(int(dt / dtstretch)), cststretch / maxi
-            return int(dt / dtstretch)
-        else:
-            print "dt, ndt , dtstretch, upperbound", dt, int(dt / dtstretch)\
-                  + 1, dt / float(int(dt / dtstretch) + 1), cststretch / maxi
-            return int(dt / dtstretch) + 1
+        # current time
+        t = simulation.time
+        self.graduv(self.velocity, self.vorticity, self.diagnostics)
+        # Compute the number of required subcycles
+        ndt, subdt = self.checkStability(dt)
+        # Call time integrator
+        print "time steps ...", dt, subdt
+        assert sum(subdt) == dt
+        for i in xrange(ndt):
+            self.vorticity.data[0], self.vorticity.data[1],\
+                self.vorticity.data[2] =\
+                self.timeIntegrator(self.vorticity.data, subdt[i], t)
+
+    def checkStability(self, dt):
+        """
+        Computes a stability condition depending on some
+        diagnostics (from GradVxW)
+        @param current time step
+        @return int,float : the number of required subcycles and
+        the subcycle time-step.
+        """
+        dt_stab = min(dt, self.cststretch / self.diagnostics[0])
+        nb_cycles = int(ceil(dt / dt_stab))
+        subdt = np.zeros((nb_cycles))
+        subdt[:] = dt_stab
+        subdt[-1] = dt - (nb_cycles - 1) * dt_stab
+        return nb_cycles, subdt
 
 if __name__ == "__main__":
     print __doc__
-    print "- Provided class : Stretching_d"
-    print Stretching_d.__doc__
+    print "- Provided class : ConservativeStretching, GradStretching"
+    print ConservativeStretching.__doc__
+
+

HySoP/hysop/operator/discrete/updateGhosts.py

@@ -117,8 +117,10 @@ class UpdateGhosts_d(DiscreteOperator):
         comm = self.topology.comm
         rank = self.topology.rank
         neighbours = self.topology.neighbours
-        exchange_dir = [d for d in xrange(self._dim)
-                        if self.topology.cutdir[d]]
+        exchange_dir = []
+        if self.topology.size > 1:
+            exchange_dir = [d for d in xrange(self._dim)
+                            if self.topology.cutdir[d]]
         i = 0
         for d in exchange_dir:
             # 1 - Fill in buffers

HySoP/hysop/operator/monitors/compute_forces.py

@@ -64,7 +64,7 @@ class Compute_forces(Monitoring):
         nbGhosts = 0
         if self.method.find('FD_order2') >= 0:
             nbGhosts = 1
-        elif self.method.find('FD_order4') >= 0:
+        elif self.method.find('FD_C4') >= 0:
             nbGhosts = 2
         else:
             nbGhosts = 2

HySoP/hysop/operator/reprojection.py

@@ -48,7 +48,7 @@ class Reprojection_criterion(Monitoring):
         nbGhosts = 0
         if self.method.find('FD_order2') >= 0:
             nbGhosts = 1
-        elif self.method.find('FD_order4') >= 0:
+        elif self.method.find('FD_C4') >= 0:
             nbGhosts = 2
         else:
             nbGhosts = 2

HySoP/hysop/operator/stretching.py

@@ -2,33 +2,40 @@
 """
 @file operator/stretching.py
 
-Stretching operator representation
-
+Computation of stretch. term in Navier-Stokes.
 
 """
-from parmepy.constants import debug
+from parmepy.constants import debug, GRADUW, CONSERVATIVE
 from parmepy.operator.continuous import Operator
-from parmepy.operator.discrete.stretching import Stretching_d
+from parmepy.operator.discrete.stretching import ConservativeStretching
+from parmepy.operator.discrete.stretching import GradStretching
+from parmepy.operator.updateGhosts import UpdateGhosts
 import numpy as np
 
 
 class Stretching(Operator):
     """
-    Stretching operator representation
     \todo write latex formulas
     """
 
     @debug
-    def __init__(self, velocity, vorticity,
-                 resolutions=None, method='FD_order4 RK3',
-                 topo=None, ghosts=None, **other_config):
+    def __init__(self, velocity, vorticity, resolutions,
+                 formulation=None, method=None, topo=None, ghosts=None,
+                 ghostUpdate=None):
         """
-        Constructor.
         Create a Stretching operator from given
         velocity and vorticity variables.
+        'FD_C4 RK3',
 
-        @param velocity ContinuousVectorField : velocity variable.
-        @param vorticity ContinuousVectorField : vorticity variable.
+        @param velocity field
+        @param vorticity field
+        @param resolutions : grid resolution of velocity and vorticity
+        @param method : solving method
+        (default = finite differences, 4th order, in space
+        and Runge-Kutta 3 in time.)
+        @param topo : a predefined topology to discretize velocity/vorticity
+        @param ghosts : number of ghosts points. Default depends on the method.
+        Autom. computed if not set.
         """
         Operator.__init__(self, [velocity, vorticity], method, topo=topo,
                           ghosts=ghosts, name="Stretching")
@@ -36,64 +43,91 @@ class Stretching(Operator):
         self.velocity = velocity
         ## vorticity variable (vector)
         self.vorticity = vorticity
+        ## Grid resolution for each variable (dictionnary)
         self.resolutions = resolutions
-        self.config = other_config
+        ## Numerical methods for time and space discretization
+        if method is None:
+            self.method = 'FD_C4_RK3'
+        else:
+            self.method = method
+        ## Formulation used for the stretching equation.
+        ## Default = conservative form.
+        if formulation is not None:
+            self.formulation = formulation
+        else:
+            self.formulation = CONSERVATIVE
+        ## Do this operator deals with ghost points update?
+        ## default = yes
+        if ghostUpdate is None:
+            self.ghostUpdate = True
+        else:
+            self.ghostUpdate = ghostUpdate
         self.input = [self.velocity, self.vorticity]
         self.output = [self.vorticity]
+        ## operator to update ghost points of the input variables
+        self.synchronize = None
 
     @debug
     def setUp(self):
-        """
-        Stretching operator discretization method.
-        Create a discrete Stretching operator from given specifications.
-
-        @param idVelocityD : Index of velocity discretisation to use.
-        @param idVorticityD : Index of vorticity discretisation to use.
-        @param propertyOp : property garantied by the stretching Op
-            (divConservation, massConservation) default : divConservation
-        @param method : the method to use :
-            space : (FD_order2, FD_order4, ..) default : FD_order4
-            and time : (Euler, RK, ..) default : RK3
-        """
         print '================ Stretching scheme ====================='
         print ' Space discretization and integr. scheme:', self.method
         print '========================================================'
         nbGhosts = 0
         if self.method.find('FD_order2') >= 0:
             nbGhosts = 1
-        elif self.method.find('FD_order4') >= 0:
+        elif self.method.find('FD_C4') >= 0:
             nbGhosts = 2
         else:
-            nbGhosts = 2
+            raise ValueError("Unknown method for stretching operator.")
 
         if self._predefinedTopo is not None:
             assert (self._predefinedTopo.ghosts >= nbGhosts).all()
             topo = self._predefinedTopo
             for v in self.variables:
                 self.discreteFields[v] = v.discretize(topo)
+
         else:
             # Variables discretization
             if self.ghosts is not None:
                 self.ghosts[self.ghosts < nbGhosts] = nbGhosts
             else:
-                self.ghosts = np.zeros((self.domain.dimension)) * nbGhosts
+                self.ghosts = np.ones((self.domain.dimension)) * nbGhosts
+
             for v in self.variables:
-                # the topology for v ...
                 topo = self.domain.getOrCreateTopology(self.domain.dimension,
                                                        self.resolutions[v],
                                                        ghosts=self.ghosts)
-
-                # ... and the corresponding discrete field
                 self.discreteFields[v] = v.discretize(topo)
 
-        self.discreteOperator = \
-            Stretching_d(self.discreteFields[self.velocity],
-                         self.discreteFields[self.vorticity],
-                         method=self.method, **self.config)
-
+        if self.formulation is CONSERVATIVE:
+            self.discreteOperator =\
+                ConservativeStretching(self.discreteFields[self.velocity],
+                                       self.discreteFields[self.vorticity],
+                                       method=self.method)
+            if self.ghostUpdate:
+                topov = self.discreteFields[self.velocity].topology
+                self.synchronize = UpdateGhosts(self.input, topo=topov)
+                self.synchronize.setUp()
+        elif self.formulation is GRADUW:
+            self.discreteOperator =\
+                GradStretching(self.discreteFields[self.velocity],
+                               self.discreteFields[self.vorticity],
+                               method=self.method)
+            if self.ghostUpdate:
+                topov = self.discreteFields[self.velocity].topology
+                self.synchronize = UpdateGhosts(self.velocity, topo=topov)
+                self.synchronize.setUp()
         self.discreteOperator.setUp()
         self._isUpToDate = True
 
+    def apply(self, simulation=None):
+        # update ghost points, if required
+        if self.synchronize is not None:
+            self.synchronize.apply()
+        # computation ...
+        self.discreteOperator.apply(simulation)
+
+
 if __name__ == "__main__":
     print __doc__
     print "- Provided class : Stretching"

HySoP/hysop/operator/tests/test_Stretching.py

 # -*- coding: utf-8 -*-
-import time
-
 import parmepy as pp
 import numpy as np
 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):
@@ -77,7 +73,7 @@ def test_Stretching():
     stretch = Stretching(velo, vorti,
                          resolutions={velo: nbElem,
                                       vorti: nbElem},
-                         method='FD_order4 RK3',
-                         propertyOp='divConservation')
+                         method='FD_C4_RK3')
     stretch.setUp()
     simulation = Simulation(0, timeStep, 20)
+    stretch.apply(simulation)

HySoP/hysop/operator/updateGhosts.py

@@ -23,7 +23,7 @@ class UpdateGhosts(Operator):
         """
         Operator.__init__(self, variables, method, topo=topo)
 
-        self.topo = topo
+        self.topo = self._predefinedTopo
 
     @debug
     def setUp(self):