I ported a Jos Stam's demo about Fluid mech to check
the difference of speed between C implementation and
Python. I think I achieved good results with Python
and there is space to improve, without to extend the
program with C routine, I mean.
--
Good hack,
Alberto Santini
(http://www.albertosantini.it/)
--------------------------- demo.py ---------------------------
""" Real-Time Fluid Dynamics for Games by Jos Stam (2003).
Parts of author's work are also protected
under U. S. patent #6,266,071 B1 [Patent].
Original paper by Jos Stam, "Real-Time Fluid Dynamics for Games".
Proceedings of the Game Developer Conference, March 2003
http://www.dgp.toronto.edu/people/st...earch/pub.html
Tested on
python 2.4
numarray 1.1.1
PyOpenGL-2.0.2.01.py2.4-numpy23
glut-3.7.6
How to use this demo:
Add densities with the right mouse button
Add velocities with the left mouse button and dragging the mouse
Toggle density/velocity display with the 'v' key
Clear the simulation by pressing the 'c' key
"""
import psyco
psyco.full()
import sys
from numarray import Float64, zeros
from solver import vel_step, dens_step
try:
from OpenGL.GLUT import *
from OpenGL.GL import *
from OpenGL.GLU import *
except:
print '''
ERROR: PyOpenGL not installed properly.
'''
sys.exit()
N = 32
size = N+2
dt = 0.1
diff = 0.0
visc = 0.0
force = 5.0
source = 100.0
dvel = False
win_x = 512
win_y = 512
omx = 0.0
omy = 0.0
mx = 0.0
my = 0.0
mouse_down = [False,False,False]
""" Start with two grids: one that contains the density values from the previous
time step and one
that will contain the new values. For each grid cell of the latter we trace the
cell's center
position backwards through the velocity field. We then linearly interpolate from
the grid of
previous density values and assign this value to the current grid cell.
"""
u = zeros((size,size), Float64) # velocity
u_prev = zeros((size,size), Float64)
v = zeros((size,size), Float64) # velocity
v_prev = zeros((size,size), Float64)
dens = zeros((size,size), Float64) # density
dens_prev = zeros((size,size), Float64)
def clear_data():
global u, v, u_prev, v_prev, dens, dens_prev, size
u[0:size,0:size] = 0.0
v[0:size,0:size] = 0.0
u_prev[0:size,0:size] = 0.0
v_prev[0:size,0:size] = 0.0
dens[0:size,0:size] = 0.0
dens_prev[0:size,0:size] = 0.0
def pre_display():
glViewport(0, 0, win_x, win_y )
glMatrixMode(GL_PROJECTION)
glLoadIdentity()
gluOrtho2D(0.0, 1.0, 0.0, 1.0)
glClearColor(0.0, 0.0, 0.0, 1.0)
glClear(GL_COLOR_BUFFER_BIT)
def post_display():
glutSwapBuffers ()
def draw_velocity():
h = 1.0/N
glColor3f(1.0, 1.0, 1.0)
glLineWidth(1.0)
glBegin(GL_LINES)
for i in range(1, N+1):
x = (i-0.5)*h;
for j in range(1, N+1):
y = (j-0.5)*h
glColor3f(1, 0, 0)
glVertex2f(x, y)
glVertex2f(x+u[i,j], y+v[i,j])
glEnd()
def draw_density():
h = 1.0/N
glBegin(GL_QUADS)
for i in range(0, N+1):
x = (i-0.5)*h
for j in range(0, N+1):
y = (j-0.5)*h
d00 = dens[i,j]
d01 = dens[i,j+1]
d10 = dens[i+1,j]
d11 = dens[i+1,j+1]
glColor3f(d00, d00, d00); glVertex2f(x, y)
glColor3f(d10, d10, d10); glVertex2f(x+h, y)
glColor3f(d11, d11, d11); glVertex2f(x+h, y+h)
glColor3f(d01, d01, d01); glVertex2f(x, y+h)
glEnd ()
def get_from_UI(d, u, v):
global omx, omy
d[0:size,0:size] = 0.0
u[0:size,0:size] = 0.0
v[0:size,0:size] = 0.0
if not mouse_down[GLUT_LEFT_BUTTON] and not mouse_down[GLUT_RIGHT_BUTTON]:
return
i = int((mx/float(win_x))*N+1)
j = int(((win_y-float(my))/float(win_y))*float(N)+1.0)
if i<1 or i>N or j<1 or j>N: return
if mouse_down[GLUT_LEFT_BUTTON]:
u[i,j] = force * (mx-omx)
v[i,j] = force * (omy-my)
if mouse_down[GLUT_RIGHT_BUTTON]:
d[i,j] = source
omx = mx
omy = my
def key_func(key, x, y):
global dvel
if key == 'c' or key == 'C':
clear_data()
if key == 'v' or key == 'V':
dvel = not dvel
def mouse_func(button, state, x, y):
global omx, omy, mx, my, mouse_down
omx = mx = x;
omy = my = y;
mouse_down[button] = (state == GLUT_DOWN)
def motion_func(x, y):
global mx, my
mx = x; my = y
def reshape_func(width, height):
global win_x, win_y
glutReshapeWindow(width, height)
win_x = width
win_y = height
def idle_func():
global dens, dens_prev, u, u_prev, v, v_prev, N, visc, dt, diff
get_from_UI(dens_prev, u_prev, v_prev)
vel_step(N, u, v, u_prev, v_prev, visc, dt)
dens_step(N, dens, dens_prev, u, v, diff, dt)
glutPostRedisplay()
def display_func():
pre_display()
if dvel: draw_velocity()
else: draw_density()
post_display()
def open_glut_window():
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE)
glutInitWindowPosition(0, 0)
glutInitWindowSize(win_x, win_y)
glutCreateWindow("Alias | wavefront (porting by Alberto Santini)")
glClearColor(0.0, 0.0, 0.0, 1.0)
glClear(GL_COLOR_BUFFER_BIT)
glutSwapBuffers()
glClear(GL_COLOR_BUFFER_BIT)
glutSwapBuffers()
pre_display ()
glutKeyboardFunc(key_func)
glutMouseFunc(mouse_func)
glutMotionFunc(motion_func)
glutReshapeFunc(reshape_func)
glutIdleFunc(idle_func)
glutDisplayFunc(display_func)
if __name__ == '__main__':
glutInit(sys.argv)
clear_data()
open_glut_window()
glutMainLoop()
---------------------------------------------------------------
--------------------------- solver.py ---------------------------
""" Real-Time Fluid Dynamics for Games by Jos Stam (2003).
Parts of author's work are also protected
under U. S. patent #6,266,071 B1 [Patent].
"""
def set_bnd(N, b, x):
"""We assume that the fluid is contained in a
box with solid walls: no flow should exit the walls. This simply means that
the horizontal
component of the velocity should be zero on the vertical walls, while the
vertical component of
the velocity should be zero on the horizontal walls. For the density and other
fields considered
in the code we simply assume continuity. The following code implements these
conditions.
"""
for i in range(1, N+1):
if b == 1: x[0,i] = -x[1,i]
else: x[0,i] = x[1,i]
if b == 1: x[N+1,i] = -x[N,i]
else: x[N+1,i] = x[N,i]
if b == 2: x[i,0] = -x[i,1]
else: x[i,0] = x[i,1]
if b == 2: x[i,N+1] = -x[i,N]
else: x[i,N+1] = x[i,N]
x[0,0] = 0.5*(x[1,0]+x[0,1])
x[0,N+1] = 0.5*(x[1,N+1]+x[0,N])
x[N+1,0] = 0.5*(x[N,0]+x[N+1,1])
x[N+1,N+1] = 0.5*(x[N,N+1]+x[N+1,N])
def lin_solve(N, b, x, x0, a, c):
for k in range(0, 20):
x[1:N+1,1:N+1] = (x0[1:N+1,1:N+1] +
a*(x[0:N,1:N+1]+x[2:N+2,1:N+1]+x[1:N+1,0:N]+x[1:N+1,2:N+2]))/c
set_bnd(N, b, x)
# Addition of forces: the density increases due to sources
def add_source(N, x, s, dt):
size = (N+2)
x[0:size,0:size] += dt*s[0:size,0:size]
# Diffusion: the density diffuses at a certain rate
def diffuse (N, b, x, x0, diff, dt):
""" The basic idea behind our method is to find the densities which when
diffused backward
in time yield the densities we started with.
The simplest iterative solver which works well in practice is Gauss-Seidel
relaxation.
"""
a = dt*diff*N*N
lin_solve(N, b, x, x0, a, 1+4*a)
# Advection: the density follows the velocity field
def advect (N, b, d, d0, u, v, dt):
""" The basic idea behind the advection step. Instead of moving the cell
centers forward in
time through the velocity field, we look for the particles which end up
exactly at
the cell centers by tracing backwards in time from the cell centers.
"""
dt0 = dt*N;
for i in range(1, N+1):
for j in range(1, N+1):
x = i-dt0*u[i,j]; y = j-dt0*v[i,j]
if x<0.5: x=0.5
if x>N+0.5: x=N+0.5
i0 = int(x); i1 = i0+1
if y<0.5: y=0.5
if y>N+0.5: y=N+0.5
j0 = int(y); j1 = j0+1
s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1
d[i,j] = s0*(t0*d0[i0,j0]+t1*d0[i0,j1])+s1*(t0*d0[i1,j0]+t1*d0[i1,j1])
set_bnd (N, b, d)
def project(N, u, v, p, div):
h = 1.0/N
div[1:N+1,1:N+1]
= -0.5*h*(u[2:N+2,1:N+1]-u[0:N,1:N+1]+v[1:N+1,2:N+2]-v[1:N+1,0:N])
p[1:N+1,1:N+1] = 0
set_bnd (N, 0, div)
set_bnd (N, 0, p)
lin_solve (N, 0, p, div, 1, 4)
u[1:N+1,1:N+1] -= 0.5*(p[2:N+2,1:N+1]-p[0:N,1:N+1])/h # ??? not in the paper
/h
v[1:N+1,1:N+1] -= 0.5*(p[1:N+1,2:N+2]-p[1:N+1,0:N])/h # ??? not in the paper
/h
set_bnd (N, 1, u)
set_bnd (N, 2, v)
# Evolving density: advection, diffusion, addition of sources
def dens_step (N, x, x0, u, v, diff, dt):
add_source(N, x, x0, dt)
x0, x = x, x0 # swap
diffuse(N, 0, x, x0, diff, dt)
x0, x = x, x0 # swap
advect(N, 0, x, x0, u, v, dt)
# Evolving velocity: self-advection, viscous diffusion, addition of forces
def vel_step (N, u, v, u0, v0, visc, dt):
add_source(N, u, u0, dt)
add_source(N, v, v0, dt);
u0, u = u, u0 # swap
diffuse(N, 1, u, u0, visc, dt)
v0, v = v, v0 # swap
diffuse(N, 2, v, v0, visc, dt)
project(N, u, v, u0, v0)
u0, u = u, u0 # swap
v0, v = v, v0 # swap
advect(N, 1, u, u0, u0, v0, dt)
advect(N, 2, v, v0, u0, v0, dt)
project(N, u, v, u0, v0)
-----------------------------------------------------------------
