Pythonのより効率的な風洞シミュレーション、numpyを使用

Question

私はPythonプログラミングコースの学校プロジェクトに取り組んでいる航空宇宙の学生です。割り当ては、Pygameとnumpyを使用しているプログラムを作成することです。私は、2次元翼上の気流をシミュレートする風洞シミュレーションを作成することに決めました。私はプログラミングの観点から計算を行うより効率的な方法があるのだろうかと思っていました。

私はここに画像を添付しています：

（定常）流れ場 は渦パネルのメソッドを使用してモデル化され、私はプログラムを説明します。基本的には、各点で速度（u、v）ベクトルが与えられるNx倍のNy点のグリッドを使用しています。次に、Pygameを使用して、これらのグリッドポイントを円としてマッピングします。その結果、これらのグリッドポイントは影響範囲に似ています。

Iは、N個の粒子を作成し、以下のように反復することによって、それらの速度を決定する：グリッド点は、次の画像における灰色の円である

粒子のリストを作成します。
グリッドリストを作成します。粒子のリスト中の各粒子について
 ：グリッドリスト内の各格子点のため

   粒子Aは、格子点の影響の領域内にあるN（XN、YN）場合：
     粒子その速度=格子点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() 

コードが流れ場自体を計算し、別のスクリプトが必要です。また、テキストファイルからデータポイントを読み取り、ウィングのジオメトリを取得します。 私はこれら2つのファイルを提供していませんが、必要に応じて追加することができます。前もって感謝します。

Answer 1

コード内のボトルネックは、衝突の可能性があります。 CollisionDetect()は、各粒子と各円の間の距離を計算します。次に、衝突が検出された場合、は円のインデックスをCircles（線形探索）と見なし、Velocitiesの同じインデックスから速度を取り出すことができます。明らかに、これは改善のために熟しています。

まず、Circleクラスのspeedx/speedy属性の速度が既にあるため、Velocitiesを削除することができます。

第2に、円が固定された位置にあるので、どの円が粒子の位置から任意の粒子に最も近いかを計算することができます。

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

Answer 2

私はあなたのコードを整理して、いくつかの変更を加えました。つまり、クラスにスコープを追加し、さらにいくつか紹介します。 Flowの詳細な知識がなければ、私はこれを完全にテストすることはできませんが、私に戻ることができれば、もう少しやり直すことができます。ここでは、「流れ場」をnumpy.meshgrid関数でシミュレートできると仮定しています。

import pygame,numpy,sys 
import pygame.locals 
import math 

class Particle: 
    def __init__(self,xpos,ypos,vx,vy): 
     self.size = numpy.array([4,4]) 
     self.radius =2 
     self.pos = numpy.array([xpos,ypos]) 
     self.speed = numpy.array([30,0]) 
     self.rectangle = numpy.hstack((self.pos,self.size)) 
    def move(self): 
     self.pos += self.speed 
     self.rectangle = numpy.hstack((self.pos,self.size)) 
    def distance(self,circle1): 
     return math.sqrt(numpy.sum((circle1.pos - self.pos)**2)) 
    def collision(self,circle1): 
     result = False 
     if self.pos[1] >430 or self.pos[1]<160: 
      result = True 
     if self.distance(circle1) <= (circle1.radius+self.radius): 
      self.speed = circle1.speed 
     return result 

class Particles: 
    def __init__(self,num_particles): 
     self.num = num_particles 
     self.particles =[] 
     self.spawn() 
    def spawn(self): 
     for i in range(self.num): 
      i=i*(300/self.num)   
      self.particles.append(Particle(0, 160+i,1,0)) 
    def move(self): 
     for particle in self.particles: 
      particle.move() 
      if particle.pos[0] >1100: 
       self.particles.remove(particle) 
     if not self.particles: self.spawn() 
    def draw(self): 
     for particle in self.particles: 
      pygame.draw.rect(Surface,(150,0,0),particle.rectangle,0) 
    def collisiondetect(self,circle1): 
     for particle in self.particles: 
      if particle.collision(circle1): 
       self.particles.remove(particle) 

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

#draw everything 
def Draw(sw,cir): 
    Surface.fill((255,255,255)) 
    cir.draw() 
    for i in range(panelcoordinates.shape[1]): 
     pygame.draw.line(Surface,(0,0,0),panelcoordinates[0,i-1],panelcoordinates[0,i],3) 
    sw.draw() 
    pygame.display.flip() 

# a grid area 
class Circle: 
    def __init__(self,xpos,ypos,vx,vy): 
     self.radius=16 
     self.pos = numpy.array([xpos,ypos]) 
     self.speed = numpy.array([vx,vy]) 

class Circles: 
    def __init__(self,columns,rows): 
     self.list = [] 
     grid_x,grid_y = numpy.meshgrid(numpy.linspace(0,1000,columns),numpy.linspace(200,400,rows)) 
     grid_u,grid_v = numpy.meshgrid(numpy.linspace(20,20,columns),numpy.linspace(-1,1,rows)) 
     for y in range(rows): 
      for x in range(columns): 
       c1= Circle(int(grid_x[y,x]),int(grid_y[y,x]),grid_u[y,x],grid_v[y,x]) 
       self.list.append(c1) 
    def draw(self): 
     for circle in self.list: 
      pygame.draw.circle(Surface,(245,245,245),circle.pos,circle.radius) 
    def detectcollision(self,parts): 
     for circle in self.list: 
      parts.collisiondetect(circle) 


if __name__ == '__main__': 
    #initialise variables 
    number_of_particles=10 
    Nx=30 
    Ny=10 

    #circles and particles 
    circles1 = Circles(Nx,Ny) 
    particles1 = Particles(number_of_particles) 

    #read the airfoil geometry 
    panel_x = numpy.array([400,425,450,500,600,500,450,425,400]) 
    panel_y = numpy.array([300,325,330,320,300,280,270,275,300]) 
    panelcoordinates= numpy.dstack((panel_x,panel_y)) 

    #initialise PyGame 
    pygame.init() 
    clock = pygame.time.Clock() 
    Surface = pygame.display.set_mode((1000,600)) 

    while True: 
     ticks = clock.tick(6) 
     GetInput() 
     circles1.detectcollision(particles1) 
     particles1.move() 
     Draw(particles1,circles1) 

私はまた、エアロフォイルを描画する際に粗い刺し傷を作りました。データファイルの座標を知ることはできません。