Ytukyg

I am making the process of moving from MatLab to Python, and am attempting a real time N-body simulation as a way of getting to grips with the language. I have completed a script that uses glumpy for the visualisation, and vectorised numpy code for the numerical integration, but it slows waaaaaay down for more than about 150 points. Essentially, I am looking for tips on how to optimize the code so that it can run in real time with at least an order of magnitude more points. Particularly, is there any way of implementing the integration within the shader code using the GPU??

I am aware of the limitations with the Euler method that I have used (and the scaling of the acceleration), and am looking for input on optimization of the software implementation, as opposed to the numerical method itself.

import numpy as np

from glumpy import app, gloo, gl



#VERTEX SHADER CODE

vertex = """

    attribute vec2 position;

    void main() {

         gl_PointSize = 5.0;

         gl_Position = vec4(position, 0.0, 1.0);



    } """



#FRAGMENT SHADER CODE

fragment = """

    void main() {

        gl_FragColor = vec4(vec3(0.5), 1.0);

    } """



#CREATE ARRAY OF POSITIONS

n_points = int(300)                         #number of points

p = np.random.uniform(-1,1,(n_points,2))    #array of positions (x,y)

v = np.zeros((n_points,2))                  #array of velocities (xdot,ydot)

m = np.ones((n_points,1))                   #array of masses

a = np.zeros((n_points,2))                  #array of accelerations

index = np.linspace(0,n_points-1,n_points)



#OPEN WINDOW AND CREATE SHADER PROGRAM

window = app.Window(512,512, color = (1,1,1,1))

points = gloo.Program(vertex,fragment,count = n_points)

points["position"] = p



@window.event

def on_draw(dt):

    global v, p, a, m, theta



    for i in range(n_points):    #for each point

        delta_x = p[index != i, 0] - p[i, 0]    #distance in x to each other point

        delta_y = p[index != i, 1] - p[i, 1]    #distance in y to each other point

        mass = m[index != i,:]    #mass of other points

        theta = np.arctan2(delta_y, delta_x)  #angle to other points

        r = (delta_x**2 + delta_y**2)*0.5 + 0.01  #distance to other points

        a[i,0] = np.sum(mass * np.cos(theta) / r**2)  #acceleration in x

        a[i,1] = np.sum(mass * np.sin(theta) / r**2)  #acceleration in y



    print(dt)

    a = a*0.00000001

    v = v + a*dt

    p = p + v*dt



    window.clear()

    points["position"] = p

    points.draw(gl.GL_POINTS)



app.run()