VB.NET: Problem converting fluid sim code from c++

[Mantu]

I found nice realtime fluid simulation code created with c++. I converted it to vb.net but i have run into trouble. You can find original c++ code here: http://www.dgp.toronto.edu/people/st...DROM_GDC03.zip . Here's my vb.net solver part of code:

    Sub set_bnd(ByVal N As Integer, ByVal b As Integer, ByVal x(,) As Single)
 
        Dim i As Integer
 
        For i = 1 To (N + 1)
 
            If b = 1 Then
                x(0, i) = -x(1, i)
            Else
                x(0, i) = x(1, i)
            End If
 
            If b = 1 Then
                x(N + 1, i) = -x(N, i)
            Else
                x(N + 1, i) = x(N, i)
            End If
 
            If b = 2 Then
                x(i, 0) = -x(i, 1)
            Else
                x(i, 0) = x(i, 1)
            End If
 
            If b = 2 Then
                x(i, N + 1) = -x(i, N)
            Else
                x(i, N + 1) = x(i, N)
            End If
        Next
 
        x(0, 0) = 0.5F * (x(1, 0) + x(0, 1))
        x(0, N + 1) = 0.5F * (x(1, N + 1) + x(0, N))
        x(N + 1, 0) = 0.5F * (x(N, 0) + x(N + 1, 1))
        x(N + 1, N + 1) = 0.5F * (x(N, N + 1) + x(N + 1, N))
    End Sub
 
    Sub lin_solve(ByVal N As Integer, ByVal b As Integer, ByVal x(,) As Single, ByVal x0(,) As Single, ByVal a As Single, ByVal c As Single)
 
        Dim i, j, k As Integer
 
        For k = 0 To 20
            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
            set_bnd(N, b, x)
        Next
    End Sub
 
    Sub add_source(ByVal N As Integer, ByVal x(,) As Single, ByVal s(,) As Single, ByVal dt As Single)
        Dim i, j As Integer
        Dim size As Integer = (N + 2)
 
        For i = 0 To size
            For j = 0 To size
                x(i, j) += dt * s(i, j)
            Next
        Next
    End Sub
 
    'Diffusion
    Sub diffuse(ByVal N As Integer, ByVal b As Integer, ByVal x(,) As Single, ByVal x0(,) As Single, ByVal diff As Single, ByVal dt As Single)
        Dim a As Single
        a = dt * diff * N * N
 
        lin_solve(N, b, x, x0, a, 1 + 4 * a)
    End Sub
 
    'Advection
    Sub advect(ByVal N As Integer, ByVal b As Integer, ByVal d(,) As Single, ByVal d0(,) As Single, ByVal u(,) As Single, ByVal v(,) As Single, ByVal dt As Single)
        Dim i, j, i0, j0, i1, j1 As Integer
        Dim x, y, s0, t0, s1, t1, dt0 As Single
 
        dt0 = dt * N
 
        For i = 1 To N
            For j = 1 To N
 
                x = i - dt0 * u(i, j)
                y = j - dt0 * v(i, j)
 
                If x < 0.5F Then x = 0.5F
                If x > (N + 0.5F) Then x = N + 0.5F
 
                i0 = CInt(x)
                i1 = i0 + 1
 
                If y < 0.5F Then y = 0.5F
                If y > (N + 0.5F) Then y = N + 0.5F
 
                j0 = CInt(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
 
        set_bnd(N, b, d)
    End Sub
 
    Sub project(ByVal N As Integer, ByVal u(,) As Single, ByVal v(,) As Single, ByVal p(,) As Single, ByVal div(,) As Single)
        Dim i, j As Integer
        Dim h As Single
 
        h = 1.0F / N
 
        For i = 1 To N + 1
            For j = 1 To N + 1
                div(i, j) = -0.5F * h * (u(i + 1, j) - u(i - 1, j) + v(i, j + 1) - v(i, j - 1))
                p(i, j) = 0
            Next
        Next
 
        set_bnd(N, 0, div)
        set_bnd(N, 0, p)
 
        lin_solve(N, 0, p, div, 1, 4)
 
        For i = 1 To N + 1
            For j = 1 To N + 1
                u(i, j) -= 0.5F * (p(i + 1, j) - p(i - 1, j)) / h
                v(i, j) -= 0.5F * (p(i, j + 1) - p(i, j - 1)) / h
            Next
        Next
 
        set_bnd(N, 1, u)
        set_bnd(N, 2, v)
    End Sub
 
    Sub dens_step(ByVal N As Integer, ByVal x(,) As Single, ByVal x0(,) As Single, ByVal u(,) As Single, ByVal v(,) As Single, ByVal diff As Single, ByVal dt As Single)
 
        Array.Clear(temp, 0, Size1D)
 
        add_source(N, x, x0, dt)
 
        'Swap arrays
        Array.Copy(x, temp, Size1D)
        Array.Copy(x, x0, Size1D)
        Array.Copy(temp, x, Size1D)
 
        diffuse(N, 0, x, x0, diff, dt)
 
        Array.Clear(temp, 0, Size1D)
 
        'Swap arrays
        Array.Copy(x, temp, Size1D)
        Array.Copy(x, x0, Size1D)
        Array.Copy(temp, x, Size1D)
 
        advect(N, 0, x, x0, u, v, dt)
    End Sub
 
    Sub vel_step(ByVal N As Integer, ByVal u(,) As Single, ByVal v(,) As Single, ByVal u0(,) As Single, ByVal v0(,) As Single, ByVal visc As Single, ByVal dt As Single)
 
        add_source(N, u, u0, dt)
        add_source(N, v, v0, dt)
 
        Array.Clear(temp, 0, Size1D)
 
        Array.Copy(u, temp, Size1D)
        Array.Copy(u, u0, Size1D)
        Array.Copy(temp, u, Size1D)
 
        diffuse(N, 1, u, u0, visc, dt)
 
        Array.Clear(temp, 0, Size1D)
 
        Array.Copy(v, temp, Size1D)
        Array.Copy(v, v0, Size1D)
        Array.Copy(temp, v, Size1D)
 
        diffuse(N, 2, v, v0, visc, dt)
 
        project(N, u, v, u0, v0)
 
        Array.Clear(temp, 0, Size1D)
 
        Array.Copy(u, temp, Size1D)
        Array.Copy(u, u0, Size1D)
        Array.Copy(temp, u, Size1D)
 
        Array.Clear(temp, 0, Size1D)
 
        Array.Copy(v, temp, Size1D)
        Array.Copy(v, v0, Size1D)
        Array.Copy(temp, v, Size1D)
 
        advect(N, 1, u, u0, u0, v0, dt)
        advect(N, 2, v, v0, u0, v0, dt)
 
        project(N, u, v, u0, v0)
    End Sub
End Class
 

Simulation seems to become very unstable with vb.net. I think problem is somewhere in vel_step, which calculates velocity components. Because when displaying velocity vectorfield its very unstable and flickering. When I tried to debug code I saw that some of the values become -1.#IND which means I'm dividing with zero or doing something else wrong. But I havent found any part where division by zero might happen and I think this is not the case.

I have checked and rechecked my code. Calculations are identical to c++. Maybe there's something wrong with .net and it cant perform this amount of calculations. :D

[Mantu]

Here's the rest of code if you want to compile it:

Imports Tao.OpenGl
Imports Tao.FreeGlut
 
Public Class NavierStokes
 
    'Resolution params
    Dim N, size, Size1D As Integer
    'Simulation params
    Dim dt, diff, visc, force, source As Single
 
    Dim win_x, win_y As Integer
 
    Dim dvel As Boolean
 
    Dim omx, omy, mouX, mouY As Integer
    Dim mouse_down() As Boolean = {False, False, False}
 
    'Using 2d arrays to store data instead of 1d like in C++ code.
    Dim u(,), u_prev(,) As Single
    Dim v(,), v_prev(,) As Single
    Dim dens(,), dens_prev(,) As Single
 
    Dim temp(,) As Single
 
    Public Sub New()
        'Give some values to params
        N = 50
        size = (N + 2)
        Size1D = size * size
        dt = 0.1F
        visc = 0.0008F 'Here's something strange c++ uses zeros here,
        diff = 0.0009F 'if you put zeros now simulation behaves strange.
        force = 5.0F
        source = 100.0F
 
        win_x = 512
        win_y = 512
 
        omx = 0
        omy = 0
        mouX = 0
        mouY = 0
 
        dvel = False
 
        'Defines how much space arrays need.
        ReDim u(size, size)
        ReDim u_prev(size, size)
        ReDim v(size, size)
        ReDim v_prev(size, size)
        ReDim dens(size, size)
        ReDim dens_prev(size, size)
        ReDim temp(size, size)
    End Sub
 
    Sub clear_data()
        Array.Clear(u, 0, Size1D)
        Array.Clear(u_prev, 0, Size1D)
        Array.Clear(v, 0, Size1D)
        Array.Clear(v_prev, 0, Size1D)
        Array.Clear(dens, 0, Size1D)
        Array.Clear(dens_prev, 0, Size1D)
    End Sub
 
    Sub pre_display()
        Gl.glViewport(0, 0, win_x, win_y)
        Gl.glMatrixMode(Gl.GL_PROJECTION)
        Gl.glLoadIdentity()
        Gl.glOrtho(0, 1, 0, 1, 0, 100)
        Gl.glClearColor(0.0F, 0.0F, 0.0F, 1.0F)
        Gl.glClear(Gl.GL_COLOR_BUFFER_BIT)
    End Sub
 
    Sub post_display()
        Glut.glutSwapBuffers()
    End Sub
 
    Sub draw_velocity()
        Dim i, j As Integer
        Dim x, y, h As Single
 
        h = 1.0F / N
 
        Gl.glColor3f(1.0F, 0.0F, 0.0F)
        Gl.glLineWidth(1.0F)
 
        Gl.glBegin(Gl.GL_LINES)
        For i = 1 To (N + 1)
 
            x = (i - 0.5F) * h
            For j = 1 To (N + 1)
 
                y = (j - 0.5F) * h
                Gl.glVertex2f(x, y)
                Gl.glVertex2f(x + u(i, j), y + v(i, j))
            Next
        Next
        Gl.glEnd()
    End Sub
 
    Sub draw_density()
        Dim i, j As Integer
        Dim x, y, h, d00, d01, d10, d11 As Single
 
        h = 1.0F / N
 
        Gl.glBegin(Gl.GL_QUADS)
        For i = 0 To (N + 1)
 
            x = (i - 0.5F) * h
            For j = 0 To (N + 1)
 
                y = (j - 0.5F) * h
                d00 = dens(i, j)
                d01 = dens(i, j + 1)
                d10 = dens(i + 1, j)
                d11 = dens(i + 1, j + 1)
 
                Gl.glColor3f(d00, d00, d00)
                Gl.glVertex2f(x, y)
                Gl.glColor3f(d10, d10, d10)
                Gl.glVertex2f(x + h, y)
                Gl.glColor3f(d11, d11, d11)
                Gl.glVertex2f(x + h, y + h)
                Gl.glColor3f(d01, d01, d01)
                Gl.glVertex2f(x, y + h)
            Next
        Next
        Gl.glEnd()
    End Sub
 
    Sub get_from_UI(ByVal d(,) As Single, ByVal u(,) As Single, ByVal v(,) As Single)
        Dim i, j As Integer
 
        Array.Clear(d, 0, Size1D)
        Array.Clear(u, 0, Size1D)
        Array.Clear(v, 0, Size1D)
 
        If Not mouse_down(Glut.GLUT_LEFT_BUTTON) And Not mouse_down(Glut.GLUT_RIGHT_BUTTON) Then
            Return
        End If
 
        i = CInt((mouX / CSng(win_x)) * N + 1)
        j = CInt(((win_y - CSng(mouY)) / CSng(win_y)) * CSng(N) + 1.0)
 
        If i < 1 Or i > N Or j < 1 Or j > N Then
            Return
        End If
 
        If mouse_down(Glut.GLUT_LEFT_BUTTON) Then
 
            u(i, j) = force * (mouX - omx)
            v(i, j) = force * (omy - mouY)
        End If
 
        If mouse_down(Glut.GLUT_RIGHT_BUTTON) Then
 
            d(i, j) = source
        End If
 
        omx = mouX
        omy = mouY
    End Sub
 
    Sub key_func(ByVal key As Byte, ByVal x As Integer, ByVal y As Integer)
 
        If key = 99 Then
            clear_data()
        End If
 
        If key = 118 Then
            dvel = Not dvel
        End If
    End Sub
 
    Sub mouse_func(ByVal button As Integer, ByVal state As Integer, ByVal x As Integer, ByVal y As Integer)
 
        mouX = x
        mouY = y
        omx = mouX
        omy = mouY
 
        mouse_down(button) = (state = Glut.GLUT_DOWN)
    End Sub
 
    Sub motion_func(ByVal x As Integer, ByVal y As Integer)
 
        mouX = x
        mouY = y
    End Sub
 
    Sub reshape_func(ByVal width As Integer, ByVal height As Integer)
 
        Glut.glutReshapeWindow(width, height)
        win_x = width
        win_y = height
    End Sub
 
    Sub idle_func()
 
        get_from_UI(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)
 
        Glut.glutPostRedisplay()
    End Sub
 
    Sub display_func()
 
        pre_display()
 
        If dvel Then
            draw_velocity()
        Else
            draw_density()
        End If
 
        post_display()
    End Sub
 
    Sub open_glut_window()
 
        Glut.glutInitDisplayMode(Glut.GLUT_RGBA Or Glut.GLUT_DOUBLE)
        Glut.glutInitWindowPosition(0, 0)
        Glut.glutInitWindowSize(win_x, win_y)
        Glut.glutCreateWindow("NavierStokes")
        Gl.glClearColor(0.0F, 0.0F, 0.0F, 1.0F)
        Gl.glClear(Gl.GL_COLOR_BUFFER_BIT)
        Glut.glutSwapBuffers()
        Gl.glClear(Gl.GL_COLOR_BUFFER_BIT)
        Glut.glutSwapBuffers()
 
        pre_display()
 
        Glut.glutKeyboardFunc(New Glut.KeyboardCallback(AddressOf key_func))
        Glut.glutMouseFunc(New Glut.MouseCallback(AddressOf mouse_func))
        Glut.glutMotionFunc(New Glut.MotionCallback(AddressOf motion_func))
        Glut.glutReshapeFunc(New Glut.ReshapeCallback(AddressOf reshape_func))
        Glut.glutIdleFunc(New Glut.IdleCallback(AddressOf idle_func))
        Glut.glutDisplayFunc(New Glut.DisplayCallback(AddressOf display_func))
    End Sub
 
    Public Shared Sub Main()
 
        Dim NavierStokes As Navier_Stokes2.NavierStokes
        NavierStokes = New Navier_Stokes2.NavierStokes
 
        Glut.glutInit()
        NavierStokes.clear_data()
        NavierStokes.open_glut_window()
        Glut.glutMainLoop()
    End Sub
 

I'm using Tao Framework for opengl and glut. You can get it here:
http://taoframework.com/

[Plater]

Hmm, I haven't looked close but it sounds like you are either overflowing your bounds, or something that you thought would be kept as a floating point is getting converted to integer. I saw singles and integers used but no real typecasting. So it is very likely you have fallen victem to it.
For example:

float f=0.0;
int a=10;
int b=3;
 
f=a/b;
 

One might assume that f=3.333 (to some percision), however, it will just be 3, you will need to typecast things out to retain the fraction values

[Mantu]

Well I checked every part where typecasting might happen. I used Cint and Csng where were typeconversions , but it didn't work. I don't know about overflowing could you tell where this could happen? Or does anyone else have other suggestions?

[Plater]

Anywhere you do math with integers.
For example

Dim h As Single
 
h = 1.0F / N
 

Here you have a floating point number 1.0F being divided by an integer, the result will be the interger portion. So if N=2 then h=0, if N=1 then f=1 and the like.
You will need to fix THOSE kind of occurances

[Mantu]

How you would fix this I tried to but it's still the same.

[Plater]

Might be overkill but should work:

Dim h As Single
 
h = (Single)(1.0F / (Single)N)
 

[Mantu]

Do you mean:

Dim h As Single
 
h = CSng(1.0F / CSng(N))
 

I tried to change every thing according to this but doesn't work, this is really driving me nuts.

[Mantu]

I found something:
http://techlists.org/archives/progra...msg04100.shtml
The same code is converted to python but why vb.net can't do it?

[Plater]

Try c#?
I find it's much easier to understand.