updates for new topology (not tested...) + pep8

HySoP/hysop/__init__.py.in

@@ -10,6 +10,7 @@ __all__ = ['Box', 'CartesianTopology']#, 'ScalarField']
 
 # MPI
 if("@USE_MPI@" is "ON"):
+    #import parmepy.mpi as mpi
     from mpi import topology
     CartesianTopology = topology.CartesianTopology
     if(mpi.main_rank == 0):

HySoP/hysop/mpi/topology.py

@@ -387,40 +387,39 @@ class SubMesh(object):
         ## Topology that creates (and owns) this mesh
         self._topology = topo
         ## Local resolution of this mesh, including ghost points
-        self._resolution = resolution
+        self.resolution = resolution
         ## index of the lowest point of this mesh
         ## (in each dir) in the global mesh
-        self._global_start = g_start
+        self.global_start = g_start
         ## index of the upper point (in each dir), global mesh
-        self._global_end = self._global_start + self._resolution - 1
+        self.global_end = self.global_start + self.resolution - 1
         ## Mesh step size
-        self._space_step_size = topo.domain.length / \
+        self.space_step_size = topo.domain.length / \
             (topo.globalMeshResolution - 1)
         ## Mesh local indices, only for "computed" points
         ## (i.e. excluding ghosts)
-        self._local_start = topo.ghosts.copy()
-        self._local_end = self._resolution.copy() - topo.ghosts[:] - 1
+        self.local_start = topo.ghosts.copy()
+        self.local_end = self.resolution.copy() - topo.ghosts[:] - 1
         ## Coordinates of the "lowest" point of this mesh (including ghost)
         ## Warning FP : this strongly depends on where is the origin
         ## of the domain, of the type of boundary conditions
         ## and if origin is on ghosts or "real" points.
-        self._origin = topo.domain.origin.copy()
-        self._origin[:] += self._space_step_size * \
-            (self._global_start[:] - topo.ghosts[:])
+        self.origin = topo.domain.origin.copy()
+        self.origin[:] += self.space_step_size * \
+            (self.global_start[:] - topo.ghosts[:])
 
     def __str__(self):
         """ Sub mesh display """
         s = 'Coords:' + str(self._topology.coords[:])
-        s += ' Sub mesh resolution:' + str(self._resolution) + '\n'
-        s += 'Starting point global indices :' + str(self._global_start) + '\n'
+        s += ' Sub mesh resolution:' + str(self.resolution) + '\n'
+        s += 'Starting point global indices :' + str(self.global_start) + '\n'
         s += 'Local indices for "computed" points (excluding ghosts) :\n'
-        s += '   start = ' + str(self._local_start)
-        s += ', end = ' + str(self._local_end) + '\n'
+        s += '   start = ' + str(self.local_start)
+        s += ', end = ' + str(self.local_end) + '\n'
         return s
 
 
 if __name__ == "__main__":
     print "This module defines the following classes:"
     print "- CartesianTopology: ", CartesianTopology.__doc__
-    print "- FFTTopology: ", FFTTopology.__doc__
     print "- LocalMesh: ", LocalMesh.__doc__

HySoP/hysop/operator/stretching.py

@@ -4,9 +4,9 @@
 
 Penalization operator representation
 """
-from .continuous import ContinuousOperator
-from .stretching_d import Stretching_d
-from ..particular_solvers.integrator.runge_kutta3 import RK3
+from continuous import ContinuousOperator
+from stretching_d import Stretching_d
+from parmepy.integrator import RK3
 
 
 class Stretching(ContinuousOperator):

HySoP/hysop/operator/stretching_d.py

 # -*- coding: utf-8 -*-
 """
-@package operator
+@package parmepy.operator.stretching_d
+
 Discrete stretching representation
 """
-from ..constants import *
-from .discrete import DiscreteOperator
-from .differentialOperator import DifferentialOperator
-from .differentialOperator_d import DifferentialOperator_d
-from .fct2op import Fct2Op
-from ..particular_solvers.integrator.euler import Euler
-from ..particular_solvers.integrator.runge_kutta2 import RK2
-from ..particular_solvers.integrator.runge_kutta3 import RK3
-from ..particular_solvers.integrator.runge_kutta4 import RK4
-from ..particular_solvers.integrator.integrator import ODESolver
-from ..tools.cpu_data_transfer import Synchronize
-import numpy as np
-import numpy.testing as npt
+
+from discrete import DiscreteOperator
+from differentialOperator_d import DifferentialOperator_d
+from fct2op import Fct2Op
+import parmepy.integrator as integrator
 import time
-import sys
 
 
 class Stretching_d(DiscreteOperator):
@@ -26,7 +18,8 @@ class Stretching_d(DiscreteOperator):
     DiscreteOperator.DiscreteOperator specialization.
     """
 
