By using this site, you agree to our updated Privacy Policy and our Terms of Use. Manage your Cookies Settings.
 428,742 Members | 1,570 Online + Ask a Question
Need help? Post your question and get tips & solutions from a community of 428,742 IT Pros & Developers. It's quick & easy.

# Real-Time Fluid Dynamics for Games...

 P: n/a 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
Share this Question
2 Replies

 P: n/a 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

 P: n/a 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 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 discussion thread is closed

Replies have been disabled for this discussion. 