diff --git a/Examples/NavierStokes3d.py b/Examples/NavierStokes3d.py
index 2ba76bae8f937a5eea3348e387685e6f6889fe6c..df497b387f102c3aadbd73ac4a369ed47b88dded 100755
--- a/Examples/NavierStokes3d.py
+++ b/Examples/NavierStokes3d.py
@@ -5,6 +5,9 @@ import numpy as np
import mpi4py.MPI as MPI
import math as m
+PARMES_REAL = pp.constants.PARMES_REAL
+ORDER = pp.constants.ORDER
+
#from numpy import linalg as LA
pi = m.pi
@@ -28,9 +31,19 @@ hx = Lx / ncells[0]
hy = Ly / ncells[1]
hz = Lz / ncells[2]
+## Obstacle
+lambd = np.array([0, 1, 10 ** 8], dtype=PARMES_REAL, order=ORDER)
+sphere = pp.Obstacle(myDomain3d, name='sphere', zlayer=0.1,
+ radius=0.1, center=[0.5, 0.5, 0.5],
+ orientation='West', porousLayer=0.05)
+
+## Post
+outputFilePrefix = './res/NS_'
+outputModulo = 1
+
# Simulation parameters
-finalTime = 1.0
timeStep = 1e-8
+finalTime = timeStep#1.0
# Fields declaration
@@ -46,18 +59,27 @@ print "FFT solver local resolution/offset: ", localres, localoffset
##topofft = poisson.getTopology()
-def computeVel():
- return np.ones((localres))
+def computeVel(x, y, z):
+ vx = 2. / np.sqrt(3) * np.sin(2. * pi / 3.) * np.sin(x) \
+ * np.cos(y) * np.cos(z)
+ vy = 2. / np.sqrt(3) * np.sin(-2. * pi / 3.) * np.cos(x) \
+ * np.sin(y) * np.cos(z)
+ vz = 0.
+ return vx, vy, vz
-def computeVort():
- return np.ones((localres))
+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)
+ return wx, wy, wz
-#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
@@ -65,15 +87,6 @@ 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)
@@ -114,21 +127,33 @@ scalesres, scalesoffset, stab_coeff = \
scales.init_advection_solver(ncells, myDomain3d.length,
topodims, order='p_O2')
+# 3 - Stretching
+stretch = pp.Stretching(vorticity, velocity)
-print "Advection solver local meshes resolution : ", scalesres
+# 4 - Penalization
+penal = pp.Penalization(velocity, vorticity, sphere, lambd)
+##topofft = poisson.getTopology()
+topofft = pp.CartesianTopology(domain=myDomain3d,
+ resolution=resolution3d, dim=3)
-# 3 - Stretching
-###stretch = pp.Stretching(vorticity, velocity)
+navierStokes = pp.Problem(topofft, [penal, stretch])
-###navierStokes = pp.Problem(topofft, [stretch])
-###printer = pp.Printer(fields=[vorticity, velocity], frequency=outputModulo,
-### outputPrefix=outputFilePrefix)
-###navierStokes.setSolver(finalTime, timeStep, solver_type='basic', io=printer)
+printer = pp.Printer(fields=[vorticity, velocity], frequency=outputModulo,
+ outputPrefix=outputFilePrefix)
+
+navierStokes.setSolver(finalTime, timeStep, solver_type='basic', io=printer)
## Problem => ParticularSover = basic.initialize()
-###navierStokes.initSolver()
+navierStokes.initSolver()
+omega_x = vorticity.discreteField[0][0]
+omega_y = vorticity.discreteField[0][1]
+omega_z = vorticity.discreteField[0][2]
+vx = velocity.discreteField[0][0]
+vy = velocity.discreteField[0][1]
+vz = velocity.discreteField[0][2]
+
## end of init ##
# Mind that scales works on arrays of size resolution - 1
@@ -140,21 +165,16 @@ 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()
+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]
nudt = 0.0001
-omega_x, omega_y, omega_z = ppfft.solve_diffusion_3d(nudt, vx, vy, vz)
+omega_x, omega_y, omega_z = \
+ ppfft.solve_diffusion_3d(nudt, vx, vy, vz, omega_x, omega_y, omega_z)
# solve poisson
-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)
## end of time loop ##
diff --git a/Examples/NavierStokes3dStretch.py b/Examples/NavierStokes3dStretch.py
deleted file mode 100644
index a8b300e86f5245044f5e47ed250c8fd4e3ace26f..0000000000000000000000000000000000000000
--- a/Examples/NavierStokes3dStretch.py
+++ /dev/null
@@ -1,142 +0,0 @@
-#!/usr/bin/python
-import parmepy as pp
-import parmepy.f2py
-import numpy as np
-import mpi4py.MPI as MPI
-import math as m
-import time
-
-#from numpy import linalg as LA
-
-pi = m.pi
-
-# Import scales and fftw solvers
-ppfft = parmepy.f2py.fftw2py
-scales = parmepy.f2py.scales2py
-
-rank = MPI.COMM_WORLD.Get_rank()
-print "Mpi process number ", rank
-
-# ----------- A 3d problem -----------
-print " ========= Start test for Navier-Stokes 3D ========="
-
-# Physical Domain description
-Lx = Ly = Lz = 2 * pi
-myDomain3d = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[0., 0., 0.])
-resolution3d = np.asarray((65, 65, 65))
-ncells = resolution3d - 1
-hx = Lx / ncells[0]
-hy = Ly / ncells[1]
-hz = Lz / ncells[2]
-
-# Simulation parameters
-finalTime = 0.05
-timeStep = 1e-3
-outputFilePrefix = './res/Stretching_'
-outputModulo = 1
-
-t0 = time.time()
-
-## Obstacle
-lambd=np.array([0, 1, 10**8],dtype = PARMES_REAL, order=ORDER)
-sphere = pp.Obstacle(box, name='sphere', zlayer=0.1, radius=0.1, center=[0.5,0.5,0.5], orientation='West', porousLayer=0.05)
-
-# Fields declaration
-
-# 1 - Poisson/diffusion solvers initialisation.
-# See poisson3d.py for a working example.
-# poisson = pp.Poisson(vorticity,velocity)
-#
-localres, localoffset = ppfft.init_fftw_solver(resolution3d,
- myDomain3d.length)
-
-print "FFT solver local resolution/offset: ", localres, localoffset
-
-##topofft = poisson.getTopology()
-topofft = pp.CartesianTopology(domain=myDomain3d, resolution=resolution3d, dim=3)
-
-
-def computeVel(x, y, z):
- vx = 2./np.sqrt(3) * np.sin(2.* pi/3.) * np.sin(x) * np.cos(y) * np.cos(z)
- vy = 2./np.sqrt(3) * np.sin(-2.* pi/3.) * np.cos(x) * np.sin(y) * np.cos(z)
- vz = 0.
- return vx , vy , vz
-
-
-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)
- return wx , wy , wz
-
-
-velocity = pp.AnalyticalField(domain=myDomain3d, formula=computeVel,
- name='Velocity', vector=True)
-vorticity = pp.AnalyticalField(domain=myDomain3d, formula=computeVort,
- name='Vorticity', vector=True)
-
-## 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.
-#scalesres, scalesoffset, stab_coeff = \
-# scales.init_advection_solver(ncells, myDomain3d.length, order='p_O2')
-
-# 3 - Stretching
-stretch = pp.Stretching(vorticity, velocity)
-
-# 4 - Penalization
-penal = pp.Penalization(velo, vorti, sphere, lambd)
-
-navierStokes = pp.Problem(topofft, [stretch,penal])
-printer = pp.Printer(fields=[vorticity, velocity], frequency=outputModulo,
- outputPrefix=outputFilePrefix)
-
-navierStokes.setSolver(finalTime, timeStep, solver_type='basic', io=printer)
-
-## Problem => ParticularSover = basic.initialize()
-navierStokes.initSolver()
-
-## 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[1][:-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]
-
-## 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)
-
-t1 = time.time()
-
-## 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)
-
-## solve poisson
-#vx, vy, vz = ppfft.solve_poisson_3d(omega_x, omega_y, omega_z)
-
-## end of time loop ##
-
-# Clean memory buffers
-#ppfft.clean_fftw_solver(myDomain3d.dimension)
-
-tf = time.time()
-
-print "\n"
-print "Total time : ", tf - t0, "sec (CPU)"
-print "Init time : ", t1 - t0, "sec (CPU)"
-print "Solving time : ", tf - t1, "sec (CPU)"