My first LBM

Lattice Boltzmann Method LBM: Theory
#1

Hello,

i’m trying to implement LBM on Visual Basic 6, visualizing the velocitiy vectors on the lattice points with directx8. It seems to work somehow, but I think there’s still something wrong in the main algorithm. Tho following code shows the streaming and collision steps, could someone overlook it and tell me if it’s wrong or right?

Commenting the code I found that my Bounce-Back-Alg. is nonsense, so I just deleted it for now…



Sub streaming()                          'STREAMING

    For i = 0 To N                           'POSITION X
        For j = 0 To N                       'POSITION Y
            For k = 0 To 8                   'DIRECTIONS 
                nextposx = i + Cx(k)         'CALCULATING THE NEXT POSITION WITH THE DIRECTION 
                nextposy = j + Cy(k)         'VECTOR  ( Cx(k), Cy(k) )
                If Check6.Value = 1 Then     'PERIODIC BOUNDARIES
                    If nextposx = -1 Then    'NEXT POSITION SHOULD BE BETWEEN 0 AND N
                        nextposx = N
                    ElseIf nextposx = N + 1 Then
                        nextposx = 0
                    End If
                    If nextposy = -1 Then
                        nextposy = N
                    ElseIf nextposy = N + 1 Then
                        nextposy = 0
                    End If
                    
                    fnew(nextposx, nextposy, k) = f(i, j, k)    'THE ACTUAL STREAMING

                Else    'NONPERIODIC BOUNDARIES, IGNORE POINTS OUTSIDE LATTICE
                    If nextposx >= 0 And nextposx <= N And nextposy >= 0 And nextposy <= N Then
                        fnew(nextposx, nextposy, k) = f(i, j, k)  'THE ACTUAL STREAMING
                    End If
                End If
            Next k
        Next j
    Next i
  
    
    For i = 0 To N                            
        For j = 0 To N
            For k = 0 To 8
                f(i, j, k) = fnew(i, j, k)   'REASSIGN THE NEW DISTRIBUTION
            Next k
        Next j
    Next i

End Sub

Sub collision()  'COLLISION

    For i = 0 To N   'POSITION X
        For j = 0 To N   'POSITION Y
            
           
            fdensity(i, j) = 0   'DENSITY CALCULATION
            For l = 0 To 8
                fdensity(i, j) = fdensity(i, j) + f(i, j, l)
            Next l
            
           
            fxvelocity(i, j) = 0   'VELOCITY CALCULATION
            fyvelocity(i, j) = 0
            For l = 0 To 8
                fxvelocity(i, j) = fxvelocity(i, j) + f(i, j, l) * Cx(l)
                fyvelocity(i, j) = fyvelocity(i, j) + f(i, j, l) * Cy(l)
            Next l
            If fdensity(i, j) <> 0 Then   'PREVENT DIV. BY ZERO
                fxvelocity(i, j) = fxvelocity(i, j) / fdensity(i, j)
                fyvelocity(i, j) = fyvelocity(i, j) / fdensity(i, j)
            Else
                fxvelocity(i, j) = 0
                fyvelocity(i, j) = 0
            End If
                                    
            For k = 0 To 8 
                                  'CALCULATE EQUILIBRIUM
                fequilibrium = w(k) * fdensity(i, j) * (1 + 3 * (Cx(k) * fxvelocity(i, j) + Cy(k) * fyvelocity(i, j)) / (c * c) + 9 / 2 * (Cx(k) * fxvelocity(i, j) + Cy(k) * fyvelocity(i, j)) / (c * c * c * c) - 3 / 2 * (fxvelocity(i, j) * fxvelocity(i, j) + fyvelocity(i, j) * fyvelocity(i, j)) / (c * c))

                                     'CALCULATE NEW DISTRIBUTION
                fnew(i, j, k) = f(i, j, k) - 1 / tau * (f(i, j, k) - fequilibrium)
            Next k
        Next j
    Next i


    For i = 0 To N
        For j = 0 To N
            For k = 0 To 8
                f(i, j, k) = fnew(i, j, k)  'REASSIGN THE NEW DISTRIBUTION
            Next k
        Next j
    Next i

End Sub

Thank you for reading until here,
Greets,
streikbrecher

#2

Hi streikbrecher,

Congratulations on your first LB code! I’ve personally never seen a LB code written in Visual Basic before and found this pretty interesting. Have you by any chance had the opportunity to measure the execution speed of the code, and scale it against other languages, especially interpreted ones?

It is not easy to find errors in an incomplete code, because of the impossibility to run it and try it out. All in all, I have however the feeling that it does the right thing, except for a minor error in the computation of the equilibrium. The equilibrium has four components, which are all summed up and multiplied by rhow:
(a) the constant term “1”
(b) the term "3 cu"
© the term “9/2 (c*u)^2”
(d) the term “-3/2 u^2”

I do suspect that there is a mistake in term © as, according to my understanding, you failed evaluating the power-two (you compute cu instead of (cu)^2).

I hope this helps. Actually, once your code runs, and if you are interested in sharing your efforts with the community, I would certainly suggest that you consider submitting it for addition in the list of codes.