python - 使用 numpy 在 Pygame 中进行更高效的风洞模拟

Question

我是一名航空学生，正在为我们的 Python 编程类(class)开展一个学校项目。作业是创建一个仅使用 Pygame 和 numpy 的程序。我决定创建一个风洞模拟来模拟二维机翼上的气流。我想知道从编程的角度来看是否有更有效的计算方法。我将解释该程序:

我在这里附上一张图片:

(稳定)流场是使用涡流板方法建模的。基本上，我使用的是 Nx 乘以 Ny 点的网格，其中每个点都有一个速度 (u,v) 向量。然后使用 Pygame，我将这些网格点映射为圆圈，这样它们就像一个影响区域。网格点是下图中的灰色圆圈:

我创建了 N 个粒子并通过如下迭代确定它们的速度:

创建一个粒子列表。
创建一个网格列表。

对于网格列表中的每个网格点:
对于粒子列表中的每个粒子:
若粒子A在网格点n(xn,yn)的影响范围内:
粒子A其速度=网格点n处的速度。

在 Pygame 中可视化所有内容。

这种基本方法是我能想到的在 Pygame 中可视化流程的唯一方法。模拟效果很好，但是如果我增加网格点的数量(增加流场的精度)，性能就会下降。我的问题是是否有更有效的方法仅使用 pygame 和 numpy 来做到这一点？

我在这里附上了代码:

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()

该代码需要另一个脚本来计算流场本身。它还从文本文件中读取数据点以获得机翼的几何形状。 我没有提供这两个文件，但如果需要我可以添加它们。提前谢谢你。

Answer 1

代码中的一个瓶颈可能是碰撞检测。 CollisionDetect() 计算每个粒子和每个圆之间的距离。然后，如果检测到碰撞，CircleCollide() 会在 Circles 中找到圆的索引(线性搜索)，以便可以从同一索引中检索速度在 Velocities 中。显然，改进的时机已经成熟。

首先，Circle 类在 speedx/speedy 属性中已经有了速度，所以可以去掉 Velocities。

其次，因为圆圈位于固定位置，您可以从粒子的位置计算出哪个圆圈最接近任何给定的粒子。

# 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

然后 CollisionDetect() 变成:

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)