Fluid simulation

[threedslider]

Hello all

I was writing for fluid simulaiton based on thttps://github.com/983/Fluid-Simulation in Purebasic but only in CPU not GPU... I think I have done from C++ to Purebasic quite good however it doesn't work very well in gl texture 2D as transparency so please checck out this C++ to compare with my source code

Here is my source code inspired from C++ :

Code: Select all

;;; 
;;; Fluid Simulation v. 0.1 
;;;

w.i = 512
h.i = 512

nx.i = 32
ny.i = 32

dt.f = 0.02
iterations.i = 5
vorticity.f = 10.0

Global texture 

Structure vec2f
  x.f
  y.f
EndStructure

Myvar.vec2f

Structure vecvalue
  Array values.vec2f(32*32)
EndStructure

Procedure idx(x.i, y.i) 
  
  Shared nx
  Shared ny
  ; wrap around
  xd = (x + nx) % nx
  yd = (y + ny) % ny
  ProcedureReturn xd + yd*nx  
  
EndProcedure
      

Structure grid
  x.i
  y.i
EndStructure

Structure gridvalue
  Array values.grid(32*32)
EndStructure


Dim old_velocity.vec2f(nx*ny)
Dim new_velocity.vec2f(nx*ny)

Dim old_density.grid(nx*ny)
Dim new_density.grid(nx*ny)

Dim pixels.grid(nx*ny)

Dim pixels_data(nx*ny)


