flow of fluid solver and volumen rendering

[D.J.Peters]

videos on youtube

download include files: solver.zip


more of this stuff: http://blog.mmacklin.com/2010/11/01/adv ... imulation/

it woks with 3D arrays const N must be a power of 2 eg. 16,32,64 ...

Velocity(0, N+1 ,N ,N) ' x forces
Velocity(1, N ,N+1 ,N) ' y forces
Velocity(2 ,N ,N, N+1) ' z forces

Bondaries(N,N,N) ' if <>0 it'a a bondary cell (the fluid can't flow in this cell)

Compression(N,N,N) ' represent the density of any fluid

inside of the sub Boundary you can setup your fluid (if it is first frame)
and you can set / modify cells of Velocity every frame of simulation.

The calculation of compression and flow are only one part of the whole story.
After every step of simulation you get as result the array Compression(N,N,N)
it's the current density of your fluid per cell.
I use a minimal ray tracer as renderer and one light for shadows.
If the ray hits a cell with compression > 0 it calculates the density of 3x3x2 neighbor cells.
If the ray hits the floor a new ray in direction of light tells the renderer if it should use the density between the floor and light as shadow.

bla bla bla :-)

Joshy

here are an example of smoke it's a god idea to compile it like:
fbc -gen gcc -O 3 smoke.bas
on win32 you can enable video capturing also (optional)

Code: Select all

#include once "solver.bi"
#include once "renderer.bi"

#ifdef __FB_WIN32__
 #ifndef __FB_64BIT__
  '#define CAPTURE_VIDEO ' on win32 you can eanble video capturing 
 #endif
#endif

#ifdef CAPTURE_VIDEO
#include once "XvidCapture.bi"
#endif

type Smoke3D
  public:
  declare constructor(byval DeltaTime as single)
  declare sub Simulate(byval blnClearBoundary as integer=1)
  declare sub Render(byval frame as integer, _
                     byval eyeXPos as single=0.0, _
                     byval eyeYPos as single=0.0, _
                     byval eyeZPos as single=0.0)
  private:
  declare sub Boundary(byval blnClearVelocity as integer)
  declare sub Project
  declare sub Advection
  const as integer N3D = 32 ' cells(32,32,32)
  const as integer N3DM = N3D-1
  as Renderer3D ptr Renderer
  as Solver3D   ptr Solver
  as P3D Velocity(2)
  as P3D Boundaries
  as P3D Compression
  as P3D Divergence
  as P3D Pressure
  as integer physicsframe
  as single DeltaTime
end type

constructor Smoke3D(byval TimeStep as single)
  DeltaTime   = TimeStep
  Renderer    = new Renderer3D(N3D)
  Solver      = new Solver3D(DeltaTime)
  ' Cell Memory Allocation
  Velocity(0) = Alloc3D(N3D+1,N3D,N3D)
  Velocity(1) = Alloc3D(N3D,N3D+1,N3D)
  Velocity(2) = Alloc3D(N3D,N3D,N3D+1)
  Compression = Alloc3D(N3D,N3D,N3D)
  Boundaries  = Alloc3D(N3D,N3D,N3D)
  Divergence  = Alloc3D(N3D,N3D,N3D)
  Pressure    = Alloc3D(N3D,N3D,N3D)
end constructor

sub Smoke3D.Project
  const as single h = 1.0/N3D
  ' Compute Divergence
  for i as integer=0 to N3DM
    dim as integer ip=i+1
    for j as integer=0 to N3DM
      dim as integer jp=j+1
      for k as integer=0 to N3DM
        dim as integer kp=k+1
        Divergence[i][j][k] = (velocity(0)[ip][j ][k ]- _
                               velocity(0)[i ][j ][k ]+ _
                               velocity(1)[i ][jp][k ]- _
                               velocity(1)[i ][j ][k ]+ _
                               velocity(2)[i ][j ][kp]- _
                               velocity(2)[i ][j ][k ]) / h
      next
    next
  next

  ' Solve Pressure
  solver->solve(Pressure, Divergence, boundaries, N3D )

  ' Subtract Pressure Gradient
  for i as integer=0 to N3DM
    dim as integer im=i-1
    for j as integer=0 to N3DM
      dim as integer jm=j-1
      for k as integer=0 to N3DM
        dim as integer km=k-1
        if( i>0) then velocity(0)[i][j][k] -= (Pressure[i][j][k]-Pressure[im][j ][k ])/h
        if( j>0) then velocity(1)[i][j][k] -= (Pressure[i][j][k]-Pressure[i ][jm][k ])/h
        if( k>0) then velocity(2)[i][j][k] -= (Pressure[i][j][k]-Pressure[i ][j ][km])/h
      next
    next
  next
