diff --git a/Examples/NavierStokes3d.py b/Examples/NavierStokes3d.py
index 4af070deb2e8401fef6c5e3f40d0c1b5591508bb..2ba76bae8f937a5eea3348e387685e6f6889fe6c 100755
--- a/Examples/NavierStokes3d.py
+++ b/Examples/NavierStokes3d.py
@@ -30,7 +30,7 @@ hz = Lz / ncells[2]
# Simulation parameters
finalTime = 1.0
-timeStep = 1e-3
+timeStep = 1e-8
# Fields declaration
@@ -54,16 +54,69 @@ def computeVort():
return np.ones((localres))
-velocity = pp.AnalyticalField(domain=myDomain3d, formula=computeVel,
- name='Velocity', vector=True)
-vorticity = pp.AnalyticalField(domain=myDomain3d, formula=computeVort,
- name='Vorticity', vector=True)
+#velocity = pp.AnalyticalField(domain=myDomain3d, formula=computeVel,
+# name='Velocity', vector=True)
+#vorticity = pp.AnalyticalField(domain=myDomain3d, formula=computeVort,
+# name='Vorticity', vector=True)
+############ REF ##############
+x = np.arange(localoffset[0], localres[0] + localoffset[0],
+ dtype='float64') * hx
+y = np.arange(localoffset[1], localres[1] + localoffset[1],
+ dtype='float64') * hy
+z = np.arange(localoffset[2], localres[2] + localoffset[2],
+ dtype='float64') * hz
+omega_x = np.zeros((localres), dtype='float64', order='Fortran')
+omega_y = np.zeros((localres), dtype='float64', order='Fortran')
+omega_z = np.zeros((localres), dtype='float64', order='Fortran')
+ref_x = np.zeros((localres), dtype='float64', order='Fortran')
+ref_y = np.zeros((localres), dtype='float64', order='Fortran')
+ref_z = np.zeros((localres), dtype='float64', order='Fortran')
+vx = np.zeros((localres), dtype='float64', order='Fortran')
+vy = np.zeros((localres), dtype='float64', order='Fortran')
+vz = np.zeros((localres), dtype='float64', order='Fortran')
+
+cden = 4 * pi ** 2 * (Ly ** 2 * Lz ** 2 + Lx ** 2 * Lz ** 2 +
+ Lx ** 2 * Ly ** 2) / (Lx ** 2 * Ly ** 2 * Lz ** 2)
+cx = 2 * pi / Lx
+cy = 2 * pi / Ly
+cz = 2 * pi / Lz
+
+# Initialize vorticity (This init should be improved or done if C or Fortran)
+## for k in range(localres[2]):
+## for j in range(localres[1]):
+## for i in range(localres[0]):
+## omega_x[i, j, k] = cden * (m.sin(cx * x[i]) *
+## m.sin(cy * y[j]) * m.cos(cz * z[k]))
+## omega_y[i, j, k] = cden * (m.cos(cx * x[i]) *
+## m.sin(cy * y[j]) * m.sin(cz * z[k]))
+## omega_z[i, j, k] = cden * (m.cos(cx * x[i]) *
+## m.cos(cy * y[j]) * m.sin(cz * z[k]))
+## ref_x[i, j, k] = -cy * (m.cos(cx * x[i]) * m.sin(cy * y[j])
+## * m.sin(cz * z[k])) - \
+## cz * (m.cos(cx * x[i]) * m.sin(cy * y[j]) * m.cos(cz * z[k]))
+
+## ref_y[i, j, k] = -cz * (m.sin(cx * x[i]) * m.sin(cy * y[j])
+## * m.sin(cz * z[k])) + \
+## cx * (m.sin(cx * x[i]) * m.cos(cy * y[j]) * m.sin(cz * z[k]))
+
+## ref_z[i, j, k] = -cx * (m.sin(cx * x[i]) * m.sin(cy * y[j])
+## * m.sin(cz * z[k])) - \
+## cy * (m.sin(cx * x[i]) * m.cos(cy * y[j]) * m.cos(cz * z[k]))
+
# 2 - Advection solver initialisation. See testScales for a working example.
# Based on scales JB solver
# Warning : fields input for scales should be of size (ncells), not localres.
+nbProcs = MPI.COMM_WORLD.Get_size()
+topodims = [1, 1, nbProcs] # fits with fftw topology
+
scalesres, scalesoffset, stab_coeff = \
- scales.init_advection_solver(ncells, myDomain3d.length, order='p_O2')
+ scales.init_advection_solver(ncells, myDomain3d.length,
+ topodims, order='p_O2')
+
+
+print "Advection solver local meshes resolution : ", scalesres
+
# 3 - Stretching
###stretch = pp.Stretching(vorticity, velocity)
@@ -80,30 +133,25 @@ scalesres, scalesoffset, stab_coeff = \
## end of init ##
# Mind that scales works on arrays of size resolution - 1
# --> pointers to subarrays of velocity/vorticity
-svx = velocity[0][:-1, :-1, :-1]
-svy = velocity.y[:-1, :-1, :-1]
-svz = velocity[2][:-1, :-1, :-1]
-somega_x = vorticity[0][:-1, :-1, :-1]
-somega_y = vorticity[1][:-1, :-1, :-1]
-somega_z = vorticity[2][:-1, :-1, :-1]
-
+#somega_x = np.zeros((localcells), dtype='float64', order='Fortran')
# solve advection
-scales.solve_advection(timeStep, svx, svy, svz, somega_x)
-scales.solve_advection(timeStep, svx, svy, svz, somega_y)
-scales.solve_advection(timeStep, svx, svy, svz, somega_z)
+omega_x = scales.solve_advection(timeStep, vx, vy, vz, omega_x)
+omega_y = scales.solve_advection(timeStep, vx, vy, vz, omega_y)
+omega_z = scales.solve_advection(timeStep, vx, vy, vz, omega_z)
# solve stretching
#navierStokes.solve()
# solve diffusion
-vx = velocity[0]
-vy = velocity[1]
-vz = velocity[2]
-omega_x = vorticity[0]
-omega_y = vorticity[1]
-omega_z = vorticity[2]
-
-omega_x, omega_y, omega_z = ppfft.solve_diffusion_3d(vx, vy, vz)
+## vx = velocity[0]
+## vy = velocity[1]
+## vz = velocity[2]
+## omega_x = vorticity[0]
+## omega_y = vorticity[1]
+## omega_z = vorticity[2]
+
+nudt = 0.0001
+omega_x, omega_y, omega_z = ppfft.solve_diffusion_3d(nudt, vx, vy, vz)
# solve poisson
vx, vy, vz = ppfft.solve_poisson_3d(omega_x, omega_y, omega_z)
diff --git a/Examples/poisson2d.py b/Examples/poisson2d.py
index 8dd8c8e937fe4c3c0d9097cb982b18f991829571..d90ccfd3f7246baf1722ae1e5286d97ec3fe6881 100755
--- a/Examples/poisson2d.py
+++ b/Examples/poisson2d.py
@@ -61,8 +61,8 @@ print "test 0", np.allclose(vely, refy1)
##########################
-localres, localoffset = ppfft.init_poisson_solver(resolution2d,
- myDomain2d.length)
+localres, localoffset = ppfft.init_fftw_solver(resolution2d,
+ myDomain2d.length)
print "topo locale :", localres, localoffset
@@ -74,14 +74,16 @@ y = np.arange(localoffset[1], localres[1] + localoffset[1],
Y, X = np.meshgrid(y, x)
omega = np.zeros((localres), order='Fortran', dtype='float64')
+vx = np.zeros((localres), order='Fortran', dtype='float64')
+vy = np.zeros((localres), order='Fortran', dtype='float64')
omega[:, :] = 4 * pi ** 2 * np.cos(2 * pi / Lx * X) * \
np.sin(2. * pi / Ly * Y) * (1. / Lx ** 2 + 1. / Ly ** 2)
refx = 2. * pi / Ly * np.cos(2 * pi / Lx * X) * np.cos(2 * pi / Ly * Y)
refy = 2. * pi / Lx * np.sin(2 * pi / Lx * X) * np.sin(2 * pi / Ly * Y)
-vx, vy = ppfft.solve_poisson_2d(omega)
+vx, vy = ppfft.solve_poisson_2d(omega, vx, vy)
print np.allclose(refx, vx)
print np.allclose(refy, vy)
-ppfft.clean_poisson_solver(myDomain2d.dimension)
+ppfft.clean_fftw_solver(myDomain2d.dimension)
diff --git a/Examples/poisson3d.py b/Examples/poisson3d.py
index a0b5c4cecbffc89d197511d0005206f3b5376cc9..3dcfdd0cf149ae3590757afee0f82c5bce7442aa 100755
--- a/Examples/poisson3d.py
+++ b/Examples/poisson3d.py
@@ -40,6 +40,9 @@ z = np.arange(localoffset[2], localres[2] + localoffset[2],
omega_x = np.zeros((localres), dtype='float64', order='Fortran')
omega_y = np.zeros((localres), dtype='float64', order='Fortran')
omega_z = np.zeros((localres), dtype='float64', order='Fortran')
+vx = np.zeros((localres), dtype='float64', order='Fortran')
+vy = np.zeros((localres), dtype='float64', order='Fortran')
+vz = np.zeros((localres), dtype='float64', order='Fortran')
ref_x = np.zeros((localres), dtype='float64', order='Fortran')
ref_y = np.zeros((localres), dtype='float64', order='Fortran')
ref_z = np.zeros((localres), dtype='float64', order='Fortran')
@@ -76,7 +79,7 @@ for k in range(localres[2]):
print 'start fftw solver ...'
start = MPI.Wtime()
-vx, vy, vz = ppfft.solve_poisson_3d(omega_x, omega_y, omega_z)
+vx, vy, vz = ppfft.solve_poisson_3d(omega_x, omega_y, omega_z, vx, vy, vz)
print "elasped time =", MPI.Wtime() - start
print np.allclose(ref_x, vx)
print np.allclose(ref_y, vy)
diff --git a/Examples/testScales.py b/Examples/testScales.py
index 049aee02a57d0608219badb5ab6488d139e9eed1..0c1f953bf02c3547530b61118c21314566d96825 100755
--- a/Examples/testScales.py
+++ b/Examples/testScales.py
@@ -20,22 +20,33 @@ nbProcs = comm.Get_size()
# Physical Domain description
Lx = Ly = Lz = 2 * pi
myDomain3d = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[0., 0., 0.])
-resolution3d = np.asarray((64, 64, 64))
+nbCells = np.asarray((16, 16, 32))
# ===== Initialize the particular solver =====
+# First choose the number of mpi process in each dir
+#topodims = MPI.Compute_dims(nbProcs, 3)
+if nbProcs > 1:
+ topodims = [1, 2, nbProcs / 2]
+else:
+ topodims = [1, 1, 1]
+
+# It must be such that nbCells is a multiple of topodims
+# and such that nbCells/topodims is a multiple of 2, 4 or 5.
localres, localoffset, stab_coeff = \
- scales.init_advection_solver(resolution3d, myDomain3d.length, order='p_O2')
+ scales.init_advection_solver(nbCells, myDomain3d.length,
+ topodims, order='p_O2')
if(rank == 0):
print "Start initialisation ..."
# ===== Create and initialize Fields =====
-scal3D = np.zeros((localres), dtype='float64', order='Fortran')
+#scal3D = np.zeros((localres), dtype='float64', order='Fortran')
+scal3D = np.zeros((localres), dtype='float64')
r02 = ((myDomain3d.length / 10.).min()) ** 2 # ??
# Coordinates of the lowest point of the current subdomain
-step = myDomain3d.length / resolution3d
+step = myDomain3d.length / nbCells
coord0 = localoffset * localres * step
vz = np.zeros((localres), dtype='float64', order='Fortran')
vx = np.zeros((localres), dtype='float64', order='Fortran') # vz.copy()
@@ -58,7 +69,7 @@ for k in range(localres[2]):
goodScal = scal3D.copy()
goodVelo = vx.copy()
-timeStep = 0.02
+timeStep = 0.01
finalTime = 2. * period
currentTime = 0.0
@@ -71,7 +82,7 @@ if(rank == 0):
# Simulation
t0 = MPI.Wtime()
-scales.solve_advection(timeStep, vx, vy, vz, scal3D)
+scal3D = scales.solve_advection(timeStep, vx, vy, vz, scal3D)
print "elapsed time on processus ", rank, ": ", MPI.Wtime() - t0