-    def __init__(self, stretch, topology=None, idVelocityD=0, idVorticityD=0, propertyOp='divConservation', method=None):
+    def __init__(self, stretch, topology=None, idVelocityD=0, idVorticityD=0,
+                 propertyOp='divConservation', method=None):
         """
         Constructor.
         Create a Stretching operator on a given continuous domain.
@@ -46,7 +39,6 @@ class Stretching_d(DiscreteOperator):
         self.method = method
         ## self.name = "stretching"
         self.topology = topology
-        self.OpSynchronize = Synchronize(self.topology)
 
     def apply(self, t, dt):
         """
@@ -56,8 +48,10 @@ class Stretching_d(DiscreteOperator):
         @param dt : time step.
 
         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).
+        @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)
         """
@@ -65,43 +59,48 @@ class Stretching_d(DiscreteOperator):
 
         if self.method is not None:
 
-            if self.propertyOp == 'divConservation' :
+            if self.propertyOp == 'divConservation':
 
                 ## Space discretisation of grad(u)
-                self.stretchOp = DifferentialOperator_d(self.vorticity, self.velocity, choice='gradV', topology=self.topology)
+                self.stretchOp = DifferentialOperator_d(
+                    self.vorticity, self.velocity, choice='gradV',
+                    topology=self.topology)
+
                 gradResult, maxgersh = self.stretchOp.apply()
 
                 ## Determination of Stretching time step (stability condition)
                 self.ndt = self.stabilityTest(dt, maxgersh[0], self.method)
-                self.dtstretch = dt/float(self.ndt)
+                self.dtstretch = dt / float(self.ndt)
 
                 print "Maxgersh", maxgersh
                 print "dtStretch, ndt", self.dtstretch, self.ndt
-                print "dtAdvec", 1./maxgersh[1]
+                print "dtAdvec", 0.5 / maxgersh[1]
 
                 ## fct2Op
-                self.fonction = Fct2Op(self.vorticity.data, gradResult,choice = 'gradV', topology=self.topology)
+                self.fonction = Fct2Op(self.vorticity.data,
+                                       gradResult, choice='gradV',
+                                       topology=self.topology)
 
-            elif self.propertyOp == 'massConservation' :
+            elif self.propertyOp == 'massConservation':
                 ## Determination of Stretching time step (stability condition)
                 self.ndt = 1
                 self.dtstretch = dt
 
                 ## fct2Op
-                self.fonction = Fct2Op(self.vorticity, self.velocity,choice = 'divWU', topology=self.topology)
+                self.fonction = Fct2Op(self.vorticity,
+                                       self.velocity, choice='divWU',
+                                       topology=self.topology)
 
             else:
-                raise ValueError("Unknown Stretching Operator property: " + str(self.propertyOp))
-
-            # Time integration 
-            for i in range (self.ndt) :
-                methodInt=self.method(self.fonction)
-                self.vorticity.data = methodInt.integrate(self.fonction.apply, t, self.dtstretch, self.vorticity.data)
-
-
-            # Ghosts synchronization
-#            self.OpSynchronize.apply([self.vorticity.data, self.velocity.data])
+                raise ValueError("Unknown Stretching Operator property: "
+                                 + str(self.propertyOp))
 
+            # Time integration
+            for i in range(self.ndt):
+                methodInt = self.method(self.fonction)
+                self.vorticity.data = methodInt.integrate(self.fonction.apply,
+                                                          t, self.dtstretch,
+                                                          self.vorticity.data)
 
             self.compute_time = time.time() - self.compute_time
             self.total_time += self.compute_time
@@ -110,25 +109,27 @@ class Stretching_d(DiscreteOperator):
 
     def stabilityTest(self, dt, maxi, method):
         dtstretch = dt
-        cststretch =0.
-        if method == Euler :
+        cststretch = 0.
+        if method == integrator.Euler:
             cststretch = 2.0
-            dtstretch = min(dtstretch, cststretch/maxi)
-        if method == RK2 :
+            dtstretch = min(dtstretch, cststretch / maxi)
+        if method == integrator.RK2:
             cststretch = 2.0
-            dtstretch = min(dtstretch, cststretch/maxi)
-        if method == RK3 :
+            dtstretch = min(dtstretch, cststretch / maxi)
+        if method == integrator.RK3:
             cststretch = 2.5127
-            dtstretch = min(dtstretch, cststretch/maxi)
-        if method == RK4 :
+            dtstretch = min(dtstretch, cststretch / maxi)
+        if method == integrator.RK4:
             cststretch = 2.7853
-            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
+            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
 
     def printComputeTime(self):
         self.timings_info[0] = "\"Stretching total\""
@@ -137,7 +138,8 @@ class Stretching_d(DiscreteOperator):
         print "Time of the last Stretching iteration :", self.compute_time
 
     def __str__(self):
-        s = "Stretching_d (DiscreteOperator). " + DiscreteOperator.__str__(self)
+        s = "Stretching_d (DiscreteOperator). " +\
+            DiscreteOperator.__str__(self)
         return s
 
 if __name__ == "__main__":

HySoP/hysop/particular_solvers/basic.py

@@ -17,6 +17,7 @@ from parmepy.tools.timer import Timer
 from parmepy.tools.printer import Printer
 from parmepy.integrator import RK3
 
+
 class ParticleSolver(Solver):
     """
     Particular solver description.
