more efficient wind tunnel simulation in Pygame, using numpy

Question

I am an aerospace student working on a school project for our python programming course. The assignment is create a program only using Pygame and numpy. I decided to create a wind tunnel simulation that simulates the airflow over a two dimensional wing. I was wondering if there is a more efficient way of doing the computation from a programming perspective. I will explain the program:

I have attached an image here:

The (steady) flow field is modeled using the vortex panel method. Basically, I am using a grid of Nx times Ny points where at each point a velocity (u,v) vector is given. Then using Pygame I map these grid points as circles, so that they resemble an area of influence. The grid points are the grey circles in the following image:

I create N particles and determine their velocities by iterating as follows:

create a list of particles.
create a grid list.

for each gridpoint in grid list:
  for each particle in list of particles:
  if particle A is within the area of influence of grid point n (xn,yn):
   particle A its velocity = velocity at grid point n.

Visualize everything in Pygame.

this basic way was the only way I could think of visualizing the flow in Pygame. The simulation works pretty well, but If I increase the number of grid points(increase the accuracy of the flow field), the performance decreases. My question is if there is a more efficient way to do this just using pygame and numpy?

I have attached the code here:

import pygame,random,sys,numpy
from Flow import Compute
from pygame.locals import *
import random, math, sys
#from PIL import Image
pygame.init()

Surface = pygame.display.set_mode((1000,600))

#read the airfoil geometry from a dat file
with open ('./resources/naca0012.dat') as file_name:
    x, y = numpy.loadtxt(file_name, dtype=float, delimiter='\t', unpack=True)    

#parameters used to describe the flow

Nx=30# 30 column grid
Ny=10#10 row grid
N=20#number of panels
alpha=0#angle of attack
u_inf=1#freestream velocity

#compute the flow field
u,v,X,Y= Compute(x,y,N,alpha,u_inf,Nx,Ny)  

#The lists used for iteration
Circles = []
Particles= []
Velocities=[]

#Scaling factors used to properly map the potential flow datapoints into Pygame
magnitude=400
vmag=30
umag=30
panel_x= numpy.multiply(x,magnitude)+315
panel_y= numpy.multiply(-y,magnitude)+308


#build the grid suited for Pygame
grid_x= numpy.multiply(X,magnitude)+300
grid_y= numpy.multiply(Y,-1*magnitude)+300

grid_u =numpy.multiply(u,umag)
grid_v =numpy.multiply(v,-vmag)
panelcoordinates=  zip(panel_x, panel_y)

# a grid area
class Circle:
    def __init__(self,xpos,ypos,vx,vy):

        self.radius=16

        self.x = xpos
        self.y = ypos
        self.speedx = 0
        self.speedy = 0

#create the grid list
for i in range(Ny):
    for s in range(Nx):
        Circles.append(Circle(int(grid_x[i][s]),int(grid_y[i][s]),grid_u[i][s],grid_v[i][s]))
        Velocities.append((grid_u[i][s],grid_v[i][s]))

#a particle
class Particle:
    def __init__(self,xpos,ypos,vx,vy):
        self.image = pygame.Surface([10, 10])
        self.image.fill((150,0,0))
        self.rect = self.image.get_rect()
        self.width=4
        self.height=4
        self.radius =2
        self.x = xpos
        self.y = ypos
        self.speedx = 30
        self.speedy = 0

#change particle velocity if collision with grid point
def CircleCollide(Circle,Particle):
    Particle.speedx = int(Velocities[Circles.index((Circle))][0])
    Particle.speedy = int(Velocities[Circles.index((Circle))][1])

#movement of particles
def Move():
    for Particle in Particles:
        Particle.x += Particle.speedx
        Particle.y += Particle.speedy
#create particle streak 
def Spawn(number_of_particles):
    for i in range(number_of_particles):
            i=i*(300/number_of_particles)        
            Particles.append(Particle(0, 160+i,1,0))

#create particles again if particles are out of wake
def Respawn(number_of_particles):
    for Particle in Particles:
        if Particle.x >1100:
            Particles.remove(Particle)
    if Particles==[]:
            Spawn(number_of_particles)

#Collsion detection using pythagoras and distance formula

def CollisionDetect():

   for Circle in Circles:  
       for Particle in Particles:
           if Particle.y >430 or Particle.y<160:
               Particles.remove(Particle)
           if math.sqrt( ((Circle.x-Particle.x)**2)  +  ((Circle.y-Particle.y)**2)  ) <= (Circle.radius+Particle.radius):       
               CircleCollide(Circle,Particle)

#draw everything
def Draw():
    Surface.fill((255,255,255))
    #Surface.blit(bg,(-300,-83))
    for Circle in Circles:
        pygame.draw.circle(Surface,(245,245,245),(Circle.x,Circle.y),Circle.radius)

    for Particle in Particles:
        pygame.draw.rect(Surface,(150,0,0),(Particle.x,Particle.y,Particle.width,Particle.height),0)


        #pygame.draw.rect(Surface,(245,245,245),(Circle.x,Circle.y,1,16),0)

    for i in range(len(panelcoordinates)-1):
        pygame.draw.line(Surface,(0,0,0),panelcoordinates[i],panelcoordinates[i+1],3)

    pygame.display.flip()


def GetInput():
    keystate = pygame.key.get_pressed()
    for event in pygame.event.get():
        if event.type == QUIT or keystate[K_ESCAPE]:
            pygame.quit();sys.exit()


def main():

    #bg = pygame.image.load("pressure.png")
    #bg = pygame.transform.scale(bg,(1600,800))
    #thesize= bg.get_rect()
    #bg= bg.convert()
    number_of_particles=10
    Spawn(number_of_particles)
    clock = pygame.time.Clock()

    while True:
        ticks = clock.tick(60)
        GetInput()
        CollisionDetect()
        Move()
        Respawn(number_of_particles)
        Draw()
if __name__ == '__main__': main()

The code requires another script that computes the flow field itself. It also reads datapoints from a textfile to get the geometry of the wing. I have not provided these two files, but I can add them if necessary. Thank you in advance.

Answer 1

One bottleneck in your code is likely collision detection. CollisionDetect() computes the distance between each particle and each circle. Then, if a collision is detected, CircleCollide() finds the index of the circle in Circles (a linear search), so that the velocities can be retrieved from the same index in Velocities. Clearly this is ripe for improvement.

First, the Circle class already has the velocities in the speedx/speedy attributes, so Velocities can be eliminated .

Second, because the circles are at fixed locations, you can calculate which circle is closest to any given particle from the position of the particle.

# You may already have these values from creating grid_x etc.
# if not, you only need to calculated them once, because the
# circles don't move
circle_spacing_x = Circles[1].x - Circles[0].x
circle_spacing_y = Circles[Nx].y - Circles[0].y

circle_first_x = Circles[0].x - circle_spacing_x / 2
circle_first_y = Circles[0].y - circle_spacing_y / 2

Then CollisionDetect() becomes:

def CollisionDetect():

    for particle in Particles:
        if particle.y >430 or particle.y<160:
           Particles.remove(particle)
           continue

        c = (particle.x - circle_first_x) // circle_spacing_x
        r = (particle.y - circle_first_y) // circle_spacing_y
        circle = Circles[r*Nx + c]

        if ((circle.x - particle.x)**2 + (circle.y - particle.y)**2
            <= (circle.radius+particle.radius)**2):       
                particle.speedx = int(circle.speedx)
                particle.speedy = int(circle.speedy)