Seite 1 von 4

Fluid Dynamics

Verfasst: 11.09.2006 16:34
von remi_meier
Ich habe auf der Website von PlasmaPong diesen Link gefunden:
http://www.dgp.toronto.edu/people/stam/ ... /GDC03.pdf
daraufhin habe ich den C-Code dort übersetzt nach PB. Das Resultat ist
hier (Maustaste drücken und bewegen):

Code: Alles auswählen

EnableExplicit

Structure FArray
  f.f[0]
EndStructure


Macro IX(i, j)
  ((i) + (n + 2) * (j))
EndMacro

Declare set_bnd(n.l, b.l, *x.FArray)
Declare project(n.l, *u.FArray, *v.FArray, *p.FArray, *div.FArray)

Procedure add_source(n.l, *x.FArray, *s.FArray, dt.f)
  Define.l i, Size = (n + 2) * (n + 2)
  
  For i = 0 To Size - 1
    *x\f[i] + dt * *s\f[i]
  Next
EndProcedure


Procedure diffuse(n.l, b.l, *x.FArray, *x0.FArray, diff.f, dt.f)
  Define.l i, j, k
  Define.f a
  a = dt * diff * n * n
  
  For k = 0 To 19
    For i = 1 To n
      For j = 1 To n
        *x\f[IX(i, j)] = (*x0\f[IX(i, j)] + a * (*x\f[IX(i-1, j)] + *x\f[IX(i+1, j)] + *x\f[IX(i, j-1)] + *x\f[IX(i, j+1)])) / (1.0 + 4.0 * a)
        ;*x\f[IX(i, j)] = *x0\f[IX(i, j)] + a * (*x0\f[IX(i-1, j)] + *x0\f[IX(i+1, j)] + *x0\f[IX(i, j-1)] + *x0\f[IX(i, j+1)] - 4 * *x0\f[IX(i, j)])
      Next
    Next
    set_bnd(n, b, *x)
  Next
EndProcedure

