## flow of fluid solver and volumen rendering

D.J.Peters
Posts: 7852
Joined: May 28, 2005 3:28

### flow of fluid solver and volumen rendering

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"#endiftype 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 DeltaTimeend typeconstructor 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 constructorsub 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  nextend subsub Smoke3D.Advection  Solver->Advect( @Velocity(0), Compression, N3D)end subsub 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 subsub 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  nextend subsub smoke3D.Simulate(byval blnClearBoundary as integer)  windowtitle "phy " & physicsframe  boundary(blnClearBoundary)  project()  advection()  physicsframe+=1end subdim as integer w=512dim as integer h=wscreenres w,h,32#ifdef CAPTURE_VIDEO' create video stream with 25 FPSvar avi = OpenVideoStream("smoke.avi",w,h,25,100)#endifdim as integer framedim 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+=1wend#ifdef CAPTURE_VIDEO' important close the vide streamif avi then CloseStream#elseBSave "tmp.bmp",0#endifWindowtitle "ok ..."sleep`
Last edited by D.J.Peters on Mar 17, 2015 17:30, edited 1 time in total.
D.J.Peters
Posts: 7852
Joined: May 28, 2005 3:28

### Re: flow of fluid solver and volumen rendering

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
Posts: 2391
Joined: May 27, 2005 4:59
Contact:

### Re: flow of fluid solver and volumen rendering

It's very cool. I thought it was real time. hehehehe... I was about to think something unthinkable. :p
Posts: 1608
Joined: May 24, 2007 22:10
Location: The Netherlands

### Re: flow of fluid solver and volumen rendering

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
Posts: 7852
Joined: May 28, 2005 3:28

### Re: flow of fluid solver and volumen rendering

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=64const as single  INVN = 1.0/N' delata time const as single  DT=0.1const as single  DTN=DT*Nconst as single  DTNN=DTN*N' N+2 rows, N+2 columsconst as integer ROWSIZE=N+2const as integer NITEMS=ROWSIZE*ROWSIZE ' X. = X*DTsub 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  nextend sub' boundaries top,bottom rows left,right columns and the four cornerssub 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 = -------------------------'               csub 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 ifend subsub 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 subsub 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 subsub 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 subsub 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 subsub 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 subsub 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 variablesdim shared as integer gShowvelocitydim shared as single gDiffusion, gViscositydim shared as single gForce, gCompdim shared as single ptr ptr gVelocityX, gVelocityY, gCompressiondim shared as single ptr ptr gVelocityXPreview, gVelocityYPreview, gCompressionPreviewdim shared as integer gWinWidth, gWinHeightdim 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#endmacrosub FreeCells  free2d(gVelocityX)  free2d(gVelocityY)  free2d(gCompression)  free2d(gVelocityXPreview)  free2d(gVelocityYPreview)  free2d(gCompressionPreview)end subsub ClearCells  clear2d(gVelocityX)  clear2d(gVelocityY)  clear2d(gCompression)  clear2d(gVelocityXPreview)  clear2d(gVelocityYPreview)  clear2d(gCompressionPreview)end subsub AllocCells  alloc2d(gVelocityX)  alloc2d(gVelocityY)  alloc2d(gCompression)  alloc2d(gVelocityXPreview)  alloc2d(gVelocityYPreview)  alloc2d(gCompressionPreview)end subsub 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 subsub 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 subsub 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 = gMouseYend subsub 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 gshowvelocityend subsub 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 subsub UpdateDisplay  glClear ( GL_COLOR_BUFFER_BIT )  if (gshowvelocity) then     DrawVelocity  else    DrawDensity  end if  flipend subsub UpdateSimulation  UpdateMouse  UpdateKeyboard  UpdateInput  UpdateVelocity(gVelocityX, gVelocityY, gVelocityXPreview, gVelocityYPreview, gviscosity)  UpdateCompression(gCompression, gCompressionPreview, gVelocityX, gVelocityY, gdiffusion)  UpdateDisplayend sub'' main'screeninfo gWinWidth,gWinHeightgWinWidth\=2gWinHeight=gWinWidth'gWinWidth=512'gWinHeight=512screenres gWinWidth,gWinHeight,32,,fb.GFX_OPENGLwindowtitle  "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)flipgDiffusion = 0.001 ' compression smoothing (how fast compression fades out)gViscosity = 0.001 ' velocity smoothing (how fast velocity comes at rest)gForce     = 5.0gComp      = 100.0AllocCellsdo  UpdateSimulation  sleep 10loop`
tinram
Posts: 88
Joined: Nov 30, 2006 13:35
Location: UK

### Re: flow of fluid solver and volumen rendering

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