Adding tests for each of the NS operators

HySoP/hysop/domain/obstacle/semi_cylinder.py

+"""
+Semi-cylinder obstacle description
+"""
+from parmepy.constants import np, PARMES_INTEGER
+from parmepy.domain.obstacle.obstacle_2D import Obstacle2D
+
+
+class SemiCylinder(Obstacle2D):
+    """
+    Discrete obstacles representation.
+    """
+
+    def __init__(self, obst, radius=0.2):
+        """
+        Creates the semi-cylinder obsctacle and y-boundaries and 
+        returns 3 arrays corresponding to the characteristic 
+        functions of three different porous media.
+
+        @param obstacle2D : Two dimensional obstacle description.
+        @param radius : Semi-cylinder radius
+        """
+        ## 2D parent obstacle
+        self.obst = obst
+        ## Radius of the semi-cylinder
+        self.radius = radius
+        ## Characteristic function (array) for y-boundaries
+        self.chiBoundary = None
+        ## Characteristic function (array) for the solid
+        self.chiSolid = None
+        ## Characteristic function (array) for the porous area
+        self.chiPorous = None
+        print 'obstacle =', self.obst.obstacleName
+
+    def setUp(self, topology):
+        """
+        Compute the characteristic functions associated 
+        to y-boundaries and semi-cylinder
+        """
+        ## Temporary arrays
+        chiBoundary_i = []
+        chiBoundary_j = []
+        chiSolid_i = []
+        chiSolid_j = []
+        chiPorous_i = []
+        chiPorous_j = []
+
+        step = topology.mesh.space_step
+        ghosts = topology.ghosts
+        local_start = topology.mesh.local_start
+        local_end = topology.mesh.local_end
+        coord_start = topology.mesh.origin + (ghosts * step)
+        coord_end = topology.mesh.end - (ghosts * step)
+        layerMin = coord_start[1] + self.obst.zlayer
+        layerMax = coord_end[1] - self.obst.zlayer
+        if not (self.obst.porousLayerThickn <= self.radius):
+            raise ValueError("Error, porous layer thickness" +
+                             "is higher than semi-cylinder radius.")
+        radiusMinuslayer = self.radius- self.obst.porousLayerThickn
+
+        print 'step, ghosts, local_start, local_end, zlayer', \
+              step, ghosts, local_start, local_end, self.obst.zlayer
+        print 'start, end, layerMin, layerMax, radiusMinuslayer', \
+              coord_start, coord_end, layerMin, layerMax, radiusMinuslayer
+
+        for j in xrange (local_start[1], local_end[1] + 1):
+            cy = coord_start[1] + j * step[1]
+            for i in xrange (local_start[0], local_end[0] + 1):
+                if (cy >= layerMax or cy <= layerMin):
+                # we are in the y-layer boundary:
+                    chiBoundary_i.append(i)
+                    chiBoundary_j.append(j)
+                else :
+                    cx = coord_start[0] + i * step[0]
+                    dist = np.sqrt((cx - self.obst.center[0]) ** 2 +
+                                   (cy - self.obst.center[1]) ** 2)
+                    if (radiusMinuslayer < dist 
+                        and dist <= self.radius + 1E-12 
+                        and cx <= self.obst.center[0]
+                        and self.obst.porousLayerThickn != 0.):
+                        # we are in the porous region of the semi-cylinder:
+                        chiPorous_i.append(i)
+                        chiPorous_j.append(j)
+                    if (dist <= radiusMinuslayer + 1E-12
+                        and cx <= self.obst.center[0]):
+                        # we are in the solid region of the semi-cylinder:
+                        chiSolid_i.append(i)
+                        chiSolid_j.append(j)
+
+        ## Characteristic function of penalized boundaries
+        chiBoundary_i = np.asarray(chiBoundary_i, dtype=PARMES_INTEGER)
+        chiBoundary_j = np.asarray(chiBoundary_j, dtype=PARMES_INTEGER)
+        self.chiBoundary = tuple([chiBoundary_i, chiBoundary_j])
+
+        ## Characteristic function of solid areas
+        chiSolid_i = np.asarray(chiSolid_i, dtype=PARMES_INTEGER)
+        chiSolid_j = np.asarray(chiSolid_j, dtype=PARMES_INTEGER)
+        self.chiSolid = tuple([chiSolid_i, chiSolid_j])
+
+        ## Characteristic function of porous areas
+        chiPorous_i = np.asarray(chiPorous_i, dtype=PARMES_INTEGER)
+        chiPorous_j = np.asarray(chiPorous_j, dtype=PARMES_INTEGER)
+        self.chiPorous = tuple([chiPorous_i, chiPorous_j])
+
+    def __str__(self):
+        """ToString method"""
+        return "SemiCylinder"
+
+if __name__ == "__main__" :
+    print "This module defines the following classe:"
+    print "SemiCylinder: ", SemiCylinder.__doc__