Procedure advect (n.l, b.l, *d.FArray, *d0.FArray, *u.FArray, *v.FArray, dt.f)
  Define.l i, j, i0, j0, i1, j1
  Define.f x, y, s0, t0, s1, t1, dt0
  dt0 = dt * n
  
  For i = 1 To n
    For j = 1 To n
      x = i - dt0 * *u\f[IX(i, j)]
      y = j - dt0 * *v\f[IX(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\f[IX(i, j)] = s0 * (t0 * *d0\f[IX(i0, j0)] + t1 * *d0\f[IX(i0, j1)]) + s1 * (t0 * *d0\f[IX(i1, j0)] + t1 * *d0\f[IX(i1, j1)])
      
    Next
  Next
  set_bnd(n, b, *d)
EndProcedure

Procedure dens_step(n.l, *x.FArray, *x0.FArray, *u.FArray, *v.FArray, diff.f, dt.f)
  add_source(n, *x, *x0, dt)
  diffuse(n, 0, *x0, *x, diff, dt )
  advect(n, 0, *x, *x0, *u, *v, dt )
EndProcedure

Procedure vel_step(n.l, *u.FArray, *v.FArray, *u0.FArray, *v0.FArray, visc.f, dt.f)
  add_source(n, *u, *u0, dt) : add_source(n, *v, *v0, dt)
  diffuse(n, 1, *u0, *u, visc, dt)
  diffuse(n, 2, *v0, *v, visc, dt)
  project(n, *u0, *v0, *u, *v)
  
  advect(n, 1, *u, *u0, *u0, *v0, dt) : advect(n, 2, *v, *v0, *u0, *v0, dt)
  project(n, *u, *v, *u0, *v0)
EndProcedure


Procedure project(n.l, *u.FArray, *v.FArray, *p.FArray, *div.FArray)
  Define.l i, j, k
  Define.f h = 1.0 / n
  
  For i = 1 To n
    For j = 1 To n
      *div\f[IX(i, j)] = -0.5 * h * (*u\f[IX(i+1, j)] - *u\f[IX(i-1, j)] + *v\f[IX(i, j+1)] - *v\f[IX(i, j-1)])
      *p\f[IX(i, j)] = 0
    Next
  Next
  set_bnd(n, 0, *div) : set_bnd(n, 0, *p)
  
  For k = 0 To 19
    For i = 1 To n
      For j = 1 To n
        *p\f[IX(i, j)] = (*div\f[IX(i, j)] + *p\f[IX(i-1, j)] + *p\f[IX(i+1, j)] + *p\f[IX(i, j-1)] + *p\f[IX(i, j+1)]) / 4
      Next
    Next
    set_bnd(n, 0, *p)
  Next
  
  For i = 1 To n
    For j = 1 To n
      *u\f[IX(i, j)] - (0.5 * (*p\f[IX(i+1, j)] - *p\f[IX(i-1, j)]) / h)
      *v\f[IX(i, j)] - (0.5 * (*p\f[IX(i, j+1)] - *p\f[IX(i, j-1)]) / h)
    Next
  Next
  set_bnd(n, 1, *u) : set_bnd(n, 2, *v)
EndProcedure


Procedure set_bnd(n.l, b.l, *x.FArray)
  Define.l i
  
  For i = 1 To n
    If b = 1
      *x\f[IX(0, i)]   = -*x\f[IX(1, i)]
      *x\f[IX(n+1, i)] = -*x\f[IX(n, i)]
    Else
      *x\f[IX(0, i)]   = *x\f[IX(1, i)]
      *x\f[IX(n+1, i)] = *x\f[IX(n, i)]
    EndIf
    
    If b = 2
      *x\f[IX(i,   0)] = -*x\f[IX(i, 1)]
      *x\f[IX(i, n+1)] = -*x\f[IX(i, n)]
    Else
      *x\f[IX(i,   0)] = *x\f[IX(i, 1)]
      *x\f[IX(i, n+1)] = *x\f[IX(i, n)]
    EndIf
  Next
  
  *x\f[IX(0,     0)] = 0.5 * (*x\f[IX(1,   0)] + *x\f[IX(0,   1)])
  *x\f[IX(0,   n+1)] = 0.5 * (*x\f[IX(1, n+1)] + *x\f[IX(0,   n)])
  *x\f[IX(n+1,   0)] = 0.5 * (*x\f[IX(n,   0)] + *x\f[IX(n+1, 1)])
  *x\f[IX(n+1, n+1)] = 0.5 * (*x\f[IX(n, n+1)] + *x\f[IX(n+1, n)])
EndProcedure

;- End Simulation
#N = 60
#SIZE = (#N + 2) * (#N + 2)


Procedure array_normalize(n.l, *a.FArray)
  Define.l i, j, Size = (n + 2) * (n + 2)
  Define.f max
  
  For i = 0 To n - 1
    For j = 0 To n - 1
      If Abs(*a\f[IX(i, j)]) > max
        max = Abs(*a\f[IX(i, j)])
      EndIf
    Next
  Next
  
  For i = 0 To Size - 1
    *a\f[i] / (max + 0.001)
  Next
EndProcedure


Procedure draw_dens(n.l, *dens.FArray)
  Define.l i, j
  Define.f max
  
  ClearScreen(0)
  ExamineKeyboard()
  If KeyboardPushed(#PB_Key_Escape)
    End
  EndIf
  
  StartDrawing(ScreenOutput())
  array_normalize(n, *dens)
  For i = 0 To n - 1
    For j = 0 To n - 1
      Box(i * 10, j * 10, 10, 10, Abs(*dens\f[IX(i, j)]) * 255)
    Next
  Next
  StopDrawing()
  
  WindowEvent()
  FlipBuffers()
EndProcedure

Procedure set_datas(n.l, *dens0.FArray, *u0.FArray, *v0.FArray)
  Define.l i, j
  Define.l z, Size = (n + 2) * (n + 2)
  
  ExamineMouse()
  i = MouseX() / 10
  j = MouseY() / 10
  
  If MouseButton(1)
    fillmemory_(*dens0, (n + 2) * (n + 2) * SizeOf(Float), 0)
    fillmemory_(*u0, (n + 2) * (n + 2) * SizeOf(Float), 0)
    fillmemory_(*v0, (n + 2) * (n + 2) * SizeOf(Float), 0)
    
    *u0\f[IX(i, j)] = 0
    *v0\f[IX(i, j)] = -99
    
    ;Beep_(1000, 50)
    For z = 0 To Size - 1
      *dens0\f[z] = 0.0005
    Next
  EndIf
EndProcedure


Procedure main()
  Define *u.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  Define *v.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  Define *u_prev.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  Define *v_prev.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  Define *dens.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  Define *dens_prev.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  
  Define.f visc = 8.0e-4, dt = 0.5, diff = 0.01
  
  Repeat
    set_datas(#N, *dens_prev, *u_prev, *v_prev)
    
    vel_step(#N, *u, *v, *u_prev, *v_prev, visc, dt)
    dens_step(#N, *dens, *dens_prev, *u, *v, diff, dt)
    
    draw_dens(#N, *dens)
    
    Delay(30)
  ForEver
EndProcedure


InitSprite()
InitMouse()
InitKeyboard()

OpenWindow(0, 100, 100, #N * 10, #N * 10, "Fluid Dynamics")
OpenWindowedScreen(WindowID(0), 0, 0, #N * 10, #N * 10, 0, 0, 0)
main()
Sieht schon nach etwas aus, jedoch wurde noch etwas getrickst. Ich denke,
die Funktion array_normalize() ist unnötig, da sie das Resultat verfälscht.
Ev. hat es also noch Fehler oder die Parameter sind nicht im richtigen
Bereich. Wenns jemand gleich sieht, danke für eine Rückmeldung :)

greetz
Remi

Verfasst: 11.09.2006 16:47
von Kaeru Gaman
yo, ganz witzig.

wenn ich die maus mit gedrücktem knopf etwas schneller bewege,
verschwindet der effekt, als hätte ich losgelassen.

muss mir die mathematik mal reinziehen, wenn ich muße hab,
ist bestimmt für einige dinge gut zu gebrauchen.

Verfasst: 11.09.2006 16:54
von remi_meier
> wenn ich die maus mit gedrücktem knopf etwas schneller bewege,
> verschwindet der effekt, als hätte ich losgelassen.
liegt daran, dass ich ein Delay() zur Verlangsamung drin hab. Mit zu
grossen Werten hat PB mit DirectX dann Probleme.

Verfasst: 12.09.2006 14:44
von remi_meier
Hab ein paar Parameter von PlasmaPong übernommen und nun hat man
Plasma :) . Ein paar Ungereimtheiten gibt es noch: Ich habe die set_bnd()
Prozedur so abgeändert, dass die Ränder alle Wellen verschlucken. Sonst
steigt die Dichte ohne grosses Zutun automatisch, was ich mir noch nicht
erklären kann. Auch ist die Dichte immer noch nicht zwischen 0 und 1,
dieses Problem habe ich allerdings mit der Sigmoiden Funktion umgangen.

Sieht schon schön aus und ist langsam recht brauchbar. Werds ev. mal
ein wenig auf SSE umbauen.

Code: Alles auswählen

; Push right and left mouse button (left: set density, right: set velocity)

EnableExplicit

Structure FArray
  f.f[0]
EndStructure


Macro IX(i, j)
  ((i) + (n + 2) * (j))
EndMacro

Declare set_bnd(n.l, b.l, *x.FArray)
Declare project(n.l, *u.FArray, *v.FArray, *p.FArray, *div.FArray)

Procedure add_source(n.l, *x.FArray, *s.FArray, dt.f)
  Define.l i, Size = (n + 2) * (n + 2)
  
  For i = 0 To Size - 1
    *x\f[i] + dt * *s\f[i]
  Next
EndProcedure


Procedure diffuse(n.l, b.l, *x.FArray, *x0.FArray, diff.f, dt.f)
  Define.l i, j, k
  Define.f a
  a = dt * diff * n * n
  
  For k = 0 To 19
    For i = 1 To n
      For j = 1 To n
        *x\f[IX(i, j)] = (*x0\f[IX(i, j)] + a * (*x\f[IX(i-1, j)] + *x\f[IX(i+1, j)] + *x\f[IX(i, j-1)] + *x\f[IX(i, j+1)])) / (1.0 + 4.0 * a)
        ;*x\f[IX(i, j)] = *x0\f[IX(i, j)] + a * (*x0\f[IX(i-1, j)] + *x0\f[IX(i+1, j)] + *x0\f[IX(i, j-1)] + *x0\f[IX(i, j+1)] - 4 * *x0\f[IX(i, j)])
      Next
    Next
    set_bnd(n, b, *x)
  Next
EndProcedure

Procedure advect (n.l, b.l, *d.FArray, *d0.FArray, *u.FArray, *v.FArray, dt.f)
  Define.l i, j, i0, j0, i1, j1
  Define.f x, y, s0, t0, s1, t1, dt0
  dt0 = dt * n
  
  For i = 1 To n
    For j = 1 To n
      x = i - dt0 * *u\f[IX(i, j)]
      y = j - dt0 * *v\f[IX(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\f[IX(i, j)] = s0 * (t0 * *d0\f[IX(i0, j0)] + t1 * *d0\f[IX(i0, j1)]) + s1 * (t0 * *d0\f[IX(i1, j0)] + t1 * *d0\f[IX(i1, j1)])
      
    Next
  Next
  set_bnd(n, b, *d)
EndProcedure

Procedure dens_step(n.l, *x.FArray, *x0.FArray, *u.FArray, *v.FArray, diff.f, dt.f)
  add_source(n, *x, *x0, dt)
  diffuse(n, 0, *x0, *x, diff, dt )
  advect(n, 0, *x, *x0, *u, *v, dt )
EndProcedure

Procedure vel_step(n.l, *u.FArray, *v.FArray, *u0.FArray, *v0.FArray, visc.f, dt.f)
  add_source(n, *u, *u0, dt) : add_source(n, *v, *v0, dt)
  diffuse(n, 1, *u0, *u, visc, dt)
  diffuse(n, 2, *v0, *v, visc, dt)
  project(n, *u0, *v0, *u, *v)
  
  advect(n, 1, *u, *u0, *u0, *v0, dt) : advect(n, 2, *v, *v0, *u0, *v0, dt)
  project(n, *u, *v, *u0, *v0)
EndProcedure


Procedure project(n.l, *u.FArray, *v.FArray, *p.FArray, *div.FArray)
  Define.l i, j, k
  Define.f h = 1.0 / n
  
  For i = 1 To n
    For j = 1 To n
      *div\f[IX(i, j)] = -0.5 * h * (*u\f[IX(i+1, j)] - *u\f[IX(i-1, j)] + *v\f[IX(i, j+1)] - *v\f[IX(i, j-1)])
      *p\f[IX(i, j)] = 0
    Next
  Next
  set_bnd(n, 0, *div) : set_bnd(n, 0, *p)
  
  For k = 0 To 19
    For i = 1 To n
      For j = 1 To n
        *p\f[IX(i, j)] = (*div\f[IX(i, j)] + *p\f[IX(i-1, j)] + *p\f[IX(i+1, j)] + *p\f[IX(i, j-1)] + *p\f[IX(i, j+1)]) / 4
      Next
    Next
    set_bnd(n, 0, *p)
  Next
  
  For i = 1 To n
    For j = 1 To n
      *u\f[IX(i, j)] - (0.5 * (*p\f[IX(i+1, j)] - *p\f[IX(i-1, j)]) / h)
      *v\f[IX(i, j)] - (0.5 * (*p\f[IX(i, j+1)] - *p\f[IX(i, j-1)]) / h)
    Next
  Next
  set_bnd(n, 1, *u) : set_bnd(n, 2, *v)
EndProcedure


Procedure set_bnd(n.l, b.l, *x.FArray)
  Define.l i
  
  For i = 1 To n
    If b = 1
      *x\f[IX(0, i)]   = 0; -*x\f[IX(1, i)]
      *x\f[IX(n+1, i)] = 0; -*x\f[IX(n, i)]
    Else
      *x\f[IX(0, i)]   = 0; *x\f[IX(1, i)]
      *x\f[IX(n+1, i)] = 0; *x\f[IX(n, i)]
    EndIf
    
    If b = 2
      *x\f[IX(i,   0)] = 0; -*x\f[IX(i, 1)]
      *x\f[IX(i, n+1)] = 0; -*x\f[IX(i, n)]
    Else
      *x\f[IX(i,   0)] = 0; *x\f[IX(i, 1)]
      *x\f[IX(i, n+1)] = 0; *x\f[IX(i, n)]
    EndIf
  Next
  
  *x\f[IX(0,     0)] = 0.5 * (*x\f[IX(1,   0)] + *x\f[IX(0,   1)])
  *x\f[IX(0,   n+1)] = 0.5 * (*x\f[IX(1, n+1)] + *x\f[IX(0,   n)])
  *x\f[IX(n+1,   0)] = 0.5 * (*x\f[IX(n,   0)] + *x\f[IX(n+1, 1)])
  *x\f[IX(n+1, n+1)] = 0.5 * (*x\f[IX(n, n+1)] + *x\f[IX(n+1, n)])
EndProcedure

;- End Simulation
#N = 100
#DSTEP = 3
#SIZE = (#N + 2) * (#N + 2)

Procedure.l color_gradient(StartColor.l, EndColor.l, f.d)
  Protected r.l, g.l, b.l
  Protected rt.l, gt.l, bt.l
  
  rt = Red(StartColor)
  gt = Green(StartColor)
  bt = Blue(StartColor)
  
  r = rt + f * (Red(EndColor) - rt)
  g = gt + f * (Green(EndColor) - gt)
  b = bt + f * (Blue(EndColor) - bt)
  ProcedureReturn = RGB(r, g, b)
EndProcedure

Procedure draw_dens(n.l, *dens.FArray)
  Define.l i, j
  Define.f max
  
  ClearScreen(0)
  ExamineKeyboard()
  If KeyboardPushed(#PB_Key_Escape)
    End
  EndIf
  
  StartDrawing(ScreenOutput())
  For i = 0 To n - 1
    For j = 0 To n - 1; $FF, $FFFF
      Box(i * #DSTEP, j * #DSTEP, #DSTEP, #DSTEP, color_gradient($55, $FFFF, 1.0 / (1 + Pow(2.71828, -5*(*dens\f[IX(i, j)] - 0.5)))))
    Next
  Next
  StopDrawing()
  
  WindowEvent()
  FlipBuffers()
EndProcedure

Procedure set_datas(n.l, *dens0.FArray, *u0.FArray, *v0.FArray)
  Define.l i, j
  Define.l z, Size = (n + 2) * (n + 2)
  
  ExamineMouse()
  i = MouseX() / #DSTEP
  j = MouseY() / #DSTEP
  
  fillmemory_(*dens0, (n + 2) * (n + 2) * SizeOf(Float), 0)
  If MouseButton(1)
    ; fillmemory_(*dens0, (n + 2) * (n + 2) * SizeOf(Float), 0)
    fillmemory_(*u0, (n + 2) * (n + 2) * SizeOf(Float), 0)
    fillmemory_(*v0, (n + 2) * (n + 2) * SizeOf(Float), 0)
    
    For z = -5 To 5
      ; *u0\f[IX(i+z, j)] = 0
      ; *v0\f[IX(i+z, j)] = -500
      *dens0\f[IX(i+z, j)] = 200
    Next
  EndIf  
  If MouseButton(2)
    
    fillmemory_(*u0, (n + 2) * (n + 2) * SizeOf(Float), 0)
    fillmemory_(*v0, (n + 2) * (n + 2) * SizeOf(Float), 0)
    
    For z = -5 To 5
      *u0\f[IX(i+z, j)] = 0
      *v0\f[IX(i+z, j)] = -500
      ; *dens0\f[IX(i+z, j)] = 200
    Next
  EndIf
EndProcedure


Procedure main()
  Define *u.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  Define *v.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  Define *u_prev.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  Define *v_prev.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  Define *dens.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  Define *dens_prev.FArray = AllocateMemory(#SIZE * SizeOf(Float))
  
  Define.f visc = 0.001, dt = 0.015, diff = 0.0018
  
  Repeat
    set_datas(#N, *dens_prev, *u_prev, *v_prev)
    
    vel_step(#N, *u, *v, *u_prev, *v_prev, visc, dt)
    dens_step(#N, *dens, *dens_prev, *u, *v, diff, dt)
    
    draw_dens(#N, *dens)
    
    Delay(0)
  ForEver
EndProcedure


InitSprite()
InitMouse()
InitKeyboard()

OpenWindow(0, 100, 100, #N * #DSTEP, #N * #DSTEP, "Fluid Dynamics")
OpenWindowedScreen(WindowID(0), 0, 0, #N * #DSTEP, #N * #DSTEP, 0, 0, 0)
main()


Verfasst: 12.09.2006 14:58
von Kaeru Gaman
bevor ich mich jetzt durch jede einzelne operation durchbeiße...

hast du nen link zur hand, von ner T&T-site wo die grundzüge der Plasma-Rechnung erklärt werden?

Verfasst: 12.09.2006 15:05
von remi_meier
Nicht direkt, da ich noch nicht versucht habe, den Code richtig zu verstehen
(die Physik und Mathematik). Im PDF wird ganz kurz die Formel angesprochen,
mit den Zusammenhängen zwischen Dichte, Kraft, etc.
Das wäre ev. noch etwas.:
http://www9.informatik.uni-erlangen.de/ ... troemungen

Ansonsten suche ich mal, wenn ich etwas mehr Zeit habe. Dann werde
ich ev. auch eine physikalisch korrekte Version hinbekommen.

Verfasst: 12.01.2007 17:14
von Vallan
Also bei mir Geht der code überhauptnicht.

ich sehe nur eine Braune Fläche :-(

Verfasst: 12.01.2007 17:39
von Andre
Vallan hat geschrieben:Also bei mir Geht der code überhauptnicht.

ich sehe nur eine Braune Fläche :-(
Kann ich so nicht ganz bestätigen: ich sehe auch die braune Fläche, aber nach "Rumdrücken" der linken + rechten Maustaste und Bewegen des (unsichtbaren) Mauscursors sehe ich schon einige Effekte, die dann langsam ausgeblendet werden. Aber alles eben in Zeitlupe...
(Selbst das Beenden mit Esc dauert etliche Sekunden, bis das Programm reagiert.)

Vielleicht hat sich jetzt bzgl. des Codes (der ja sicherlich für PB4.00 geschrieben wurde) jetzt etwas bei PB 4.02 geändert?


//Edit: mit ausgeschaltetem Debugger gehts natürlich schneller... :oops:
Aber eben trotzdem immer noch sehr langsam. /:->

Verfasst: 12.01.2007 18:01
von remi_meier
Jetzt tut doch nicht so, muss ja auch viel berechnet werden :mrgreen:

Verfasst: 12.01.2007 20:35
von Kekskiller
Ich sehe auch nur braun. Nicht mal das "Gedrücke", wie so schön gesagt funzt. Ein Berechnungsbalken wäre für sowas ev. ganz sinnvoll.. nur mal zum Testen, damit man weiß wo hier die Performance hingeht, bzw. wie lange es dauert, wie man was sieht.........