whole sequential Navier-Stokes problem

Examples/NavierStokes3d_linked.py

+#!/usr/bin/python
+import parmepy as pp
+import parmepy.f2py
+import numpy as np
+import mpi4py.MPI as MPI
+import math as m
+
+PARMES_REAL = pp.constants.PARMES_REAL
+ORDER = pp.constants.ORDER
+
+#from numpy import linalg as LA
+
+pi = m.pi
+
+## Import scales and fftw solvers
+ppfft = parmepy.f2py.fftw2py
+scales = parmepy.f2py.scales2py
+
+rank = MPI.COMM_WORLD.Get_rank()
+print "Mpi process number ", rank
+
+## ----------- A 3d problem -----------
+print " ========= Start test for Navier-Stokes 3D ========="
+
+## Physical Domain description
+Lx = Ly = Lz = 2 * pi
+myDomain3d = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[0., 0., 0.])
+resolution3d = np.asarray((129, 129, 129))
+ncells = resolution3d - 1
+hx = Lx / ncells[0]
+hy = Ly / ncells[1]
+hz = Lz / ncells[2]
+
+## Obstacle
+lambd = np.array([0, 10, 10E8], dtype=PARMES_REAL, order=ORDER)
+sphere = pp.Obstacle(myDomain3d, name='sphere', zlayer=0.1,
+                     radius=0.5, center=[Lx/2., Ly/2., Lz/2.],
+                     orientation='West', porousLayer=0.)
+
+## Post
+outputFilePrefix = './res/NS_'
+outputModulo = 1
+
+## Simulation parameters
+timeStep = 1e-8
+finalTime = timeStep#1.0
+
+
+## Topologies definitions
+# 1---- FFT topology + diffusion/poisson solvers initialization -----------
+localres, localoffset = ppfft.init_fftw_solver(resolution3d,
+                                               myDomain3d.length)
+print "FFT solver local resolution/offset: ", localres, localoffset
+##topofft = poisson.getTopology()
+
+# ---- Cartesian topology -----
+topoCart = pp.CartesianTopology(domain=myDomain3d,
+                               resolution=resolution3d, dim=3, ghosts=[2,2,2])
+
+# ---- JB's topology ----------
+nbProcs = MPI.COMM_WORLD.Get_size()
+topodims = [1, 1, nbProcs]  # fits with fftw topology
+
+
+## Fields declaration
+def computeVel(x, y, z):
+#    vx = 2. / np.sqrt(3) * np.sin(2. * pi / 3.) * np.sin(x) \
+#        * np.cos(y) * np.cos(z)
+#    vy = 2. / np.sqrt(3) * np.sin(-2. * pi / 3.) * np.cos(x) \
+#        * np.sin(y) * np.cos(z)
+#    vz = 0.
+    uinf = 1.
+    vx = 0.
+    module = (x-Lx/2)**2 + (y-Ly/2)**2 + (z-Lz/2)**2
+    if (module >= (0.5)**2):
+        vx = uinf * (1-(0.5)**2/module)
+    vy = 0.
+    vz = 0.
+    return vx, vy, vz
+
+
+def computeVort(x, y, z):
+#    wx = - np.cos(x) * np.sin(y) * np.sin(z)
+#    wy = - np.sin(x) * np.cos(y) * np.sin(z)
+#    wz = 2. * np.sin(x) * np.sin(y) * np.cos(z)
+    wx = 0.
+    wy = 0.
+    wz = 0.
+    return wx, wy, wz
+
+def computeDensity(x, y, z):
+    if (x>=0):
+        rho = 1000.
+    else : 
+        rho = 1.
+    return rho
+
+velocity = pp.AnalyticalField(domain=myDomain3d, formula=computeVel,
+                              name='Velocity', vector=True)
+vorticity = pp.AnalyticalField(domain=myDomain3d, formula=computeVort,
+                               name='Vorticity', vector=True)
+#density = pp.AnalyticalField(domain=myDomain3d, formula=computeDensity,
+#                             name='Density', vector=False)
+
+############ REF ##############
+x = np.arange(localoffset[0], localres[0] + localoffset[0],
+              dtype='float64') * hx
+y = np.arange(localoffset[1], localres[1] + localoffset[1],
+              dtype='float64') * hy
+z = np.arange(localoffset[2], localres[2] + localoffset[2],
+              dtype='float64') * hz
+
+cden = 4 * pi ** 2 * (Ly ** 2 * Lz ** 2 + Lx ** 2 * Lz ** 2 +
+                      Lx ** 2 * Ly ** 2) / (Lx ** 2 * Ly ** 2 * Lz ** 2)
+cx = 2 * pi / Lx
+cy = 2 * pi / Ly
+cz = 2 * pi / Lz
+
+
+## 2 - Advection solver initialization. See testScales for a working example.
+# Based on scales JB solver
+# Warning : fields input for scales should be of size (ncells), not localres.
+
+scalesres, scalesoffset, stab_coeff = \
+    scales.init_advection_solver(ncells, myDomain3d.length,
+                                 topodims, order='p_O2')
+
+advecVort = pp.Advec_scales(stab_coeff, velocity, vorticity)
+#advecDensity = pp.Advec_scales(stab_coeff, velocity, density) 
+
+# ----> Change TOPO from Scales to Cartesian
+
+## 3 - Stretching
+stretch = pp.Stretching(vorticity, velocity)
+
+# ----> Change TOPO from Cartesian to FFT
+
+## 4 - Diffusion
+diffusion = pp.Diffusion(velocity, vorticity, viscosity=0.001)
+
+## 4 - Poisson
+poisson = pp.Poisson(velocity, vorticity)
+
+# ----> Change TOPO from FFT to Cartesian
+
+## 6 - Penalization
+penal = pp.Penalization(velocity, vorticity, sphere, lambd)
+
+# ----> Change TOPO from Cartesian to Scales
+
+## Define the problem to solve
+navierStokes = pp.Problem(topoCart, [advecVort, stretch, diffusion, poisson, penal])
+
+printer = pp.Printer(fields=[vorticity, velocity], frequency=outputModulo,
+                     outputPrefix=outputFilePrefix)
+
+navierStokes.setSolver(finalTime, timeStep, solver_type='basic', io=printer)
+
+## Problem => ParticularSover = basic.initialize()
+navierStokes.initSolver()
+
+penal.apply(0., timeStep)
+
+## solve stretching
+navierStokes.solve()
+
+print 'vorticity', vorticity.discreteField[0][0] , vorticity.discreteField[0][1] ,  vorticity.discreteField[0][2]
+
+## end of time loop ##
+
+# Clean memory buffers
+ppfft.clean_fftw_solver(myDomain3d.dimension)