HySoP/hysop/operator/tests/test_Stretching.py

+# -*- coding: utf-8 -*-
+import time
+
+import parmepy as pp
+import numpy as np
+import numpy.testing as npt
+from parmepy.fields.analytical import AnalyticalField
+from parmepy.operator.stretching import Stretching
+from parmepy.problem.navier_stokes import NSProblem
+
+
+def computeVel(x, y, z):
+    amodul = np.cos(np.pi * 1. / 3.)
+    pix = np.pi * x
+    piy = np.pi * y
+    piz = np.pi * z
+    pi2x = 2. * pix
+    pi2y = 2. * piy
+    pi2z = 2. * piz
+    vx = 2. * np.sin(pix) * np.sin(pix) \
+         * np.sin(pi2y) * np.sin(pi2z) * amodul
+    vy = - np.sin(pi2x) * np.sin(piy) \
+         * np.sin(piy) * np.sin(pi2z) * amodul
+    vz = - np.sin(pi2x) * np.sin(piz) \
+         * np.sin(piz) * np.sin(pi2y) * amodul
+    return vx, vy, vz
+
+def computeVort(x, y, z):
+    amodul = np.cos(np.pi * 1. / 3.)
+    pix = np.pi * x
+    piy = np.pi * y
+    piz = np.pi * z
+    pi2x = 2. * pix
+    pi2y = 2. * piy
+    pi2z = 2. * piz
+    wx = 2.* np.pi * np.sin(pi2x) * amodul*(- np.cos(pi2y) \
+         * np.sin(piz) * np.sin(piz) \
+         + np.sin(piy) * np.sin(piy) * np.cos(pi2z))
+    wy = 2.* np.pi * np.sin(pi2y) * amodul*(2. * np.cos(pi2z) \
+         * np.sin(pix) * np.sin(pix) \
+         + np.sin(piz) * np.sin(piz) * np.cos(pi2x))
+    wz = -2.* np.pi * np.sin(pi2z) * amodul *(np.cos(pi2x) \
+         * np.sin(piy) * np.sin(piy) \
+         + np.sin(pix) * np.sin(pix) * np.cos(pi2y))
+    return wx, wy, wz
+
+def test_Stretching():
+    # Parameters
+    nb = 33
+    dim = 3
+    boxLength = [1., 1., 1.]
+    boxMin = [0., 0., 0.]
+    nbElem = [nb, nb, nb]
+
+    timeStep = 0.05
+    finalTime = timeStep
+
+    t0 = time.time()
+
+    ## Domain
+    box = pp.Box(dim, length=boxLength, origin=boxMin)
+
+    ## Fields
+    velo = AnalyticalField(domain=box, formula=computeVel, 
+                           name='Velocity', vector=True)
+    vorti = AnalyticalField(domain=box, formula=computeVort, 
+                            name='Vorticity', vector=True)
+
+    ## Operators
+    stretch = Stretching(velo, vorti,
+                         resolutions={velo: nbElem,
+                                      vorti: nbElem},
+                         method='FD_order4 RK3',
+                         propertyOp='divConservation'
+                        )
+
+    ##Problem
+    pb = NSProblem([stretch])
+
+    ## Setting solver to Problem
+    pb.setUp(finalTime, timeStep)
+
+    t1 = time.time()
+
+    ## Solve problem
+    timings = pb.solve()
+
+    tf = time.time()
+
+    print "\n"
+    print "Total time : ", tf - t0, "sec (CPU)"
+    print "Init time : ", t1 - t0, "sec (CPU)"
+    print "Solving time : ", tf - t1, "sec (CPU)"
+
+
+if __name__ == "__main__":
+    test_Stretching()

