diff --git a/Examples/NavierStokes3d_linked.py b/Examples/NavierStokes3d_linked.py
index dc6d93126938eafcf4f144842e4340dad9168193..3e2c541ead155bc0c6111c283e5012f9715fe33b 100755
--- a/Examples/NavierStokes3d_linked.py
+++ b/Examples/NavierStokes3d_linked.py
@@ -22,10 +22,30 @@ print "Mpi process number ", rank
## ----------- A 3d problem -----------
print " ========= Start test for Navier-Stokes 3D ========="
+## Opening/reading input data file
+inputData = open("input.dat", "r")
+
+Ox = np.float64(inputData.readline().rstrip('\n\r'))
+Oy = np.float64(inputData.readline().rstrip('\n\r'))
+Oz = np.float64(inputData.readline().rstrip('\n\r'))
+Lx = np.float64(inputData.readline().rstrip('\n\r'))
+Ly = np.float64(inputData.readline().rstrip('\n\r'))
+Lz = np.float64(inputData.readline().rstrip('\n\r'))
+resX = np.uint32(inputData.readline().rstrip('\n\r'))
+resY = np.uint32(inputData.readline().rstrip('\n\r'))
+resZ = np.uint32(inputData.readline().rstrip('\n\r'))
+timeStep = np.float64(inputData.readline().rstrip('\n\r'))
+finalTime = np.float64(inputData.readline().rstrip('\n\r'))
+outputModulo = np.uint32(inputData.readline().rstrip('\n\r'))
+outputFilePrefix = inputData.readline().rstrip('\n\r')
+nbGhosts = np.uint32(inputData.readline().rstrip('\n\r'))
+visco = np.float64(inputData.readline().rstrip('\n\r'))
+
+inputData.close()
+
## Physical Domain description
-Lx = Ly = Lz = 2 * pi
-myDomain3d = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[0., 0., 0.])
-resolution3d = np.asarray((129, 129, 129))
+myDomain3d = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[Ox, Oy, Oz])
+resolution3d = np.asarray((resX, resY, resZ))
ncells = resolution3d - 1
hx = Lx / ncells[0]
hy = Ly / ncells[1]
@@ -33,17 +53,17 @@ hz = Lz / ncells[2]
## Obstacle
lambd = np.array([0, 10, 10E8], dtype=PARMES_REAL, order=ORDER)
-sphere = pp.Obstacle(myDomain3d, name='sphere', zlayer=0.1,
- radius=0.5, center=[Lx/2., Ly/2., Lz/2.],
+sphere = pp.Obstacle(myDomain3d, name='sphere', zlayer=0.05,
+ radius=0.5, center=[0., 0., 0.],
orientation='West', porousLayer=0.)
-## Post
-outputFilePrefix = './res/NS_'
-outputModulo = 1
+### Post
+#outputFilePrefix = './res/NS_'
+#outputModulo = 50
-## Simulation parameters
-timeStep = 1e-8
-finalTime = timeStep#1.0
+### Simulation parameters
+#timeStep = 1.#1e-1
+#finalTime = 10000.#20.*timeStep
## Topologies definitions
@@ -55,7 +75,7 @@ print "FFT solver local resolution/offset: ", localres, localoffset
# ---- Cartesian topology -----
topoCart = pp.CartesianTopology(domain=myDomain3d,
- resolution=resolution3d, dim=3, ghosts=[2,2,2])
+ resolution=resolution3d, dim=3, ghosts=[nbGhosts,nbGhosts,nbGhosts])
# ---- JB's topology ----------
nbProcs = MPI.COMM_WORLD.Get_size()
@@ -69,11 +89,16 @@ def computeVel(x, y, z):
# vy = 2. / np.sqrt(3) * np.sin(-2. * pi / 3.) * np.cos(x) \
# * np.sin(y) * np.cos(z)
# vz = 0.
- uinf = 1.
+#------------------
+# uinf = 1.
+# vx = 0.
+# module = (x)**2 + (y)**2 + (z)**2
+# if (module >= (0.5)**2):
+# vx = uinf * (1-(0.5)**2/module)
+# vy = 0.
+# vz = 0.
+#-----------------
vx = 0.
- module = (x-Lx/2)**2 + (y-Ly/2)**2 + (z-Lz/2)**2
- if (module >= (0.5)**2):
- vx = uinf * (1-(0.5)**2/module)
vy = 0.
vz = 0.
return vx, vy, vz
@@ -83,9 +108,28 @@ def computeVort(x, y, z):
# wx = - np.cos(x) * np.sin(y) * np.sin(z)
# wy = - np.sin(x) * np.cos(y) * np.sin(z)
# wz = 2. * np.sin(x) * np.sin(y) * np.cos(z)
+#-----------------
+# wx = 0.
+# wy = 0.
+# wz = 0.
+ xc = 6./2.
+ yc = 6./2.
+ zc = 6./6.
+ R = 1.5
+ sigma = R / 3.
+ Gamma = 0.0075
+ dist = m.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) * m.exp(-(s2 / sigma**2))
+ wx = - wTheta * sinTheta
+ wy = wTheta * cosTheta
+ wz = 0.
return wx, wy, wz
def computeDensity(x, y, z):
@@ -125,18 +169,19 @@ scalesres, scalesoffset, stab_coeff = \
scales.init_advection_solver(ncells, myDomain3d.length,
topodims, order='p_O2')
+print 'advection stability coeff', stab_coeff
advecVort = pp.Advec_scales(stab_coeff, velocity, vorticity)
#advecDensity = pp.Advec_scales(stab_coeff, velocity, density)
# ----> Change TOPO from Scales to Cartesian
## 3 - Stretching
-stretch = pp.Stretching(vorticity, velocity)
+stretch = pp.Stretching(velocity, vorticity)
# ----> Change TOPO from Cartesian to FFT
## 4 - Diffusion
-diffusion = pp.Diffusion(velocity, vorticity, viscosity=0.001)
+diffusion = pp.Diffusion(velocity, vorticity, viscosity=visco)
## 4 - Poisson
poisson = pp.Poisson(velocity, vorticity)
@@ -149,7 +194,8 @@ penal = pp.Penalization(velocity, vorticity, sphere, lambd)
# ----> Change TOPO from Cartesian to Scales
## Define the problem to solve
-navierStokes = pp.Problem(topoCart, [advecVort, stretch, diffusion, poisson, penal])
+navierStokes = pp.Problem(topoCart, [poisson, diffusion, advecVort, stretch])
+#navierStokes = pp.Problem(topoCart, [advecVort, stretch, poisson])
printer = pp.Printer(fields=[vorticity, velocity], frequency=outputModulo,
outputPrefix=outputFilePrefix)
@@ -159,12 +205,13 @@ navierStokes.setSolver(finalTime, timeStep, solver_type='basic', io=printer)
## Problem => ParticularSover = basic.initialize()
navierStokes.initSolver()
-penal.apply(0., timeStep)
+#penal.apply(0., timeStep)
+#poisson.apply(0., timeStep)
## solve stretching
navierStokes.solve()
-print 'vorticity', vorticity.discreteField[0][0] , vorticity.discreteField[0][1] , vorticity.discreteField[0][2]
+#print 'vorticity', vorticity.discreteField[0][0] , vorticity.discreteField[0][1] , vorticity.discreteField[0][2]
## end of time loop ##
diff --git a/Examples/input.dat b/Examples/input.dat
new file mode 100644
index 0000000000000000000000000000000000000000..c7ac525a062f61dd44c09dd892ddcc6b37c39073
--- /dev/null
+++ b/Examples/input.dat
@@ -0,0 +1,15 @@
+0.0
+0.0
+0.0
+6.0
+6.0
+6.0
+257
+257
+257
+1.
+10000.
+50
+./res/NS_
+2
+0.000001
diff --git a/Examples/oarscript.sh b/Examples/oarscript.sh
new file mode 100755
index 0000000000000000000000000000000000000000..fc4fea944eea0568ba0c23afd010699104612b15
--- /dev/null
+++ b/Examples/oarscript.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+#OAR -n NS_256
+#OAR -l /core=1,walltime=80:00:00
+
+# OAR cree un repertoire temporaire
+export TMPDIR=/scratch/$USER/oar.$OAR_JOB_NAME
+
+# Liste des noeuds (nom des noeuds repete autant de fois que nb de coeurs specifies)
+export NODES=`awk -v ORS=, '{print}' $OAR_FILE_NODES|sed 's/,$//'`
+# Chargement des modules necessaires
+module add CMake-2.8.9 fftw-3.3.1-gnu python-2.7.2
+# Lancement du calcul
+mkdir -p $TMPDIR
+mkdir -p $TMPDIR/res
+cp input.dat $TMPDIR
+
+export EXE_DIR=`pwd`
+cd $TMPDIR
+ipython $EXE_DIR/NavierStokes3d_linked.py
diff --git a/HySoP/hysop/operator/differentialOperator_d.py b/HySoP/hysop/operator/differentialOperator_d.py
index 0d2f90135e0c802a8d2f6011176e013679eb6f71..6cf36f79c6e3a2930b7e72f82119664e227fc2b2 100755
--- a/HySoP/hysop/operator/differentialOperator_d.py
+++ b/HySoP/hysop/operator/differentialOperator_d.py
@@ -186,18 +186,21 @@ class DifferentialOperator_d(DiscreteOperator):
## Fourth order scheme
# X components of temp and result
+# temp1 = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER)
temp1 = self.field2[0][...] * 0.
temp1[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[0][ind0-2,ind1a:ind1b,ind2a:ind2b] -
8.0 * self.field2[0][ind0-1,ind1a:ind1b,ind2a:ind2b] +\
8.0 * self.field2[0][ind0+1,ind1a:ind1b,ind2a:ind2b] -\
1.0 * self.field2[0][ind0+2,ind1a:ind1b,ind2a:ind2b]) / (12. * self.meshSize[0])
+# temp2 = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER)
temp2 = self.field2[0][...] * 0.
temp2[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[0][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
8.0 * self.field2[0][ind0a:ind0b,ind1-1,ind2a:ind2b] +\
8.0 * self.field2[0][ind0a:ind0b,ind1+1,ind2a:ind2b] -\
1.0 * self.field2[0][ind0a:ind0b,ind1+2,ind2a:ind2b]) / (12. * self.meshSize[1])
-
+
+# temp3 = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER)
temp3 = self.field2[0][...] * 0.
temp3[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
8.0 * self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2-1] +\
@@ -208,48 +211,54 @@ class DifferentialOperator_d(DiscreteOperator):
maxadv1= np.max(abs(temp1))
# Y components of temp and result
+# temp4 = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER)
temp4 = self.field2[1][...] * 0.
temp4[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[1][ind0-2,ind1a:ind1b,ind2a:ind2b] -
8.0 * self.field2[1][ind0-1,ind1a:ind1b,ind2a:ind2b] +\
8.0 * self.field2[1][ind0+1,ind1a:ind1b,ind2a:ind2b] -\
1.0 * self.field2[1][ind0+2,ind1a:ind1b,ind2a:ind2b]) / (12. * self.meshSize[0])
+# temp5 = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER)
temp5 = self.field2[1][...] * 0.
temp5[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[1][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
8.0 * self.field2[1][ind0a:ind0b,ind1-1,ind2a:ind2b] +\
8.0 * self.field2[1][ind0a:ind0b,ind1+1,ind2a:ind2b] -\
1.0 * self.field2[1][ind0a:ind0b,ind1+2,ind2a:ind2b]) / (12. * self.meshSize[1])
-
+
+# temp6 = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER)
temp6 = self.field2[1][...] * 0.
temp6[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
8.0 * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2-1] +\
8.0 * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2+1] -\
1.0 * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2+2]) / (12. * self.meshSize[2])
- maxstr2= np.max(abs(temp1)+abs(temp2)+abs(temp3))
- maxadv2= np.max(abs(temp2))
+ maxstr2= np.max(abs(temp4)+abs(temp5)+abs(temp6))
+ maxadv2= np.max(abs(temp5))
# Z components of temp and result
+# temp7 = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER)
temp7 = self.field2[2][...] * 0.
temp7[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[2][ind0-2,ind1a:ind1b,ind2a:ind2b] -
8.0 * self.field2[2][ind0-1,ind1a:ind1b,ind2a:ind2b] +\
8.0 * self.field2[2][ind0+1,ind1a:ind1b,ind2a:ind2b] -\
1.0 * self.field2[2][ind0+2,ind1a:ind1b,ind2a:ind2b]) / (12. * self.meshSize[0])
+# temp8 = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER)
temp8 = self.field2[2][...] * 0.
temp8[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[2][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
8.0 * self.field2[2][ind0a:ind0b,ind1-1,ind2a:ind2b] +\
8.0 * self.field2[2][ind0a:ind0b,ind1+1,ind2a:ind2b] -\
1.0 * self.field2[2][ind0a:ind0b,ind1+2,ind2a:ind2b]) / (12. * self.meshSize[1])
-
+
+# temp9 = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER)
temp9 = self.field2[2][...] * 0.
temp9[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
8.0 * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2-1] +\
8.0 * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2+1] -\
1.0 * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2+2]) / (12. * self.meshSize[2])
- maxstr3= np.max(abs(temp1)+abs(temp2)+abs(temp3))
- maxadv3= np.max(abs(temp3))
+ maxstr3= np.max(abs(temp7)+abs(temp8)+abs(temp9))
+ maxadv3= np.max(abs(temp9))
maxgersh[0] = max(maxstr1,maxstr2,maxstr3)
maxgersh[1] = max(maxadv1,maxadv2,maxadv3)
result = np.concatenate((np.array([temp1, temp2, temp3]),np.array([temp4, temp5, temp6]),np.array([temp7, temp8, temp9])))
diff --git a/HySoP/hysop/operator/diffusion_d.py b/HySoP/hysop/operator/diffusion_d.py
index 6c6a5d5e68bef84709d6e44edc457268c41fabdc..3fd674034eb6d8db75d1166550ea966639779ec8 100644
--- a/HySoP/hysop/operator/diffusion_d.py
+++ b/HySoP/hysop/operator/diffusion_d.py
@@ -55,7 +55,6 @@ class Diffusion_d(DiscreteOperator):
for i in xrange (self.dim):
vorticityNoG[i] = self.vorticity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
velocityNoG[i] = self.velocity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
- print 'shape vortiNoG, veloNoG', np.asarray(vorticityNoG).shape, np.asarray(velocityNoG).shape
# Vorticity diffusion
vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...] = \
@@ -65,7 +64,6 @@ class Diffusion_d(DiscreteOperator):
for i in xrange (self.dim):
self.vorticity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = vorticityNoG[i][...]
- print 'shape vorti, velo', np.asarray(self.vorticity.data).shape, np.asarray(self.velocity.data).shape
# print 'shape vorti, velo, diffusion', np.asarray(self.vorticity.data).shape, np.asarray(self.velocity.data).shape
self.compute_time = time.time() - self.compute_time
diff --git a/HySoP/hysop/operator/poisson_d.py b/HySoP/hysop/operator/poisson_d.py
index c537c5c75b206e290c3f9a6991e65a4fcfea99bf..eec04199c86d612eb5bbaa95a784abc2fe30aad9 100644
--- a/HySoP/hysop/operator/poisson_d.py
+++ b/HySoP/hysop/operator/poisson_d.py
@@ -30,6 +30,7 @@ class Poisson_d(DiscreteOperator):
self.ghosts = topology.ghosts
self.resolution = topology.mesh.resolution
self.dim = topology.dim
+ self.coords = topology.mesh.coords
# self.ppfft = parmepy.f2py.fftw2py
self.compute_time = 0.
@@ -54,7 +55,6 @@ class Poisson_d(DiscreteOperator):
for i in xrange (self.dim):
vorticityNoG[i] = self.vorticity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
velocityNoG[i] = self.velocity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
- print 'shape vortiNoG, veloNoG', np.asarray(vorticityNoG).shape, np.asarray(velocityNoG).shape
# Recover velocity field from vorticity field
velocityNoG[0][...], velocityNoG[1][...], velocityNoG[2][...] = \
@@ -63,13 +63,59 @@ class Poisson_d(DiscreteOperator):
for i in xrange (self.dim):
self.velocity[i][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = velocityNoG[i][...]
- print 'shape vorti, velo', np.asarray(self.vorticity.data).shape, np.asarray(self.velocity.data).shape
-# print 'shape vorti, velo, poisson', np.asarray(self.vorticity.data).shape, np.asarray(self.velocity.data).shape
+# #======== 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/(4. * np.pi)**2
+
+# #======== 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] / (4. * np.pi)**2
+
+# 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/(4. * np.pi)**2
+
+# #======== 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] / (4. * np.pi)**2
+
+# 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/(4. * np.pi)**2
+
+ # PENSER AU MPI_REDUCE !!
+
+
self.compute_time = time.time() - self.compute_time
self.total_time += self.compute_time
return self.compute_time
-# return result
def printComputeTime(self):
self.timings_info[0] = "\"Poisson solver total\""
diff --git a/HySoP/src/fftw/fft3d.f90 b/HySoP/src/fftw/fft3d.f90
index e9802ad32492774bf7e143c9605cfcb014e3dba0..5f8ed18c8ed221af6bcce34c34ac9b8e4e47b63e 100755
--- a/HySoP/src/fftw/fft3d.f90
+++ b/HySoP/src/fftw/fft3d.f90
@@ -573,7 +573,8 @@ contains
do j = 1,local_resolution(c_Y)
do k = 1, fft_resolution(c_Z)
do i = 1, local_resolution(c_X)
- coeff = Icmplx/(1. + nudt * kx(i)**2+ky(j)**2+kz(k)**2)
+ !!coeff = Icmplx/(1. + nudt * kx(i)**2+ky(j)**2+kz(k)**2)
+ coeff = Icmplx/(1. + nudt * (kx(i)**2+ky(j)**2+kz(k)**2))
buffer1 = dataout1(i,k,j)
buffer2 = dataout2(i,k,j)
dataout1(i,k,j) = coeff*(ky(j)*dataout3(i,k,j)-kz(k)*dataout2(i,k,j))