end sub

sub Smoke3D.Advection
  Solver->Advect( @Velocity(0), Compression, N3D)
end sub

sub Smoke3D.Render(byval frame as integer, _
                   byval eyeXPos as single, _
                   byval eyeYPos as single, _
                   byval eyeZPos as single)
  Renderer->Render(frame,compression,eyeXPos,eyeYPos,eyeZPos)
end sub


sub Smoke3D.Boundary(byval blnClearVelocity as integer)
  if blnClearVelocity then

    for i as integer=0 to N3D
      for j as integer=0 to N3DM
        for k as integer=0 to N3DM
          ' x-velocity left and right wall 
          if (i=0 or i=N3D) then Velocity(0)[i][j][k] = 0.0
        next
      next
    next

    for i as integer=0 to N3DM
      for j as integer=0 to N3D
        for k as integer=0 to N3DM
          ' y-velocity bottom and top wall 
          if (j=0 or j=N3D) then Velocity(1)[i][j][k] = 0.0
        next
      next
    next

    for i as integer=0 to N3DM
      for j as integer=0 to N3DM
        for k as integer=0 to N3D
          ' z-velocity near and far wall
          if(k=0 or k=N3D) then Velocity(2)[i][j][k] = 0.0
        next
      next
    next
  end if
  
  
  '##############
  '# init smoke #
  '##############
  ' if first frame set some cells with initial Smoke
  if( physicsframe < 1 ) then
    for x as integer=1 to N3DM
      for y as integer=N3D\2 to N3D\2
        for z as integer =1 to N3DM
          'Velocity(1)[x][y][z] = 0
          Compression[x][y][z] = 1
        next
      next
    next
  end if

  ' Set cells marked as boundariy to 0
  for i as integer=0 to N3DM
    dim as integer ip=i+1
    for j as integer=0 to N3DM
      dim as integer jp=j+1
      for k as integer=0 to N3DM
        if( boundaries[i][j][k] )=1 then
          dim as integer kp=k+1
          Compression[i ][j ][k ] = 0.0
          Velocity(0)[i ][j ][k ] = 0.0 
          Velocity(1)[i ][j ][k ] = 0.0
          Velocity(2)[i ][j ][k ] = 0.0
          Velocity(0)[ip][j ][k ] = 0.0 
          Velocity(1)[i ][jp][k ] = 0.0
          Velocity(2)[i ][j ][kp] = 0.0
        end if
      next
    next
  next

  ' Give some external forces
  for x as integer=1 to N3DM
    for y as integer=1 to N3DM
      for z as integer=1 to N3DM
        if boundaries[x][y][z] then continue for
        dim as single c=Compression[x][y][z]
        'if c<0.1 then
        '  Compression[x][y][z]=c*0.99
        'end if
        if x<N3D\2 then 
          ' left side +y velocty
          Velocity(1)[x][y][z] += .01
        else  
          ' right side -y velocty
          Velocity(1)[x][y][z] -= .01
        end if  
      next
    next
  next
end sub

sub smoke3D.Simulate(byval blnClearBoundary as integer)
  windowtitle "phy " & physicsframe
  boundary(blnClearBoundary)
  project()
  advection()
  physicsframe+=1
end sub

dim as integer w=512
dim as integer h=w
screenres w,h,32

#ifdef CAPTURE_VIDEO
' create video stream with 25 FPS
var avi = OpenVideoStream("smoke.avi",w,h,25,100)
#endif

dim as integer frame
dim as Smoke3D Smoke = Smoke3D(0.1)
while Inkey()="" andalso frame<24*30
  Smoke.Simulate(0)
  #ifdef CAPTURE_VIDEO
  ' render every frame
  Smoke.Render frame
  ' capture frame
  if avi then WriteRGBScreen
  #else
  ' no video render only some frames
  if frame mod 8=0 then
    Smoke.Render frame
  end if 
  #endif
  frame+=1
wend
#ifdef CAPTURE_VIDEO
' important close the vide stream
if avi then CloseStream
#else
BSave "tmp.bmp",0
#endif
Windowtitle "ok ..."
sleep

[D.J.Peters]

some kind of smoothing is done via spline and a cubic function I found this C code

Code: Select all