HySoP/hysop/operator/tests/test_advec_scales.py

+# -*- coding: utf-8 -*-
+import time
+
+import parmepy as pp
+import numpy as np
+import numpy.testing as npt
+from parmepy.fields.analytical import AnalyticalField
+from parmepy.operator.advection import Advection
+from parmepy.problem.navier_stokes import NSProblem
+from parmepy.operator.monitors.printer import Printer
+from math import sqrt, pi, exp
+
+
+def computeVel(x, y, z):
+    vx = 0.
+    vy = 0.
+    vz = 0.1
+    return vx, vy, vz
+
+def computeVort(x, y, z):
+    xc = 1. / 2.
+    yc = 1. / 2.
+    zc = 1. / 4.
+    R = 0.2
+    sigma = R / 2.
+    Gamma = 0.0075
+    dist = sqrt((x-xc) ** 2 + (y-yc) ** 2)
+    s2 = (z - zc) ** 2 + (dist - R) ** 2
+    wx = 0.
+    wy = 0.
+    wz = 0.
+    if (dist != 0.):
+        cosTheta = (x - xc) / dist
+        sinTheta = (y - yc) / dist
+        wTheta = Gamma / (pi * sigma ** 2) * \
+                 exp(-(s2 / sigma ** 2))
+        wx = - wTheta * sinTheta
+        wy = wTheta * cosTheta
+        wz = 0.
+    return wx, wy, wz
+
+def test_Advection_scales():
+    # Parameters
+    nb = 65
+    dim = 3
+    boxLength = [1., 1., 1.]
+    boxMin = [0., 0., 0.]
+    nbElem = [nb, nb, nb]
+
+    timeStep = 0.01
+    finalTime = 150 * timeStep
+    outputFilePrefix = './scales_'
+    outputModulo = 10
+
+    t0 = time.time()
+
+    ## Domain
+    box = pp.Box(dim, length=boxLength, origin=boxMin)
+
+    ## Fields
+    velo = AnalyticalField(domain=box, formula=computeVel, 
+                           name='Velocity', vector=True)
+    vorti = AnalyticalField(domain=box, formula=computeVort, 
+                            name='Vorticity', vector=True)
+
+    ## Advection (with scales) operator
+    advec = Advection(velo, vorti,
+                      resolutions={velo: nbElem,
+                                   vorti: nbElem},
+                      method='scales_p_02'
+                     )
+
+    ## Problem
+    pb = NSProblem([advec],
+                    monitors=[Printer(fields=[vorti],
+                              frequency=outputModulo,
+                              outputPrefix=outputFilePrefix)])
+
+    ## Setting solver to Problem
+    pb.setUp(finalTime, timeStep)
+
+    t1 = time.time()
+
+    ## Solve problem
+    timings = pb.solve()
+
+    tf = time.time()
+
+    print "\n"
+    print "Total time : ", tf - t0, "sec (CPU)"
+    print "Init time : ", t1 - t0, "sec (CPU)"
+    print "Solving time : ", tf - t1, "sec (CPU)"
+
+
+if __name__ == "__main__":
+    test_Advection_scales()

HySoP/hysop/operator/tests/test_diff_poisson.py