HySoP/hysop/operator/advec_scales.py

+# -*- coding: utf-8 -*-
+"""
+@package parmepy.operator Advec_scales
+Continous Advection operator using Scales code (Jean-Baptiste)
+"""
+from .continuous import ContinuousOperator
+from .advec_scales_d import Advec_scales_d
+
+
+class Advec_scales(ContinuousOperator):
+    """
+    Advection scales operator representation
+    """
+
+    def __init__(self, stab_coeff=None, velocity=None, advectedField=None):
+        """
+        Constructor.
+        Create an Advection operator using scales code (JB) from given velocity variables.
+
+        @param velocity ContinuousVectorField : velocity variable.
+        @param advectedField ContinuousVector/ScalarField : variable to advect.
+        """
+        ContinuousOperator.__init__(self)
+        self.stab_coeff = stab_coeff
+        self.velocity = velocity
+        self.advectedField = advectedField
+        self.addVariable([velocity, advectedField])
+
+    def discretize(self, idVelocityD=0, idAdvectedFieldD=0, topology=None):
+        """
+        Advection operator discretization method.
+        Create a discrete advection operator from given specifications.
+
+        @param idVelocityD : Index of velocity discretisation to use for advection.
+        @param idAdvectedFieldD : Index of field discretisation to advect.
+        @param topology : topology of the input fields
+        """
+        if self.discreteOperator is None:
+            self.discreteOperator = Advec_scales_d(self, idVelocityD, idAdvectedFieldD, self.stab_coeff, topology)
+
+    def __str__(self):
+        """ToString method"""
+        s = " Advection scales operator (ContinuousOperator)"
+        if self.discreteOperator is not None:
+            s += "with the following discretization:\n"
+            s += str(self.discreteOperator)
+        else:
+            s += ". Not discretised"
+        return s + "\n"
+
+if __name__ == "__main__":
+    print __doc__
+    print "- Provided class : Advec_scales"
+    print Advec_scales.__doc__