static FLOAT spline_cubic(const FLOAT a[4], FLOAT x) {
	int i, j;
	FLOAT alpha[4], l[4], mu[4], z[4];
	FLOAT b[4], c[4], d[4];
	for(i = 1; i < 3; i++) {
		alpha[i] = 3.0 * (a[i+1] - a[i]) - 3.0 * (a[i] - a[i-1]);
	}

	l[0] = 1.0;
	mu[0] = 0.0;
	z[0] = 0.0;
	for(i = 1; i < 3; i++) {
		l[i] = 4.0 - mu[i-1];
		mu[i] = 1.0 / l[i];
		z[i] = (alpha[i] - z[i-1]) / l[i];
	}
	l[3] = 1.0;
	z[3] = 0.0;
	c[3] = 0.0;
	for(j = 2; 0 <= j; j--) {
		c[j] = z[j] - mu[j] * c[j+1];
		b[j] = a[j+1] - a[j] - (c[j+1] + 2.0 * c[j]) / 3.0;
		d[j] = (c[j+1] - c[j]) / 3.0;
	}

	FLOAT minv = min(a[1],a[2]);
	FLOAT maxv = max(a[2],a[1]);
	return min(maxv,max(minv,(a[1] + b[1] * x + c[1] * x * x + d[1] * x * x * x )));
}

this part of code are used very often and is slow
so I optimized it like a machine :-)
Do you believe this FreeBASIC code (without any loop and local arrays) does the same job 10 times faster ?

Code: Select all

function Solver3D.Cubic(byval a as const single ptr,byval x as single) as single
  dim as single xx=x*x
  dim as single xxx=xx*x
  dim as single p213=3.0 * (a[2] - a[1])
  dim as single a1 = p213 - 3.0 * (a[1] - a[0])
  dim as single a2 = 3.0 * (a[3] - a[2]) - p213
  dim as single z1 = a1 * 0.25
  dim as single z2 = (a2 - z1) * 0.266666666
  dim as single c1 = z1 - 0.25 * z2
  dim as single b1 = a[2] - a[1] - (z2 + 2*c1)*0.333333333
  dim as single d1 = (z2 - c1) * 0.333333333
  dim as single minv = min(a[1],a[2])
  dim as single maxv = max(a[2],a[1])
  dim as single mami = a[1] + b1*x + c1*xx + d1*xxx 
  return min(maxv,max(minv,mami))
end function

[Dr_D]

It's very cool. I thought it was real time. hehehehe... I was about to think something unthinkable. :p

[badidea]

Nice, I always wanted to do something like this, but seemed so complicated.

Is this a Navier-Stokes equations solver, or a more simple simulation? How accurate is this for simulating e.g. water flow in a pipe?

[D.J.Peters]

Here are flow of fluid inside a 2D grid only.
right mouse = add compression -> density
left mouse = add forces -> velocity
[v] = toggle show velocity/compression
[q] or [esc] = quit

Joshy

Code: Select all

' based on paper: "Real-Time Fluid Dynamics for Games" by Jos Stam.
' optimized version by d.j.peters aka. Joshy

#include once "fbgfx.bi"
#include once "GL/gl.bi"

' number of 2D Cells[N][N]
const as integer N=64
const as single  INVN = 1.0/N
' delata time 
const as single  DT=0.1
const as single  DTN=DT*N
const as single  DTNN=DTN*N

' N+2 rows, N+2 colums
const as integer ROWSIZE=N+2
const as integer NITEMS=ROWSIZE*ROWSIZE
 
' X. = X*DT
sub AddSource(byval Value as single ptr ptr,byval ValuePreview as single ptr ptr)
  dim as integer i,j
  for i=0 to ROWSIZE-1
    for j=0 to ROWSIZE-1
      Value[i][j] += ValuePreview[i][j]*DT
    next
  next
end sub