+# -*- coding: utf-8 -*-
+import time
+
+import parmepy as pp
+import numpy as np
+from parmepy.constants import PARMES_REAL, ORDER
+from parmepy.fields.analytical import AnalyticalField
+from parmepy.operator.poisson import Poisson
+from parmepy.operator.diffusion import Diffusion
+from parmepy.problem.navier_stokes import NSProblem
+from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
+from math import sqrt, pi, exp
+
+
+def computeVel(x, y, z):
+    vx = 0.
+    vy = 0.
+    vz = 0.
+    return vx, vy, vz
+
+def computeVort(x, y, z):
+    xc = 1. / 2.
+    yc = 1. / 2.
+    zc = 1. / 4.
+    R = 0.2
+    sigma = R / 2.
+    Gamma = 0.0075
+    dist = sqrt((x-xc) ** 2 + (y-yc) ** 2)
+    s2 = (z - zc) ** 2 + (dist - R) ** 2
+    wx = 0.
+    wy = 0.
+    wz = 0.
+    if (dist != 0.):
+        cosTheta = (x - xc) / dist
+        sinTheta = (y - yc) / dist
+        wTheta = Gamma / (pi * sigma ** 2) * \
+                 exp(-(s2 / sigma ** 2))
+        wx = - wTheta * sinTheta
+        wy = wTheta * cosTheta
+        wz = 0.
+    return wx, wy, wz
+
+def test_Diff_Poisson():
+    # Parameters
+    nb = 65
+    dim = 3
+    boxLength = [1., 1., 1.]
+    boxMin = [0., 0., 0.]
+    nbElem = [nb, nb, nb]
+
+    timeStep = 0.01
+    finalTime = 150 * timeStep
+    outputFilePrefix = './Energies'
+    outputModulo = 10
+
+    t0 = time.time()
+
+    ## Domain
+    box = pp.Box(dim, length=boxLength, origin=boxMin)
+
+    ## Fields
+    velo = AnalyticalField(domain=box, formula=computeVel, 
+                           name='Velocity', vector=True)
+    vorti = AnalyticalField(domain=box, formula=computeVort, 
+                            name='Vorticity', vector=True)
+
+    ## FFT Diffusion operators and FFT Poisson solver
+    diffusion = Diffusion(velo, vorti,
+                          resolutions={velo: nbElem,
+                                       vorti: nbElem},
+                          method='',
+                          viscosity=0.002, 
+                          with_curl=False
+                         )
+
+    poisson = Poisson(velo, vorti,
+                      resolutions={velo: nbElem,
+                                   vorti: nbElem},
+                      method='',
+                      domain=box
+                     )
+
+    ## Problem
+    pb = NSProblem([diffusion, poisson],
+                    monitors=[Energy_enstrophy(
+                                fields=[velo, vorti],
+                                frequency=outputModulo,
+                                outputPrefix=outputFilePrefix)])
+
+
+    ## Setting solver to Problem
+    pb.setUp(finalTime, timeStep)
+
+    t1 = time.time()
+
+    ## Solve problem
+    timings = pb.solve()
+
+    tf = time.time()
+
+    print "\n"
+    print "Total time : ", tf - t0, "sec (CPU)"
+    print "Init time : ", t1 - t0, "sec (CPU)"
+    print "Solving time : ", tf - t1, "sec (CPU)"
+
+
+if __name__ == "__main__":
+    test_Diff_Poisson()

HySoP/hysop/operator/tests/test_penalization_2D.py

+# -*- coding: utf-8 -*-
+import time
+import parmepy as pp
+import numpy as np
+from parmepy.constants import PARMES_REAL, ORDER
+from parmepy.domain.obstacle.obstacle_2D import Obstacle2D
+from parmepy.fields.analytical import AnalyticalField
+from parmepy.operator.penalization import Penalization
+from parmepy.problem.navier_stokes import NSProblem
+from parmepy.operator.monitors.printer import Printer
+
+
+def computeVel(x, y, ):
+    vx = 1.
+    vy = 1.
+    return vx, vy
+
+def computeVort(x, y):
+    wx = 0.
+    wy = 0.
+    return wx, wy
+
+def test_Penal2D():
+    # Parameters
+    nb = 101
+    dim = 2
+    boxLength = [1., 1.]
+    boxMin = [0., 0.]
+    nbElem = [nb, nb]
+
+    timeStep = 0.05
+    finalTime = timeStep
+    outputFilePrefix = './penal2D_'
+    outputModulo = 1
+
+    t0 = time.time()
+
+    ## Domain
+    box = pp.Box(dim, length=boxLength, origin=boxMin)
+
+    ## Obstacle
+    lambd = np.array([10.E8, 10.], 
+                     dtype=PARMES_REAL, order=ORDER)
+    cylinder = Obstacle2D(box, center=[0.5, 0.5], 
+                          zlayer=0.02, 
+                          porousLayerThickn=0.05, 
+                          porousLayerConfig='', 
+                          obstacleName='semi_cylinder', 
+                          radius=0.2)
+
+    ## Fields
+    velo = AnalyticalField(domain=box, formula=computeVel, 
+                           name='Velocity', vector=True)
+    vorti = AnalyticalField(domain=box, formula=computeVort, 
+                            name='Vorticity', vector=True)
+
+    ## Operators
+    penal = Penalization(velo, vorti,
+                         obstacle=cylinder,
+                         lambd=lambd,
+                         resolutions={velo: nbElem,
+                                      vorti: nbElem},
+                         method='',
+                         )
+
+    ##Problem
+    pb = NSProblem([penal],
+                    monitors=[Printer(fields=[velo],
+                              frequency=outputModulo,
+                              outputPrefix=outputFilePrefix)])
+
+    ## Setting solver to Problem
+    pb.setUp(finalTime, timeStep)
+
+    t1 = time.time()
+
+    ## Solve problem
+    timings = pb.solve()
+
+    tf = time.time()
+
+
+    print "\n"
+    print "Total time : ", tf - t0, "sec (CPU)"
+    print "Init time : ", t1 - t0, "sec (CPU)"
+    print "Solving time : ", tf - t1, "sec (CPU)"
+    print ""
+#    return pb.timings_info[0] + "\"Solving time\" \"Init time\" \"Total time\"", pb.timings_info[1] + str(tf - t1) + " " + str(t1 - t0) + " " + str(tf - t0)
+
+if __name__ == "__main__":
+    test_Penal2D()

