From fe9aff9bdb08ad02b07d67a7efc86cf81ab6eee1 Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Fri, 15 Mar 2013 09:04:39 +0000
Subject: [PATCH] Adding tests for each of the NS operators
---
HySoP/hysop/domain/obstacle/semi_cylinder.py | 110 ++++++++++++++++++
HySoP/hysop/operator/tests/test_Stretching.py | 97 +++++++++++++++
.../hysop/operator/tests/test_advec_scales.py | 96 +++++++++++++++
.../hysop/operator/tests/test_diff_poisson.py | 108 +++++++++++++++++
.../operator/tests/test_penalization_2D.py | 91 +++++++++++++++
.../operator/tests/test_penalization_3D.py | 95 +++++++++++++++
6 files changed, 597 insertions(+)
create mode 100644 HySoP/hysop/domain/obstacle/semi_cylinder.py
create mode 100755 HySoP/hysop/operator/tests/test_Stretching.py
create mode 100755 HySoP/hysop/operator/tests/test_advec_scales.py
create mode 100755 HySoP/hysop/operator/tests/test_diff_poisson.py
create mode 100644 HySoP/hysop/operator/tests/test_penalization_2D.py
create mode 100644 HySoP/hysop/operator/tests/test_penalization_3D.py
diff --git a/HySoP/hysop/domain/obstacle/semi_cylinder.py b/HySoP/hysop/domain/obstacle/semi_cylinder.py
new file mode 100644
index 000000000..c5bca7414
--- /dev/null
+++ b/HySoP/hysop/domain/obstacle/semi_cylinder.py
@@ -0,0 +1,110 @@
+"""
+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__
diff --git a/HySoP/hysop/operator/tests/test_Stretching.py b/HySoP/hysop/operator/tests/test_Stretching.py
new file mode 100755
index 000000000..62912496f
--- /dev/null
+++ b/HySoP/hysop/operator/tests/test_Stretching.py
@@ -0,0 +1,97 @@
+# -*- 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()
diff --git a/HySoP/hysop/operator/tests/test_advec_scales.py b/HySoP/hysop/operator/tests/test_advec_scales.py
new file mode 100755
index 000000000..dbd2f71e3
--- /dev/null
+++ b/HySoP/hysop/operator/tests/test_advec_scales.py
@@ -0,0 +1,96 @@
+# -*- 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()
diff --git a/HySoP/hysop/operator/tests/test_diff_poisson.py b/HySoP/hysop/operator/tests/test_diff_poisson.py
new file mode 100755
index 000000000..7ecebea58
--- /dev/null
+++ b/HySoP/hysop/operator/tests/test_diff_poisson.py
@@ -0,0 +1,108 @@
+# -*- 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()
diff --git a/HySoP/hysop/operator/tests/test_penalization_2D.py b/HySoP/hysop/operator/tests/test_penalization_2D.py
new file mode 100644
index 000000000..b5e86515a
--- /dev/null
+++ b/HySoP/hysop/operator/tests/test_penalization_2D.py
@@ -0,0 +1,91 @@
+# -*- 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()
diff --git a/HySoP/hysop/operator/tests/test_penalization_3D.py b/HySoP/hysop/operator/tests/test_penalization_3D.py
new file mode 100644
index 000000000..59159c214
--- /dev/null
+++ b/HySoP/hysop/operator/tests/test_penalization_3D.py
@@ -0,0 +1,95 @@
+# -*- 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()
--
GitLab