' boundaries top,bottom rows left,right columns and the four corners
sub SetBoundaries(byval Bound as integer,byval Cells as single ptr ptr)
  dim as integer i
  ' rows collumns
  if Bound=0 then
    for i=1 to N
      Cells[0  ][i] = Cells[1][i] ' top
      Cells[N+1][i] = Cells[N][i] ' bottom
    next
    for i=1 to N
      Cells[i][0  ] = Cells[i][1] ' left
      Cells[i][N+1] = Cells[i][N] ' right
    next
  elseif Bound=1 then
    for i=1 to N
      Cells[0  ][i] =-Cells[1][i] ' -top
      Cells[N+1][i] =-Cells[N][i] ' -bottom
    next
  elseif Bound=2 then
    for i=1 to N
      Cells[i][0  ] =-Cells[i][1] ' -left
      Cells[i][N+1] =-Cells[i][N] ' -right
    next
  end if
  ' the four corners
  Cells[0  ][0  ] = 0.5*(Cells[1][0  ]+Cells[0  ][1])
  Cells[0  ][N+1] = 0.5*(Cells[1][N+1]+Cells[0  ][N])
  Cells[N+1][0  ] = 0.5*(Cells[N][0  ]+Cells[N+1][1])
  Cells[N+1][N+1] = 0.5*(Cells[N][N+1]+Cells[N+1][N])
end sub

'     (X + a*(neighbors cells))
' X = -------------------------
'               c
sub Smooth(byval Bound as integer,byval x as single ptr ptr, byval x0 as single ptr ptr,byval a as single,byval c as single)
  dim as integer i,ip,im,j,jm,jp, k
  c=1.0/c
  if a=1.0 then
    for k=1 to 20
      im=0
      for i=1 to N
        ip=i+1:jm=0
        for j=1 to N
          jp=j+1
          x[i][j] = (x0[i][j] + (x[im][j] + _ 
                                 x[ip][j] + _
                                 x[i][jm] + _
                                 x[i][jp])    )*c
          jm=j
        next
        im=i
      next
      SetBoundaries(Bound,x )
    next
  elseif a=0 then
    for k=1 to 20
      im=0
      for i=1 to N
        ip=i+1:jm=0
        for j=1 to N
          jp=j+1
          x[i][j] = x0[i][j]*c
          jm=j
        next
        im=i
      next
      SetBoundaries(Bound,x )
    next
  else
    for k=1 to 20
      im=0
      for i=1 to N
        ip=i+1:jm=0
        for j=1 to N
          jp=j+1
          x[i][j] = (x0[i][j] + a * (x[im][j] + _ 
                                     x[ip][j] + _
                                     x[i][jm] + _
                                     x[i][jp])    )*c
          jm=j
        next
        im=i
      next
      SetBoundaries(Bound,x )
    next
  end if
end sub

sub Diffuse(byval b as integer, byval x as single ptr ptr, byval x0 as single ptr ptr,byval diffusion as single)
  dim as single a=DTNN*diffusion
  Smooth (b, x, x0,a, 1+4*a)
end sub

sub Advect(byval b as integer,byval dOut as single ptr ptr,byval dIn as single ptr ptr,byval u as single ptr ptr,byval v as single ptr ptr)
  dim as integer i, j, i0, j0, i1, j1
  dim as single x, y, s0, t0, s1, t1
  for i=1 to N
    for j=1 to N
      x = i-DTN*u[i][j] 
      y = j-DTN*v[i][j]
      if (x<0.5f) then 
        x=0.5f 
      elseif (x>N+0.5f) then 
        x=N+0.5f 
      end if 
      if (y<0.5f) then 
        y=0.5f 
      elseif (y>N+0.5f) then 
        y=N+0.5f
      end if
      i0=int(x) : i1=i0+1
      s1=x - i0 : s0=1-s1
      j0=int(y) : j1=j0+1
      t1=y - j0 : t0=1-t1
      dOut[i][j] = s0*( t0*dIn[i0][j0]+ t1*dIn[i0][j1]) _
              + s1*( t0*dIn[i1][j0]+ t1*dIn[i1][j1])
    next
  next
  SetBoundaries(b,dOut)
end sub

sub Project(byval u as single ptr ptr,byval v as single ptr ptr,byval p as single ptr ptr,byval div as single ptr ptr)
  dim as integer i, j
  for i=1 to N
    dim as integer ip=i+1
    dim as integer im=i-1
    for j=1 to N
      div[i][j] = -0.5f*(u[ip][j]-u[im][j]+v[i][j+1]-v[i][j-1])/N
      p[i][j] = 0
    next
  next
  SetBoundaries(0, div )
  SetBoundaries(0, p )
  Smooth(0, p, div, 1, 4 )

  for i=1 to N
    dim as integer ip=i+1
    dim as integer im=i-1
    for j=1 to N
      u[i][j] -= 0.5f*N*(p[ip][j]-p[im][j])
      v[i][j] -= 0.5f*N*(p[i][j+1]-p[i][j-1])
    next
  next
  SetBoundaries (1, u )
  SetBoundaries (2, v )
