diff --git a/Examples/NS_bluff_bodies.py b/Examples/NS_bluff_bodies.py
index 7c7fbe1b6ffd502ed3b7fc96cb1ec7aec9f23fb2..4fa4de483862eaf62c3b6dadaaa708c446e348c9 100755
--- a/Examples/NS_bluff_bodies.py
+++ b/Examples/NS_bluff_bodies.py
@@ -43,13 +43,13 @@ sin = np.sin
## Domain
#box = pp.Box(dim, length=[1., 1., 1.], origin=[0., 0., 0.])
-box = pp.Box(dim, length=[2.0 * pi, pi, pi], origin=[0., 0., 0.])
+box = pp.Box(dim, length=[4.0 * pi, 2.0 * pi, 2.0 * pi], origin=[-pi, -pi, -pi])
#box = pp.Box(dim, length=[12., 10., 10.], origin=[-4., -5., -5.])
#box = pp.Box(dim, length=[2.0 * pi, 2.0 * pi, 2.0 * pi])
## Global resolution
#nbElem = [nb] * dim
-nbElem = [65, 33, 33]
+nbElem = [129, 65, 65]
# Upstream flow velocity
uinf = 1.0
@@ -68,6 +68,9 @@ def computeVel(res, x, y, z, t):
res[0][...] = - ((uinf * x) / module) * (1 - (RADIUS ** 3 / module ** 3))
res[1][...] = ((uinf * y) / module) * (1 + (RADIUS ** 3 / (2. * module ** 2)))
res[2][...] = 0.
+# res[0][...] = uinf
+# res[1][...] = 0.
+# res[2][...] = 0.
return res
# Function to compute initial vorticity
@@ -152,15 +155,6 @@ poisson = Poisson(velo, vorti,
vorti: nbElem},
projection=PROJ)
-correction = VelocityCorrection(velo, vorti,
- resolutions={velo: nbElem,
- vorti: nbElem},
- uinf=uinf)
-
-penal = Penalization(velo, [sphere],
- coeff=[1e8],
- resolutions={velo: nbElem})
-
curl = Curl(velo, vorti,
resolutions={velo: nbElem,
vorti: nbElem},
@@ -182,7 +176,7 @@ distrCurlAdv = Redistribute([vorti, velo], curl, advec)
## Simulation with fixed time step
simu = Simulation(tinit=0.0,
- tend=5., timeStep=0.005,
+ tend=10., timeStep=0.005,
iterMax=1000000)
## Define the problem to solve
@@ -196,20 +190,17 @@ simu = Simulation(tinit=0.0,
# distrStrDiff,
# diffusion,
# poisson],
-## correction,
-## penal],
# simulation=simu, dumpFreq=-1)
pb = NSProblem(operators=[curl,
- distrAdvStr_velo, # Redistribute while advec
+ distrCurlAdv,
+ distrCurlStr_velo, # Redistribute while advec
advec, # <--------------------:
distrAdvStr_vorti,
stretch,
distrStrDiff,
diffusion,
- poisson,
- correction,
- penal],
+ poisson],
simulation=simu, dumpFreq=-1)
## Setting solver to Problem (only operators for computational tasks)
@@ -218,6 +209,18 @@ pb.pre_setUp()
## Topology without ghost points
topofft = poisson.discreteFields[poisson.vorticity].topology
+## Add operators related to the problem and depending on FFT topology
+correction = VelocityCorrection(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ uinf=uinf,
+ topo=topofft)
+
+penal = Penalization(velo, [sphere],
+ coeff=[1e8],
+ topo=topofft,
+ resolutions={velo: nbElem})
+
## Diagnostics/monitors related to the problem
forces = Forces(velo, vorti, topo,
@@ -226,7 +229,7 @@ forces = Forces(velo, vorti, topo,
frequency=1, prefix='./res/Noca_sphere.dat')
printer = Printer(variables=[velo, vorti],
- topo=topo,
+ topo=topofft,
frequency=100,
prefix='./res/NS_sphere',
formattype=VTK)
@@ -234,14 +237,18 @@ printer = Printer(variables=[velo, vorti],
energy = Energy_enstrophy(velo, vorti,
topo=topofft,
viscosity=VISCOSITY,
- isNormalized=False,
+ isNormalized=True,
frequency=OUTPUT_FREQ,
prefix=FILENAME)
## Add monitors and setting solver to Problem
vorti.setTopoInit(topofft)
velo.setTopoInit(topofft)
-pb.addMonitors([energy, printer])
+pb.addMonitors([correction, penal, energy, printer])
+penal.discretize()
+correction.discretize()
+correction.setUp()
+penal.setUp()
pb.setUp()
penal.apply(simu)
diff --git a/Examples/dataNS_bb.py b/Examples/dataNS_bb.py
index adacf14feffdbaf2ea4d5b7e01953a199fd1c220..e7783e24c32b401da71a320848754d8ff0f5af30 100644
--- a/Examples/dataNS_bb.py
+++ b/Examples/dataNS_bb.py
@@ -50,11 +50,11 @@ PROJ = [False, 10]
## pi constant
pi = np.pi
# Obstacle description
-OBST_Ox = np.pi / 2.0
-OBST_Oy = np.pi / 2.0
-OBST_Oz = np.pi / 2.0
+OBST_Ox = 0.
+OBST_Oy = 0.
+OBST_Oz = 0.
RADIUS = 0.5
# post treatment ...
OUTPUT_FREQ = 1
-FILENAME = './res/energy_NS_potFlow_res2_2.dat'
+FILENAME = './res/energy_NS_sphere_FFTcurl_8cpu.dat'
diff --git a/HySoP/hysop/operator/discrete/velocity_correction.py b/HySoP/hysop/operator/discrete/velocity_correction.py
index e262b7ceb6199b37eca25e12c9e94648a2f7baee..f36f88ad470c7eb2afa5449169c1fd95aabba301 100755
--- a/HySoP/hysop/operator/discrete/velocity_correction.py
+++ b/HySoP/hysop/operator/discrete/velocity_correction.py
@@ -50,16 +50,15 @@ class VelocityCorrection_D(DiscreteOperator):
def setUp(self):
spaceStep = self.velocity.topology.mesh.space_step
- self.subLength = self.velocity.topology.mesh.end - \
- self.velocity.topology.mesh.origin
- self.coeff_mean = np.prod(spaceStep) / np.prod(self.subLength)
+ lengths = self.velocity.topology.domain.length
+ self.coeff_mean = np.prod(spaceStep) / np.prod(lengths)
self.x0_coord = self.velocity.topology.domain.origin[XDIR]
self.x_coord = self.velocity.topology.mesh.coords[XDIR]
if self.case is '2D':
- self.surf = self.subLength[YDIR]
+ self.surf = lengths[YDIR]
self.dvol = spaceStep[YDIR]
elif self.case is '3D':
- self.surf = self.subLength[YDIR] * self.subLength[ZDIR]
+ self.surf = lengths[YDIR] * lengths[ZDIR]
self.dvol = spaceStep[YDIR] * spaceStep[ZDIR]
else:
raise ValueError("invalid problem dimension")
@@ -69,11 +68,12 @@ class VelocityCorrection_D(DiscreteOperator):
@timed_function
def apply(self, simulation):
- # Calling for requirements completion
+ ## Calling for requirements completion
DiscreteOperator.apply(self, simulation)
ctime = MPI.Wtime()
- # Computation of the flowrates evaluated from current (ie non corrected) velocity
+ ## Computation of the flowrates evaluated from
+ ## current (ie non corrected) velocity
local_comput_flowrates = npw.zeros((3))
for i in xrange(self.velocity.nbComponents):
local_comput_flowrates[i] = np.sum(self.velocity[i][...])
@@ -85,26 +85,45 @@ class VelocityCorrection_D(DiscreteOperator):
for i in xrange(self.velocity.nbComponents):
self.comput_flowrates[i] = recvbuff[i] * self.dvol
- # Correction of the X-velocity component
+ ## Correction of the X-velocity component
self.velocity[XDIR][...] += self.uinf - \
self.comput_flowrates[XDIR] / self.surf
- # Correction of the other velocity components
+ ## Correction of the other velocity components
if self.case is '2D':
- vort_mean = np.sum(self.vorticity[0][...]) * self.coeff_mean
+ # Computation of vorticity mean (in space)
+ local_vort_mean = np.sum(self.vorticity[0][...])
+ recvbuff = 0.0
+ self.velocity.topology.topo.Allreduce(local_vort_mean,
+ recvbuff)
+ vort_mean = recvbuff * self.coeff_mean
+
+ # Correction of the Y-velocity component
self.velocity[YDIR][...] += vort_mean * \
(self.x_coord - self.x0_coord) - \
- self.comput_flowrates[YDIR] / self.surf
+ self.comput_flowrates[YDIR] \
+ / self.surf
elif self.case is '3D':
- vort_Y_mean = np.sum(self.vorticity[YDIR][...]) * self.coeff_mean
- vort_Z_mean = np.sum(self.vorticity[ZDIR][...]) * self.coeff_mean
- self.velocity[YDIR][...] += vort_Z_mean * \
+ # Computation of Y and Z-vorticity means (in space)
+ local_vort_Y_mean = np.sum(self.vorticity[YDIR][...])
+ local_vort_Z_mean = np.sum(self.vorticity[ZDIR][...])
+
+ sendbuff = npw.zeros((2))
+ recvbuff = npw.zeros((2))
+ sendbuff[:] = [local_vort_Y_mean, local_vort_Z_mean]
+ self.velocity.topology.topo.Allreduce(sendbuff, recvbuff)
+ vort_mean = recvbuff * self.coeff_mean
+
+ # Correction of the Y and Z-velocity components
+ self.velocity[YDIR][...] += vort_mean[1] * \
(self.x_coord - self.x0_coord) - \
- self.comput_flowrates[YDIR] / self.surf
- self.velocity[ZDIR][...] += vort_Y_mean * \
+ self.comput_flowrates[YDIR] \
+ / self.surf
+ self.velocity[ZDIR][...] += vort_mean[0] * \
(self.x_coord - self.x0_coord) - \
- self.comput_flowrates[ZDIR] / self.surf
+ self.comput_flowrates[ZDIR] \
+ / self.surf
else:
raise ValueError("invalid problem dimension")