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)

-----------------------------------------------------------------