end sub

sub SmoothCompression(byval x as single ptr ptr, byval x0 as single ptr ptr,byval diffusion as single)
  Smooth (0, x, x0,0, 1+4*diffusion)
end sub

sub UpdateCompression(byval Compression as single ptr ptr, byval CompressionPreview as single ptr ptr, _
                      byval xVelocity   as single ptr ptr, byval yVelocity          as single ptr ptr, _
                      byval Diffusion   as single)
  AddSource(Compression,CompressionPreview)
  SmoothCompression(CompressionPreview,Compression, Diffusion)
  Advect   (0, Compression,CompressionPreview, xVelocity, yVelocity)
end sub

sub UpdateVelocity(byval xVelocity        as single ptr ptr,byval yVelocity        as single ptr ptr, _
                   byval xVelocityPreview as single ptr ptr,byval yVelocityPreview as single ptr ptr, _
                   byval Viscosity        as single)
  AddSource(xVelocity, xVelocityPreview)
  AddSource(yVelocity ,yVelocityPreview)
  Diffuse(1,xVelocityPreview, xVelocity, Viscosity)
  Diffuse(2,yVelocityPreview, yVelocity, Viscosity)
  Project(xVelocityPreview, yVelocityPreview, xVelocity, yVelocity)
  Advect(1, xVelocity, xVelocityPreview, xVelocityPreview, yVelocityPreview) 
  Advect(2, yVelocity, yVelocityPreview, xVelocityPreview, yVelocityPreview)
  Project(xVelocity, yVelocity, xVelocityPreview, yVelocityPreview)
end sub




' global variables

dim shared as integer gShowvelocity
dim shared as single gDiffusion, gViscosity
dim shared as single gForce, gComp

dim shared as single ptr ptr gVelocityX, gVelocityY, gCompression
dim shared as single ptr ptr gVelocityXPreview, gVelocityYPreview, gCompressionPreview
dim shared as integer gWinWidth, gWinHeight
dim shared as integer gMouseButtons(2)
dim shared as integer gMouseXOld, gMouseYOld, gMouseX, gMouseY

' free/clear/allocate 2D cells
#macro free2d(XX)
  if XX then
    for ii as integer=0 to ROWSIZE-1
      if XX[ii] then deallocate XX[ii]
    next
    deallocate XX
  end if
#endmacro

#macro clear2d(XX)
  if XX then
    for ii as integer=0 to ROWSIZE-1
      for jj as integer=0 to ROWSIZE-1
        XX[ii][jj]=0
      next
    next
  end if
#endmacro

#macro alloc2d(XX) 
  XX = allocate(ROWSIZE*sizeof(single ptr) )
  for ii as integer=0 to ROWSIZE-1
    XX[ii] = callocate(ROWSIZE*sizeof(single) )
  next
#endmacro

sub FreeCells
  free2d(gVelocityX)
  free2d(gVelocityY)
  free2d(gCompression)
  free2d(gVelocityXPreview)
  free2d(gVelocityYPreview)
  free2d(gCompressionPreview)
end sub

sub ClearCells
  clear2d(gVelocityX)
  clear2d(gVelocityY)
  clear2d(gCompression)
  clear2d(gVelocityXPreview)
  clear2d(gVelocityYPreview)
  clear2d(gCompressionPreview)
end sub

sub AllocCells
  alloc2d(gVelocityX)
  alloc2d(gVelocityY)
  alloc2d(gCompression)
  alloc2d(gVelocityXPreview)
  alloc2d(gVelocityYPreview)
  alloc2d(gCompressionPreview)
end sub

sub DrawVelocity
  dim as integer i,j
  dim as single x,y,u,v,r,g,b
  glBegin ( GL_LINES )
  glColor3f (1,1,1)
  for i = 1 to N
    x = (i-0.5)*INVN
    for j = 1 to N
      y = (j-0.5)*INVN
      u=gVelocityX[i][j]*5
      v=gVelocityY[i][j]*5
      
      glVertex2f ( x, y )
      glVertex2f ( x+u, y+v)
    next
  next
  glEnd ()

end sub

