473,231 Members | 1,439 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,231 software developers and data experts.

Real-Time Fluid Dynamics for Games...

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)
-----------------------------------------------------------------
Jul 18 '05 #1
2 3671
Your two email addresses bouce emails back, so I post a shortened
version of my comment here.
I haven't installed:
PyOpenGL-2.0.2.01.py2.4-numpy23
glut-3.7.6
Therefore at the moment I cannot try your interesting code.
What's the speed of this Python code on your computer?
I'd like to see a screenshoot of the running Python Program...

Some people are doing in Python some things that require lots of
computations, like:
http://www.joachim-bauch.de/projects/python/pytrace

Bye,
Bearophile

Jul 18 '05 #2
You can find some screenshot in the Stam's original paper.
I didn't do any serious benchmark. I compared the speed
of C version with the Python one. It seems enough good.

I advice you to download from Stam's site paper and C code,
compile the C code and verify yourself the results.

I think the solver of Navier-Stokes equations is a piece
of cake(remember, it's patented): one page of code or less. :)

I don't like the use of global in the callback functions
of OpenGL.

--
Regards,
Alberto Santini
<be************@lycos.com> ha scritto nel messaggio
news:11*********************@z14g2000cwz.googlegro ups.com...
Your two email addresses bouce emails back, so I post a shortened
version of my comment here.
I haven't installed:
PyOpenGL-2.0.2.01.py2.4-numpy23
glut-3.7.6
Therefore at the moment I cannot try your interesting code.
What's the speed of this Python code on your computer?
I'd like to see a screenshoot of the running Python Program...

Some people are doing in Python some things that require lots of
computations, like:
http://www.joachim-bauch.de/projects/python/pytrace

Bye,
Bearophile

Jul 18 '05 #3

This thread has been closed and replies have been disabled. Please start a new discussion.

Similar topics

4
by: Aaron W. West | last post by:
Timings... sometimes there are almost too many ways to do the same thing. The only significant findings I see from all the below timings is: 1) Integer math is generally fastest, naturally....
2
by: cwdjr | last post by:
Real One has a page to copy on their site that detects if the browser of a viewer of a page has Real One installed. The page is located at...
4
by: Allan Adler | last post by:
I'm trying to reinstall RedHat 7.1 Linux on a PC that was disabled when I tried to upgrade from RH7.1 to RH9. This is presenting lots of unexpected difficulties. Apart from wanting to keep the old...
10
by: Pavils Jurjans | last post by:
Hallo, It is know issue that due to the fact that computer has to store the real numbers in limited set of bytes, thus causing a minor imprecision from the decimal value that likely was stored....
17
by: David Scemama | last post by:
Hi, I'm writing a program using VB.NET that needs to communicate with a DOS Pascal program than cannot be modified. The communication channel is through some file databases, and I have a huge...
5
by: Henry Wu | last post by:
Hi, now that in .NET everything is on millimeter, I was wondering how can one convert Pixel to Millimeter and any user's screen/monitor. I saw the following code on how to convert pixel to...
0
by: support | last post by:
Veteran Real Estate Investor Shares some of his best Insider Secrets for successful investments! www.RealEstateBeginners.ws Have you ever wondered about investing in real estate? Maybe one...
12
by: Raymond Hettinger | last post by:
I am evaluating a request for an alternate version of itertools.izip() that has a None fill-in feature like the built-in map function: >>> map(None, 'abc', '12345') # demonstrate map's None...
16
by: DirtyHarry | last post by:
Good day everyone. This sounds like a stupid question, but I became just curious yesterday, and I looked up several textbooks. However, no textbooks on computer language (that I have ) mentioned...
2
by: Tim | last post by:
Folks, Can anyone thow some clarifying light on the following? I have come across a column with the same name and same data contents defined on different tables, on some the column is defined...
0
by: VivesProcSPL | last post by:
Obviously, one of the original purposes of SQL is to make data query processing easy. The language uses many English-like terms and syntax in an effort to make it easy to learn, particularly for...
3
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 3 Jan 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). For other local times, please check World Time Buddy In...
0
by: abbasky | last post by:
### Vandf component communication method one: data sharing ​ Vandf components can achieve data exchange through data sharing, state sharing, events, and other methods. Vandf's data exchange method...
0
by: stefan129 | last post by:
Hey forum members, I'm exploring options for SSL certificates for multiple domains. Has anyone had experience with multi-domain SSL certificates? Any recommendations on reliable providers or specific...
1
by: davi5007 | last post by:
Hi, Basically, I am trying to automate a field named TraceabilityNo into a web page from an access form. I've got the serial held in the variable strSearchString. How can I get this into the...
0
by: MeoLessi9 | last post by:
I have VirtualBox installed on Windows 11 and now I would like to install Kali on a virtual machine. However, on the official website, I see two options: "Installer images" and "Virtual machines"....
0
by: DolphinDB | last post by:
The formulas of 101 quantitative trading alphas used by WorldQuant were presented in the paper 101 Formulaic Alphas. However, some formulas are complex, leading to challenges in calculation. Take...
0
by: Aftab Ahmad | last post by:
Hello Experts! I have written a code in MS Access for a cmd called "WhatsApp Message" to open WhatsApp using that very code but the problem is that it gives a popup message everytime I clicked on...
0
by: Aftab Ahmad | last post by:
So, I have written a code for a cmd called "Send WhatsApp Message" to open and send WhatsApp messaage. The code is given below. Dim IE As Object Set IE =...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.