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
[code=vb]
Dim h As Single

h = 1.0F / N
[/code]
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:
[code=vb]
Dim h As Single

h = (Single)(1.0F / (Single)N)
[/code]

[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/programming/pythonlist/2005-02/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.