sub DrawDensity
  dim as integer i,ip,j,jp
  dim as single x,xh,y,yh,c
  glClear(GL_COLOR_BUFFER_BIT)
  
  glBegin ( GL_QUADS )
  x=0:xh=INVN
  for i=0 to N
    ip=i+1 : y=0: yh=INVN
    for j=0 to N
      jp=j+1
      c = gCompression[i ][j ]:glColor3f(c,c,c):glVertex2f(x , y )
      c = gCompression[i ][jp]:glColor3f(c,c,c):glVertex2f(x , yh)
      c = gCompression[ip][jp]:glColor3f(c,c,c):glVertex2f(xh, yh)
      c = gCompression[ip][j ]:glColor3f(c,c,c):glVertex2f(xh,  y)
      y=yh:yh+=INVN
    next
    x=xh:xh+=INVN
  next
  glEnd ()
end sub

sub UpdateInput
  clear2d(gVelocityXPreview)
  clear2d(gVelocityYPreview)
  clear2d(gCompressionPreview)

  if ( gMouseButtons(0)=0 andalso gMouseButtons(1)=0) then return

  if gMouseX<1 then
    gMouseX=1
  elseif gMouseX>gWinWidth-1 then
    gMouseX=gWinWidth-1
  end if

  if gMouseY<1 then
    gMouseY=1
  elseif gMouseY>gWinHeight-1 then
    gMouseY=gWinHeight-1
  end if

  dim as integer i = gMouseX/gWinWidth  * (ROWSIZE-1)
  dim as integer j = (gWinHeight-gMouseY)/gWinHeight * (ROWSIZE-1)

  if ( gMouseButtons(0) ) then ' left button set force
    gVelocityXPreview[i][j] = gForce * (gMouseX-gMouseXOld)
    gVelocityYPreview[i][j] = gForce * (gMouseYOld-gMouseY)
  end if

  if ( gMouseButtons(1) ) then ' right button set compression
    gCompressionPreview[i][j] = gComp
  end if
  gMouseXOld = gMouseX
  gMouseYOld = gMouseY
end sub




sub UpdateKeyboard
  dim as string key=inkey()
  if len(key)<1 then return
  if key="c" then ClearCells
  if key="q" or asc(key)=27 then FreeCells:end
  if key="v" then gshowvelocity= not gshowvelocity
end sub

sub UpdateMouse
  static as integer omb=0
  dim as integer nx,ny,nb
  if getmouse(nx,ny,,nb)=0 then
    if omb=nb then
      ' mouse motion
      gMouseX=nx
      gMouseY=ny
    else
      ' a mouse button change
      gMouseButtons(0) = iif(nb and 1,1,0)
      gMouseButtons(1) = iif(nb and 2,1,0) 
      gMouseXOld=nx
      gMouseYOld=ny
      gMouseX=nx
      gMouseY=ny
      omb = nb
    end if
  end if 
end sub

sub UpdateDisplay
  glClear ( GL_COLOR_BUFFER_BIT )
  if (gshowvelocity) then 
    DrawVelocity
  else
    DrawDensity
  end if
  flip
end sub

sub UpdateSimulation
  UpdateMouse
  UpdateKeyboard
  UpdateInput
  UpdateVelocity(gVelocityX, gVelocityY, gVelocityXPreview, gVelocityYPreview, gviscosity)
  UpdateCompression(gCompression, gCompressionPreview, gVelocityX, gVelocityY, gdiffusion)
  UpdateDisplay
end sub


'
' main
'
screeninfo gWinWidth,gWinHeight
gWinWidth\=2
gWinHeight=gWinWidth
'gWinWidth=512
'gWinHeight=512
screenres gWinWidth,gWinHeight,32,,fb.GFX_OPENGL
windowtitle  "RM=compression LM=forces [v]=show velocity [ESC]=quit"
glViewport ( 0, 0, gWinWidth, gWinHeight )
glMatrixMode ( GL_PROJECTION )
glLoadIdentity ()
glOrtho(0,1, 0,1, -1,1)
glClearColor(0,0,0,1)
glClear ( GL_COLOR_BUFFER_BIT )
glDisable(GL_DEPTH_TEST)
flip


gDiffusion = 0.001 ' compression smoothing (how fast compression fades out)
gViscosity = 0.001 ' velocity smoothing (how fast velocity comes at rest)
gForce     = 5.0
gComp      = 100.0
AllocCells

do
  UpdateSimulation
  sleep 10
loop

[tinram]

Really good Joshy.

The stirring smoothness and speed are impressive. OpenGL is a good move.

(I believe Duke4e did a Jos Stam conversion in the forum a few years ago, which was good, but pixelated.)

Processing, usually fast at graphics, can struggle with this sort of thing:
http://www.openprocessing.org/sketch/197