@@ -46,7 +47,8 @@ class ParticleSolver(Solver):
         Solver.__init__(self, problem)
         ## Problem to solve
         self.problem = problem
-        ## Advection operator of the problem. Advection need special treatment in particle methods.
+        ## Advection operator of the problem. Advection need special
+        ## treatment in particle methods.
         self.advection = None
         self.stretch = None
         self.penal = None
@@ -87,19 +89,23 @@ class ParticleSolver(Solver):
             #if isinstance(op, Velocity)
                 #self.velocity = op
         #if self.advection is None:
-        #    raise ValueError("Particular Solver : Cannot create from a problem with no Transport operator")
+        #  raise ValueError("Particular Solver :
+        # Cannot create from a problem with no Transport operator")
 
     def initialize(self):
         """
         Solver initialisation.
         Initialize a particle method solver regarding the problem.
         """
-        self.p_position = ContinuousField(domain=self.problem.domains[0], name='ParticlePosition')
-        self.p_scalar = ContinuousField(domain=self.problem.domains[0], name='ParticleScalar')
+        self.p_position = ContinuousField(domain=self.problem.domains[0],
+                                          name='ParticlePosition')
+        self.p_scalar = ContinuousField(domain=self.problem.domains[0],
+                                        name='ParticleScalar')
         self.problem.addVariable([self.p_position, self.p_scalar])
         ## Discretise domains
+        ## Note FP : this is probably useless???
         for d in self.problem.domains:
-            d.discretize(self.problem.topology.resolution)
+            d.discretize(self.problem.topology.globalMeshResolution)
         ## Discretise fields and initialize
         for v in self.problem.variables:
             v.discretize(self.problem.topology)
@@ -110,9 +116,12 @@ class ParticleSolver(Solver):
         ## Discretise operators
         for op in self.problem.operators:
             if op is self.advection:
-                op.discretize(result_position=self.p_position, result_scalar=self.p_scalar, method=self.ODESolver)
+                op.discretize(result_position=self.p_position,
+                              result_scalar=self.p_scalar,
+                              method=self.ODESolver)
             if op is self.stretch:
-                op.discretize(topology=self.problem.topology, method=self.ODESolver)
+                op.discretize(topology=self.problem.topology,
+                              method=self.ODESolver)
             if op is self.penal:
                 op.obstacle.discretize(self.problem.topology)
                 op.discretize()
@@ -124,22 +133,27 @@ class ParticleSolver(Solver):
                 op.discretize(topology=self.problem.topology)
             ## Velocity is only program on gpu for instance
             #elif op is self.velocity:
-            #    op.discretize(result_position=self.p_position, result_scalar=self.p_scalar, method=self.ODESolver)
+            #    op.discretize(result_position=self.p_position,
+            # result_scalar=self.p_scalar, method=self.ODESolver)
             else:
                 op.discretize()
         if self.advection is not None and not self.advection.include_remesh:
             ## Create remeshing operator
-            self.remeshing = Remeshing(partPositions=self.advection.discreteOperator.res_position,
-                                    partScalar=self.advection.discreteOperator.res_scalar,
-                                    resscal=self.advection.discreteOperator.scalar,
-                                    method=self.RemeshingMethod)
+            self.remeshing = Remeshing(
+                partPositions=self.advection.discreteOperator.res_position,
+                partScalar=self.advection.discreteOperator.res_scalar,
+                resscal=self.advection.discreteOperator.scalar,
+                method=self.RemeshingMethod)
+
         for op in self.problem.operators:
             if op.needSplitting:
                 index = self.problem.operators.index(op)
                 if op is self.advection and not self.advection.include_remesh:
-                    self.problem.operators[index] = Splitting([op, self.remeshing], dim=self.problem.topology.dim)
+                    self.problem.operators[index] = Splitting(
+                        [op, self.remeshing], dim=self.problem.topology.dim)
                 else:
-                    self.problem.operators[index] = Splitting([op], dim=self.problem.topology.dim)
+                    self.problem.operators[index] = Splitting(
+                        [op], dim=self.problem.topology.dim)
         self.isInitialized = True
 
     def end(self):
@@ -147,10 +161,10 @@ class ParticleSolver(Solver):
 
     def __str__(self):
         """ToString method"""
-        s = " Particular solver "
+        s = " Particle solver "
         return s
 
 if __name__ == "__main__":
     print __doc__
-    print "- Provided class : ParticularSolver"
-    print ParticularSolver.__doc__
+    print "- Provided class : ParticleSolver"
+    print ParticleSolver.__doc__