HySoP/hysop/operator/tests/test_penalization_3D.py

+# -*- coding: utf-8 -*-
+import time
+import parmepy as pp
+import numpy as np
+from parmepy.constants import PARMES_REAL, ORDER
+from parmepy.domain.obstacle.obstacle_3D import Obstacle3D
+from parmepy.fields.analytical import AnalyticalField
+from parmepy.operator.penalization import Penalization
+from parmepy.problem.navier_stokes import NSProblem
+from parmepy.operator.monitors.printer import Printer
+
+
+def computeVel(x, y, z):
+    vx = 1.
+    vy = 1.
+    vz = 1.
+    return vx, vy, vz
+
+def computeVort(x, y, z):
+    wx = 0.
+    wy = 0.
+    wz = 0.
+    return wx, wy, wz
+
+def test_Penal3D():
+    ## Parameters
+    nb = 129
+    dim = 3
+    boxLength = [1., 1., 1.]
+    boxMin = [0., 0., 0.]
+    nbElem = [nb, nb, nb]
+
+    timeStep = 0.05
+    finalTime = timeStep
+    outputFilePrefix = './penal3D_'
+    outputModulo = 1
+
+    t0 = time.time()
+
+    ## Domain
+    box = pp.Box(dim, length=boxLength, origin=boxMin)
+
+    ## Obstacle
+    obstacle = Obstacle3D(box, center=[0.5, 0.5, 0.5], 
+                        zlayer=0.02, 
+                        porousLayerThickn=0.1, 
+                        porousLayerConfig='', 
+                        obstacleName='hemisphere', 
+                        radius=0.3)
+
+    ## Penalization parameter lambda
+    lambd = np.array([10.E8, 10.], 
+                     dtype=PARMES_REAL, order=ORDER)
+
+    ## Fields
+    velo = AnalyticalField(domain=box, formula=computeVel, 
+                           name='Velocity', vector=True)
+    vorti = AnalyticalField(domain=box, formula=computeVort, 
+                            name='Vorticity', vector=True)
+
+    ## Penalization operator
+    penal = Penalization(velo, vorti,
+                         obstacle=obstacle,
+                         lambd=lambd,
+                         resolutions={velo: nbElem,
+                                      vorti: nbElem},
+                         method='',
+                         )
+
+    ## Problem
+    pb = NSProblem([penal],
+                    monitors=[Printer(fields=[velo],
+                              frequency=outputModulo,
+                              outputPrefix=outputFilePrefix)])
+
+    ## Setting solver to Problem
+    pb.setUp(finalTime, timeStep)
+
+    t1 = time.time()
+
+    ## Solve problem
+    timings = pb.solve()
+
+    tf = time.time()
+
+
+    print "\n"
+    print "Total time : ", tf - t0, "sec (CPU)"
+    print "Init time : ", t1 - t0, "sec (CPU)"
+    print "Solving time : ", tf - t1, "sec (CPU)"
+    print ""
+#    return pb.timings_info[0] + "\"Solving time\" \"Init time\" \"Total time\"", pb.timings_info[1] + str(tf - t1) + " " + str(t1 - t0) + " " + str(tf - t0)
+
+if __name__ == "__main__":
+    test_Penal3D()