HySoP/hysop/operator/advec_scales_d.py

+# -*- coding: utf-8 -*-
+"""
+@package operator
+Discrete Advection operator using scales code (Jean-Baptiste)
+"""
+from ..f2py import scales2py
+import mpi4py.MPI as MPI
+from .discrete import DiscreteOperator
+from ..constants import *
+import numpy as np
+import time
+import sys
+
+class Advec_scales_d(DiscreteOperator):
+    """
+    Operator representation.
+    """
+
+    def __init__(self, advec, idVelocityD=0, idAdvectedFieldD=0, stab_coeff=None, topology=None):
+        """
+        Constructor.
+        @param stab_coeff : stability coefficient
+        @param velocity : descretized velocity to use.
+        @param advectedField : discretized field to advect.
+        """
+        DiscreteOperator.__init__(self)
+        self.advec = advec
+        self.velocity = advec.velocity.discreteField[idVelocityD]
+        self.advectedField = advec.advectedField.discreteField[idAdvectedFieldD]
+        self.stab_coeff = stab_coeff
+        self.topology = topology
+        self.ghosts = topology.ghosts
+        self.resolution = topology.mesh.resolution
+        self.dim = topology.dim
+#        self.scales = parmepy.f2py.scales2py
+        self.compute_time = 0.
+
+    def apply(self, t, dt):
+        """
+        Apply operator.
+        """
+        self.compute_time = time.time()
+
+        ind0a = self.ghosts[0]
+        ind0b=self.resolution[0]-self.ghosts[0]
+        ind1a = self.ghosts[1]
+        ind1b=self.resolution[1]-self.ghosts[1]
+        ind2a = self.ghosts[2]
+        ind2b=self.resolution[2]-self.ghosts[2]
+        ind0=np.arange(ind0a,ind0b)
+        ind1=np.arange(ind1a,ind1b)
+        ind2=np.arange(ind2a,ind2b)
+
+        ndt = 1
+        if (dt % self.stab_coeff == 0):
+            ndt = int(dt / self.stab_coeff)
+        else :
+            ndt = int(dt / self.stab_coeff) + 1
+
+        dtadvec = dt/float(ndt)
+        for i in range (ndt) :
+            print 'advection cycle n°', i
+            if self.advectedField.vector :
+
+                advecFieldNoG = [np.zeros((self.resolution - 2 * self.ghosts), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dim)]
+                velocityNoG = [np.zeros((self.resolution - 2 * self.ghosts), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dim)]
+                for i in xrange (self.dim):
+                    advecFieldNoG[i] = self.advectedField[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
+                    velocityNoG[i] = self.velocity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
+                print 'shape vortiNoG, veloNoG', np.asarray(advecFieldNoG).shape, np.asarray(velocityNoG).shape
+
+                # Vector field advection
+                advecFieldNoG[0][...] = scales2py.solve_advection(dtadvec, velocityNoG[0][...], \
+                                             velocityNoG[1][...], velocityNoG[2][...], advecFieldNoG[0][...])
+                advecFieldNoG[1][...] = scales2py.solve_advection(dtadvec, velocityNoG[0][...], \
+                                             velocityNoG[1][...], velocityNoG[2][...], advecFieldNoG[1][...])
+                advecFieldNoG[2][...] = scales2py.solve_advection(dtadvec, velocityNoG[0][...], \
+                                             velocityNoG[1][...], velocityNoG[2][...], advecFieldNoG[2][...])
+
+                for i in xrange (self.dim):
+                    self.advectedField[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = advecFieldNoG[i][...]
+                print 'shape vorti, velo', np.asarray(self.advectedField.data).shape, np.asarray(self.velocity.data).shape
+
+            else :
+
+                advecFieldNoG = self.advectedField[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
+                velocityNoG = self.velocity[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
+                # Scalar field advection
+                advecFieldNoG[...] = scales2py.solve_advection(dtadvec, velocityNoG[0][...], \
+                                          velocityNoG[1][...], velocityNoG[2][...], advecFieldNoG[...])
+
+                self.advectedField[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = advecFieldNoG[...]
+
+        self.compute_time = time.time() - self.compute_time
+        self.total_time += self.compute_time
+        return self.compute_time
+
+    def printComputeTime(self):
+        self.timings_info[0] = "\"Advection calculation total\""
+        self.timings_info[1] = str(self.total_time)
+#        print "Advection scales total time : ", self.total_time
+        print "Time of the last advection calculation loop :", self.compute_time 
+
+    def __str__(self):
+        s = "Advec_scales_d (DiscreteOperator). " + DiscreteOperator.__str__(self)
+        return s
+
+if __name__ == "__main__":
+    print __doc__
+    print "- Provided class : Advec_scales_d"
+    print Advec_scales_d.__doc__

HySoP/hysop/operator/diffusion.py

+# -*- coding: utf-8 -*-
+"""
+@package parmepy.operator Diffusion
+Continous Diffusion operator (F2PY)
+"""
+from .continuous import ContinuousOperator
+from .diffusion_d import Diffusion_d
+
+
+class Diffusion(ContinuousOperator):
+    """
+    Diffusion operator representation using FFT
+    """
+
+    def __init__(self, velocity=None, vorticity=None, viscosity=0.1):
+        """
+        Constructor.
+        Create a Diffusion operator using FFT.
+
+        @param velocity ContinuousVectorField : velocity variable.
+        @param vorticity ContinuousVectorField : vorticity variable.
+        @param viscosity : viscosity of the considered medium.
+        """
+        ContinuousOperator.__init__(self)
+        self.velocity = velocity
+        self.vorticity = vorticity
+        self.viscosity = viscosity
+        self.addVariable([velocity, vorticity])
+
+    def discretize(self, idVelocityD=0, idVorticityD=0, topology=None):
+        """
+        Diffusion operator discretization method.
+        Create a discrete Diffusion operator from given specifications.
+
+        @param idVelocityD : Index of velocity discretisation to use.
+        @param idVelocityD : Index of vorticity discretisation to update.
+        """
+        if self.discreteOperator is None:
+            self.discreteOperator = Diffusion_d(self, idVelocityD, idVorticityD, self.viscosity, topology)
+
+    def __str__(self):
+        """ToString method"""
+        s = " Diffusion operator (ContinuousOperator)"
+        if self.discreteOperator is not None:
+            s += "with the following discretization:\n"
+            s += str(self.discreteOperator)
+        else:
+            s += ". Not discretised"
+        return s + "\n"
+
+if __name__ == "__main__":
+    print __doc__
+    print "- Provided class : Diffusion"
+    print Diffusion.__doc__

HySoP/hysop/operator/diffusion_d.py

+# -*- coding: utf-8 -*-
+"""
+@package operator
+Discrete Diffusion operator using FFT
+"""
+from ..f2py import fftw2py
+import mpi4py.MPI as MPI
+from .discrete import DiscreteOperator
+from ..constants import *
+import numpy as np
+import time
+import sys
+
+class Diffusion_d(DiscreteOperator):
+    """
+    Operator representation.
+    """
+
+    def __init__(self, diff, idVelocityD=0, idVorticityD=0, viscosity=0.1, topology=None):
+        """
+        Constructor.
+        @param velocity : descretized velocity to use.
+        @param vorticity : discretized vorticity to update.
+        """
+        DiscreteOperator.__init__(self)
+        self.diff = diff
+        self.velocity = diff.velocity.discreteField[idVelocityD]
+        self.vorticity = diff.vorticity.discreteField[idVorticityD]
+        self.viscosity = viscosity
+        self.topology = topology
+        self.ghosts = topology.ghosts
+        self.resolution = topology.mesh.resolution
+        self.dim = topology.dim
+#        self.ppfft = parmepy.f2py.fftw2py
+        self.compute_time = 0.
+
+    def apply(self, t, dt):
+        """
+        Apply operator.
+        """
+        self.compute_time = time.time()
+
+        ind0a = self.ghosts[0]
+        ind0b=self.resolution[0]-self.ghosts[0]
+        ind1a = self.ghosts[1]
+        ind1b=self.resolution[1]-self.ghosts[1]
+        ind2a = self.ghosts[2]
+        ind2b=self.resolution[2]-self.ghosts[2]
+        ind0=np.arange(ind0a,ind0b)
+        ind1=np.arange(ind1a,ind1b)
+        ind2=np.arange(ind2a,ind2b)
+
+        vorticityNoG = [np.zeros((self.resolution - 2 * self.ghosts), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dim)]
+        velocityNoG = [np.zeros((self.resolution - 2 * self.ghosts), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dim)]
+        for i in xrange (self.dim):
+            vorticityNoG[i] = self.vorticity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
+            velocityNoG[i] = self.velocity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
+        print 'shape vortiNoG, veloNoG', np.asarray(vorticityNoG).shape, np.asarray(velocityNoG).shape
+
+        # Vorticity diffusion
+        vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...] = \
+            fftw2py.solve_diffusion_3d(self.viscosity * dt, 
+                                          velocityNoG[0][...], velocityNoG[1][...], velocityNoG[2][...], 
+                                          vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...])
+
+        for i in xrange (self.dim):
+            self.vorticity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = vorticityNoG[i][...]
+        print 'shape vorti, velo', np.asarray(self.vorticity.data).shape, np.asarray(self.velocity.data).shape
+
+#        print 'shape vorti, velo, diffusion', np.asarray(self.vorticity.data).shape, np.asarray(self.velocity.data).shape
+        self.compute_time = time.time() - self.compute_time
+        self.total_time += self.compute_time
+        return self.compute_time
+#        return result
+
+    def printComputeTime(self):
+        self.timings_info[0] = "\"Diffusion calculation total\""
+        self.timings_info[1] = str(self.total_time)
+#        print "Diffusion total time : ", self.total_time
+        print "Time of the last diffusion calculation loop :", self.compute_time 
+
+    def __str__(self):
+        s = "Diffusion_d (DiscreteOperator). " + DiscreteOperator.__str__(self)
+        return s
+
+if __name__ == "__main__":
+    print __doc__
+    print "- Provided class : Diffusion_d"
+    print Diffusion_d.__doc__

HySoP/hysop/operator/poisson.py

+# -*- coding: utf-8 -*-
+"""
+@package parmepy.operator Poisson solver
+Continous Poisson solver operator (F2PY)
+"""
+from .continuous import ContinuousOperator
+from .poisson_d import Poisson_d
+
+
+class Poisson(ContinuousOperator):
+    """
+    Poisson solver operator representation using FFT
+    """
+
+    def __init__(self, velocity=None, vorticity=None):
+        """
+        Constructor.
+        Create a Poisson solver operator using FFT.
+
+        @param velocity ContinuousVectorField : velocity variable.
+        @param vorticity ContinuousVectorField : vorticity variable.
+        """
+        ContinuousOperator.__init__(self)
+        self.velocity = velocity
+        self.vorticity = vorticity
+        self.addVariable([velocity, vorticity])
+
+    def discretize(self, idVelocityD=0, idVorticityD=0, topology=None):
+        """
+        Poisson solver operator discretization method.
+        Create a discrete Poisson solver operator from given specifications.
+
+        @param idVelocityD : Index of velocity discretisation to update from vorticity.
+        @param idVelocityD : Index of vorticity.
+        @param topology : topology of the input fields
+        """
+        if self.discreteOperator is None:
+            self.discreteOperator = Poisson_d(self, idVelocityD, idVorticityD, topology)
+
+    def __str__(self):
+        """ToString method"""
+        s = " Poisson solver operator (ContinuousOperator)"
+        if self.discreteOperator is not None:
+            s += "with the following discretization:\n"
+            s += str(self.discreteOperator)
+        else:
+            s += ". Not discretised"
+        return s + "\n"
+
+if __name__ == "__main__":
+    print __doc__
+    print "- Provided class : Poisson"
+    print Poisson.__doc__

HySoP/hysop/operator/poisson_d.py

+# -*- coding: utf-8 -*-
+"""
+@package operator
+Discrete Poisson solver operator using FFT
+"""
+from ..f2py import fftw2py
+import mpi4py.MPI as MPI
+from .discrete import DiscreteOperator
+from ..constants import *
+import numpy as np
+import time
+import sys
+
+class Poisson_d(DiscreteOperator):
+    """
+    Operator representation.
+    """
+
+    def __init__(self, poisson, idVelocityD=0, idVorticityD=0, topology=None):
+        """
+        Constructor.
+        @param velocity : descretized velocity to update from vorticity.
+        @param vorticity
+        """
+        DiscreteOperator.__init__(self)
+        self.poisson = poisson
+        self.velocity = poisson.velocity.discreteField[idVelocityD]
+        self.vorticity = poisson.vorticity.discreteField[idVorticityD]
+        self.topology = topology
+        self.ghosts = topology.ghosts
+        self.resolution = topology.mesh.resolution
+        self.dim = topology.dim
+#        self.ppfft = parmepy.f2py.fftw2py
+        self.compute_time = 0.
+
+    def apply(self, t, dt):
+        """
+        Apply operator.
+        """
+        self.compute_time = time.time()
+
+        ind0a = self.ghosts[0]
+        ind0b=self.resolution[0]-self.ghosts[0]
+        ind1a = self.ghosts[1]
+        ind1b=self.resolution[1]-self.ghosts[1]
+        ind2a = self.ghosts[2]
+        ind2b=self.resolution[2]-self.ghosts[2]
+        ind0=np.arange(ind0a,ind0b)
+        ind1=np.arange(ind1a,ind1b)
+        ind2=np.arange(ind2a,ind2b)
+
+        vorticityNoG = [np.zeros((self.resolution - 2 * self.ghosts), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dim)]
+        velocityNoG = [np.zeros((self.resolution - 2 * self.ghosts), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dim)]
+        for i in xrange (self.dim):
+            vorticityNoG[i] = self.vorticity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
+            velocityNoG[i] = self.velocity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
+        print 'shape vortiNoG, veloNoG', np.asarray(vorticityNoG).shape, np.asarray(velocityNoG).shape
+
+        # Recover velocity field from vorticity field
+        velocityNoG[0][...], velocityNoG[1][...], velocityNoG[2][...] = \
+            fftw2py.solve_poisson_3d(vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...],
+                                       velocityNoG[0][...], velocityNoG[1][...], velocityNoG[2][...])
+
+        for i in xrange (self.dim):
+            self.velocity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = velocityNoG[i][...]
+        print 'shape vorti, velo', np.asarray(self.vorticity.data).shape, np.asarray(self.velocity.data).shape
+
+#        print 'shape vorti, velo, poisson', np.asarray(self.vorticity.data).shape, np.asarray(self.velocity.data).shape
+        self.compute_time = time.time() - self.compute_time
+        self.total_time += self.compute_time
+        return self.compute_time
+#        return result
+
+    def printComputeTime(self):
+        self.timings_info[0] = "\"Poisson solver total\""
+        self.timings_info[1] = str(self.total_time)
+#        print "Poisson total time : ", self.total_time
+        print "Time of the last Poisson solver loop :", self.compute_time 
+
+    def __str__(self):
+        s = "Poisson_d (DiscreteOperator). " + DiscreteOperator.__str__(self)
+        return s
+
+if __name__ == "__main__":
+    print __doc__
+    print "- Provided class : Poisson_d"
+    print Poisson_d.__doc__