From af37b87fcaa2b33a6e79c27be05c5f6a752ae815 Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Fri, 15 Mar 2013 09:07:30 +0000
Subject: [PATCH] Operators and obstacles updates
---
HySoP/hysop/domain/obstacle/cylinder.py | 20 +-
HySoP/hysop/domain/obstacle/hemisphere.py | 28 +-
HySoP/hysop/domain/obstacle/obstacle_2D.py | 19 +-
HySoP/hysop/domain/obstacle/obstacle_3D.py | 20 +-
HySoP/hysop/domain/obstacle/sphere.py | 24 +-
HySoP/hysop/numerics/differentialOperator.py | 36 +-
HySoP/hysop/numerics/fct2op.py | 8 +-
HySoP/hysop/numerics/integrators/euler.py | 2 +-
HySoP/hysop/operator/diffusion.py | 24 +-
HySoP/hysop/operator/diffusion_fft.py | 132 +++---
HySoP/hysop/operator/energy_enstrophy.py | 61 ++-
.../hysop/operator/monitors/compute_forces.py | 431 +++++++++++-------
HySoP/hysop/operator/penalization.py | 1 +
HySoP/hysop/operator/penalization_d.py | 1 +
HySoP/hysop/operator/poisson.py | 18 +
HySoP/hysop/operator/poisson_fft.py | 170 +++----
HySoP/hysop/operator/scales_advection.py | 13 +-
HySoP/hysop/operator/stretching.py | 1 +
HySoP/hysop/operator/stretching_d.py | 5 +-
19 files changed, 609 insertions(+), 405 deletions(-)
diff --git a/HySoP/hysop/domain/obstacle/cylinder.py b/HySoP/hysop/domain/obstacle/cylinder.py
index 8f7779cf6..b17b7a9c9 100644
--- a/HySoP/hysop/domain/obstacle/cylinder.py
+++ b/HySoP/hysop/domain/obstacle/cylinder.py
@@ -19,7 +19,8 @@ class Cylinder(Obstacle2D):
@param obstacle2D : Two dimensional obstacle description.
@param radius : Cylinder radius
"""
- Obstacle2D.__init__(self, obst)
+ ## 2D parent obstacle
+ self.obst = obst
## Radius of the cylinder
self.radius = radius
## Characteristic function (array) for y-boundaries
@@ -28,6 +29,7 @@ class Cylinder(Obstacle2D):
self.chiSolid = None
## Characteristic function (array) for the porous area
self.chiPorous = None
+ print 'obstacle =', self.obst.obstacleName
def setUp(self, topology):
"""
@@ -48,15 +50,15 @@ class Cylinder(Obstacle2D):
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.zlayer
- layerMax = coord_end[1] - self.zlayer
- if not (self.porousLayerThickn <= self.radius):
+ 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 cylinder radius.")
- radiusMinuslayer = self.radius- self.porousLayerThickn
+ radiusMinuslayer = self.radius- self.obst.porousLayerThickn
print 'step, ghosts, local_start, local_end, zlayer', \
- step, ghosts, local_start, local_end, self.zlayer
+ step, ghosts, local_start, local_end, self.obst.zlayer
print 'start, end, layerMin, layerMax, radiusMinuslayer', \
coord_start, coord_end, layerMin, layerMax, radiusMinuslayer
@@ -69,11 +71,11 @@ class Cylinder(Obstacle2D):
chiBoundary_j.append(j)
else :
cx = coord_start[0] + i * step[0]
- dist = np.sqrt((cx - self.center[0]) ** 2 +
- (cy - self.center[1]) ** 2)
+ 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 self.porousLayerThickn != 0.):
+ and self.obst.porousLayerThickn != 0.):
# we are in the porous region of the cylinder:
chiPorous_i.append(i)
chiPorous_j.append(j)
diff --git a/HySoP/hysop/domain/obstacle/hemisphere.py b/HySoP/hysop/domain/obstacle/hemisphere.py
index dcf8b986c..1c6e4aa31 100644
--- a/HySoP/hysop/domain/obstacle/hemisphere.py
+++ b/HySoP/hysop/domain/obstacle/hemisphere.py
@@ -16,10 +16,11 @@ class Hemisphere(Obstacle3D):
returns 3 arrays corresponding to the characteristic
functions of three different porous media.
- @param obstacle3D : Three dimensional obstacle description.
+ @param obst : Three dimensional parent obstacle.
@param radius : Hemisphere radius
"""
- Obstacle3D.__init__(self, obst)
+ ## 3D parent obstacle
+ self.obst = obst
## Radius of the hemisphere
self.radius = radius
## Characteristic function (array) for z-boundaries
@@ -28,6 +29,7 @@ class Hemisphere(Obstacle3D):
self.chiSolid = None
## Characteristic function (array) for the porous area
self.chiPorous = None
+ print 'obstacle =', self.obst.obstacleName
def setUp(self, topology):
"""
@@ -51,15 +53,15 @@ class Hemisphere(Obstacle3D):
local_end = topology.mesh.local_end
coord_start = topology.mesh.origin + (ghosts * step)
coord_end = topology.mesh.end - (ghosts * step)
- layerMin = coord_start[2] + self.zlayer
- layerMax = coord_end[2] - self.zlayer
- if not (self.porousLayerThickn <= self.radius):
+ layerMin = coord_start[2] + self.obst.zlayer
+ layerMax = coord_end[2] - self.obst.zlayer
+ if not (self.obst.porousLayerThickn <= self.radius):
raise ValueError("Error, porous layer thickness" +
"is higher than hemisphere radius.")
- radiusMinuslayer = self.radius- self.porousLayerThickn
+ radiusMinuslayer = self.radius- self.obst.porousLayerThickn
print 'step, ghosts, local_start, local_end, zlayer', \
- step, ghosts, local_start, local_end, self.zlayer
+ step, ghosts, local_start, local_end, self.obst.zlayer
print 'start, end, layerMin, layerMax, radiusMinuslayer', \
coord_start, coord_end, layerMin, layerMax, radiusMinuslayer
@@ -75,19 +77,19 @@ class Hemisphere(Obstacle3D):
chiBoundary_k.append(k)
else :
cx = coord_start[0] + i * step[0]
- dist = np.sqrt((cx - self.center[0]) ** 2 +
- (cy - self.center[1]) ** 2 +
- (cz - self.center[2]) ** 2)
+ dist = np.sqrt((cx - self.obst.center[0]) ** 2 +
+ (cy - self.obst.center[1]) ** 2 +
+ (cz - self.obst.center[2]) ** 2)
if (radiusMinuslayer < dist
and dist <= self.radius + 1E-12
- and cx <= self.center[0]
- and self.porousLayerThickn != 0.):
+ and cx <= self.obst.center[0]
+ and self.obst.porousLayerThickn != 0.):
# we are in the porous region of the hemisphere:
chiPorous_i.append(i)
chiPorous_j.append(j)
chiPorous_k.append(k)
if (dist <= radiusMinuslayer + 1E-12
- and cx <= self.center[0]):
+ and cx <= self.obst.center[0]):
# we are in the solid region of the hemisphere:
chiSolid_i.append(i)
chiSolid_j.append(j)
diff --git a/HySoP/hysop/domain/obstacle/obstacle_2D.py b/HySoP/hysop/domain/obstacle/obstacle_2D.py
index dc6565e1b..0d75d475d 100644
--- a/HySoP/hysop/domain/obstacle/obstacle_2D.py
+++ b/HySoP/hysop/domain/obstacle/obstacle_2D.py
@@ -41,23 +41,22 @@ class Obstacle2D(Obstacle):
self.porousLayerConfig = porousLayerConfig
## Obstacle name
self.obstacleName = obstacleName
- print 'obstacle name', self.obstacleName
## Other cbstacle configurations
self.config = otherConfig
-# self.setUp()
def setUp(self):
if self.definedObst is None:
- if self.obstacleName.find('cylinder') >= 0:
+ if self.obstacleName.find('semi_cylinder') >= 0:
+ from parmepy.domain.obstacle.semi_cylinder import SemiCylinder
+ self.definedObst = SemiCylinder(self, **self.config)
+ elif self.obstacleName.find('cylinder') >= 0:
from parmepy.domain.obstacle.cylinder import Cylinder
self.definedObst = Cylinder(self, **self.config)
- if self.obstacleName.find('half_cylinder') >= 0:
- from parmepy.domain.obstacle.halfCylinder import HalfCylinder
- self.definedObst = HalfCylinder(self, **self.config)
-# else:
-# print "Using default 2D obstacle : cylinder"
-# from parmepy.domain.obstacle.Cylinder import Cylinder
-# self.definedObst = Cylinder(self, **self.config)
+ else:
+ print "Unknown obstacle '", self.obstacleName, \
+ "' --> Using default 2D obstacle : cylinder"
+ from parmepy.domain.obstacle.Cylinder import Cylinder
+ self.definedObst = Cylinder(self, **self.config)
if __name__ == "__main__" :
print __doc__
diff --git a/HySoP/hysop/domain/obstacle/obstacle_3D.py b/HySoP/hysop/domain/obstacle/obstacle_3D.py
index 320fa4996..dcf2c0858 100644
--- a/HySoP/hysop/domain/obstacle/obstacle_3D.py
+++ b/HySoP/hysop/domain/obstacle/obstacle_3D.py
@@ -13,7 +13,7 @@ class Obstacle3D(Obstacle):
def __init__(self, domain,
center=[0.5, 0.5, 0.5],
zlayer=0.02,
- porousLayerThickn=0.05,
+ porousLayerThickn=0.,
porousLayerConfig='',
obstacleName='',
**otherConfig):
@@ -41,23 +41,23 @@ class Obstacle3D(Obstacle):
self.porousLayerConfig = porousLayerConfig
## Obstacle name
self.obstacleName = obstacleName
- print 'obstacle name', self.obstacleName
## Other cbstacle configurations
self.config = otherConfig
-# self.setUp()
def setUp(self):
if self.definedObst is None:
- if self.obstacleName.find('sphere') >= 0:
- from parmepy.domain.obstacle.sphere import Sphere
- self.definedObst = Sphere(self, **self.config)
if self.obstacleName.find('hemisphere') >= 0:
from parmepy.domain.obstacle.hemisphere import Hemisphere
self.definedObst = Hemisphere(self, **self.config)
-# else:
-# print "Using default 3D obstacle : sphere"
-# from parmepy.domain.obstacle.sphere import Sphere
-# self.definedObst = Sphere(self, **self.config)
+ elif self.obstacleName.find('sphere') >= 0:
+ from parmepy.domain.obstacle.sphere import Sphere
+ self.definedObst = Sphere(self, **self.config)
+ else:
+ print "Unknown obstacle '", self.obstacleName, \
+ "' --> Using default 3D obstacle : sphere"
+ from parmepy.domain.obstacle.sphere import Sphere
+ self.obstacleName = 'sphere'
+ self.definedObst = Sphere(self, **self.config)
if __name__ == "__main__" :
print __doc__
diff --git a/HySoP/hysop/domain/obstacle/sphere.py b/HySoP/hysop/domain/obstacle/sphere.py
index d56728631..32a725b09 100644
--- a/HySoP/hysop/domain/obstacle/sphere.py
+++ b/HySoP/hysop/domain/obstacle/sphere.py
@@ -16,10 +16,11 @@ class Sphere(Obstacle3D):
returns 3 arrays corresponding to the characteristic
functions of three different porous media.
- @param obstacle3D : Three dimensional obstacle description.
+ @param obst : Three dimensional parent obstacle.
@param radius : Sphere radius
"""
- Obstacle3D.__init__(self, obst)
+ ## 3D parent obstacle
+ self.obst = obst
## Radius of the sphere
self.radius = radius
## Characteristic function (array) for z-boundaries
@@ -28,6 +29,7 @@ class Sphere(Obstacle3D):
self.chiSolid = None
## Characteristic function (array) for the porous area
self.chiPorous = None
+ print 'obstacle =', self.obst.obstacleName
def setUp(self, topology):
"""
@@ -51,15 +53,15 @@ class Sphere(Obstacle3D):
local_end = topology.mesh.local_end
coord_start = topology.mesh.origin + (ghosts * step)
coord_end = topology.mesh.end - (ghosts * step)
- layerMin = coord_start[2] + self.zlayer
- layerMax = coord_end[2] - self.zlayer
- if not (self.porousLayerThickn <= self.radius):
+ layerMin = coord_start[2] + self.obst.zlayer
+ layerMax = coord_end[2] - self.obst.zlayer
+ if not (self.obst.porousLayerThickn <= self.radius):
raise ValueError("Error, porous layer thickness" +
"is higher than sphere radius.")
- radiusMinuslayer = self.radius- self.porousLayerThickn
+ radiusMinuslayer = self.radius- self.obst.porousLayerThickn
print 'step, ghosts, local_start, local_end, zlayer', \
- step, ghosts, local_start, local_end, self.zlayer
+ step, ghosts, local_start, local_end, self.obst.zlayer
print 'start, end, layerMin, layerMax, radiusMinuslayer', \
coord_start, coord_end, layerMin, layerMax, radiusMinuslayer
@@ -75,12 +77,12 @@ class Sphere(Obstacle3D):
chiBoundary_k.append(k)
else :
cx = coord_start[0] + i * step[0]
- dist = np.sqrt((cx - self.center[0]) ** 2 +
- (cy - self.center[1]) ** 2 +
- (cz - self.center[2]) ** 2)
+ dist = np.sqrt((cx - self.obst.center[0]) ** 2 +
+ (cy - self.obst.center[1]) ** 2 +
+ (cz - self.obst.center[2]) ** 2)
if (radiusMinuslayer < dist
and dist <= self.radius + 1E-12
- and self.porousLayerThickn != 0.):
+ and self.obst.porousLayerThickn != 0.):
# we are in the porous region of the sphere:
chiPorous_i.append(i)
chiPorous_j.append(j)
diff --git a/HySoP/hysop/numerics/differentialOperator.py b/HySoP/hysop/numerics/differentialOperator.py
index 569364713..6227dd0f9 100755
--- a/HySoP/hysop/numerics/differentialOperator.py
+++ b/HySoP/hysop/numerics/differentialOperator.py
@@ -116,17 +116,23 @@ class DifferentialOperator(NumMethod):
else:
# X components of temp and result
- temp1 = (1.0 * self.field1[0][ind0-2,ind1a:ind1b,ind2a:ind2b] * self.field2[0][ind0-2,ind1a:ind1b,ind2a:ind2b] -
+ temp1 = self.field2[0][...] * 0.
+ temp1[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = \
+ (1.0 * self.field1[0][ind0-2,ind1a:ind1b,ind2a:ind2b] * self.field2[0][ind0-2,ind1a:ind1b,ind2a:ind2b] -
8.0 * self.field1[0][ind0-1,ind1a:ind1b,ind2a:ind2b] * self.field2[0][ind0-1,ind1a:ind1b,ind2a:ind2b] +\
8.0 * self.field1[0][ind0+1,ind1a:ind1b,ind2a:ind2b] * self.field2[0][ind0+1,ind1a:ind1b,ind2a:ind2b] -\
1.0 * self.field1[0][ind0+2,ind1a:ind1b,ind2a:ind2b] * self.field2[0][ind0+2,ind1a:ind1b,ind2a:ind2b]) / (12. * self.meshSize[0])
- temp2 = (1.0 * self.field1[1][ind0a:ind0b,ind1-2,ind2a:ind2b] * self.field2[0][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
+ temp2 = self.field2[0][...] * 0.
+ temp2[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = \
+ (1.0 * self.field1[1][ind0a:ind0b,ind1-2,ind2a:ind2b] * self.field2[0][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
8.0 * self.field1[1][ind0a:ind0b,ind1-1,ind2a:ind2b] * self.field2[0][ind0a:ind0b,ind1-1,ind2a:ind2b] +\
8.0 * self.field1[1][ind0a:ind0b,ind1+1,ind2a:ind2b] * self.field2[0][ind0a:ind0b,ind1+1,ind2a:ind2b] -\
1.0 * self.field1[1][ind0a:ind0b,ind1+2,ind2a:ind2b] * self.field2[0][ind0a:ind0b,ind1+2,ind2a:ind2b]) / (12. * self.meshSize[1])
- temp3 = (1.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2-2] * self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
+ temp3 = self.field2[0][...] * 0.
+ temp3[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = \
+ (1.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2-2] * self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
8.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2-1] * self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2-1] +\
8.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2+1] * self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2+1] -\
1.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2+2] * self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2+2]) / (12. * self.meshSize[2])
@@ -134,17 +140,23 @@ class DifferentialOperator(NumMethod):
tmp1 = np.array([temp1 + temp2 + temp3])
# Y components of temp and result
- temp1 = (1.0 * self.field1[0][ind0-2,ind1a:ind1b,ind2a:ind2b] * self.field2[1][ind0-2,ind1a:ind1b,ind2a:ind2b] -
+ temp1 = self.field2[1][...] * 0.
+ temp1[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = \
+ (1.0 * self.field1[0][ind0-2,ind1a:ind1b,ind2a:ind2b] * self.field2[1][ind0-2,ind1a:ind1b,ind2a:ind2b] -
8.0 * self.field1[0][ind0-1,ind1a:ind1b,ind2a:ind2b] * self.field2[1][ind0-1,ind1a:ind1b,ind2a:ind2b] +\
8.0 * self.field1[0][ind0+1,ind1a:ind1b,ind2a:ind2b] * self.field2[1][ind0+1,ind1a:ind1b,ind2a:ind2b] -\
1.0 * self.field1[0][ind0+2,ind1a:ind1b,ind2a:ind2b] * self.field2[1][ind0+2,ind1a:ind1b,ind2a:ind2b]) / (12. * self.meshSize[0])
- temp2 = (1.0 * self.field1[1][ind0a:ind0b,ind1-2,ind2a:ind2b] * self.field2[1][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
+ temp2 = self.field2[1][...] * 0.
+ temp2[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = \
+ (1.0 * self.field1[1][ind0a:ind0b,ind1-2,ind2a:ind2b] * self.field2[1][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
8.0 * self.field1[1][ind0a:ind0b,ind1-1,ind2a:ind2b] * self.field2[1][ind0a:ind0b,ind1-1,ind2a:ind2b] +\
8.0 * self.field1[1][ind0a:ind0b,ind1+1,ind2a:ind2b] * self.field2[1][ind0a:ind0b,ind1+1,ind2a:ind2b] -\
1.0 * self.field1[1][ind0a:ind0b,ind1+2,ind2a:ind2b] * self.field2[1][ind0a:ind0b,ind1+2,ind2a:ind2b]) / (12. * self.meshSize[1])
- temp3 = (1.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2-2] * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
+ temp3 = self.field2[1][...] * 0.
+ temp3[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = \
+ (1.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2-2] * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
8.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2-1] * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2-1] +\
8.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2+1] * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2+1] -\
1.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2+2] * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2+2]) / (12. * self.meshSize[2])
@@ -152,17 +164,23 @@ class DifferentialOperator(NumMethod):
tmp2 = np.array([temp1 + temp2 + temp3])
# Z components of temp and result
- temp1 = (1.0 * self.field1[0][ind0-2,ind1a:ind1b,ind2a:ind2b] * self.field2[2][ind0-2,ind1a:ind1b,ind2a:ind2b] -
+ temp1 = self.field2[2][...] * 0.
+ temp1[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = \
+ (1.0 * self.field1[0][ind0-2,ind1a:ind1b,ind2a:ind2b] * self.field2[2][ind0-2,ind1a:ind1b,ind2a:ind2b] -
8.0 * self.field1[0][ind0-1,ind1a:ind1b,ind2a:ind2b] * self.field2[2][ind0-1,ind1a:ind1b,ind2a:ind2b] +\
8.0 * self.field1[0][ind0+1,ind1a:ind1b,ind2a:ind2b] * self.field2[2][ind0+1,ind1a:ind1b,ind2a:ind2b] -\
1.0 * self.field1[0][ind0+2,ind1a:ind1b,ind2a:ind2b] * self.field2[2][ind0+2,ind1a:ind1b,ind2a:ind2b]) / (12. * self.meshSize[0])
- temp2 = (1.0 * self.field1[1][ind0a:ind0b,ind1-2,ind2a:ind2b] * self.field2[2][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
+ temp2 = self.field2[2][...] * 0.
+ temp2[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = \
+ (1.0 * self.field1[1][ind0a:ind0b,ind1-2,ind2a:ind2b] * self.field2[2][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
8.0 * self.field1[1][ind0a:ind0b,ind1-1,ind2a:ind2b] * self.field2[2][ind0a:ind0b,ind1-1,ind2a:ind2b] +\
8.0 * self.field1[1][ind0a:ind0b,ind1+1,ind2a:ind2b] * self.field2[2][ind0a:ind0b,ind1+1,ind2a:ind2b] -\
1.0 * self.field1[1][ind0a:ind0b,ind1+2,ind2a:ind2b] * self.field2[2][ind0a:ind0b,ind1+2,ind2a:ind2b]) / (12. * self.meshSize[1])
- temp3 = (1.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2-2] * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
+ temp3 = self.field2[2][...] * 0.
+ temp3[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = \
+ (1.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2-2] * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
8.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2-1] * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2-1] +\
8.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2+1] * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2+1] -\
1.0 * self.field1[2][ind0a:ind0b,ind1a:ind1b,ind2+2] * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2+2]) / (12. * self.meshSize[2])
diff --git a/HySoP/hysop/numerics/fct2op.py b/HySoP/hysop/numerics/fct2op.py
index f6390def0..a11dda9a0 100644
--- a/HySoP/hysop/numerics/fct2op.py
+++ b/HySoP/hysop/numerics/fct2op.py
@@ -77,13 +77,13 @@ class Fct2Op(object):
return result
if self.choice == 'hessianU' : # field2 = velocity field
- gradU, maxgersh = DifferentialOperator(self.field1, self.field2, choice='gradV', topology=self.topology).__call__()
+ gradU, maxgersh = DifferentialOperator(self.field1, self.field2, self.method, choice='gradV').__call__()
GradUline1 = np.concatenate(([gradU[0]],[gradU[1]],[gradU[2]]))
GradUline2 = np.concatenate(([gradU[3]],[gradU[4]],[gradU[5]]))
GradUline3 = np.concatenate(([gradU[6]],[gradU[7]],[gradU[8]]))
- result1, maxgersh = DifferentialOperator(self.field1, GradUline1, choice='gradV', topology=self.topology).__call__()
- result2, maxgersh = DifferentialOperator(self.field1, GradUline2, choice='gradV', topology=self.topology).__call__()
- result3, maxgersh = DifferentialOperator(self.field1, GradUline3, choice='gradV', topology=self.topology).__call__()
+ result1, maxgersh = DifferentialOperator(self.field1, GradUline1, self.method, choice='gradV').__call__()
+ result2, maxgersh = DifferentialOperator(self.field1, GradUline2, self.method, choice='gradV').__call__()
+ result3, maxgersh = DifferentialOperator(self.field1, GradUline3, self.method, choice='gradV').__call__()
return result1, result2, result3
else :
print 'choice', choice , 'in Fct2Op is not defined'
diff --git a/HySoP/hysop/numerics/integrators/euler.py b/HySoP/hysop/numerics/integrators/euler.py
index fffd7f60c..fab93edd1 100755
--- a/HySoP/hysop/numerics/integrators/euler.py
+++ b/HySoP/hysop/numerics/integrators/euler.py
@@ -28,7 +28,7 @@ class Euler(ODESolver):
ODESolver.__init__(self, f)
self.name = 'euler'
- def integrate(self, f, t, dt, y):
+ def __call__(self, f, t, dt, y):
"""
Integration step for Euler Method.
up = f[t, u(t,x,y,z)].
diff --git a/HySoP/hysop/operator/diffusion.py b/HySoP/hysop/operator/diffusion.py
index 1bd91d4c2..98aaafb33 100644
--- a/HySoP/hysop/operator/diffusion.py
+++ b/HySoP/hysop/operator/diffusion.py
@@ -4,7 +4,10 @@
Continous Diffusion operator (F2PY)
"""
from parmepy.operator.continuous import Operator
+from parmepy.f2py import fftw2py
+from parmepy.mpi.topology import Cartesian
from diffusion_fft import DiffusionFFT
+from parmepy.constants import debug
class Diffusion(Operator):
@@ -13,7 +16,7 @@ class Diffusion(Operator):
"""
@debug
- def __init__(self, velocity=None, vorticity=None, viscosity=0.1,
+ def __init__(self, velocity=None, vorticity=None,
resolutions=None, method='',
**other_config):
"""
@@ -27,7 +30,6 @@ class Diffusion(Operator):
Operator.__init__(self, [velocity, vorticity])
self.velocity = velocity
self.vorticity = vorticity
- self.viscosity = viscosity
self.resolutions = resolutions
self.method = method
self.config = other_config
@@ -37,13 +39,23 @@ class Diffusion(Operator):
"""
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.
"""
Operator.setUp(self)
+ print self.method
+ localres, localoffset = fftw2py.init_fftw_solver(
+ self.resolutions[self.vorticity],
+ self.domain.length)
+ topodims = self.resolutions[self.vorticity] / localres
+ print 'topodims DIFFUSION', topodims
+ #variables discretization
+ for v in self.variables:
+ topoId, topo = self.domain.addTopology(
+ Cartesian.withResolution(
+ domain=self.domain, topoResolution=topodims,
+ globalMeshResolution=self.resolutions[v]))
+ df, dfId = v.discretize(topo)
+ self.discreteFieldId[v] = dfId
self.discreteOperator = DiffusionFFT(self, self.method,
- self.viscosity,
**self.config)
self.discreteOperator.setUp()
diff --git a/HySoP/hysop/operator/diffusion_fft.py b/HySoP/hysop/operator/diffusion_fft.py
index d121a58e8..b13a06f65 100644
--- a/HySoP/hysop/operator/diffusion_fft.py
+++ b/HySoP/hysop/operator/diffusion_fft.py
@@ -5,7 +5,7 @@ Discrete Diffusion operator using FFT
"""
from parmepy.f2py import fftw2py
from parmepy.operator.discrete import DiscreteOperator
-from parmepy.constants import np, ORDER, PARMES_REAL
+from parmepy.constants import np, ORDER, PARMES_REAL, debug
import time
@@ -16,29 +16,30 @@ class DiffusionFFT(DiscreteOperator):
@debug
def __init__(self, diff, method='',
- viscosity=0.1):
+ viscosity=0.1, with_curl=False):
"""
Constructor.
@param velocity : descretized velocity to use.
@param vorticity : discretized vorticity to update.
"""
- DiscreteOperator.__init__(self)
+ DiscreteOperator.__init__(self, diff)
+ ## Parent Continous Op.
self.diff = diff
- self.velocity = self.diff.velocity.getDiscreteField(
- self.diff.resolutions[self.diff.velocity])
- self.vorticity = self.diff.vorticity.getDiscreteField(
- self.diff.resolutions[self.diff.vorticity])
+ ## Velocity.
+ self.velocity = self.diff.velocity.discreteField[
+ self.diff.discreteFieldId[self.diff.velocity]]
+ ## Vorticity.
+ self.vorticity = self.diff.vorticity.discreteField[
+ self.diff.discreteFieldId[self.diff.vorticity]]
+ ## Viscosity.
self.viscosity = viscosity
+ ## Boolean : pure vort diffusion or curl(u) + vort diffusion.
+ self.with_curl = with_curl
self.compute_time = 0.
@debug
- def setUp(self, method='', topology=None):
- ## La topologie se recupere depuis les variables discretes
- self.topology = self.velocity.topology
- self.topology = self.vorticity.topology
- self.ghosts = self.topology.ghosts
- self.resolution = self.topology.mesh.resolution
- self.dim = self.topology.dim
+ def setUp(self):
+ pass
@debug
def apply(self, t, dt, ite):
@@ -47,53 +48,76 @@ class DiffusionFFT(DiscreteOperator):
"""
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]
-
- 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]
-
- # Curl + Vorticity diffusion
+ if self.with_curl:
+
+ self.vorticity.data[0], \
+ self.vorticity.data[1], \
+ self.vorticity.data[2] = \
+ fftw2py.solve_curl_diffusion_3d(self.viscosity * dt,
+ self.velocity.data[0],
+ self.velocity.data[1],
+ self.velocity.data[2],
+ self.vorticity.data[0],
+ self.vorticity.data[1],
+ self.vorticity.data[2])
+
+ else:
+ self.vorticity.data[0], \
+ self.vorticity.data[1], \
+ self.vorticity.data[2] = \
+ fftw2py.solve_diffusion_3d(self.viscosity * dt,
+ self.vorticity.data[0],
+ self.vorticity.data[1],
+ self.vorticity.data[2])
+
+
+# ind0a = self.topology.mesh.local_start[0]
+# ind0b = self.topology.mesh.local_end[0] + 1
+# ind1a = self.topology.mesh.local_start[1]
+# ind1b = self.topology.mesh.local_end[1] + 1
+# ind2a = self.topology.mesh.local_start[2]
+# ind2b = self.topology.mesh.local_end[2] + 1
+
+# 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]
+
+# # Curl + Vorticity diffusion
+## vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...] = \
+## fftw2py.solve_curl_diffusion_3d(self.viscosity * dt,
+## velocityNoG[0][...],
+## velocityNoG[1][...],
+## velocityNoG[2][...],
+## vorticityNoG[0][...],
+## vorticityNoG[1][...],
+## vorticityNoG[2][...])
+
+# # Pure Vorticity diffusion
# vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...] = \
-# fftw2py.solve_curl_diffusion_3d(self.viscosity * dt,
-# velocityNoG[0][...],
-# velocityNoG[1][...],
-# velocityNoG[2][...],
-# vorticityNoG[0][...],
-# vorticityNoG[1][...],
-# vorticityNoG[2][...])
-
- # Pure Vorticity diffusion
- vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...] = \
- fftw2py.solve_diffusion_3d(self.viscosity * dt,
- vorticityNoG[0][...],
- vorticityNoG[1][...],
- vorticityNoG[2][...])
-
- for i in xrange(self.dim):
- self.vorticity[i][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = \
- vorticityNoG[i][...]
+# fftw2py.solve_diffusion_3d(self.viscosity * dt,
+# vorticityNoG[0][...],
+# vorticityNoG[1][...],
+# vorticityNoG[2][...])
+
+# for i in xrange(self.dim):
+# self.vorticity[i][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = \
+# vorticityNoG[i][...]
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] = "\"Diffusion calculation total\""
- self.timings_info[1] = str(self.total_time)
+# self.timings_info[0] = "\"Diffusion calculation total\""
+# self.timings_info[1] = str(self.total_time)
print "Time of the last diffusion calculation loop :", \
self.compute_time
diff --git a/HySoP/hysop/operator/energy_enstrophy.py b/HySoP/hysop/operator/energy_enstrophy.py
index 07342c153..2dfba8505 100644
--- a/HySoP/hysop/operator/energy_enstrophy.py
+++ b/HySoP/hysop/operator/energy_enstrophy.py
@@ -6,6 +6,7 @@ Compute Energy and Enstrophy
from parmepy.constants import np
import mpi4py.MPI as MPI
from parmepy.operator.monitors.monitoring import Monitoring
+from parmepy.mpi import main_comm, main_rank, MPI
import time
@@ -14,27 +15,28 @@ class Energy_enstrophy(Monitoring):
Compute the forces according the Noca s formula
"""
- def __init__(self, topology=None, velocity=None, vorticity=None,
+ def __init__(self, fields=[],
frequency=None, outputPrefix=None):
"""
Constructor.
- @param topology : Local topology
- @param velocity : Continous Velocity field
- @param vorticity : Continous Vorticity field
+ @param fields : List of fields
@param frequency : output file producing frequency.
@param outputPrefix : file name prefix, contains relative path.
"""
- Monitoring.__init__(self, frequency)
- self.topology = topology
- self.velocity = velocity
- self.vorticity = vorticity
+ Monitoring.__init__(self, fields, frequency)
self.outputPrefix = outputPrefix
- self.dim = topology.dim
- self.resolution = topology.mesh.resolution
- self.step = topology.mesh.size
- self.ghosts = topology.ghosts
+ self.velocity = None
+ self.vorticity = None
+ for f in fields:
+ if f.name == 'Velocity':
+ self.velocity = f
+ if f.name == 'Vorticity':
+ self.vorticity = f
+ if self.velocity == None or self.vorticity == None:
+ raise ValueError("Missing velocity or vorticity field")
+
self.bufferVelo = 0.
- if (self.topology.rank == 0):
+ if (main_rank == 0):
self.f = open(self.outputPrefix, 'w')
def apply(self, t, dt, ite):
@@ -44,16 +46,17 @@ class Energy_enstrophy(Monitoring):
"""
velo = self.velocity.discreteField[0]
vort = self.vorticity.discreteField[0]
- ind0a = self.ghosts[0]
- ind0b = self.resolution[0] - self.ghosts[0]
- if (self.dim > 1):
- ind1a = self.ghosts[1]
- ind1b = self.resolution[1] - self.ghosts[1]
- if (self.dim > 2):
- ind2a = self.ghosts[2]
- ind2b = self.resolution[2] - self.ghosts[2]
- dvol = np.prod(self.step)
+ ind0a = velo.topology.mesh.local_start[0]
+ ind0b = velo.topology.mesh.local_end[0] + 1
+ if (velo.dimension > 1):
+ ind1a = velo.topology.mesh.local_start[1]
+ ind1b = velo.topology.mesh.local_end[1] + 1
+ if (velo.dimension > 2):
+ ind2a = velo.topology.mesh.local_start[2]
+ ind2b = velo.topology.mesh.local_end[2] + 1
+
+ dvol = np.prod(velo.topology.mesh.space_step)
squareNormVel = np.sum(velo[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] ** 2 + \
velo[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] ** 2 + \
@@ -61,6 +64,17 @@ class Energy_enstrophy(Monitoring):
energy = 0.5 * squareNormVel
+ ind0a = vort.topology.mesh.local_start[0]
+ ind0b = vort.topology.mesh.local_end[0] + 1
+ if (vort.dimension > 1):
+ ind1a = vort.topology.mesh.local_start[1]
+ ind1b = vort.topology.mesh.local_end[1] + 1
+ if (vort.dimension > 2):
+ ind2a = vort.topology.mesh.local_start[2]
+ ind2b = vort.topology.mesh.local_end[2] + 1
+
+ dvol = np.prod(vort.topology.mesh.space_step)
+
enstrophy = np.sum(vort[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] ** 2 + \
vort[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] ** 2 + \
vort[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] ** 2) * dvol
@@ -69,8 +83,7 @@ class Energy_enstrophy(Monitoring):
self.bufferVelo = squareNormVel
- if ((self.topology.rank == 0) and (ite % self.frequency == 0) and (t > dt)):
- # print time and forces values in the following order : time, cX, cY, cZ
+ if ((main_rank == 0) and (ite % self.freq == 0) and (t > dt)):
self.f.write("%s %s %s %s\n" % (t, energy, enstrophy, recoveredViscosity))
def __str__(self):
diff --git a/HySoP/hysop/operator/monitors/compute_forces.py b/HySoP/hysop/operator/monitors/compute_forces.py
index 16236415d..de0ca75b0 100644
--- a/HySoP/hysop/operator/monitors/compute_forces.py
+++ b/HySoP/hysop/operator/monitors/compute_forces.py
@@ -4,10 +4,11 @@
Compute forces
"""
from parmepy.constants import np, PARMES_REAL, PARMES_INTEGER, PI
+from parmepy.mpi.topology import Cartesian
from parmepy.operator.fct2op import Fct2Op
from parmepy.operator.differentialOperator_d import DifferentialOperator_d
from parmepy.operator.monitors.monitoring import Monitoring
-from parmepy.mpi import main_comm, MPI
+from parmepy.mpi import main_comm, main_rank, MPI
import time
@@ -16,9 +17,9 @@ class Compute_forces(Monitoring):
Compute the forces according the Noca s formula
"""
- def __init__(self, velocity=None, vorticity=None, topology=None,
+ def __init__(self, velocity=None, vorticity=None,
obstacle=None, boxMin= None, boxMax=None, Reynolds=None,
- frequency=None, outputPrefix=None):
+ method='', frequency=None, outputPrefix=None):
"""
Constructor.
@param velocity : Continuous velocity field
@@ -27,34 +28,76 @@ class Compute_forces(Monitoring):
@param obstacle : Continuous obstacle
@param boxMin : Global minimum coordinates of the control volume
@param boxMax : Global maximum coordinates of the control volume
+ @param Reynolds : Reynolds number value
+ @param method : numerical method for Differential operators discretizations
@param frequency : output file producing frequency.
@param outputPrefix : file name prefix, contains relative path.
"""
Monitoring.__init__(self, frequency)
self.velocity = velocity
self.vorticity = vorticity
- self.topology = topology
self.obstacle = obstacle
self.boxMin = boxMin
self.boxMax = boxMax
self.Re = Reynolds
+ self.method = method
self.outputPrefix = outputPrefix
- self.dim = topology.dim
- self.res = topology.mesh.resolution
- self.step = topology.mesh.size
- self.coordMin = topology.mesh.origin + (topology.ghosts* self.step)
- self.coordMax = topology.mesh.end * self.step + topology.domain.origin - 2.* (topology.ghosts* self.step)
- self.ghosts = topology.ghosts
- self.periods = topology.periods
-# print 'res, dx, cMin, cMax, bmin, bMax', self.res, self.step, self.coordMin, self.coordMax, self.boxMin, self.boxMax
+
+ def setUp(self):
+ """
+ Transport operator discretization method.
+ Create an discrete Transport operator from given specifications.
+
+ """
+ nbGhosts = 0
+ if self.method.find('FD_order2') >= 0:
+ nbGhosts = 1
+ elif self.method.find('FD_order4') >= 0:
+ nbGhosts = 2
+ else:
+ nbGhosts = 2
+
+ # velocity discretization
+ topoId, topo = self.velocity.domain.addTopology(
+ Cartesian(domain=self.domain, dim=self.domain.dimension,
+ globalMeshResolution=self.resolutions[self.velocity],
+ ghosts=[nbGhosts, nbGhosts, nbGhosts]))
+ df, dfId = self.velocity.discretize(topo)
+ self.discreteFieldId[self.velocity] = dfId
+ self.velocity = self.velocity.dicreteField[dfIf]
+
+ # vorticity discretization
+ topoId, topo = self.vorticity.domain.addTopology(
+ Cartesian(domain=self.domain, dim=self.domain.dimension,
+ globalMeshResolution=self.resolutions[self.vorticity]))
+ df, dfId = self.vorticity.discretize(topo)
+ self.discreteFieldId[self.vorticity] = dfId
+ self.vorticity = self.vorticity.dicreteField[dfIf]
+
self.compute_control_box()
- if (self.topology.rank == 0):
+
+ if (main_rank == 0):
self.f = open(self.outputPrefix, 'w')
def compute_control_box(self) :
"""
- Compute indicator functions for the control box (including the obstacle)
+ Compute indicator functions for the control box
+ (including the obstacle)
+
+ The computation of index is based on velocity topology
"""
+ self.dim = self.velocity.topology.dim
+ self.res = self.velocity.topology.mesh.resolution
+ self.step = self.velocity.topology.mesh.space_step
+ self.ghosts = self.velocity.topology.ghosts
+ self.periods = self.velocity.topology.periods
+ self.coord_start = self.velocity.topology.mesh.origin + \
+ (self.ghosts * self.step)
+ self.coord_end = self.velocity.topology.mesh.end - \
+ (self.ghosts * self.step)
+
+# print 'res, dx, cMin, cMax, bmin, bMax', self.res, self.step, self.coord_start, self.coord_end, self.boxMin, self.boxMax
+
# normal vectors definitions
self.normal = np.zeros([self.dim*2, self.dim], dtype=PARMES_REAL)
self.Up = 0
@@ -73,18 +116,18 @@ class Compute_forces(Monitoring):
# control box definition
tmp_ind = np.zeros([self.dim, 2], dtype=PARMES_INTEGER)
- distMin = self.boxMin - self.coordMin
- distMax = self.boxMax - self.coordMax
+ distMin = self.boxMin - self.coord_start
+ distMax = self.boxMax - self.coord_end
for i in xrange(self.dim):
if (distMin[i]>=0. and distMax[i]<=0.):
# the control volume is included inside the local domain
tmp_ind[i][0] = distMin[i] // self.step[i] # Down, West, South
tmp_ind[i][1] = self.res[i] -2. * self.ghosts[i] + distMax[i] // self.step[i] -1# Up, East, North
- if (distMin[i]>=0. and self.boxMin[i]<= self.coordMax[i] and distMax[i]>=0.):
+ if (distMin[i]>=0. and self.boxMin[i]<= self.coord_end[i] and distMax[i]>=0.):
# only min corner of control volume is included inside the local domain
tmp_ind[i][0] = distMin[i] // self.step[i]
tmp_ind[i][1] = self.res[i] -2. * self.ghosts[i] - 1
- if (distMin[i]<=0. and self.boxMax[i]>= self.coordMin[i] and distMax[i]<=0.) :
+ if (distMin[i]<=0. and self.boxMax[i]>= self.coord_start[i] and distMax[i]<=0.) :
# only max corner of control volume is included inside the local domain
tmp_ind[i][1] = self.res[i] -2.*self.ghosts[i] + distMax[i] // self.step[i] -1
if (distMin[i]<=0. and distMax[i]>=0.):
@@ -93,7 +136,7 @@ class Compute_forces(Monitoring):
# otherwise, there is no volume control inside the local domain
- ind = np.zeros(self.dim*2, dtype=PARMES_INTEGER)
+ ind = np.zeros(self.dim * 2, dtype=PARMES_INTEGER)
ind[self.Down] = tmp_ind[0][0]
ind[self.Up] = tmp_ind[0][1]
ind[self.West] = tmp_ind[1][0]
@@ -119,17 +162,17 @@ class Compute_forces(Monitoring):
# print 'nb points plan Down', nbPoints[self.Down]
# if there is no volume control inside the local domain => nbPoints=0 => forces are not computed
- if (distMin[0]<0. and distMax[0]<0. and self.boxMax[0]< self.coordMin[0]):
+ if (distMin[0]<0. and distMax[0]<0. and self.boxMax[0]< self.coord_start[0]):
nbPoints[self.Down] = 0
- if (distMin[0]>0. and distMax[0]>0. and self.boxMin[0]> self.coordMax[0]):
+ if (distMin[0]>0. and distMax[0]>0. and self.boxMin[0]> self.coord_end[0]):
nbPoints[self.Up] = 0
- if (distMin[1]<0. and distMax[1]<0. and self.boxMax[1]< self.coordMin[1]):
+ if (distMin[1]<0. and distMax[1]<0. and self.boxMax[1]< self.coord_start[1]):
nbPoints[self.West] = 0
- if (distMin[1]>0. and distMax[1]>0. and self.boxMin[1]> self.coordMax[1]):
+ if (distMin[1]>0. and distMax[1]>0. and self.boxMin[1]> self.coord_end[1]):
nbPoints[self.East] = 0
- if (distMin[2]<0. and distMax[2]<0. and self.boxMax[2]< self.coordMin[2]):
+ if (distMin[2]<0. and distMax[2]<0. and self.boxMax[2]< self.coord_start[2]):
nbPoints[self.South] = 0
- if (distMin[2]>0. and distMax[2]>0. and self.boxMin[2]> self.coordMax[2]):
+ if (distMin[2]>0. and distMax[2]>0. and self.boxMin[2]> self.coord_end[2]):
nbPoints[self.North] = 0
boxDim = np.zeros(self.dim, dtype=PARMES_INTEGER)
@@ -140,76 +183,142 @@ class Compute_forces(Monitoring):
# print 'nb points box', nbPointsBox
- tmp_box = np.zeros([nbPointsBox,self.dim], dtype=PARMES_INTEGER)
+# tmp_box = np.zeros([nbPointsBox,self.dim], dtype=PARMES_INTEGER)
+ chi_box_i=[]
+ chi_box_j=[]
+ chi_box_k=[]
coords = np.zeros(self.dim, dtype=PARMES_REAL)
count_box = 0
if(boxDim.all() >0):
for k in xrange (ind[self.South],ind[self.North]+1): ## enlever les boucles !!
- coords[2] = self.coordMin[2] + k * self.step[2]
+ coords[2] = self.coord_start[2] + k * self.step[2]
for j in xrange (ind[self.West],ind[self.East]+1):
- coords[1] = self.coordMin[1] + j * self.step[1]
+ coords[1] = self.coord_start[1] + j * self.step[1]
for i in xrange (ind[self.Down],ind[self.Up]+1):
- coords[0] = self.coordMin[0] + i * self.step[0]
+ coords[0] = self.coord_start[0] + i * self.step[0]
dist = np.vdot(coords - self.obstacle.center,coords - self.obstacle.center) - self.obstacle.radius ** 2
if(dist >=0.0): # We are on or outside the sphere ...
- tmp_box[count_box] = [i, j, k]
- count_box = count_box + 1
+# tmp_box[count_box] = [i, j, k]
+# count_box = count_box + 1
+ chi_box_i.append(i)
+ chi_box_j.append(j)
+ chi_box_k.append(k)
# print 'tmp_box', tmp_box, tmp_box.shape, count_box
# local indices of points in the control box which are outside the obstacle
- self.chi_box = np.zeros([count_box,self.dim], dtype=PARMES_INTEGER)
- self.chi_box = tmp_box[0:count_box,:] ### PBM !!!
+# self.chi_box = np.zeros([count_box,self.dim], dtype=PARMES_INTEGER)
+# self.chi_box = tmp_box[0:count_box,:] ### PBM !!!
+
+ chi_box_i = np.asarray(chi_box_i, dtype=PARMES_INTEGER)
+ chi_box_j = np.asarray(chi_box_j, dtype=PARMES_INTEGER)
+ chi_box_k = np.asarray(chi_box_k, dtype=PARMES_INTEGER)
+ self.chi_box = tuple([chi_box_i, chi_box_j, chi_box_k])
- self.chi_up = np.zeros([nbPoints[self.Up],self.dim], dtype=PARMES_INTEGER)
- self.chi_down = np.zeros([nbPoints[self.Down],self.dim], dtype=PARMES_INTEGER)
- self.chi_east = np.zeros([nbPoints[self.East],self.dim], dtype=PARMES_INTEGER)
- self.chi_west = np.zeros([nbPoints[self.West],self.dim], dtype=PARMES_INTEGER)
- self.chi_north = np.zeros([nbPoints[self.North],self.dim], dtype=PARMES_INTEGER)
- self.chi_south = np.zeros([nbPoints[self.South],self.dim], dtype=PARMES_INTEGER)
+# self.chi_up = np.zeros([nbPoints[self.Up],self.dim], dtype=PARMES_INTEGER)
+# self.chi_down = np.zeros([nbPoints[self.Down],self.dim], dtype=PARMES_INTEGER)
+# self.chi_east = np.zeros([nbPoints[self.East],self.dim], dtype=PARMES_INTEGER)
+# self.chi_west = np.zeros([nbPoints[self.West],self.dim], dtype=PARMES_INTEGER)
+# self.chi_north = np.zeros([nbPoints[self.North],self.dim], dtype=PARMES_INTEGER)
+# self.chi_south = np.zeros([nbPoints[self.South],self.dim], dtype=PARMES_INTEGER)
# local indices of points located on the surfaces of the control box
- if(nbPoints[self.South]!=0.): # South boundary is in the domain
- count = 0
+# if(nbPoints[self.South]!=0.): # South boundary is in the domain
+# count = 0
+# for j in xrange (ind[self.West],ind[self.East] + 1):
+# for i in xrange (ind[self.Down],ind[self.Up] + 1):
+# self.chi_south[count] = [i, j, ind[self.South]]
+# count = count + 1
+
+ # South boundary is in the domain
+ chi_south_i=[]
+ chi_south_j=[]
+ chi_south_k=[]
+ if(nbPoints[self.South]!=0.):
for j in xrange (ind[self.West],ind[self.East] + 1):
for i in xrange (ind[self.Down],ind[self.Up] + 1):
- self.chi_south[count] = [i, j, ind[self.South]]
- count = count + 1
-
- if(nbPoints[self.East]!=0.): # East boundary is in the domain
- count = 0
+ chi_south_i.append(i)
+ chi_south_j.append(j)
+ chi_south_k.append(ind[self.South])
+ chi_south_i = np.asarray(chi_south_i, dtype=PARMES_INTEGER)
+ chi_south_j = np.asarray(chi_south_j, dtype=PARMES_INTEGER)
+ chi_south_k = np.asarray(chi_south_k, dtype=PARMES_INTEGER)
+ self.chi_south = tuple([chi_south_i, chi_south_j, chi_south_k])
+
+ # East boundary is in the domain
+ chi_east_i=[]
+ chi_east_j=[]
+ chi_east_k=[]
+ if(nbPoints[self.East]!=0.):
for k in xrange (ind[self.South],ind[self.North] + 1):
for i in xrange (ind[self.Down],ind[self.Up] + 1):
- self.chi_east[count] = [i, ind[self.East], k]
- count = count + 1
-
- if(nbPoints[self.Down]!=0.): # Down boundary is in the domain
- count = 0
+ chi_east_i.append(i)
+ chi_east_j.append(ind[self.East])
+ chi_east_k.append(k)
+ chi_east_i = np.asarray(chi_east_i, dtype=PARMES_INTEGER)
+ chi_east_j = np.asarray(chi_east_j, dtype=PARMES_INTEGER)
+ chi_east_k = np.asarray(chi_east_k, dtype=PARMES_INTEGER)
+ self.chi_east = tuple([chi_east_i, chi_east_j, chi_east_k])
+
+ # Down boundary is in the domain
+ chi_down_i=[]
+ chi_down_j=[]
+ chi_down_k=[]
+ if(nbPoints[self.Down]!=0.):
for k in xrange (ind[self.South],ind[self.North] + 1):
for j in xrange (ind[self.West],ind[self.East] + 1):
- self.chi_down[count] = [ind[self.Down], j, k]
- count = count + 1
-
- if(nbPoints[self.North]!=0.): # North boundary is in the domain
- count = 0
+ chi_down_i.append(ind[self.Down])
+ chi_down_j.append(j)
+ chi_down_k.append(k)
+ chi_down_i = np.asarray(chi_down_i, dtype=PARMES_INTEGER)
+ chi_down_j = np.asarray(chi_down_j, dtype=PARMES_INTEGER)
+ chi_down_k = np.asarray(chi_down_k, dtype=PARMES_INTEGER)
+ self.chi_east = tuple([chi_down_i, chi_down_j, chi_down_k])
+
+ # North boundary is in the domain
+ chi_north_i=[]
+ chi_north_j=[]
+ chi_north_k=[]
+ if(nbPoints[self.North]!=0.):
for j in xrange (ind[self.West],ind[self.East] + 1):
for i in xrange (ind[self.Down],ind[self.Up] + 1):
- self.chi_north[count] = [i, j, ind[self.North]]
- count = count + 1
-
- if(nbPoints[self.West]!=0.): # West boundary is in the domain
- count = 0
+ chi_north_i.append(i)
+ chi_north_j.append(j)
+ chi_north_k.append(ind[self.North])
+ chi_north_i = np.asarray(chi_north_i, dtype=PARMES_INTEGER)
+ chi_north_j = np.asarray(chi_north_j, dtype=PARMES_INTEGER)
+ chi_north_k = np.asarray(chi_north_k, dtype=PARMES_INTEGER)
+ self.chi_north = tuple([chi_north_i, chi_north_j, chi_north_k])
+
+ # West boundary is in the domain
+ chi_west_i=[]
+ chi_west_j=[]
+ chi_west_k=[]
+ if(nbPoints[self.West]!=0.):
for k in xrange (ind[self.South],ind[self.North] + 1):
for i in xrange (ind[self.Down],ind[self.Up] + 1):
- self.chi_west[count] = [i, ind[self.West], k]
- count = count + 1
-
- if(nbPoints[self.Up]!=0.): # Up boundary is in the domain
- count = 0
+ chi_west_i.append(i)
+ chi_west_j.append(ind[self.West])
+ chi_west_k.append(k)
+ chi_west_i = np.asarray(chi_west_i, dtype=PARMES_INTEGER)
+ chi_west_j = np.asarray(chi_west_j, dtype=PARMES_INTEGER)
+ chi_west_k = np.asarray(chi_west_k, dtype=PARMES_INTEGER)
+ self.chi_west = tuple([chi_west_i, chi_west_j, chi_west_k])
+
+ # Up boundary is in the domain
+ chi_up_i=[]
+ chi_up_j=[]
+ chi_up_k=[]
+ if(nbPoints[self.Up]!=0.):
for k in xrange (ind[self.South],ind[self.North] + 1):
for j in xrange (ind[self.West],ind[self.East] + 1):
- self.chi_up[count] = [ind[self.Up], j, k]
- count = count + 1
+ chi_up_i.append(ind[self.Up])
+ chi_up_j.append(j)
+ chi_up_k.append(k)
+ chi_up_i = np.asarray(chi_up_i, dtype=PARMES_INTEGER)
+ chi_up_j = np.asarray(chi_up_j, dtype=PARMES_INTEGER)
+ chi_up_k = np.asarray(chi_up_k, dtype=PARMES_INTEGER)
+ self.chi_up = tuple([chi_up_i, chi_up_j, chi_up_k])
self.bufferForce = 0.
# Compute coef used to determine drag coefficient : cD = coef.Fx = 2.Fx/rho.uinf**2.D
@@ -221,13 +330,13 @@ class Compute_forces(Monitoring):
Computation of the drag according to the "impulsion" formula presented by
Noca et. al in Journal of Fluids and Structures, 1999
"""
- velo = self.velocity.discreteField[0]
- vort = self.vorticity.discreteField[0]
+ velo = self.velocity
+ vort = self.vorticity
localForce = np.zeros(self.dim, dtype=PARMES_REAL)
force = np.zeros(self.dim, dtype=PARMES_REAL)
# computation of hessian matrices (d2_u1, d2_u2 and d2_u3) and jacobian matrix (d_u)
- hessian1, hessian2, hessian3 = Fct2Op(vort, velo, choice = 'hessianU', topology=self.topology).apply(None,None)
- self.jacobian, maxgersh = DifferentialOperator_d(vort, velo, choice='gradV', topology=self.topology).apply()
+ hessian1, hessian2, hessian3 = Fct2Op(velo, velo, self.method, choice = 'hessianU').apply(None,None)
+ self.jacobian, maxgersh = DifferentialOperator(velo, velo, self.method, choice='gradV').apply()
# computation of the components of Div(T) where T = 1/Re * (Nabla(U) + Nabla(U)T)
self.a = 1. / self.Re * (2. * hessian1[0] + hessian1[4] + hessian2[1] + hessian1[8] + hessian3[2])
self.b = 1. / self.Re * (hessian2[0] + hessian1[1] + 2. * hessian2[4] + hessian2[8] + hessian3[5])
@@ -253,26 +362,23 @@ class Compute_forces(Monitoring):
comm = main_comm
comm.Reduce([localForce, MPI.DOUBLE], [force, MPI.DOUBLE], op=MPI.SUM, root=0)
- if ((self.topology.rank == 0) and (ite % self.frequency == 0) and (t > dt)):
# print time and forces values in the following order : time, cX, cY, cZ
+ if ((main_rank == 0) and (ite % self.frequency == 0) and (t > dt)):
self.f.write("%s %s %s %s\n" % (t, force[0], force[1], force[2]))
- def integrateOnBox(self, force, vort, chi, dvol, dt):
+ def integrateOnBox(self, force, vort, chi, dvol, dt):
+ ## ATTENTION : pour l'instant ne marche que si les topo de vort et
+ ## velo ont la meme global mesh resolution !!!
"""
Return -1/(dim-1)*d/dt int_over_control_box coord X vorticity
"""
fact = - dvol / ((self.dim - 1) * dt)
int1 = np.zeros(self.dim, dtype=PARMES_REAL)
# For all points in the box
- for ind in xrange (chi.shape[1]):
- # adjust indices for vector fields with ghost points
- i = chi[ind][0] + self.ghosts[0]
- j = chi[ind][1] + self.ghosts[1]
- k = chi[ind][2] + self.ghosts[2]
- # coordinates of the current point
- coords = self.coordMin + chi[ind] * self.step
- # part1 of the force
- int1 = int1 + np.vdot(coords,vort[i,j,k])
+ # coordinates of the current point
+ coords = self.coord_start + chi * self.step
+ # part1 of the force
+ int1 = int1 + np.vdot(coords, vort[chi[0], chi[1], chi[2]])
force = force + fact * (int1 - self.bufferForce) # 1st order time integration
self.bufferForce = int1 # Save for next time step ...
@@ -280,6 +386,8 @@ class Compute_forces(Monitoring):
return force
def integrateOnSurface(self, force, velo, vort, chi, NormalVec, Re, dsurf):
+ ## ATTENTION : pour l'instant ne marche que si les topo de vort et
+ ## velo ont la meme global mesh resolution !!!
"""
Compute integrals on surface to calculate forces acting on the body.
See (2.1) of Noca 1999 or (52) of Plouhmans 2002
@@ -292,76 +400,93 @@ class Compute_forces(Monitoring):
X_n_DivT = np.zeros(self.dim, dtype=PARMES_REAL)
n_T = np.zeros(self.dim, dtype=PARMES_REAL)
# For each point of the current plane
- for ind in xrange (chi.shape[0]):
- # adjust indices for vector fields with ghost points
- i = chi[ind][0] + self.ghosts[0]
- j = chi[ind][1] + self.ghosts[1]
- k = chi[ind][2] + self.ghosts[2]
- # part1 = 1/2(velocity.velocity)n - (n.velocity)velocity - 1/(dim-1)((n.velocity)(coord X vorticity) + 1/(dim-1)(n.vorticity)(coord X velocity)
-
- # (velocity.velocity)
- u_u = np.vdot(velo[i,j,k],velo[i,j,k])
- # normal.velocity
- n_u = np.vdot(velo[i,j,k],NormalVec[:])
- # normal.vorticity
- n_w = np.vdot(vort[i,j,k],NormalVec[:])
- # coordinates of the current point
- coords = self.coordMin[:] + chi[ind] * self.step[:]
- # part1 of the force
- int1 = int1 + 0.5 * u_u * NormalVec[:] - n_u * velo[i,j,k] \
- - fact * n_u * np.cross(coords,vort[i,j,k]) \
- + fact * n_w * np.cross(coords,velo[i,j,k])
-
- # part 2 of the force : part 2 = 1/(dim-1)(X^(n^Div(T)) + n.T
- if (chi.all() == self.chi_down.all()): # normalVec = (-1,0,0)
- X_n_DivT[0] = - coords[1] * self.b[i,j,k] - coords[2] * self.c[i,j,k]
- X_n_DivT[1] = coords[0] * self.b[i,j,k]
- X_n_DivT[2] = coords[0] * self.c[i,j,k]
- n_T[0] = - 1. / self.Re * 2 * self.jacobian[0][i,j,k]
- n_T[1] = - 1. / self.Re * (self.jacobian[3][i,j,k] + self.jacobian[1][i,j,k])
- n_T[2] = - 1. / self.Re * (self.jacobian[6][i,j,k] + self.jacobian[2][i,j,k])
-
- if (chi.all() == self.chi_up.all()): # normalVec = (1,0,0)
- X_n_DivT[0] = coords[1] * self.b[i,j,k] + coords[2] * self.c[i,j,k]
- X_n_DivT[1] = - coords[0] * self.b[i,j,k]
- X_n_DivT[2] = - coords[0] * self.c[i,j,k]
- n_T[0] = 1. / self.Re * 2 * self.jacobian[0][i,j,k]
- n_T[1] = 1. / self.Re * (self.jacobian[3][i,j,k] + self.jacobian[1][i,j,k])
- n_T[2] = 1. / self.Re * (self.jacobian[6][i,j,k] + self.jacobian[2][i,j,k])
-
- if (chi.all() == self.chi_west.all()): # normalVec = (0,-1,0)
- X_n_DivT[0] = coords[1] * self.a[i,j,k]
- X_n_DivT[1] = - coords[0] * self.a[i,j,k] - coords[2] * self.c[i,j,k]
- X_n_DivT[2] = coords[1] * self.c[i,j,k]
- n_T[0] = - 1. / self.Re * (self.jacobian[1][i,j,k] + self.jacobian[3][i,j,k])
- n_T[1] = - 1. / self.Re * 2 * self.jacobian[4][i,j,k]
- n_T[2] = - 1. / self.Re * (self.jacobian[7][i,j,k] + self.jacobian[5][i,j,k])
-
- if (chi.all() == self.chi_east.all()): # normalVect = (0,1,0)
- X_n_DivT[0] = - coords[1] * self.a[i,j,k]
- X_n_DivT[1] = coords[0] * self.a[i,j,k] + coords[2] * self.c[i,j,k]
- X_n_DivT[2] = - coords[1] * self.c[i,j,k]
- n_T[0] = 1. / self.Re * (self.jacobian[1][i,j,k] + self.jacobian[3][i,j,k])
- n_T[1] = 1. / self.Re * 2 * self.jacobian[4][i,j,k]
- n_T[2] = 1. / self.Re * (self.jacobian[7][i,j,k] + self.jacobian[5][i,j,k])
-
- if (chi.all() == self.chi_south.all()): # normalVec = (0,0,-1)
- X_n_DivT[0] = coords[2] * self.a[i,j,k]
- X_n_DivT[1] = coords[2] * self.b[i,j,k]
- X_n_DivT[2] = - coords[0] * self.a[i,j,k] - coords[1] * self.b[i,j,k]
- n_T[0] = - 1. / self.Re * (self.jacobian[2][i,j,k] + self.jacobian[6][i,j,k])
- n_T[1] = - 1. / self.Re * (self.jacobian[5][i,j,k] + self.jacobian[7][i,j,k])
- n_T[2] = - 1. / self.Re * 2 * self.jacobian[8][i,j,k]
-
- if (chi.all() == self.chi_north.all()): # normalVec = (0,0,1)
- X_n_DivT[0] = - coords[2] * self.a[i,j,k]
- X_n_DivT[1] = - coords[2] * self.b[i,j,k]
- X_n_DivT[2] = coords[0] * self.a[i,j,k] + coords[1] * self.b[i,j,k]
- n_T[0] = 1. / self.Re * (self.jacobian[2][i,j,k] + self.jacobian[6][i,j,k])
- n_T[1] = 1. / self.Re * (self.jacobian[5][i,j,k] + self.jacobian[7][i,j,k])
- n_T[2] = 1. / self.Re * 2 * self.jacobian[8][i,j,k]
-
- int2 = int2 + fact * X_n_DivT + n_T
+ # part1 = 1/2(velocity.velocity)n - (n.velocity)velocity - 1/(dim-1)((n.velocity)(coord X vorticity) + 1/(dim-1)(n.vorticity)(coord X velocity)
+
+ # (velocity.velocity)
+ u_u = np.vdot(velo[chi[0], chi[1], chi[2]],
+ velo[chi[0], chi[1], chi[2]])
+ # normal.velocity
+ n_u = np.vdot(velo[chi[0], chi[1], chi[2]],
+ NormalVec[:])
+ # normal.vorticity
+ n_w = np.vdot(vort[chi[0], chi[1], chi[2]],
+ NormalVec[:])
+ # coordinates of the current point
+ coords = self.coord_start[:] + chi * self.step[:]
+ # part1 of the force
+ int1 = int1 + 0.5 * u_u * NormalVec[:] \
+ - n_u * velo[chi[0], chi[1], chi[2]] \
+ - fact * n_u * np.cross(coords, vort[chi[0], chi[1], chi[2]]) \
+ + fact * n_w * np.cross(coords, velo[chi[0], chi[1], chi[2]])
+
+ # part 2 of the force : part 2 = 1/(dim-1)(X^(n^Div(T)) + n.T
+ if (chi.all() == self.chi_down.all()): # normalVec = (-1,0,0)
+ X_n_DivT[0] = - coords[1] * self.b[chi[0], chi[1], chi[2]] \
+ - coords[2] * self.c[chi[0], chi[1], chi[2]]
+ X_n_DivT[1] = coords[0] * self.b[chi[0], chi[1], chi[2]]
+ X_n_DivT[2] = coords[0] * self.c[chi[0], chi[1], chi[2]]
+ n_T[0] = - 1. / self.Re * 2 * self.jacobian[0][chi[0], chi[1], chi[2]]
+ n_T[1] = - 1. / self.Re * (self.jacobian[3][chi[0], chi[1], chi[2]] \
+ + self.jacobian[1][chi[0], chi[1], chi[2]])
+ n_T[2] = - 1. / self.Re * (self.jacobian[6][chi[0], chi[1], chi[2]] \
+ + self.jacobian[2][chi[0], chi[1], chi[2]])
+
+ if (chi.all() == self.chi_up.all()): # normalVec = (1,0,0)
+ X_n_DivT[0] = coords[1] * self.b[chi[0], chi[1], chi[2]] \
+ + coords[2] * self.c[chi[0], chi[1], chi[2]]
+ X_n_DivT[1] = - coords[0] * self.b[chi[0], chi[1], chi[2]]
+ X_n_DivT[2] = - coords[0] * self.c[chi[0], chi[1], chi[2]]
+ n_T[0] = 1. / self.Re * 2 * self.jacobian[0][chi[0], chi[1], chi[2]]
+ n_T[1] = 1. / self.Re * (self.jacobian[3][chi[0], chi[1], chi[2]] \
+ + self.jacobian[1][chi[0], chi[1], chi[2]])
+ n_T[2] = 1. / self.Re * (self.jacobian[6][chi[0], chi[1], chi[2]] \
+ + self.jacobian[2][chi[0], chi[1], chi[2]])
+
+ if (chi.all() == self.chi_west.all()): # normalVec = (0,-1,0)
+ X_n_DivT[0] = coords[1] * self.a[chi[0], chi[1], chi[2]]
+ X_n_DivT[1] = - coords[0] * self.a[chi[0], chi[1], chi[2]] \
+ - coords[2] * self.c[chi[0], chi[1], chi[2]]
+ X_n_DivT[2] = coords[1] * self.c[chi[0], chi[1], chi[2]]
+ n_T[0] = - 1. / self.Re * (self.jacobian[1][chi[0], chi[1], chi[2]] \
+ + self.jacobian[3][chi[0], chi[1], chi[2]])
+ n_T[1] = - 1. / self.Re * 2 * self.jacobian[4][chi[0], chi[1], chi[2]]
+ n_T[2] = - 1. / self.Re * (self.jacobian[7][chi[0], chi[1], chi[2]] \
+ + self.jacobian[5][chi[0], chi[1], chi[2]])
+
+ if (chi.all() == self.chi_east.all()): # normalVect = (0,1,0)
+ X_n_DivT[0] = - coords[1] * self.a[chi[0], chi[1], chi[2]]
+ X_n_DivT[1] = coords[0] * self.a[chi[0], chi[1], chi[2]] \
+ + coords[2] * self.c[chi[0], chi[1], chi[2]]
+ X_n_DivT[2] = - coords[1] * self.c[chi[0], chi[1], chi[2]]
+ n_T[0] = 1. / self.Re * (self.jacobian[1][chi[0], chi[1], chi[2]] \
+ + self.jacobian[3][chi[0], chi[1], chi[2]])
+ n_T[1] = 1. / self.Re * 2 * self.jacobian[4][chi[0], chi[1], chi[2]]
+ n_T[2] = 1. / self.Re * (self.jacobian[7][chi[0], chi[1], chi[2]] \
+ + self.jacobian[5][chi[0], chi[1], chi[2]])
+
+ if (chi.all() == self.chi_south.all()): # normalVec = (0,0,-1)
+ X_n_DivT[0] = coords[2] * self.a[chi[0], chi[1], chi[2]]
+ X_n_DivT[1] = coords[2] * self.b[chi[0], chi[1], chi[2]]
+ X_n_DivT[2] = - coords[0] * self.a[chi[0], chi[1], chi[2]] \
+ - coords[1] * self.b[chi[0], chi[1], chi[2]]
+ n_T[0] = - 1. / self.Re * (self.jacobian[2][chi[0], chi[1], chi[2]] \
+ + self.jacobian[6][chi[0], chi[1], chi[2]])
+ n_T[1] = - 1. / self.Re * (self.jacobian[5][chi[0], chi[1], chi[2]] \
+ + self.jacobian[7][chi[0], chi[1], chi[2]])
+ n_T[2] = - 1. / self.Re * 2 * self.jacobian[8][chi[0], chi[1], chi[2]]
+
+ if (chi.all() == self.chi_north.all()): # normalVec = (0,0,1)
+ X_n_DivT[0] = - coords[2] * self.a[chi[0], chi[1], chi[2]]
+ X_n_DivT[1] = - coords[2] * self.b[chi[0], chi[1], chi[2]]
+ X_n_DivT[2] = coords[0] * self.a[chi[0], chi[1], chi[2]] \
+ + coords[1] * self.b[chi[0], chi[1], chi[2]]
+ n_T[0] = 1. / self.Re * (self.jacobian[2][chi[0], chi[1], chi[2]] \
+ + self.jacobian[6][chi[0], chi[1], chi[2]])
+ n_T[1] = 1. / self.Re * (self.jacobian[5][chi[0], chi[1], chi[2]] \
+ + self.jacobian[7][chi[0], chi[1], chi[2]])
+ n_T[2] = 1. / self.Re * 2 * self.jacobian[8][chi[0], chi[1], chi[2]]
+
+ int2 = int2 + fact * X_n_DivT + n_T
# Product with element of surface and sum to the total (input) force
result = force + (int1 + int2) * dsurf
diff --git a/HySoP/hysop/operator/penalization.py b/HySoP/hysop/operator/penalization.py
index 0e8178747..65c79157d 100644
--- a/HySoP/hysop/operator/penalization.py
+++ b/HySoP/hysop/operator/penalization.py
@@ -4,6 +4,7 @@
Penalization operator representation
"""
+from parmepy.constants import debug
from parmepy.operator.continuous import Operator
from parmepy.mpi.topology import Cartesian
from parmepy.operator.penalization_d import Penalization_d
diff --git a/HySoP/hysop/operator/penalization_d.py b/HySoP/hysop/operator/penalization_d.py
index 60c162de5..05a507e09 100644
--- a/HySoP/hysop/operator/penalization_d.py
+++ b/HySoP/hysop/operator/penalization_d.py
@@ -3,6 +3,7 @@
@package operator
Discrete penalization representation
"""
+from parmepy.constants import debug
from parmepy.operator.discrete import DiscreteOperator
from parmepy.numerics.differentialOperator import DifferentialOperator
from parmepy.constants import np
diff --git a/HySoP/hysop/operator/poisson.py b/HySoP/hysop/operator/poisson.py
index e70f8d075..56a5d8d10 100644
--- a/HySoP/hysop/operator/poisson.py
+++ b/HySoP/hysop/operator/poisson.py
@@ -4,7 +4,10 @@
Continous Poisson solver operator (F2PY)
"""
from parmepy.operator.continuous import Operator
+from parmepy.f2py import fftw2py
+from parmepy.mpi.topology import Cartesian
from parmepy.operator.poisson_fft import PoissonFFT
+from parmepy.constants import debug
class Poisson(Operator):
@@ -41,6 +44,21 @@ class Poisson(Operator):
@param topology : topology of the input fields
"""
Operator.setUp(self)
+ print self.method
+ localres, localoffset = fftw2py.init_fftw_solver(
+ self.resolutions[self.vorticity],
+ self.domain.length)
+ print 'localres, localoffset', localres, localoffset
+ topodims = self.resolutions[self.vorticity] / localres
+ print 'topodims POISSON', topodims
+ #variables discretization
+ for v in self.variables:
+ topoId, topo = self.domain.addTopology(
+ Cartesian.withResolution(
+ domain=self.domain, topoResolution=topodims,
+ globalMeshResolution=self.resolutions[v]))
+ df, dfId = v.discretize(topo)
+ self.discreteFieldId[v] = dfId
self.discreteOperator = PoissonFFT(self, method=self.method,
**self.config)
self.discreteOperator.setUp()
diff --git a/HySoP/hysop/operator/poisson_fft.py b/HySoP/hysop/operator/poisson_fft.py
index 866c3eb72..d3d85b8ac 100644
--- a/HySoP/hysop/operator/poisson_fft.py
+++ b/HySoP/hysop/operator/poisson_fft.py
@@ -5,7 +5,7 @@ Discrete Poisson solver operator using FFT
"""
from parmepy.f2py import fftw2py
from parmepy.operator.discrete import DiscreteOperator
-from parmepy.constants import np, ORDER, PARMES_REAL
+from parmepy.constants import np, ORDER, PARMES_REAL, debug
import time
@@ -15,30 +15,33 @@ class PoissonFFT(DiscreteOperator):
"""
@debug
- def __init__(self, poisson, method=''):
+ def __init__(self, poisson, method='', domain=None):
"""
Constructor.
@param velocity : descretized velocity to update from vorticity.
@param vorticity
"""
- DiscreteOperator.__init__(self)
+ DiscreteOperator.__init__(self, poisson)
+ ## Parent Continous Op.
self.poisson = poisson
- self.velocity = self.poisson.velocity.getDiscreteField(
- self.poisson.resolutions[self.poisson.velocity])
- self.vorticity = self.poisson.vorticity.getDiscreteField(
- self.poisson.resolutions[self.poisson.vorticity])
- self.compute_time = 0.
+ ## Velocity.
+ self.velocity = self.poisson.velocity.discreteField[
+ self.poisson.discreteFieldId[self.poisson.velocity]]
+ ## Vorticity.
+ self.vorticity = self.poisson.vorticity.discreteField[
+ self.poisson.discreteFieldId[self.poisson.vorticity]]
self.method = method
+ self.domain = domain
+ self.compute_time = 0.
@debug
def setUp(self):
- #on recupere la topologie par les variables discretes
+ # topology is recovered from discrete velocity field topology
self.topology = self.velocity.topology
- self.topology = self.vorticity.topology
- self.ghosts = self.topology.ghosts
- self.resolution = self.topology.mesh.resolution
- self.dim = self.topology.dim
self.coords = self.topology.mesh.coords
+ self.step = self.topology.mesh.space_step
+ self.local_start = self.topology.mesh.local_start
+ self.local_end = self.topology.mesh.local_end
@debug
def apply(self, t, dt, ite):
@@ -47,90 +50,69 @@ class PoissonFFT(DiscreteOperator):
"""
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]
-
- 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]
-
- # 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][...])
+ #=== RECOVER VELOCITY FIELD FROM VORTICITY FIELD ===
+ self.velocity.data[0],\
+ self.velocity.data[1],\
+ self.velocity.data[2] = \
+ fftw2py.solve_poisson_3d(self.vorticity.data[0],
+ self.vorticity.data[1],
+ self.vorticity.data[2],
+ self.velocity.data[0],
+ self.velocity.data[1],
+ self.velocity.data[2])
- for i in xrange(self.dim):
- self.velocity[i][ind0a:ind0b,
- ind1a:ind1b, ind2a:ind2b] = velocityNoG[i][...]
- #======== VELOCITY CORRECTION =====================
+ #============== VELOCITY CORRECTION ================
- surf = (12.)**2
+ surf = (self.domain.length[1]-self.domain.origin[1]) * \
+ (self.domain.length[2]-self.domain.origin[2])
- #======== CORRECTION VEL X : u_x ==================
+ #============= CORRECTION VEL X : u_x ==============
# computation of the current flow rate: int_y int_z uX dydz
- debX = 0.
- debX = np.sum(self.velocity[0][0,ind1a:ind1b, ind2a:ind2b])
- debX = debX * self.topology.mesh.size[1] * self.topology.mesh.size[2]
-
- self.velocity[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = \
- self.velocity[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] \
- + 1. - debX / surf
-
- #======== CORRECTION VEL Y : u_y ==================
-
- # computation of the current flow rate: int_y int_z uY dydz
- debY = 0.
- debY = np.sum(self.velocity[1][0,ind1a:ind1b, ind2a:ind2b])
- debY = debY * self.topology.mesh.size[1] * self.topology.mesh.size[2]
-
- # computation of omega_z circulation
- omZbar = 0.
- omZbar = np.sum(self.vorticity[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b])
- omZbar = omZbar * self.topology.mesh.size[1] * self.topology.mesh.size[2] / surf
-
- self.velocity[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = \
- self.velocity[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] \
- + omZbar * self.coords[0][ind0a:ind0b] \
- - omZbar * self.coords[0][ind0a] \
- - debY / surf
-
- #======== CORRECTION VEL Z : u_z ==================
-
- # computation of the current flow rate: int_y int_z uZ dydz
- debZ = 0.
- debZ = np.sum(self.velocity[2][0,ind1a:ind1b, ind2a:ind2b])
- debZ = debZ * self.topology.mesh.size[1] * self.topology.mesh.size[2]
-
- # computation of omega_y circulation
- omYbar = 0.
- omYbar = np.sum(self.vorticity[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b])
- omYbar = omYbar * self.topology.mesh.size[1] * self.topology.mesh.size[2] / surf
-
- self.velocity[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = \
- self.velocity[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] \
- - omYbar * self.coords[0][ind0a:ind0b] \
- + omYbar * self.coords[0][ind0a] \
- - debZ / surf
+# debX = 0.
+# debX = np.sum(self.velocity.data[0][0, :, :])
+# debX = debX * self.step[1] * self.step[2]
+
+# self.velocity.data[0] = \
+# self.velocity.data[0] + 1. - debX / surf
+
+# #============= CORRECTION VEL Y : u_y ==============
+
+# # computation of the current flow rate: int_y int_z uY dydz
+# debY = 0.
+# debY = np.sum(self.velocity.data[1][0, :, :])
+# debY = debY * self.step[1] * self.step[2]
+
+# # computation of omega_z circulation
+# omZbar = 0.
+# omZbar = np.sum(self.vorticity.data[2])
+# omZbar = omZbar * self.topology.mesh.space_step_size[1] * \
+# self.topology.mesh.space_step_size[2] / surf
+
+# self.velocity.data[1] = \
+# self.velocity.data[1] \
+# + omZbar * self.coords[0][self.local_start:self.local_end] \
+# - omZbar * self.coords[0][self.local_start] \
+# - debY / surf
+
+# #============= CORRECTION VEL Z : u_z ==============
+
+# # computation of the current flow rate: int_y int_z uZ dydz
+# debZ = 0.
+# debZ = np.sum(self.velocity.data[2][0, :, :])
+# debZ = debZ * self.step[1] * self.step[2]
+
+# # computation of omega_y circulation
+# omYbar = 0.
+# omYbar = np.sum(self.vorticity.data[1])
+# omYbar = omYbar *self.step[1] * self.step[2] / surf
+
+# self.velocity.data[2] = \
+# self.velocity.data[2] \
+# - omYbar * self.coords[0][self.local_start:self.local_end] \
+# + omYbar * self.coords[0][self.local_start] \
+# - debZ / surf
# PENSER AU MPI_REDUCE !!
@@ -140,9 +122,9 @@ class PoissonFFT(DiscreteOperator):
return self.compute_time
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
+# 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):
diff --git a/HySoP/hysop/operator/scales_advection.py b/HySoP/hysop/operator/scales_advection.py
index fdc70b306..d541c5e83 100644
--- a/HySoP/hysop/operator/scales_advection.py
+++ b/HySoP/hysop/operator/scales_advection.py
@@ -41,11 +41,12 @@ class ScalesAdvection(DiscreteOperator):
@debug
def setUp(self):
+ pass
#self.stab_coeff = stab_coeff
- self.topology = self.velocity.topology
- self.ghosts = self.topology.ghosts
- self.resolution = self.topology.mesh.resolution
- self.dim = self.topology.dim
+# self.topology = self.velocity.topology
+# self.ghosts = self.topology.ghosts
+# self.resolution = self.topology.mesh.resolution
+# self.dim = self.topology.dim
@debug
def apply(self, t, dt, ite):
@@ -55,12 +56,12 @@ class ScalesAdvection(DiscreteOperator):
self.compute_time = time.time()
for adF in self.advectedFields:
if isinstance(adF, VectorField):
- for d in xrange(self.dim):
+ for d in xrange(adF.dimension):
adF[d][...] = scales2py.solve_advection(
dt, self.velocity.data[0],
self.velocity.data[1],
self.velocity.data[2],
- adF.data[d])
+ adF[d][...])
else:
adF[...] = scales2py.solve_advection(
dt, self.velocity.data[0],
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index 928f11249..2f5b24a43 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -6,6 +6,7 @@ Stretching operator representation
"""
+from parmepy.constants import debug
from parmepy.operator.continuous import Operator
from parmepy.mpi.topology import Cartesian
from parmepy.operator.stretching_d import Stretching_d
diff --git a/HySoP/hysop/operator/stretching_d.py b/HySoP/hysop/operator/stretching_d.py
index c243d1b5d..23b58188c 100755
--- a/HySoP/hysop/operator/stretching_d.py
+++ b/HySoP/hysop/operator/stretching_d.py
@@ -5,6 +5,8 @@
Discrete stretching representation
"""
+from parmepy.constants import debug
+from numpy import linalg as LA
from parmepy.operator.discrete import DiscreteOperator
from parmepy.numerics.differentialOperator import DifferentialOperator
from parmepy.numerics.fct2op import Fct2Op
@@ -125,12 +127,13 @@ class Stretching_d(DiscreteOperator):
choice='gradV',
topology=self.vorticity.topology)
- # Time integration
+ ## Time integration
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)
+# print 'norm omega', LA.norm(self.vorticity.data)
self.compute_time = time.time() - self.compute_time
self.total_time += self.compute_time
--
GitLab