Fluid 2D with the method of Navier Stokes

threedslider · Post by **threedslider** » Tue Aug 12, 2025 7:25 pm

Hi

I share you my code for Fluid 2D with Navier Stokes (not really true but it works

)

; --- Paramètres ---
#N = 128
#Scale = 8
#Size = #N * #Scale
#Dt = 0.1
#Diff = 0.0001
#Visc = 0.0001
#Iterations = 10

Structure rond
  x.i
  y.i
  rayon.i
EndStructure


Myvar.rond

myvar\x = 10
myvar\y = 10
myvar\rayon = 5

; --- Tableaux globaux ---
Global Dim u.f(#N+1, #N+1)
Global Dim v.f(#N+1, #N+1)
Global Dim u_prev.f(#N+1, #N+1)
Global Dim v_prev.f(#N+1, #N+1)
Global Dim dens.f(#N+1, #N+1)
Global Dim dens_prev.f(#N+1, #N+1)


; --- Macro ---
Macro SetIf(cond, var, trueVal, falseVal)
  If cond
    var = trueVal
  Else
    var = falseVal
  EndIf
EndMacro


; --- Fonctions ---
Procedure CopyField(Array dst.f(2), Array src.f(2))
  For i = 0 To #N
    For j = 0 To #N
      dst(i, j) = src(i, j)
    Next
  Next
EndProcedure

Procedure AddSource(Array x.f(2), Array s.f(2), dt.f)
  For i = 0 To #N
    For j = 0 To #N
      x(i,j) + s(i,j) * dt
    Next
  Next
EndProcedure

Procedure SetBound(b, Array x.f(2))
   For i = 1 To #N-1
    SetIf(b = 1, x(0, i),  -x(1, i),  x(1, i))
    SetIf(b = 1, x(#N, i), -x(#N-1, i), x(#N-1, i))
    SetIf(b = 2, x(i, 0),  -x(i, 1),  x(i, 1))
    SetIf(b = 2, x(i, #N), -x(i, #N-1), x(i, #N-1))
  Next
  
  x(0,0)      = 0.5 * (x(1,0) + x(0,1))
  x(0,#N)     = 0.5 * (x(1,#N) + x(0,#N-1))
  x(#N,0)     = 0.5 * (x(#N-1,0) + x(#N,1))
  x(#N,#N)    = 0.5 * (x(#N-1,#N) + x(#N,#N-1))
EndProcedure

Procedure LinSolve(b, Array x.f(2), Array x0.f(2), a.f, c.f)
  For k = 0 To #Iterations
    For i = 1 To #N-1
      For j = 1 To #N-1
        x(i,j) = (x0(i,j) + a * (x(i-1,j) + x(i+1,j) + x(i,j-1) + x(i,j+1))) / c
      Next
    Next
    SetBound(b, x())
  Next
EndProcedure

Procedure Diffuse(b, Array x.f(2), Array x0.f(2), diff.f, dt.f)
  a.f = dt * diff * #N * #N
  LinSolve(b, x(), x0(), a, 1 + 4 * a)
EndProcedure

Procedure Advect(b, Array d.f(2), Array d0.f(2), Array u.f(2), Array v.f(2), dt.f)
  Protected i, j, i0, j0, i1, j1
  Protected x.f, y.f, s0.f, t0.f, s1.f, t1.f
  Protected dt0.f = dt * #N

  For i = 1 To #N-1
    For j = 1 To #N-1
      x = i - dt0 * u(i,j)
      y = j - dt0 * v(i,j)
      If x < 0.5 : x = 0.5 : EndIf
      If x > #N + 0.5 : x = #N + 0.5 : EndIf
      i0 = Int(x) : i1 = i0 + 1
      If y < 0.5 : y = 0.5 : EndIf
      If y > #N + 0.5 : y = #N + 0.5 : EndIf
      j0 = Int(y) : j1 = j0 + 1
      s1 = x - i0 : s0 = 1 - s1
      t1 = y - j0 : t0 = 1 - t1
      d(i,j) = s0 * (t0 * d0(i0,j0) + t1 * d0(i0,j1)) + s1 * (t0 * d0(i1,j0) + t1 * d0(i1,j1))
    Next
  Next
  SetBound(b, d())
EndProcedure

Procedure Project(Array u.f(2), Array v.f(2), Array p.f(2), Array div.f(2))
  For i = 1 To #N-1
    For j = 1 To #N-1
      div(i,j) = -0.5 * (u(i+1,j) - u(i-1,j) + v(i,j+1) - v(i,j-1)) / #N
      p(i,j) = 0
    Next
  Next
  SetBound(0, div()) : SetBound(0, p())
  LinSolve(0, p(), div(), 1, 4)
  For i = 1 To #N-1
    For j = 1 To #N-1
      u(i,j) - 0.5 * #N * (p(i+1,j) - p(i-1,j))
      v(i,j) - 0.5 * #N * (p(i,j+1) - p(i,j-1))
    Next
  Next
  SetBound(1, u()) : SetBound(2, v())
EndProcedure

Procedure VelocityStep(Array u.f(2), Array v.f(2), Array u0.f(2), Array v0.f(2), visc.f, dt.f)
  AddSource(u(), u0(), dt)
  AddSource(v(), v0(), dt)
  CopyField(u0(), u()) : Diffuse(1, u(), u0(), visc, dt)
  CopyField(v0(), v()) : Diffuse(2, v(), v0(), visc, dt)
  Project(u(), v(), u0(), v0())
  CopyField(u0(), u())
  CopyField(v0(), v())
  Advect(1, u(), u0(), u0(), v0(), dt)
  Advect(2, v(), v0(), u0(), v0(), dt)
  Project(u(), v(), u0(), v0())
EndProcedure

Procedure DensityStep(Array x.f(2), Array x0.f(2), Array u.f(2), Array v.f(2), diff.f, dt.f)
  ;Debug x(0,1)
  AddSource(x(), x0(), dt)
  CopyField(x0(), x()) : Diffuse(0, x(), x0(), diff, dt)
  CopyField(x0(), x()) : Advect(0, x(), x0(), u(), v(), dt)
EndProcedure

Procedure DrawDensity()
  Protected i, j, col
  For i = 0 To #N-1
    For j = 0 To #N-1
      col = Int(dens(i,j) * 255)      
      If col > 0
        If col > 255 : col = 255 : EndIf
        Box(i * #Scale, j * #Scale, #Scale, #Scale, RGB(col, col, col))
      EndIf
    Next
  Next
EndProcedure

InitMouse()
InitSprite()

; --- Programme principal ---
OpenWindow(0, 0, 0, #Size/DesktopResolutionX(), #Size/DesktopResolutionY(), "Fluide 2D - PureBasic", #PB_Window_SystemMenu| #PB_Window_ScreenCentered)
OpenWindowedScreen(WindowID(0), 0, 0, #Size, #Size, 0, 0, 0)

Global running.b = #True

Repeat 
  Event = WindowEvent()
 
  ExamineMouse()
  
  For i = 0 To #N
    For j = 0 To #N
      ;u(i,j) = 0
      ;v(i,j) = 0
      ;dens(i,j) = 0
      dens_prev(i,j) = 0
      u_prev(i,j) = 0
      v_prev(i,j) = 0
    Next
  Next

  If MouseButton(#PB_MouseButton_Left)
    mx = Int(MouseX() / #Scale)
    my = Int(MouseY() / #Scale)
    If mx >= 1 And mx < #N-1 And my >= 1 And my < #N-1
      
      dens_prev(mx,my) = 100
      
      u_prev(mx,my) = 5
      v_prev(mx,my) = 5
      
     ; Debug "Inject at " + Str(mx) + "," + Str(my) + " dens = " + StrF(dens_prev(mx, my))
    EndIf
  EndIf

  VelocityStep(u(), v(), u_prev(), v_prev(), #Visc, #Dt)
  DensityStep(dens(), dens_prev(), u(), v(), #Diff, #Dt)

  ClearScreen(RGB(0.0,0.0,0.0))
  StartDrawing(ScreenOutput())
  DrawDensity()
  StopDrawing()
  FlipBuffers()

  
Until event = #PB_Event_CloseWindow

Warning it is quite unstable but work though

Enjoy and happy coding !

SPH · Post by **SPH** » Tue Aug 12, 2025 7:58 pm

No way

miso · Post by **miso** » Tue Aug 12, 2025 8:41 pm

Cool

Thanks for sharing.
I also made a test, where I replaced the 2d drawcalls with sprite draws to test the speed. This asks for c backend + optimization and debugger off.
Will quit if escape is pressed.

Code: Select all

If #PB_Compiler_Backend <> #PB_Backend_C Or #PB_Compiler_Optimizer<>1 Or #PB_Compiler_Debugger = 1
  OpenConsole("Warning")
  PrintN("Please compile with C backend, optimization on and without debugger")
  PrintN("Press Enter to quit")
  Input()
  End
EndIf


; --- Paramètres ---
#N = 128
#Scale = 8
#Size = #N * #Scale
#Dt = 0.1
#Diff = 0.0001
#Visc = 0.0001
#Iterations = 10
#SPRITE_ID = 1

Structure rond
  x.i
  y.i
  rayon.i
EndStructure


Myvar.rond

myvar\x = 10
myvar\y = 10
myvar\rayon = 5

; --- Tableaux globaux ---
Global Dim u.f(#N+1, #N+1)
Global Dim v.f(#N+1, #N+1)
Global Dim u_prev.f(#N+1, #N+1)
Global Dim v_prev.f(#N+1, #N+1)
Global Dim dens.f(#N+1, #N+1)
Global Dim dens_prev.f(#N+1, #N+1)


; --- Macro ---
Macro SetIf(cond, var, trueVal, falseVal)
  If cond
    var = trueVal
  Else
    var = falseVal
  EndIf
EndMacro


; --- Fonctions ---
Procedure CopyField(Array dst.f(2), Array src.f(2))
  For i = 0 To #N
    For j = 0 To #N
      dst(i, j) = src(i, j)
    Next
  Next
EndProcedure

Procedure AddSource(Array x.f(2), Array s.f(2), dt.f)
  For i = 0 To #N
    For j = 0 To #N
      x(i,j) + s(i,j) * dt
    Next
  Next
EndProcedure

Procedure SetBound(b, Array x.f(2))
   For i = 1 To #N-1
    SetIf(b = 1, x(0, i),  -x(1, i),  x(1, i))
    SetIf(b = 1, x(#N, i), -x(#N-1, i), x(#N-1, i))
    SetIf(b = 2, x(i, 0),  -x(i, 1),  x(i, 1))
    SetIf(b = 2, x(i, #N), -x(i, #N-1), x(i, #N-1))
  Next
  
  x(0,0)      = 0.5 * (x(1,0) + x(0,1))
  x(0,#N)     = 0.5 * (x(1,#N) + x(0,#N-1))
  x(#N,0)     = 0.5 * (x(#N-1,0) + x(#N,1))
  x(#N,#N)    = 0.5 * (x(#N-1,#N) + x(#N,#N-1))
EndProcedure

Procedure LinSolve(b, Array x.f(2), Array x0.f(2), a.f, c.f)
  For k = 0 To #Iterations
    For i = 1 To #N-1
      For j = 1 To #N-1
        x(i,j) = (x0(i,j) + a * (x(i-1,j) + x(i+1,j) + x(i,j-1) + x(i,j+1))) / c
      Next
    Next
    SetBound(b, x())
  Next
EndProcedure

Procedure Diffuse(b, Array x.f(2), Array x0.f(2), diff.f, dt.f)
  a.f = dt * diff * #N * #N
  LinSolve(b, x(), x0(), a, 1 + 4 * a)
EndProcedure

Procedure Advect(b, Array d.f(2), Array d0.f(2), Array u.f(2), Array v.f(2), dt.f)
  Protected i, j, i0, j0, i1, j1
  Protected x.f, y.f, s0.f, t0.f, s1.f, t1.f
  Protected dt0.f = dt * #N

  For i = 1 To #N-1
    For j = 1 To #N-1
      x = i - dt0 * u(i,j)
      y = j - dt0 * v(i,j)
      If x < 0.5 : x = 0.5 : EndIf
      If x > #N + 0.5 : x = #N + 0.5 : EndIf
      i0 = Int(x) : i1 = i0 + 1
      If y < 0.5 : y = 0.5 : EndIf
      If y > #N + 0.5 : y = #N + 0.5 : EndIf
      j0 = Int(y) : j1 = j0 + 1
      s1 = x - i0 : s0 = 1 - s1
      t1 = y - j0 : t0 = 1 - t1
      d(i,j) = s0 * (t0 * d0(i0,j0) + t1 * d0(i0,j1)) + s1 * (t0 * d0(i1,j0) + t1 * d0(i1,j1))
    Next
  Next
  SetBound(b, d())
EndProcedure

Procedure Project(Array u.f(2), Array v.f(2), Array p.f(2), Array div.f(2))
  For i = 1 To #N-1
    For j = 1 To #N-1
      div(i,j) = -0.5 * (u(i+1,j) - u(i-1,j) + v(i,j+1) - v(i,j-1)) / #N
      p(i,j) = 0
    Next
  Next
  SetBound(0, div()) : SetBound(0, p())
  LinSolve(0, p(), div(), 1, 4)
  For i = 1 To #N-1
    For j = 1 To #N-1
      u(i,j) - 0.5 * #N * (p(i+1,j) - p(i-1,j))
      v(i,j) - 0.5 * #N * (p(i,j+1) - p(i,j-1))
    Next
  Next
  SetBound(1, u()) : SetBound(2, v())
EndProcedure

Procedure VelocityStep(Array u.f(2), Array v.f(2), Array u0.f(2), Array v0.f(2), visc.f, dt.f)
  AddSource(u(), u0(), dt)
  AddSource(v(), v0(), dt)
  CopyField(u0(), u()) : Diffuse(1, u(), u0(), visc, dt)
  CopyField(v0(), v()) : Diffuse(2, v(), v0(), visc, dt)
  Project(u(), v(), u0(), v0())
  CopyField(u0(), u())
  CopyField(v0(), v())
  Advect(1, u(), u0(), u0(), v0(), dt)
  Advect(2, v(), v0(), u0(), v0(), dt)
  Project(u(), v(), u0(), v0())
EndProcedure

Procedure DensityStep(Array x.f(2), Array x0.f(2), Array u.f(2), Array v.f(2), diff.f, dt.f)
  ;Debug x(0,1)
  AddSource(x(), x0(), dt)
  CopyField(x0(), x()) : Diffuse(0, x(), x0(), diff, dt)
  CopyField(x0(), x()) : Advect(0, x(), x0(), u(), v(), dt)
EndProcedure

Procedure DrawDensity()
  Protected i, j, col
  For i = 0 To #N-1
    For j = 0 To #N-1
      col = Int(dens(i,j) * 255)      
      If col > 0
        If col > 255 : col = 255 : EndIf
          DisplayTransparentSprite(#SPRITE_ID,i*#scale,j*#scale,255,RGBA(col,col,col,255))
        ;Box(i * #Scale, j * #Scale, #Scale, #Scale, RGB(col, col, col))
      EndIf
    Next
  Next
EndProcedure

InitKeyboard()
InitMouse()
InitSprite()

; --- Programme principal ---
OpenWindow(0, 0, 0, #Size/DesktopResolutionX(), #Size/DesktopResolutionY(), "Fluide 2D - PureBasic", #PB_Window_SystemMenu| #PB_Window_ScreenCentered)
OpenWindowedScreen(WindowID(0), 0, 0, #Size, #Size, 0, 0, 0)

CreateSprite(#SPRITE_ID,1,1,#PB_Sprite_AlphaBlending)
StartDrawing(SpriteOutput(#SPRITE_ID))
DrawingMode(#PB_2DDrawing_AllChannels)
Box(0,0,OutputWidth(),OutputHeight(),RGBA(255,255,255,255))
StopDrawing()
ZoomSprite(#SPRITE_ID,#Scale,#Scale)

Global running.b = #True

Repeat 
  Repeat
    Event = WindowEvent()
    If event = #PB_Event_CloseWindow : End : EndIf
  Until event = 0
 
  ExamineMouse()
  ExamineKeyboard()
  
  For i = 0 To #N
    For j = 0 To #N
      ;u(i,j) = 0
      ;v(i,j) = 0
      ;dens(i,j) = 0
      dens_prev(i,j) = 0
      u_prev(i,j) = 0
      v_prev(i,j) = 0
    Next
  Next

  If MouseButton(#PB_MouseButton_Left)
    mx = Int(MouseX() / #Scale)
    my = Int(MouseY() / #Scale)
    If mx >= 1 And mx < #N-1 And my >= 1 And my < #N-1
      
      dens_prev(mx,my) = 100
      
      u_prev(mx,my) = 5
      v_prev(mx,my) = 5
      
     ; Debug "Inject at " + Str(mx) + "," + Str(my) + " dens = " + StrF(dens_prev(mx, my))
    EndIf
  EndIf

  VelocityStep(u(), v(), u_prev(), v_prev(), #Visc, #Dt)
  DensityStep(dens(), dens_prev(), u(), v(), #Diff, #Dt)

  ClearScreen(RGB(0.0,0.0,0.0))

  DrawDensity()
  DisplayTransparentSprite(#SPRITE_ID,MouseX(),MouseY())
  FlipBuffers()

  
Until KeyboardPushed(#PB_Key_Escape)

threedslider · Post by **threedslider** » Wed Aug 13, 2025 9:16 am

@miso : Wow ! Thanks a lot ! It is more stable than mine

Thank you for sharing

Caronte3D · Post by **Caronte3D** » Wed Aug 13, 2025 11:00 am

Very nice!

PureBasic Forums - English

Fluid 2D with the method of Navier Stokes

Fluid 2D with the method of Navier Stokes

Re: Fluid 2D with the method of Navier Stokes

Re: Fluid 2D with the method of Navier Stokes

Re: Fluid 2D with the method of Navier Stokes

Re: Fluid 2D with the method of Navier Stokes