Procedure draw(Array var.vec2f(1), n.i, mode.i)
  glEnableClientState_(#GL_VERTEX_ARRAY)
  glVertexPointer_(2, #GL_FLOAT, 0, var())
  glDrawArrays_(mode, 0, n)
EndProcedure



Procedure draw_circle(x.f, y.f, r.f, n.i = 100)
  Dim pos.vec2f(n)
  
  MyVar.vec2f
  
  For i = 1 To n
    angle.f = 2.0*3.14159*i/n
    MyVar\x = x
    MyVar\y = y
    pos(i)\x = MyVar\x + r*Cos(angle)
    pos(i)\y = MyVar\y + r*Sin(angle)
  Next
  
  draw(pos(), n, #GL_LINE_LOOP)
  
  FreeArray(pos())
  
EndProcedure
  

Procedure init()
  
  Shared nx
  Shared ny
  
  Shared texture
  
  Shared old_density()
  Shared old_velocity()
  
  Shared Myvar.vec2f
  
  Myvar\x = 0.0
  Myvar\y = 0.0
  
  For y = 0 To ny-1
    For n = 0 To nx-1
      
      old_density(idx(n,y))\x = 0
      old_density(idx(n,y))\y = 0
      old_velocity(idx(n,y))\x = Myvar\x
      old_velocity(idx(n,y))\y = Myvar\y
        
    Next
  Next
  

    glEnable_(#GL_TEXTURE_2D)
    glEnable_(#GL_BLEND)
    glBlendFunc_(#GL_SRC_ALPHA, #GL_ONE_MINUS_SRC_ALPHA)

    glGenTextures_(1, @texture)
    glBindTexture_(#GL_TEXTURE_2D, texture)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MIN_FILTER, #GL_LINEAR)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MAG_FILTER, #GL_LINEAR)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_S, #GL_REPEAT)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_T, #GL_REPEAT)
    glTexImage2D_(#GL_TEXTURE_2D, 0, #GL_RGBA, nx, ny, 0, #GL_RGBA, #GL_UNSIGNED_BYTE, #NUL)
      
EndProcedure

Procedure.f lerp(a.f, b.f, c.f)
  ProcedureReturn (1.0 - c)*a + c*b
EndProcedure

Procedure.f interpolationX_d(Array g.grid(1), *var.vec2f, ax.i, ay.i)
  ix.i = *var\x 
  iy.i = *var\y
  
  ux.f = *var\x - ix
  uy.f = *var\y - iy
  
 ProcedureReturn lerp(lerp(g(idx(ax+0,ay+0))\x, g(idx(ax+1,ay+0))\x, ux), lerp(g(idx(ax+0,ay+1))\x, g(idx(ax+1,ay+1))\x + 0, ux), uy)
  
EndProcedure

Procedure.f interpolationY_d(Array g.grid(1), *var.vec2f, ax.i, ay.i)
  ix.i = *var\x 
  iy.i = *var\y
  
  ux.f = *var\x - ix
  uy.f = *var\y - iy
  
 ProcedureReturn lerp(lerp(g(idx(ax+0,ay+0))\y, g(idx(ax+1,ay+0))\y, ux), lerp(g(idx(ax+0,ay+1))\y, g(idx(ax+1,ay+1))\y + 0, ux), uy)
  
EndProcedure

Procedure.f interpolationX_v(Array g.vec2f(1), *var.vec2f, ax.i, ay.i)
  ix.i = *var\x 
  iy.i = *var\y
  
  ux.f = *var\x - ix
  uy.f = *var\y - iy
  
 ProcedureReturn lerp(lerp(g(idx(ax+0,ay+0))\x, g(idx(ax+1,ay+0))\x, ux), lerp(g(idx(ax+0,ay+1))\x, g(idx(ax+1,ay+1))\x + 0, ux), uy)
  
EndProcedure

Procedure.f interpolationY_v(Array g.vec2f(1), *var.vec2f, ax.i, ay.i)
  ix.i = *var\x 
  iy.i = *var\y
  
  ux.f = *var\x - ix
  uy.f = *var\y - iy
  
 ProcedureReturn lerp(lerp(g(idx(ax+0,ay+0))\y, g(idx(ax+1,ay+0))\y, ux), lerp(g(idx(ax+0,ay+1))\y, g(idx(ax+1,ay+1))\y + 0, ux), uy)
  
 EndProcedure

Procedure swap_od_2_nd_x(Array od.grid(1), Array nd.grid(1), ax.i, ay.i)
  Swap od(idx(ax,ay))\x, nd(idx(ax,ay))\x
EndProcedure

Procedure swap_od_2_nd_y(Array od.grid(1), Array nd.grid(1), ax.i, ay.i)
  Swap od(idx(ax,ay))\y, nd(idx(ax,ay))\y
EndProcedure

Procedure swap_ov_2_nv_x(Array ov.vec2f(1), Array nv.vec2f(1), ax.i, ay.i)
  Swap ov(idx(ax,ay))\x, nv(idx(ax,ay))\x
EndProcedure

Procedure swap_ov_2_nv_y(Array ov.vec2f(1), Array nv.vec2f(1), ax.i, ay.i)
  Swap ov(idx(ax,ay))\y, nv(idx(ax,ay))\y
EndProcedure

 
 
 Procedure advect_density()
   Shared nx
   Shared ny
   
   Shared old_density()
   Shared new_density()
   Shared old_velocity()
   
   Shared Myvar.vec2f
   
  
  
   pos.vec2f
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       pos\x = Myvar\x - dt* old_velocity(idx(x,y))\x
       pos\y = Myvar\y - dt* old_velocity(idx(x,y))\y
       
       new_density(idx(x,y))\x = interpolationX_d(old_density(), @pos, x,y)
       new_density(idx(x,y))\y = interpolationY_d(old_density(), @pos, y,y)
     Next
   Next
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       swap_od_2_nd_x(old_density(), new_density(), x,y)
       swap_od_2_nd_y(old_density(), new_density(), x,y)
     Next 
   Next   
 EndProcedure
 
 
  Procedure advect_velocity()
   Shared nx
   Shared ny
   
   Shared old_velocity()
   Shared old_density()
   
   Shared new_density()
   Shared new_velocity()
   
   Shared Myvar.vec2f
   
   pos.vec2f
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       pos\x = Myvar\x - dt* old_velocity(idx(x,y))\x
       pos\y = Myvar\y - dt* old_velocity(idx(x,y))\y
       
       new_velocity(idx(x,y))\x = interpolationX_v(old_velocity(), @pos, x,y)
       new_velocity(idx(x,y))\y = interpolationY_v(old_velocity(), @pos, x,y)
     Next
   Next
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       swap_ov_2_nv_x(old_velocity(), new_velocity(), x,y)
       swap_ov_2_nv_y(old_velocity(), new_velocity(), x,y)
     Next 
   Next   
 EndProcedure
 
 
 Procedure diffuse_density()
   Shared dt
   Shared nx
   Shared ny
   
   Shared old_density()
   Shared new_density()
   
   sumX = 0
   sumY = 0
   
   diffusion.f = dt * 100.01
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       
       sumX = diffusion * ( old_density(idx(x-1,y+0))\x + old_density(idx(x+1,y+0))\x + old_density(idx(x+0,y-1))\x + old_density(idx(x+0,y+1))\x ) + old_density(idx(x+0,y+0))\x 
       sumY = diffusion * ( old_density(idx(x-1,y+0))\y + old_density(idx(x+1,y+0))\y + old_density(idx(x+0,y-1))\y + old_density(idx(x+0,y+1))\y ) + old_density(idx(x+0,y+0))\y 
       
       new_density(idx(x,y))\x = 1.0 / (1.0 + 4.0 * viscosity) * sumX
       new_density(idx(x,y))\y = 1.0 / (1.0 + 4.0 * viscosity) * sumY
     Next
   Next
   
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       swap_od_2_nd_x(old_density(), new_density(), x,y)
       swap_od_2_nd_y(old_density(), new_density(), x,y)
     Next
   Next
   
 EndProcedure
 
 
  Procedure diffuse_velocity()
   Shared dt
   Shared nx
   Shared ny
   
   Shared old_velocity()
   Shared new_velocity()
   
   sum.vec2f
   
   viscosity.f = dt * 0.000001
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       
       sum\x = viscosity * ( old_velocity(idx(x-1,y+0))\x +  old_velocity(idx(x+1,y+0))\x +  old_velocity(idx(x+0,y-1))\x +  old_velocity(idx(x+0,y+1))\x ) +  old_velocity(idx(x+0,y+0))\x 
       sum\y = viscosity * ( old_velocity(idx(x-1,y+0))\y +  old_velocity(idx(x+1,y+0))\y +  old_velocity(idx(x+0,y-1))\y +  old_velocity(idx(x+0,y+1))\y ) +  old_velocity(idx(x+0,y+0))\y 
       
       new_velocity(idx(x,y))\x = 1.0 / (1.0 + 4.0 * viscosity) * sum\x
       new_velocity(idx(x,y))\y = 1.0 / (1.0 + 4.0 * viscosity) * sum\y
     Next
   Next
   
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       swap_ov_2_nv_x(old_velocity(), new_velocity(), x,y)
       swap_ov_2_nv_y(old_velocity(), new_velocity(), x,y)
       
     Next
   Next
   
 EndProcedure    
 
 
 Procedure project_velocity()
   
   Shared nx
   Shared ny
   
   Shared iterations
   
   Shared old_velocity()
   
   Dim p.grid(nx*ny)
   Dim p2.grid(nx*ny)
   Dim div.grid(nx*ny)
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       dx.f = old_velocity(idx(x+1,y+0))\x - old_velocity(idx(x-1,y+0))\x
       dy.f = old_velocity(idx(x+0,y+1))\y - old_velocity(idx(x+0,y-1))\y
       
       div(idx(x,y))\x = dx + dy
       div(idx(x,y))\y = dx + dy
       
       p(idx(x,y))\x = 0.0
       p(idx(x,y))\y = 0.0
     Next
   Next
   
   
   For k = 1 To iterations
     
     For y = 0 To ny-1
       For x = 0 To nx-1
         sumX.f = -div(idx(x,y))\x + p(idx(x+1,y+0))\x + p(idx(x-1,y+0))\x + p(idx(x+0,y+1))\x + p(idx(x+0,y-1))\x 
         sumY.f = -div(idx(x,y))\y + p(idx(x+1,y+0))\y + p(idx(x-1,y+0))\y + p(idx(x+0,y+1))\y + p(idx(x+0,y-1))\y
         
         p2(idx(x,y))\x = 0.25 * sumX
         p2(idx(x,y))\y = 0.25 * sumY
       Next
     Next
     
     For y = 0 To ny-1
       For x = 0 To nx-1
         swap_od_2_nd_x(p(), p2(), x,y)
         swap_od_2_nd_y(p(), p2(), x,y)
       Next
     Next
     
   Next
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       old_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x - 0.5 * p(idx(x+1,y+0))\x - p(idx(x-1,y+0))\x
       old_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y - 0.5 * p(idx(x+0,y+1))\y - p(idx(x+0,y-1))\y
     Next
   Next
   
   
   FreeArray(p())
   FreeArray(p2())
   FreeArray(div())
 EndProcedure
 
 
 Procedure.f curl(ax.i, ay.i)
   
   Shared old_velocity()
   
   ProcedureReturn old_velocity(idx(ax,ay+1))\x + old_velocity(idx(ax,ay-1))\x + old_velocity(idx(ax-1,ay))\y + old_velocity(idx(ax+1,ay))\y
   
 EndProcedure
 
 Procedure.f dot(*A_v.vec2f, *B_v.vec2f)
   ProcedureReturn *A_v\x * *B_v\x + *B_v\y * *B_v\y
 EndProcedure
  
 Procedure.f length(*var.vec2f)
   ProcedureReturn Sqr(dot(var, var))
 EndProcedure
 
 Procedure vorticity_confinement()
   
   Shared nx
   Shared ny
   
   Shared new_velocity()
   Shared old_velocity()
   
   Dim abs_curl.grid(nx*ny)
   
   Shared vorticity
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       abs_curl(idx(x,y))\x = Abs(curl(x,y))
       abs_curl(idx(x,y))\y = Abs(curl(x,y))
     Next
   Next
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       direction.vec2f
       
       direction\x = abs_curl(idx(x+0,y-1))\x - abs_curl(idx(x+0,y+1))\x
       direction\y = abs_curl(idx(x+1,y+0))\y - abs_curl(idx(x-1,y+0))\y
       
       direction\x = vorticity / (direction\x + 1e-5) * direction\x
       direction\y = vorticity / (direction\y + 1e-5) * direction\y
 
       If x < nx/2
         direction\x = direction\x * 0.0
         direction\y = direction\y * 0.0
       EndIf
       
       new_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x + dt * curl(x, y) * direction\x
       new_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y + dt * curl(x, y) * direction\y
     Next
   Next
   
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       swap_ov_2_nv_x(old_velocity(), new_velocity(), x,y)
       swap_ov_2_nv_y(old_velocity(), old_velocity(), x,y)
     Next
   Next
   
FreeArray(abs_curl())   
   
EndProcedure


Procedure.f clamp(x.f, a.f, b.f)
  If x < a
    ProcedureReturn a
  ElseIf x > b
    ProcedureReturn b
  Else
    ProcedureReturn x
  EndIf
EndProcedure
  
 
 Procedure.f smoothstep(a.f, b.f, u.f)
   t.f = clamp((u-a) / (b-a), 0.0, 1.0)
   
   ProcedureReturn t*t * (3.0 - 2.0 * t)
 EndProcedure
 
   
 
 
 Procedure add_density(px.i, py.i, r.i = 10, value.f = 0.5)
   
   Shared old_density()
   
   For y = -r To r
     For x = -r To r
       d.f = Sqr(x*x + y*y)
       u.f = smoothstep((r), 0.0, d)
       old_density(idx(x,y))\x = old_density(idx(x,y))\x + u*value
       old_density(idx(x,y))\y = old_density(idx(x,y))\y + u*value
     Next
   Next
 EndProcedure
 
 Procedure randf(a.f, b.f)
    u.f = Random(1.0)
    ProcedureReturn lerp(a, b, u)
EndProcedure
 
 Procedure fluid_simulation_step()
   
   Shared nx
   Shared ny
   
   Shared dt
   
   Shared old_velocity()
   Shared old_density()
   
   Shared Myvar.vec2f
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       If (x > nx*0.5) 
         Continue
       EndIf
       

        r.f = 10.0
        old_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x + Random(RandomSeed(r)) 
        old_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y + Random(RandomSeed(r))        
      Next
    Next
    

    ; dense regions rise up
    For y = 0 To ny-1
     For x = 0 To nx-1
       old_velocity(idx(x,y))\x = ( old_velocity(idx(x,y) )\x * 20.0 - 5.0)*dt
       old_velocity(idx(x,y))\y = ( old_velocity(idx(x,y) )\y * 20.0 - 5.0)*dt
      Next
    Next
    

    ;add_density(mouseX, mouseY, 10, 0.5)

    ; fast movement is dampened
    For y = 0 To ny-1
     For x = 0 To nx-1
       old_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x * 0.999
       old_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y * 0.999
      Next
    Next
    

    ; fade away
    For y = 0 To ny-1
     For x = 0 To nx-1
       old_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x  * 0.99
       old_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y  * 0.99
      Next
    Next
    

    add_density(nx*0.25, 30)
    add_density(nx*0.75, 30)
    
    vorticity_confinement()
  
    ; diffuse_velocity()
    ; project_velocity()
    advect_velocity()
  
    project_velocity()
  

    ; diffuse_density()
    advect_density()
    
    Myvar\x = 0.0
    Myvar\y = 0.0
   
    ; zero out stuff at bottom
    For y = 0 To ny-1
     For x = 0 To nx-1
        If (y < 10)
          old_density(idx(x,y))\x = 0.0
          old_density(idx(x,y))\y = 0.0
          old_velocity(idx(x,y))\x = Myvar\x
          old_velocity(idx(x,y))\y = Myvar\y
          EndIf
        Next
      Next
      
          
EndProcedure
    
    
Procedure My_Clamp(x.i, a.i, b.i)
  If x < a
    ProcedureReturn a 
  ElseIf x > b 
    ProcedureReturn b 
  Else
    ProcedureReturn x
  EndIf
EndProcedure


 Procedure rgba32(r.i, g.i, b.i, a.i)
    r = My_Clamp(r, 0, 255)
    g = My_Clamp(g, 0, 255)
    b = My_Clamp(b, 0, 255)
    a = My_Clamp(a, 0, 255)
    ProcedureReturn (a << 24) | (b << 16) | (g << 8) | r
EndProcedure

Procedure My_RGBA(r.f, g.f,  b.f,  a.f)
    ProcedureReturn rgba32(r*256, g*256, b*256, a*256)
EndProcedure

 
 Procedure on_frame()
   
   Shared nx
   Shared ny
   
   Shared old_density()
   Shared pixels_data()
   
   ;Shared pixels()
   
    glClearColor_(0.0, 0.0, 0.0, 0.0)
    glClear_(#GL_COLOR_BUFFER_BIT | #GL_DEPTH_BUFFER_BIT)

    fluid_simulation_step()

  
    ; density field To pixels
    For y = 0 To ny-1
      For x = 0 To nx-1
        fdx.f = old_density(idx(x,y))\x
        fdy.f = old_density(idx(x,y))\y
        fd = Log((fdx+fdy)*0.25 + 1.0)
        f3.f = fd*fd*fd
        r.f = 1.5*fd
        g.f = 1.5*f3
        b.f = f3*f3
      
        pixels_data(idx(x,y)) = RGBA(r, g, b, 1.0)
      Next 
    Next
    
 
      

        ; upload pixels To texture
        glBindTexture_(#GL_TEXTURE_2D, texture)
        glTexSubImage2D_(#GL_TEXTURE_2D, 0, 0, 0, nx, ny, #GL_RGBA, #GL_UNSIGNED_BYTE, pixels_data() )
  


    ; draw texture
    glBegin_(#GL_QUADS)
    ;glColor4f_(1.0, 0.0, 0.0, 0.5)
    glTexCoord2f_(0.0, 0.0) 
    glVertex2f_(-1.0, -1.0)
    glTexCoord2f_(1.0, 0.0) 
    glVertex2f_(1.0, -1.0)
    glTexCoord2f_(1.0, 1.0) 
    glVertex2f_(1.0, 1.0)
    glTexCoord2f_(0.0, 1.0) 
    glVertex2f_(-1.0, 1.0)
    glEnd_()

  EndProcedure
  
  
  
InitKeyboard()
InitSprite()

OpenWindow(0, 0, 0, 800, 600, "OpenGL : Test for Fire Simulation ", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
OpenWindowedScreen(WindowID(0), 0, 0, 800, 600)


glMatrixMode_(#GL_PROJECTION)
glLoadIdentity_() 
gluPerspective_(45.0, 800/600, 1.0, 60.0)
glMatrixMode_(#GL_MODELVIEW)
; glShadeModel_(#GL_SMOOTH)
; glEnable_(#GL_LIGHT0)
; glEnable_(#GL_LIGHTING)
; glEnable_(#GL_COLOR)
; glEnable_(#GL_TEXTURE)
; glEnable_(#GL_STENCIL)
glDisable_(#GL_DEPTH_TEST)


;glEnable_(#GL_CULL_FACE)   
glViewport_(0, 0, 800, 600)


glTranslatef_(0.0, 0.0, -2.0)

Repeat
;glClearColor_(0.0, 0.0, 0.0, 0.0) ; background color
 ; glClear_(#GL_COLOR_BUFFER_BIT | #GL_DEPTH_BUFFER_BIT)
  Repeat
    event = WindowEvent()
    Select Event
      Case #PB_Event_CloseWindow
        quit = 1
     EndSelect
   Until event = 0
   
   init()
   
   on_frame()
   
 
   
 FlipBuffers()
 
 ExamineKeyboard()
 
Until KeyboardPushed(#PB_Key_Escape) Or quit = 1

FreeArray(old_density())
FreeArray(new_density())
FreeArray(old_velocity())
FreeArray(new_velocity())
FreeArray(pixels())
FreeArray(pixels_data())
  

What is the problem with all that ? So I need your help about on how to fix that ? Thanks.

Most important thing is my source code in Purebasic are well writing this way ?? Yes or No ?

[juergenkulow]

@threedslider
What do these lines of code in the procedure on_frame() do?

Code: Select all

    ; density field To pixels
    For y = 0 To ny-1
      For x = 0 To nx-1
        fdx.f = old_density(idx(x,y))\x
        fdy.f = old_density(idx(x,y))\y
        If 0.0<>fdx
          Debug fdx
        EndIf   
        fd = Log((fdx+fdy)*0.25 + 1.0)
        f3.f = fd*fd*fd
        r.f = 1.5*fd
        g.f = 1.5*f3
        b.f = f3*f3
      
        pixels_data(idx(x,y)) = RGBA(r, g, b, 1.0)
        If $1000000<> pixels_data(idx(x,y))
          Debug Hex(pixels_data(idx(x,y))) 
        EndIf   
      Next 
    Next

[threedslider]

@juergenkulow
Thank you, honestly I have not fully tested this my fluid simulation, the problem I think here are on number 134 :

Code: Select all

Procedure.f interpolationX_d(Array g.grid(1), *var.vec2f, ax.i, ay.i)
  ix.i = *var\x 
  iy.i = *var\y
  
  ux.f = *var\x - ix
  uy.f = *var\y - iy
  
 ProcedureReturn lerp(lerp(g(idx(ax+0,ay+0))\x, g(idx(ax+1,ay+0))\x, ux), lerp(g(idx(ax+0,ay+1))\x, g(idx(ax+1,ay+1))\x + 0, ux), uy)
  
EndProcedure

Meaning nothing cause it make zero all and similar to interpolationY(), I don't know how to make ix (variable) for function as floor() equivalent in C for Purebasic... As showing in C++ from github of Fluid simulaiton
So in number 196 :

Code: Select all

 Procedure advect_density()
   Shared nx
   Shared ny
   
   Shared old_density()
   Shared new_density()
   Shared old_velocity()
   
   Shared Myvar.vec2f
   
  
  
   pos.vec2f
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       pos\x = Myvar\x - dt* old_velocity(idx(x,y))\x
       pos\y = Myvar\y - dt* old_velocity(idx(x,y))\y
       
       ;Debug old_velocity(idx(x,y))\x
       
       new_density(idx(x,y))\x = interpolationX_d(old_density(), @pos, x,y)
       new_density(idx(x,y))\y = interpolationY_d(old_density(), @pos, y,y)
       
       Debug new_density(idx(x,y))\x
     Next
   Next
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       swap_od_2_nd_x(old_density(), new_density(), x,y)
       swap_od_2_nd_y(old_density(), new_density(), x,y)
       
       ;Debug old_density(idx(x,y))\x
       
     Next 
   Next   
 EndProcedure

It causes all zeros...

Another here is my random doesn't work well in number 508 : ( fluid_simulation_step() )

Code: Select all

r.f = 10.0
        old_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x + Random(RandomSeed(r)) 
        old_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y + Random(RandomSeed(r))     

I would like to add this function randf() in number 503 :

Code: Select all

 Procedure randf(a.f, b.f)
    u.f = Random(1.0)
    ProcedureReturn lerp(a, b, u)
EndProcedure

But by adding this randf is an error for me

[pjay]

I can't give you a solution, but I can advise you to add:

Code: Select all

EnableExplicit

at the top of your code to highlight your missed definitions (DT, Viscosity etc,.)

Floor = Round(#Num,#PB_Round_Down)

Random() returns an int, try:

Code: Select all

Procedure randf(a.f, b.f)
   ProcedureReturn lerp(a, b, Random(1000000) / 1000000.0)
EndProcedure

I've written a number of these GPU > CPU type solutions; for me to get any decent performance out of them required a lot of optimizations (i.e. multi-threading and dynamic resolution).

[pjay]

There are lots of other little bugs, but you're close.

Code: Select all

r.f = 10.0
        old_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x + Random(RandomSeed(r)) 
        old_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y + Random(RandomSeed(r))     

You need to remove RandomSeed call.

The defined Grid structure needs to be float type, not int type.

Code: Select all

pixels_data(idx(x,y)) = RGBA(r, g, b, 1.0)

... is wrong as RGBA() expects values in the 0 to 255 range (use your My_RGBA() function instead, easy fix).

You're calling init() every frame, but should just be the once.

[threedslider]

@pjay

Thanks to point me that ! But still not working... I have fixed all thank by your tips Even if randf works I have disabled with "EnableExpllicit" thought and it works now without warnings...

Code: Select all

;;; 
;;; Fluid Simulation v. 0.1 
;;;

;EnableExplicit

Define w.i, h.i, nx.i, ny.i, dt.f, iterations.i, vorticity.f
       
w = 512
h = 512

nx = 32
ny = 32

dt = 0.02
iterations = 5
vorticity = 10.0

Global texture

Structure vec2f
  x.f
  y.f
EndStructure

Structure vecvalue
  Array values.vec2f(32*32)
EndStructure

Procedure idx(x.i, y.i) 
  
  Shared nx
  Shared ny
  ; wrap around
  
  Define xd.i, yd.i
  
  xd = (x + nx) % nx
  yd = (y + ny) % ny
  ProcedureReturn xd + yd*nx  
  
EndProcedure
      

Structure grid
  x.f
  y.f
EndStructure

Structure gridvalue
  Array values.grid(32*32)
EndStructure


Dim old_velocity.vec2f(nx*ny)
Dim new_velocity.vec2f(nx*ny)

Dim old_density.grid(nx*ny)
Dim new_density.grid(nx*ny)

Dim pixels.grid(nx*ny)

Dim pixels_data(nx*ny)


Procedure draw(Array var.vec2f(1), n.i, mode.i)
  glEnableClientState_(#GL_VERTEX_ARRAY)
  glVertexPointer_(2, #GL_FLOAT, 0, var())
  glDrawArrays_(mode, 0, n)
EndProcedure



Procedure draw_circle(x.f, y.f, r.f, n.i = 100)
  Dim pos.vec2f(n)
  
  Define Myvar.vec2f
  
  For i = 1 To n
    angle.f = 2.0*3.14159*i/n
    MyVar\x = x
    MyVar\y = y
    pos(i)\x = MyVar\x + r*Cos(angle)
    pos(i)\y = MyVar\y + r*Sin(angle)
  Next
  
  draw(pos(), n, #GL_LINE_LOOP)
  
  FreeArray(pos())
  
EndProcedure
  

Procedure init()
  
  Shared nx
  Shared ny
  
  Shared texture
  
  Shared old_density()
  Shared old_velocity()
  
  Shared Myvar.vec2f
  
  Myvar\x = 0.0
  Myvar\y = 0.0
  
  For y = 0 To ny-1
    For n = 0 To nx-1
      
      old_density(idx(n,y))\x = 0
      old_density(idx(n,y))\y = 0
      old_velocity(idx(n,y))\x = Myvar\x
      old_velocity(idx(n,y))\y = Myvar\y
        
    Next
  Next
  

    glEnable_(#GL_TEXTURE_2D)
    glEnable_(#GL_BLEND)
    glBlendFunc_(#GL_SRC_ALPHA, #GL_ONE_MINUS_SRC_ALPHA)

    glGenTextures_(1, @texture)
    glBindTexture_(#GL_TEXTURE_2D, texture)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MIN_FILTER, #GL_LINEAR)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MAG_FILTER, #GL_LINEAR)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_S, #GL_REPEAT)
    glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_T, #GL_REPEAT)
    glTexImage2D_(#GL_TEXTURE_2D, 0, #GL_RGBA, nx, ny, 0, #GL_RGBA, #GL_UNSIGNED_BYTE, #NUL)
    
;     If texture
;       Debug "work !"
;     Else
;       Debug "not working !"
;     EndIf
    
      
EndProcedure

Procedure.f lerp(a.f, b.f, c.f)
  ProcedureReturn (1.0 - c)*a + c*b
EndProcedure

Procedure.f interpolationX_d(Array g.grid(1), *var.vec2f, ax.i, ay.i)
  ix.i = Round(*var\x,#PB_Round_Down) 
  iy.i = Round(*var\y,#PB_Round_Down)
  
  ux.f = *var\x - ix
  uy.f = *var\y - iy
  
 ProcedureReturn lerp(lerp(g(idx(ax+0,ay+0))\x, g(idx(ax+1,ay+0))\x, ux), lerp(g(idx(ax+0,ay+1))\x, g(idx(ax+1,ay+1))\x + 0, ux), uy)
  
EndProcedure

Procedure.f interpolationY_d(Array g.grid(1), *var.vec2f, ax.i, ay.i)
  ix.i = Round(*var\x,#PB_Round_Down)
  iy.i = Round(*var\y,#PB_Round_Down)
  
  ux.f = *var\x - ix
  uy.f = *var\y - iy
  
 ProcedureReturn lerp(lerp(g(idx(ax+0,ay+0))\y, g(idx(ax+1,ay+0))\y, ux), lerp(g(idx(ax+0,ay+1))\y, g(idx(ax+1,ay+1))\y + 0, ux), uy)
  
EndProcedure

Procedure.f interpolationX_v(Array g.vec2f(1), *var.vec2f, ax.i, ay.i)
  ix.i = Round(*var\x,#PB_Round_Down)
  iy.i = Round(*var\y,#PB_Round_Down)
  
  ux.f = *var\x - ix
  uy.f = *var\y - iy
  
 ProcedureReturn lerp(lerp(g(idx(ax+0,ay+0))\x, g(idx(ax+1,ay+0))\x, ux), lerp(g(idx(ax+0,ay+1))\x, g(idx(ax+1,ay+1))\x + 0, ux), uy)
  
EndProcedure

Procedure.f interpolationY_v(Array g.vec2f(1), *var.vec2f, ax.i, ay.i)
  ix.i = Round(*var\x,#PB_Round_Down)
  iy.i = Round(*var\y,#PB_Round_Down)
  
  ux.f = *var\x - ix
  uy.f = *var\y - iy
  
 ProcedureReturn lerp(lerp(g(idx(ax+0,ay+0))\y, g(idx(ax+1,ay+0))\y, ux), lerp(g(idx(ax+0,ay+1))\y, g(idx(ax+1,ay+1))\y + 0, ux), uy)
  
 EndProcedure

Procedure swap_od_2_nd_x(Array od.grid(1), Array nd.grid(1), ax.i, ay.i)
  Swap od(idx(ax,ay))\x, nd(idx(ax,ay))\x
EndProcedure

Procedure swap_od_2_nd_y(Array od.grid(1), Array nd.grid(1), ax.i, ay.i)
  Swap od(idx(ax,ay))\y, nd(idx(ax,ay))\y
EndProcedure

Procedure swap_ov_2_nv_x(Array ov.vec2f(1), Array nv.vec2f(1), ax.i, ay.i)
  Swap ov(idx(ax,ay))\x, nv(idx(ax,ay))\x
EndProcedure

Procedure swap_ov_2_nv_y(Array ov.vec2f(1), Array nv.vec2f(1), ax.i, ay.i)
  Swap ov(idx(ax,ay))\y, nv(idx(ax,ay))\y
EndProcedure

 
 
 Procedure advect_density()
   Shared nx
   Shared ny
   
   Shared old_density()
   Shared new_density()
   Shared old_velocity()
   
   Shared Myvar.vec2f
   
  
  
   pos.vec2f
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       pos\x = Myvar\x - dt* old_velocity(idx(x,y))\x
       pos\y = Myvar\y - dt* old_velocity(idx(x,y))\y
       
       ;Debug old_velocity(idx(x,y))\x
       
       new_density(idx(x,y))\x = interpolationX_d(old_density(), @pos, x,y)
       new_density(idx(x,y))\y = interpolationY_d(old_density(), @pos, y,y)
       
       ;Debug new_density(idx(x,y))\x
     Next
   Next
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       swap_od_2_nd_x(old_density(), new_density(), x,y)
       swap_od_2_nd_y(old_density(), new_density(), x,y)
       
       ;Debug old_density(idx(x,y))\x
       
     Next 
   Next   
 EndProcedure
 
 
  Procedure advect_velocity()
   Shared nx
   Shared ny
   
   Shared old_velocity()
   Shared old_density()
   
   Shared new_density()
   Shared new_velocity()
   
   Shared Myvar.vec2f
   
   pos.vec2f
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       pos\x = Myvar\x - dt* old_velocity(idx(x,y))\x
       pos\y = Myvar\y - dt* old_velocity(idx(x,y))\y
       
       new_velocity(idx(x,y))\x = interpolationX_v(old_velocity(), @pos, x,y)
       new_velocity(idx(x,y))\y = interpolationY_v(old_velocity(), @pos, x,y)
     Next
   Next
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       swap_ov_2_nv_x(old_velocity(), new_velocity(), x,y)
       swap_ov_2_nv_y(old_velocity(), new_velocity(), x,y)
     Next 
   Next   
 EndProcedure
 
 
 Procedure diffuse_density()
   Shared dt
   Shared nx
   Shared ny
   
   Shared old_density()
   Shared new_density()
   
   sumX = 0
   sumY = 0
   
   diffusion.f = dt * 100.01
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       
       sumX = diffusion * ( old_density(idx(x-1,y+0))\x + old_density(idx(x+1,y+0))\x + old_density(idx(x+0,y-1))\x + old_density(idx(x+0,y+1))\x ) + old_density(idx(x+0,y+0))\x 
       sumY = diffusion * ( old_density(idx(x-1,y+0))\y + old_density(idx(x+1,y+0))\y + old_density(idx(x+0,y-1))\y + old_density(idx(x+0,y+1))\y ) + old_density(idx(x+0,y+0))\y 
       
       new_density(idx(x,y))\x = 1.0 / (1.0 + 4.0 * viscosity) * sumX
       new_density(idx(x,y))\y = 1.0 / (1.0 + 4.0 * viscosity) * sumY
     Next
   Next
   
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       swap_od_2_nd_x(old_density(), new_density(), x,y)
       swap_od_2_nd_y(old_density(), new_density(), x,y)
     Next
   Next
   
 EndProcedure
 
 
  Procedure diffuse_velocity()
   Shared dt
   Shared nx
   Shared ny
   
   Shared old_velocity()
   Shared new_velocity()
   
   sum.vec2f
   
   viscosity.f = dt * 0.000001
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       
       sum\x = viscosity * ( old_velocity(idx(x-1,y+0))\x +  old_velocity(idx(x+1,y+0))\x +  old_velocity(idx(x+0,y-1))\x +  old_velocity(idx(x+0,y+1))\x ) +  old_velocity(idx(x+0,y+0))\x 
       sum\y = viscosity * ( old_velocity(idx(x-1,y+0))\y +  old_velocity(idx(x+1,y+0))\y +  old_velocity(idx(x+0,y-1))\y +  old_velocity(idx(x+0,y+1))\y ) +  old_velocity(idx(x+0,y+0))\y 
       
       new_velocity(idx(x,y))\x = 1.0 / (1.0 + 4.0 * viscosity) * sum\x
       new_velocity(idx(x,y))\y = 1.0 / (1.0 + 4.0 * viscosity) * sum\y
     Next
   Next
   
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       swap_ov_2_nv_x(old_velocity(), new_velocity(), x,y)
       swap_ov_2_nv_y(old_velocity(), new_velocity(), x,y)
       
     Next
   Next
   
 EndProcedure    
 
 
 Procedure project_velocity()
   
   Shared nx
   Shared ny
   
   Shared iterations
   
   Shared old_velocity()
   
   Dim p.grid(nx*ny)
   Dim p2.grid(nx*ny)
   Dim div.grid(nx*ny)
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       dx.f = old_velocity(idx(x+1,y+0))\x - old_velocity(idx(x-1,y+0))\x
       dy.f = old_velocity(idx(x+0,y+1))\y - old_velocity(idx(x+0,y-1))\y
       
       div(idx(x,y))\x = dx + dy
       div(idx(x,y))\y = dx + dy
       
       p(idx(x,y))\x = 0.0
       p(idx(x,y))\y = 0.0
     Next
   Next
   
   
   For k = 1 To iterations
     
     For y = 0 To ny-1
       For x = 0 To nx-1
         sumX.f = -div(idx(x,y))\x + p(idx(x+1,y+0))\x + p(idx(x-1,y+0))\x + p(idx(x+0,y+1))\x + p(idx(x+0,y-1))\x 
         sumY.f = -div(idx(x,y))\y + p(idx(x+1,y+0))\y + p(idx(x-1,y+0))\y + p(idx(x+0,y+1))\y + p(idx(x+0,y-1))\y
         
         p2(idx(x,y))\x = 0.25 * sumX
         p2(idx(x,y))\y = 0.25 * sumY
       Next
     Next
     
     For y = 0 To ny-1
       For x = 0 To nx-1
         swap_od_2_nd_x(p(), p2(), x,y)
         swap_od_2_nd_y(p(), p2(), x,y)
       Next
     Next
     
   Next
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       old_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x - 0.5 * p(idx(x+1,y+0))\x - p(idx(x-1,y+0))\x
       old_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y - 0.5 * p(idx(x+0,y+1))\y - p(idx(x+0,y-1))\y
     Next
   Next
   
   
   FreeArray(p())
   FreeArray(p2())
   FreeArray(div())
 EndProcedure
 
 
 Procedure.f curl(ax.i, ay.i)
   
   Shared old_velocity()
   
   ProcedureReturn old_velocity(idx(ax,ay+1))\x + old_velocity(idx(ax,ay-1))\x + old_velocity(idx(ax-1,ay))\y + old_velocity(idx(ax+1,ay))\y
   
 EndProcedure
 
 Procedure.f dot(*A_v.vec2f, *B_v.vec2f)
   ProcedureReturn *A_v\x * *B_v\x + *B_v\y * *B_v\y
 EndProcedure
  
 Procedure.f length(*var.vec2f)
   ProcedureReturn Sqr(dot(var, var))
 EndProcedure
 
 Procedure vorticity_confinement()
   
   Shared nx
   Shared ny
   
   Shared new_velocity()
   Shared old_velocity()
   
   Dim abs_curl.grid(nx*ny)
   
   Shared vorticity
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       abs_curl(idx(x,y))\x = Abs(curl(x,y))
       abs_curl(idx(x,y))\y = Abs(curl(x,y))
     Next
   Next
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       direction.vec2f
       
       direction\x = abs_curl(idx(x+0,y-1))\x - abs_curl(idx(x+0,y+1))\x
       direction\y = abs_curl(idx(x+1,y+0))\y - abs_curl(idx(x-1,y+0))\y
       
       direction\x = vorticity / (direction\x + 1e-5) * direction\x
       direction\y = vorticity / (direction\y + 1e-5) * direction\y
 
       If x < nx/2
         direction\x = direction\x * 0.0
         direction\y = direction\y * 0.0
       EndIf
       
       new_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x + dt * curl(x, y) * direction\x
       new_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y + dt * curl(x, y) * direction\y
     Next
   Next
   
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       swap_ov_2_nv_x(old_velocity(), new_velocity(), x,y)
       swap_ov_2_nv_y(old_velocity(), old_velocity(), x,y)
     Next
   Next
   
FreeArray(abs_curl())   
   
EndProcedure


Procedure.f clamp(x.f, a.f, b.f)
  If x < a
    ProcedureReturn a
  ElseIf x > b
    ProcedureReturn b
  Else
    ProcedureReturn x
  EndIf
EndProcedure
  
 
 Procedure.f smoothstep(a.f, b.f, u.f)
   t.f = clamp((u-a) / (b-a), 0.0, 1.0)
   
   ProcedureReturn t*t * (3.0 - 2.0 * t)
 EndProcedure
 
   
 
 
 Procedure add_density(px.i, py.i, r.i = 10, value.f = 0.5)
   
   Shared old_density()
   
   For y = -r To r
     For x = -r To r
       d.f = Sqr(x*x + y*y)
       u.f = smoothstep((r), 0.0, d)
       old_density(idx(x,y))\x = old_density(idx(x,y))\x + u*value
       old_density(idx(x,y))\y = old_density(idx(x,y))\y + u*value
     Next
   Next
 EndProcedure
 
 Procedure randf(a.f, b.f)
    ProcedureReturn lerp(a, b, Random(1000000) / 1000000.0)
EndProcedure
 
 Procedure fluid_simulation_step()
   
   Shared nx
   Shared ny
   
   Shared dt
   
   Shared old_velocity()
   Shared old_density()
   
   Shared Myvar.vec2f
   
   For y = 0 To ny-1
     For x = 0 To nx-1
       If (x > nx*0.5) 
         Continue
       EndIf
       

        r.f = 10.0
        old_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x + randf(-r, r)
        old_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y + randf(-r, r)    
      Next
    Next
    

    ; dense regions rise up
    For y = 0 To ny-1
     For x = 0 To nx-1
       old_velocity(idx(x,y))\x = ( old_velocity(idx(x,y) )\x * 20.0 - 5.0)*dt
       old_velocity(idx(x,y))\y = ( old_velocity(idx(x,y) )\y * 20.0 - 5.0)*dt
      Next
    Next
    

    ;add_density(mouseX, mouseY, 10, 0.5)

    ; fast movement is dampened
    For y = 0 To ny-1
     For x = 0 To nx-1
       old_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x * 0.999
       old_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y * 0.999
      Next
    Next
    

    ; fade away
    For y = 0 To ny-1
     For x = 0 To nx-1
       old_velocity(idx(x,y))\x = old_velocity(idx(x,y))\x  * 0.99
       old_velocity(idx(x,y))\y = old_velocity(idx(x,y))\y  * 0.99
      Next
    Next
    

    add_density(nx*0.25, 30)
    add_density(nx*0.75, 30)
    
    vorticity_confinement()
  
    ; diffuse_velocity()
    ; project_velocity()
    advect_velocity()
  
    project_velocity()
  

    ; diffuse_density()
    advect_density()
    
    Myvar\x = 0.0
    Myvar\y = 0.0
   
    ; zero out stuff at bottom
    For y = 0 To ny-1
     For x = 0 To nx-1
        If (y < 10)
          old_density(idx(x,y))\x = 0.0
          old_density(idx(x,y))\y = 0.0
          old_velocity(idx(x,y))\x = Myvar\x
          old_velocity(idx(x,y))\y = Myvar\y
          EndIf
        Next
      Next
      
          
EndProcedure
    
    
Procedure My_Clamp(x.i, a.i, b.i)
  If x < a
    ProcedureReturn a 
  ElseIf x > b 
    ProcedureReturn b 
  Else
    ProcedureReturn x
  EndIf
EndProcedure


 Procedure rgba32(r.i, g.i, b.i, a.i)
    r = My_Clamp(r, 0, 255)
    g = My_Clamp(g, 0, 255)
    b = My_Clamp(b, 0, 255)
    a = My_Clamp(a, 0, 255)
    ProcedureReturn (a << 24) | (b << 16) | (g << 8) | r
EndProcedure

Procedure My_RGBA(r.f, g.f,  b.f,  a.f)
    ProcedureReturn rgba32(r*256, g*256, b*256, a*256)
EndProcedure

 
 Procedure on_frame()
   
   Shared nx
   Shared ny
   
   Shared old_density()
   Shared pixels_data()
   
   ;Shared pixels()
   
    glClearColor_(0.0, 0.0, 0.0, 0.0)
    glClear_(#GL_COLOR_BUFFER_BIT | #GL_DEPTH_BUFFER_BIT)

    fluid_simulation_step()

  
    ; density field To pixels
    For y = 0 To ny-1
      For x = 0 To nx-1
        fdx.f = old_density(idx(x,y))\x
        
        ;Debug fdx
        
        fdy.f = old_density(idx(x,y))\y
        fd = Log((fdx+fdy)*0.25 + 1.0)
        f3.f = fd*fd*fd
        r.f = 1.5*fd
        g.f = 1.5*f3
        b.f = f3*f3
      
        pixels_data(idx(x,y)) = My_RGBA(r, g, b, 1.0)
        ;Debug pixels_data(idx(x,y))
      Next 
    Next


        ; upload pixels To texture
        glBindTexture_(#GL_TEXTURE_2D, texture)
        glTexSubImage2D_(#GL_TEXTURE_2D, 0, 0, 0, nx, ny, #GL_RGBA, #GL_UNSIGNED_BYTE, pixels_data() )
        
    ; draw texture
    glBegin_(#GL_QUADS)
    ;glColor4f_(1.0, 0.0, 0.0, 0.5)
    glTexCoord2f_(0.0, 0.0) 
    glVertex2f_(-1.0, -1.0)
    glTexCoord2f_(1.0, 0.0) 
    glVertex2f_(1.0, -1.0)
    glTexCoord2f_(1.0, 1.0) 
    glVertex2f_(1.0, 1.0)
    glTexCoord2f_(0.0, 1.0) 
    glVertex2f_(-1.0, 1.0)
    glEnd_()

  EndProcedure
  
  
  
InitKeyboard()
InitSprite()

OpenWindow(0, 0, 0, 800, 600, "OpenGL : Test for Fire Simulation ", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
OpenWindowedScreen(WindowID(0), 0, 0, 800, 600)


glMatrixMode_(#GL_PROJECTION)
glLoadIdentity_() 
gluPerspective_(45.0, 800/600, 1.0, 60.0)
glMatrixMode_(#GL_MODELVIEW)
glLoadIdentity_() 
glMatrixMode_(#GL_MODELVIEW)

; glShadeModel_(#GL_SMOOTH)
; glEnable_(#GL_LIGHT0)
; glEnable_(#GL_LIGHTING)
; glEnable_(#GL_COLOR)
; glEnable_(#GL_TEXTURE)
; glEnable_(#GL_STENCIL)
glDisable_(#GL_DEPTH_TEST)


;glEnable_(#GL_CULL_FACE)   
glViewport_(0, 0, 800, 600)


glTranslatef_(0.0, 0.0, -2.0)

init()


Repeat
  Repeat
    event = WindowEvent()
    Select Event
      Case #PB_Event_CloseWindow
        quit = 1
     EndSelect
   Until event = 0
   
  
   
  on_frame() 
   
 
   
 FlipBuffers()
 ExamineKeyboard()
 
Until KeyboardPushed(#PB_Key_Escape) Or quit = 1

FreeArray(old_density())
FreeArray(new_density())
FreeArray(old_velocity())
FreeArray(new_velocity())
FreeArray(pixels())
FreeArray(pixels_data())
  

It seems to work now but none of display to the fluid simulation....

Maybe it needs for OpenGL gadgets for working ?

[pjay]

Re-enable EnableExplicit and remove *all* issues from there, you still have a number of them. It will highlight that you're trying to share several non-global variables & issues where you're not sharing variables where you should be.

It'll also highlight this code as a something which needs fixing: (although you're not actually using it...)

Code: Select all

 Procedure.f length(*var.vec2f)
   ProcedureReturn Sqr(dot(var, var))
 EndProcedure

If you haven't done so already, I recommend getting a simple app working that draws a green quad on an OpenGLGadget that is using Ortho (2d) projection. Once you've done that you can apply that knowledge to this code.

[SMaag]

There is no StartDrawing() and StopDrawing().
I don't know exactly if you need StartDrawing with this OpenGl stuff.

Where do you specify the output for drawing?

[threedslider]

Re-enable EnableExplicit and remove *all* issues from there, you still have a number of them. It will highlight that you're trying to share several non-global variables & issues where you're not sharing variables where you should be.

It'll also highlight this code as a something which needs fixing: (although you're not actually using it...)

Code: Select all

 Procedure.f length(*var.vec2f)
   ProcedureReturn Sqr(dot(var, var))
 EndProcedure

Strange I have not this isssue of hightlighting, mine it compiles fine ... Maybe for me I am missing something to activate the transparency. And of course i don't use the length thought.

If you haven't done so already, I recommend getting a simple app working that draws a green quad on an OpenGLGadget that is using Ortho (2d) projection. Once you've done that you can apply that knowledge to this code.

I will make without OpenGLGadget to see if it works, I will let you to know it.

[threedslider]

There is no StartDrawing() and StopDrawing().
I don't know exactly if you need StartDrawing with this OpenGl stuff.

Where do you specify the output for drawing?

It draws from "on_frame() function in last down page. OpenGL doesn't need to be draw under to startdrawing/stopdrawing.

[threedslider]

@pjay

I have done with openglgadget and not working at all for my fluid simulation
Maybe I may be wrong to all my code !

@mpz if you have the time, you could help me please how i fix that all ? Thanks.

I need to review all my code thought ....

[pjay]

Hi threeslider,

I got it working with your original code but it wasn't very fluid (excuse the pun!), so I re-wrote it in a style I was more comfortable with & to the point where I thought it was working OK.

Code: Select all

; fluid dynamics (flame) - PB conversion of https://github.com/983/Fluid-Simulation  - Phil James 2023

EnableExplicit

CompilerIf #PB_Compiler_OS=#PB_OS_Linux
  Structure point
    x.i : y.i
  EndStructure
CompilerEndIf

Structure vec2f
  x.f : y.f
EndStructure

;#nx = 256 : #ny = 256 : #Upscale = 4 ; ~7ms on 5800h
;#nx = 160 : #ny = 120 : #Upscale = 8 ; ~2ms on 5800h
;#nx = 320 : #ny = 160 : #Upscale = 6 ; ~5ms on 5800h
;#nx = 512 : #ny = 256 : #Upscale = 4 ; ~12ms on 5800h
#nx = 320 : #ny = 400 : #Upscale = 2 ; ~12ms on 5800h

#nx1 = #nx - 1 : #ny1 = #ny - 1 : #nx2 = #nx - 2 : #ny2 = #ny - 2
Global dt.f = 0.02, iterations = 5, vorticity.f = 8.0, mouse.point, lmb, w.i = #nx * #Upscale, h = #ny * #Upscale, texture
Global v1size = #nx * #ny * SizeOf(float), v2size = v1size * 2, onframe

Global Dim old_velocity.vec2f(#ny1,#nx1), Dim new_velocity.vec2f(#ny1,#nx1)
Global Dim old_density.f(#ny1, #nx1), Dim new_density.f(#ny1, #nx1)
Global Dim pixels.l(#ny1,#nx1)

Macro LengthM(x,y) : Sqr((x*x) + (y*y)) : EndMacro
Macro LerpM(a, b, t) : ( a + t * (b - a) ) : EndMacro
Macro Log2M(v) : (Log(v) / 0.6931471805599453094) : EndMacro ; not sure if this is correct Log2...
Macro RandomFloatM(a, b) : LerpM(a, b, Random(1000000) / 1000000.0) : EndMacro
 
Procedure Init()
  Protected x, y
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_density(y,x) = 0.0 : old_velocity(y,x)\x = 0 : old_velocity(y,x)\y = 0
  Next : Next

  OpenWindow(0,0,0,w,h,"Fluid-Simulation (Flame)",#PB_Window_ScreenCentered|#PB_Window_SystemMenu)
  OpenGLGadget(0,0,0,w,h,#PB_OpenGL_NoFlipSynchronization)
  glOrtho_(0,1,1,0,0,1)
  glEnable_(#GL_TEXTURE_2D);
  
  glGenTextures_(1, @texture);
  glBindTexture_(#GL_TEXTURE_2D, texture);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MIN_FILTER, #GL_LINEAR);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MAG_FILTER, #GL_LINEAR);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_S, #GL_REPEAT)    ;
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_T, #GL_REPEAT)    ;
  glTexImage2D_(#GL_TEXTURE_2D, 0, #GL_RGBA, #nx, #ny, 0, #GL_RGBA, #GL_UNSIGNED_BYTE, #Null);
EndProcedure
Procedure.f interpolate(Array grid.f(2), *p.vec2f, *d.float)
  Protected ix = Round(*p\x,#PB_Round_Down), iy = Round(*p\y,#PB_Round_Down), ux.f = *p\x - ix, uy.f = *p\y - iy
  *d\f = LerpM(LerpM(grid(iy,ix),grid(iy,ix + 1), ux),LerpM(grid(iy+1,ix),grid(iy+1,ix+1), ux), uy)
EndProcedure
Procedure.f interpolateVec2(Array grid.vec2f(2), *p.vec2f, *d.vec2f)
  Protected ix = Round(*p\x,#PB_Round_Down), iy = Round(*p\y,#PB_Round_Down), ux.f = *p\x - ix, uy.f = *p\y - iy
  *d\x = LerpM(LerpM(grid(iy,ix)\x,grid(iy,ix+1)\x, ux), LerpM(grid(iy+1,ix)\x,grid(iy+1,ix+1)\x, ux), uy)
  *d\y = LerpM(LerpM(grid(iy,ix)\y,grid(iy,ix+1)\y, ux), LerpM(grid(iy+1,ix)\y,grid(iy+1,ix+1)\y, ux), uy)
EndProcedure

Procedure advect_density()
  Protected pos.vec2f, x, y
  For y = 1 To #ny2 : For x = 1 To #nx2
      pos\x = x - dt * old_velocity(y,x)\x
      pos\y = y - dt * old_velocity(y,x)\y
      If pos\y < 1 : pos\y = 1 : ElseIf pos\y > #ny2 : pos\y = #ny2 : EndIf 
      If pos\x < 1 : pos\x = 1 : ElseIf pos\x > #nx2 : pos\x = #nx2 : EndIf
      interpolate(old_density(), @pos, @new_density(y,x));
  Next : Next
  CopyMemory(@new_density(),@old_density(),v1size)
EndProcedure
Procedure advect_velocity()
  Protected pos.vec2f, x, y
  For y = 1 To #ny2 : For x = 1 To #nx2
      pos\x = x - dt * old_velocity(y,x)\x;
      pos\y = y - dt * old_velocity(y,x)\y;
      If pos\y < 1 : pos\y = 1 : ElseIf pos\y > #ny2 : pos\y = #ny2 : EndIf 
      If pos\x < 1 : pos\x = 1 : ElseIf pos\x > #nx2 : pos\x = #nx2 : EndIf
      interpolateVec2(old_velocity(), @pos,@new_velocity(y,x))
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure diffuse_density()
  Protected pos.vec2f, x, y, diffusion.f = dt * 100.01, sum.f
  For y = 1 To #ny2 : For x = 1 To #nx2
      sum = diffusion * (old_density(y, x - 1) + old_density(y, x + 1) + old_density(y - 1, x) + old_density(y + 1, x)) + old_density(y,x);
      new_density(y, x) = 1.0 / (1.0 + 4.0 * diffusion) * sum 
  Next : Next
  CopyMemory(@new_density(),@old_density(),v1size)
EndProcedure
Procedure diffuse_velocity()
  Protected viscosity = dt*0.000001, x, y, sum.vec2f, gotnan.a  = #False
  For y = 1 To #ny2 : For x = 1 To #nx2
      sum\x = viscosity * (old_velocity(y, x - 1)\x + old_velocity(y, x + 1)\x + old_velocity(y - 1, x)\x + old_velocity(y + 1, x)\x) + old_velocity(y, x)\x
      sum\y = viscosity * (old_velocity(y, x - 1)\y + old_velocity(y, x + 1)\y + old_velocity(y - 1, x)\y + old_velocity(y + 1, x)\y) + old_velocity(y, x)\y
      new_velocity(y,x)\x = 1.0 / (1.0 + 4.0 * viscosity) * sum\x
      new_velocity(y,x)\y = 1.0 / (1.0 + 4.0 * viscosity) * sum\y
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure project_velocity()
  Protected x, y, dx.f, dy.f, k, sum.f
  Static Dim p.f(#ny,#nx), Dim p2.f(#ny,#nx), Dim div.f(#ny,#nx)

  For y = 1 To #ny2 : For x = 1 To #nx2
      dx = old_velocity(y    , x + 1)\x - old_velocity(y    ,x - 1)\x
      dy = old_velocity(y + 1, x    )\y - old_velocity(y - 1,x    )\y
      div(y,x) = dx + dy : p(y,x) = 0.0
  Next : Next
  For k = 1 To iterations
    For y = 1 To #ny2 : For x = 1 To #nx2
        sum = -div(y,x) + p(y ,x + 1) + p(y ,x - 1) + p(y + 1,x) + p(y - 1,x)
        p2(y,x) = 0.25 * sum
    Next : Next
    CopyMemory(@p2(),@p(),v1size)
  Next
  For y = 1 To #ny2 : For x = 1 To #nx2
      old_velocity(y,x)\x - (0.5 * (p(y    ,x + 1) - p(y    , x - 1)))
      old_velocity(y,x)\y - (0.5 * (p(y + 1,x  ) - p(y - 1, x )))
  Next : Next
EndProcedure
Procedure.f curl(y,x)
  ProcedureReturn old_velocity(y + 1,x)\x - old_velocity(y - 1,x)\x + old_velocity(y,x - 1)\y - old_velocity(y, x + 1)\y;
EndProcedure
Procedure vorticity_confinement()
  Protected direction.vec2f, length.f, x, y;
  Static Dim abs_curl.f(#ny1,#nx1)
  For y = 1 To #ny2 : For x = 1 To #nx2
      abs_curl(y,x) = Abs(curl(y,x));
  Next : Next
  
  For y = 1 To #ny2 : For x = 1 To #nx2
      direction\x = abs_curl(y - 1, x    ) - abs_curl(y + 1, x    );
      direction\y = abs_curl(y    , x + 1) - abs_curl(y    , x - 1);
      length = lengthM(direction\x,direction\y)
      direction\x = vorticity / (length + 1e-5) * direction\x;
      direction\y = vorticity / (length + 1e-5) * direction\y;
      new_velocity(y,x)\x = old_velocity(y,x)\x + dt * curl(y,x) * direction\x;
      new_velocity(y,x)\y = old_velocity(y,x)\y + dt * curl(y,x) * direction\y;
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure.f Clamp(v.f,Low.f = 0.0, High.f = 1.0)
  If v < Low : ProcedureReturn Low : ElseIf v > High : ProcedureReturn High : EndIf
  ProcedureReturn v
EndProcedure
Procedure.f SmoothStep(Low.f, High.f, v.f)
  V = Clamp((v - low) / (high - low));
  ProcedureReturn v * v * (3.0 - 2.0 * v);
EndProcedure

Procedure add_density(px, py, r = 10, value.f = 0.5)
  Protected x, y, u.f, rf.f = r, x2, y2
  For y = -r To r : For x = -r To  r
      y2 = py + y : x2 = px + x
      If x2 < 0 : x2 = 0 : ElseIf x2 > #nx1 : x2 = #nx1 : EndIf
      If y2 < 0 : y2 = 0 : ElseIf y2 > #ny1 : y2 = #ny1 : EndIf
      u = SmoothStep(rf, 0.0, Sqr(x*x + y*y));
      old_density(y2, x2) + u * value
  Next : Next
  
EndProcedure

Procedure fluid_simulation_step()
  Protected x,y, r.f = 10.0
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_velocity(y,x)\x + RandomFloatM(-r, r) : old_velocity(y,x)\y + RandomFloatM(-r, r);
  Next : Next
  
  ;// dense regions rise up
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_velocity(y,x)\y + (old_density(y,x) * 20.0 - 5.0) * dt;
  Next : Next
 
  ;// fast movement is dampened
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_velocity(y,x)\x * 0.999 : old_velocity(y,x)\y * 0.999;
  Next : Next
  
  ;// fade away
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_density(y,x) * 0.99;
  Next : Next

  add_density(#nx*0.25, #ny*0.1, 10,0.5);
  add_density(#nx*0.75, #ny*0.1, 10,0.5);
  
  If mouse\x > 0 And mouse\y > 0 And mouse\x < #nx And mouse\y < #ny1 : add_density(Mouse\x, #ny - Mouse\y, #nx * 0.025,0.5) : EndIf
  
  vorticity_confinement()
  advect_velocity()
  project_velocity()
  diffuse_velocity()
  diffuse_density();
  advect_density()
  
  ;/ wipe out borders as prevents crashes
  For y = 0 To #ny1 
      old_density(y,0) = 0.0 : old_velocity(y,0)\x = 0 : old_velocity(y,0)\y = 0      
      old_density(y,1) = 0.0 : old_velocity(y,1)\x = 0 : old_velocity(y,1)\y = 0    
      old_density(y,#nx1) = 0.0 : old_velocity(y,#nx1)\x = 0 : old_velocity(y,#nx1)\y = 0      
      old_density(y,#nx2) = 0.0 : old_velocity(y,#nx2)\x = 0 : old_velocity(y,#nx2)\y = 0      
  Next
  For x = 0 To #nx1 
      old_density(0,x) = 0.0 : old_velocity(0,x)\x = 0 : old_velocity(0,x)\y = 0      
      old_density(#ny1,x) = 0.0 : old_velocity(#ny1,x)\x = 0 : old_velocity(#ny1,x)\y = 0      
      old_density(1,x) = 0.0 : old_velocity(1,x)\x = 0 : old_velocity(1,x)\y = 0      
      old_density(#ny2,x) = 0.0 : old_velocity(#ny2,x)\x = 0 : old_velocity(#ny2,x)\y = 0      
  Next
EndProcedure

Procedure.l RGBAgl(r.f, g.f, b.f, a.f)
  r = clamp(r * 255, 0, 255) : g = clamp(g * 255, 0, 255) : b = clamp(b * 255, 0, 255) : a = clamp(a * 255, 0, 255)
  ProcedureReturn (Int(a) << 24) | (Int(b) << 16) | (Int(g) << 8) | Int(r);
EndProcedure
Procedure Render()
  ;// density field To pixels
  Protected x, y, f.f, f3.f, r.f, g.f, b.f
  For y = 0 To #ny1 : For x = 0 To #nx1
      f = old_density(y,x);
      f = log2m( f * 0.25 + 1.0);
      f3 = f*f*f
      r = 1.5 * f : g = 1.5 * f3 : b = (f3 *f3)
      pixels(y,x) = RGBAgl(r, g, b, 1.0);
  Next : Next
  
  ;// upload pixels() array To texture
  glTexSubImage2D_(#GL_TEXTURE_2D, 0, 0, 0, #nx, #ny, #GL_RGBA, #GL_UNSIGNED_BYTE, @pixels());
  
  ;// draw texture
  glBegin_(#GL_QUADS);
  glTexCoord2f_(0.0, 1.0): glVertex2f_(0  , 0)
  glTexCoord2f_(1.0, 1.0): glVertex2f_(1.0, 0)
  glTexCoord2f_(1.0, 0.0): glVertex2f_(1.0, 1.0)
  glTexCoord2f_(0.0, 0.0): glVertex2f_(0  , 1.0)
  glEnd_()
  
  SetGadgetAttribute(0,#PB_OpenGL_FlipBuffers,#True)
EndProcedure
Init()

Define Event, Exit, time.i

Repeat
  Repeat
    Event = WindowEvent()
    Select Event
      Case #PB_Event_CloseWindow : Exit = #True
      Case #PB_Event_Gadget
        Select EventType()
          Case #PB_EventType_MouseMove
            Mouse\x = (GetGadgetAttribute(0,#PB_OpenGL_MouseX) / (w / #nx)) / DesktopResolutionX()
            Mouse\y = (GetGadgetAttribute(0,#PB_OpenGL_MouseY) / (h / #ny)) / DesktopResolutionY()
          Case #PB_EventType_LeftButtonDown
            LMB = #True
            add_density(Mouse\x, #ny - Mouse\y, #nx * 0.2,8.0); burst
          Case #PB_EventType_LeftButtonUp
            LMB = #False
        EndSelect
    EndSelect
  Until Event = 0

  time = ElapsedMilliseconds()
  fluid_simulation_step();
  Render() : onframe + 1
  time = ElapsedMilliseconds() - time
  If onframe = 60 : onframe = 1
    SetWindowTitle(0,"Fluid dynamics w/vorticity - Phil James 2023 - Process time: "+ Str(time)+"ms.")
  EndIf
  
Until Exit = #True
glDeleteTextures_(1,@texture)

[juergenkulow]

Optimized with C backend compiler and optimize compiler option:

and

Code: Select all

DisableDebugger
CompilerIf #PB_Compiler_OS=#PB_OS_Linux
  Structure point
    x.i
    y.i
  EndStructure
CompilerEndIf

Does the image match the video?

[pf shadoko]

great !

[pjay]

Thanks for the Linux tip, added point structure as suggested. It's a good showcase to see the difference the C compiler with optimations enabled can make.

What I posted had been modified a bit from the base code - Is below code more similar to the video?

Code: Select all

; fluid dynamics (flame) - PB conversion of https://github.com/983/Fluid-Simulation  - Phil James 2023

EnableExplicit

CompilerIf #PB_Compiler_OS<>#PB_OS_Windows
  Structure point
    x.i : y.i
  EndStructure
CompilerEndIf

Structure vec2f
  x.f : y.f
EndStructure

;#nx = 256 : #ny = 256 : #Upscale = 4 ; ~7ms on 5800h
;#nx = 160 : #ny = 120 : #Upscale = 8 ; ~2ms on 5800h
;#nx = 320 : #ny = 160 : #Upscale = 6 ; ~5ms on 5800h
;#nx = 512 : #ny = 256 : #Upscale = 4 ; ~12ms on 5800h
#nx = 320 : #ny = 400 : #Upscale = 2 ; ~12ms on 5800h

#nx1 = #nx - 1 : #ny1 = #ny - 1 : #nx2 = #nx - 2 : #ny2 = #ny - 2
Global dt.f = 0.02, iterations = 5, vorticity.f = 10.0, mouse.point, lmb, w.i = #nx * #Upscale, h = #ny * #Upscale, texture
Global v1size = #nx * #ny * SizeOf(float), v2size = v1size * 2, onframe

Global Dim old_velocity.vec2f(#ny1,#nx1), Dim new_velocity.vec2f(#ny1,#nx1)
Global Dim old_density.f(#ny1, #nx1), Dim new_density.f(#ny1, #nx1)
Global Dim pixels.l(#ny1,#nx1)

Macro LengthM(x,y) : Sqr((x*x) + (y*y)) : EndMacro
Macro LerpM(a, b, t) : ( a + t * (b - a) ) : EndMacro
Macro Log2M(v) : (Log(v) / 0.6931471805599453094) : EndMacro ; not sure if this is correct Log2...
Macro RandomFloatM(a, b) : LerpM(a, b, Random(1000000) / 1000000.0) : EndMacro
 
Procedure Init()
  Protected x, y
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_density(y,x) = 0.0 : old_velocity(y,x)\x = 0 : old_velocity(y,x)\y = 0
  Next : Next

  OpenWindow(0,0,0,w,h,"Fluid-Simulation (Flame)",#PB_Window_ScreenCentered|#PB_Window_SystemMenu)
  OpenGLGadget(0,0,0,w,h);,#PB_OpenGL_NoFlipSynchronization)
  glOrtho_(0,1,1,0,0,1)
  glEnable_(#GL_TEXTURE_2D);
  
  glGenTextures_(1, @texture);
  glBindTexture_(#GL_TEXTURE_2D, texture);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MIN_FILTER, #GL_LINEAR);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_MAG_FILTER, #GL_LINEAR);
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_S, #GL_REPEAT)    ;
  glTexParameteri_(#GL_TEXTURE_2D, #GL_TEXTURE_WRAP_T, #GL_REPEAT)    ;
  glTexImage2D_(#GL_TEXTURE_2D, 0, #GL_RGBA, #nx, #ny, 0, #GL_RGBA, #GL_UNSIGNED_BYTE, #Null);
EndProcedure
Procedure.f interpolate(Array grid.f(2), *p.vec2f, *d.float)
  Protected ix = Round(*p\x,#PB_Round_Down), iy = Round(*p\y,#PB_Round_Down), ux.f = *p\x - ix, uy.f = *p\y - iy
  *d\f = LerpM(LerpM(grid(iy,ix),grid(iy,ix + 1), ux),LerpM(grid(iy+1,ix),grid(iy+1,ix+1), ux), uy)
EndProcedure
Procedure.f interpolateVec2(Array grid.vec2f(2), *p.vec2f, *d.vec2f)
  Protected ix = Round(*p\x,#PB_Round_Down), iy = Round(*p\y,#PB_Round_Down), ux.f = *p\x - ix, uy.f = *p\y - iy
  *d\x = LerpM(LerpM(grid(iy,ix)\x,grid(iy,ix+1)\x, ux), LerpM(grid(iy+1,ix)\x,grid(iy+1,ix+1)\x, ux), uy)
  *d\y = LerpM(LerpM(grid(iy,ix)\y,grid(iy,ix+1)\y, ux), LerpM(grid(iy+1,ix)\y,grid(iy+1,ix+1)\y, ux), uy)
EndProcedure

Procedure advect_density()
  Protected pos.vec2f, x, y
  For y = 1 To #ny2 : For x = 1 To #nx2
      pos\x = x - dt * old_velocity(y,x)\x
      pos\y = y - dt * old_velocity(y,x)\y
      If pos\y < 1 : pos\y = 1 : ElseIf pos\y > #ny2 : pos\y = #ny2 : EndIf 
      If pos\x < 1 : pos\x = 1 : ElseIf pos\x > #nx2 : pos\x = #nx2 : EndIf
      interpolate(old_density(), @pos, @new_density(y,x));
  Next : Next
  CopyMemory(@new_density(),@old_density(),v1size)
EndProcedure
Procedure advect_velocity()
  Protected pos.vec2f, x, y
  For y = 1 To #ny2 : For x = 1 To #nx2
      pos\x = x - dt * old_velocity(y,x)\x;
      pos\y = y - dt * old_velocity(y,x)\y;
      If pos\y < 1 : pos\y = 1 : ElseIf pos\y > #ny2 : pos\y = #ny2 : EndIf 
      If pos\x < 1 : pos\x = 1 : ElseIf pos\x > #nx2 : pos\x = #nx2 : EndIf
      interpolateVec2(old_velocity(), @pos,@new_velocity(y,x))
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure diffuse_density()
  Protected pos.vec2f, x, y, diffusion.f = dt * 100.01, sum.f
  For y = 1 To #ny2 : For x = 1 To #nx2
      sum = diffusion * (old_density(y, x - 1) + old_density(y, x + 1) + old_density(y - 1, x) + old_density(y + 1, x)) + old_density(y,x);
      new_density(y, x) = 1.0 / (1.0 + 4.0 * diffusion) * sum 
  Next : Next
  CopyMemory(@new_density(),@old_density(),v1size)
EndProcedure
Procedure diffuse_velocity()
  Protected viscosity = dt*0.000001, x, y, sum.vec2f
  For y = 1 To #ny2 : For x = 1 To #nx2
      sum\x = viscosity * (old_velocity(y, x - 1)\x + old_velocity(y, x + 1)\x + old_velocity(y - 1, x)\x + old_velocity(y + 1, x)\x) + old_velocity(y, x)\x
      sum\y = viscosity * (old_velocity(y, x - 1)\y + old_velocity(y, x + 1)\y + old_velocity(y - 1, x)\y + old_velocity(y + 1, x)\y) + old_velocity(y, x)\y
      new_velocity(y,x)\x = 1.0 / (1.0 + 4.0 * viscosity) * sum\x
      new_velocity(y,x)\y = 1.0 / (1.0 + 4.0 * viscosity) * sum\y
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure project_velocity()
  Protected x, y, dx.f, dy.f, k, sum.f
  Static Dim p.f(#ny,#nx), Dim p2.f(#ny,#nx), Dim div.f(#ny,#nx)

  For y = 1 To #ny2 : For x = 1 To #nx2
      dx = old_velocity(y    , x + 1)\x - old_velocity(y    ,x - 1)\x
      dy = old_velocity(y + 1, x    )\y - old_velocity(y - 1,x    )\y
      div(y,x) = dx + dy : p(y,x) = 0.0
  Next : Next
  For k = 1 To iterations
    For y = 1 To #ny2 : For x = 1 To #nx2
        sum = -div(y,x) + p(y ,x + 1) + p(y ,x - 1) + p(y + 1,x) + p(y - 1,x)
        p2(y,x) = 0.25 * sum
    Next : Next
    CopyMemory(@p2(),@p(),v1size)
  Next
  For y = 1 To #ny2 : For x = 1 To #nx2
      old_velocity(y,x)\x - (0.5 * (p(y    ,x + 1) - p(y    , x - 1)))
      old_velocity(y,x)\y - (0.5 * (p(y + 1,x  ) - p(y - 1, x )))
  Next : Next
EndProcedure
Procedure.f curl(y,x)
  ProcedureReturn old_velocity(y + 1,x)\x - old_velocity(y - 1,x)\x + old_velocity(y,x - 1)\y - old_velocity(y, x + 1)\y;
EndProcedure
Procedure vorticity_confinement()
  Protected direction.vec2f, length.f, x, y;
  Static Dim abs_curl.f(#ny1,#nx1)
  For y = 1 To #ny2 : For x = 1 To #nx2
      abs_curl(y,x) = Abs(curl(y,x));
  Next : Next
  
  For y = 1 To #ny2 : For x = 1 To #nx2
      direction\x = abs_curl(y - 1, x    ) - abs_curl(y + 1, x    );
      direction\y = abs_curl(y    , x + 1) - abs_curl(y    , x - 1);
      length = lengthM(direction\x,direction\y)
      direction\x = vorticity / (length + 1e-5) * direction\x;
      direction\y = vorticity / (length + 1e-5) * direction\y;
      
      If x < #nx * 0.5 : direction\x = 0.0 : direction\y = 0.0 : EndIf
      
      new_velocity(y,x)\x = old_velocity(y,x)\x + dt * curl(y,x) * direction\x;
      new_velocity(y,x)\y = old_velocity(y,x)\y + dt * curl(y,x) * direction\y;
  Next : Next
  CopyMemory(@new_velocity(),@old_velocity(),v2size)
EndProcedure
Procedure.f Clamp(v.f,Low.f = 0.0, High.f = 1.0)
  If v < Low : ProcedureReturn Low : ElseIf v > High : ProcedureReturn High : EndIf
  ProcedureReturn v
EndProcedure
Procedure.f SmoothStep(Low.f, High.f, v.f)
  V = Clamp((v - low) / (high - low));
  ProcedureReturn v * v * (3.0 - 2.0 * v);
EndProcedure

Procedure add_density(px, py, r = 10, value.f = 0.5)
  Protected x, y, u.f, rf.f = r, x2, y2
  For y = -r To r : For x = -r To  r
      y2 = py + y : x2 = px + x
      If x2 < 0 : x2 = 0 : ElseIf x2 > #nx1 : x2 = #nx1 : EndIf
      If y2 < 0 : y2 = 0 : ElseIf y2 > #ny1 : y2 = #ny1 : EndIf
      u = SmoothStep(rf, 0.0, Sqr(x*x + y*y));
      old_density(y2, x2) + u * value
  Next : Next
  
EndProcedure

Procedure fluid_simulation_step()
  Protected x,y, r.f = 10.0 ; r is resolution dependant
  
  For y = 0 To #ny1 : For x = 0 To #nx1
      If x > #nx * 0.5
        old_velocity(y,x)\x + RandomFloatM(-r, r)
        old_velocity(y,x)\y + RandomFloatM(-r, r)
      EndIf
  Next : Next
  
  ;// dense regions rise up
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_velocity(y,x)\y + (old_density(y,x) * 20.0 - 5.0) * dt;
  Next : Next
 
  ;// fast movement is dampened
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_velocity(y,x)\x * 0.999 : old_velocity(y,x)\y * 0.999;
  Next : Next
  
  ;// fade away
  For y = 0 To #ny1 : For x = 0 To #nx1
      old_density(y,x) * 0.99;
  Next : Next

  add_density(#nx*0.25, #ny*0.1, 10.0, 0.5);
  add_density(#nx*0.75, #ny*0.1, 10.0, 0.5);
  
  If mouse\x > 0 And mouse\y > 0 And mouse\x < #nx And mouse\y < #ny1 : add_density(Mouse\x, #ny - Mouse\y, #nx * 0.025,0.5) : EndIf
  
  vorticity_confinement()
  ;diffuse_velocity()
  ;project_velocity()
  advect_velocity()
  project_velocity()
  ;diffuse_density();
  advect_density()

  ;/ wipe out borders as prevents crashes
  For y = 0 To #ny1 
      old_density(y,0) = 0.0 : old_velocity(y,0)\x = 0 : old_velocity(y,0)\y = 0      
      old_density(y,1) = 0.0 : old_velocity(y,1)\x = 0 : old_velocity(y,1)\y = 0    
      old_density(y,#nx1) = 0.0 : old_velocity(y,#nx1)\x = 0 : old_velocity(y,#nx1)\y = 0      
      old_density(y,#nx2) = 0.0 : old_velocity(y,#nx2)\x = 0 : old_velocity(y,#nx2)\y = 0      
  Next
  For x = 0 To #nx1 
      old_density(0,x) = 0.0 : old_velocity(0,x)\x = 0 : old_velocity(0,x)\y = 0      
      old_density(#ny1,x) = 0.0 : old_velocity(#ny1,x)\x = 0 : old_velocity(#ny1,x)\y = 0      
      old_density(1,x) = 0.0 : old_velocity(1,x)\x = 0 : old_velocity(1,x)\y = 0      
      old_density(#ny2,x) = 0.0 : old_velocity(#ny2,x)\x = 0 : old_velocity(#ny2,x)\y = 0      
  Next
EndProcedure

Procedure.l RGBAgl(r.f, g.f, b.f, a.f)
  r = clamp(r * 255, 0, 255) : g = clamp(g * 255, 0, 255) : b = clamp(b * 255, 0, 255) : a = clamp(a * 255, 0, 255)
  ProcedureReturn (Int(a) << 24) | (Int(b) << 16) | (Int(g) << 8) | Int(r);
EndProcedure
Procedure Render()
  ;// density field To pixels
  Protected x, y, f.f, f3.f, r.f, g.f, b.f
  For y = 0 To #ny1 : For x = 0 To #nx1
      f = old_density(y,x);
      f = log2m( f * 0.25 + 1.0);
      f3 = f*f*f
      r = 1.5 * f : g = 1.5 * f3 : b = (f3 *f3)
      pixels(y,x) = RGBAgl(r, g, b, 1.0);
  Next : Next
  
  ;// upload pixels() array To texture
  glTexSubImage2D_(#GL_TEXTURE_2D, 0, 0, 0, #nx, #ny, #GL_RGBA, #GL_UNSIGNED_BYTE, @pixels());
  
  ;// draw texture
  glBegin_(#GL_QUADS);
  glTexCoord2f_(0.0, 1.0): glVertex2f_(0  , 0)
  glTexCoord2f_(1.0, 1.0): glVertex2f_(1.0, 0)
  glTexCoord2f_(1.0, 0.0): glVertex2f_(1.0, 1.0)
  glTexCoord2f_(0.0, 0.0): glVertex2f_(0  , 1.0)
  glEnd_()
  
  SetGadgetAttribute(0,#PB_OpenGL_FlipBuffers,#True)
EndProcedure
Init()

Define Event, Exit, time.i

Repeat
  Repeat
    Event = WindowEvent()
    Select Event
      Case #PB_Event_CloseWindow : Exit = #True
      Case #PB_Event_Gadget
        Select EventType()
          Case #PB_EventType_MouseMove
            Mouse\x = (GetGadgetAttribute(0,#PB_OpenGL_MouseX) / (w / #nx)) / DesktopResolutionX()
            Mouse\y = (GetGadgetAttribute(0,#PB_OpenGL_MouseY) / (h / #ny)) / DesktopResolutionY()
          Case #PB_EventType_LeftButtonDown
            LMB = #True
            add_density(Mouse\x, #ny - Mouse\y, #nx * 0.2,8.0); burst
          Case #PB_EventType_LeftButtonUp
            LMB = #False
        EndSelect
    EndSelect
  Until Event = 0

  time = ElapsedMilliseconds()
  fluid_simulation_step();
  Render() : onframe + 1
  time = ElapsedMilliseconds() - time
  If onframe = 60 : onframe = 1
    SetWindowTitle(0,"Fluid dynamics w/vorticity - Phil James 2023 - Process time: "+ Str(time)+"ms.")
  EndIf
  
Until Exit = #True
glDeleteTextures_(1,@texture)

Edited to add cross-platform hints from juergenkulow & Shardik