Skip to content
Snippets Groups Projects
Commit 92ed116f authored by Franck Pérignon's avatar Franck Pérignon
Browse files

update NS (stretch+penal+poisson ...)

parent 2d793b78
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Examples/NavierStokes3d.py
+ 55
− 35
View file @ 92ed116f
......@@ -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 ##
......
Examples/NavierStokes3dStretch.py deleted 100644 → 0
+ 0
− 142
View file @ 2d793b78
#!/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)"
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment