2016-03-02 89 views
1

我使用vispy来生成三维海面。但它不是很光滑。我不知道应该改进哪一部分。有谁能告诉我吗?下面是代码:如何让我的python动画流畅?

from numpy import linspace,zeros,sin,pi,exp,sqrt 
from vispy import app,scene 
import sys 
from vispy.util.filter import gaussian_filter 

def I(x, y): 
    return exp(-(x-Lx/2.0)**2/2.0 -(y-Ly/2.0)**2/2.0) 

def f(x, y, t): 
    return sin(2*x*y*pi*t/Lx) #defined by myself 

def bc(x, y, t): 
    return 0.0 

def solver0(I, f, c, bc, Lx, Ly, nx, ny, dt, t, user_action=None): 
    dx = Lx/float(nx) 
    dy = Ly/float(ny) 
    x = linspace(0, Lx, nx+1) #grid points in x dir 
    y = linspace(0, Ly, ny+1) #grid points in y dir 
    if dt <= 0:    #max time step? 
     dt = (1/float(c))*(1/sqrt(1/dx**2 + 1/dy**2)) 
    Cx2 = (c*dt/dx)**2 
    Cy2 = (c*dt/dy)**2 #help variables 
    dt2 = dt**2 

    up = zeros((nx+1,ny+1)) #solution array 
    u = up.copy()   #solution at t-dt 
    um = up.copy()   #solution at t-2*dt 

    #set initial condition: 
    #t =0.0 
    for i in range(0,nx): 
     for j in range(0,ny): 
      u[i,j] = I(x[i], y[j]) 
    for i in range(1,nx-1): 
     for j in range(1,ny-1): 
      um[i,j] = u[i,j] + \ 
         0.5*Cx2*(u[i-1,j] - 2*u[i,j] + u[i+1,j]) + \ 
         0.5*Cy2*(u[i,j-1] - 2*u[i,j] + u[i,j+1]) + \ 
         dt2*f(x[i], y[j], t) 
    #boundary values of um (equals t=dt when du/dt=0) 
    i = 0 
    for j in range(0,ny): um[i,j] = bc(x[i], y[j], t+dt) 
    j = 0 
    for i in range(0,nx): um[i,j] = bc(x[i], y[j], t+dt) 
    i = nx 
    for j in range(0,ny): um[i,j] = bc(x[i], y[j], t+dt) 
    j = ny 
    for i in range(0,nx): um[i,j] = bc(x[i], y[j], t+dt) 

    if user_action is not None: 
     user_action(u, x, y, t) #allow user to plot etc. 

    while t <= tstop: 
     t_old = t 
     t += dt 

     #update all inner points: 
     for i in range(1,nx-1): 
      for j in range(1,ny-1): 
       up[i,j] = -um[i,j] + 2*u[i,j] + \ 
          Cx2*(u[i-1,j] - 2*u[i,j] + u[i+1,j]) + \ 
          Cy2*(u[i,j-1] - 2*u[i,j] + u[i,j+1]) + \ 
          dt2*f(x[i], y[j], t_old) 

     #insert boundary conditions: 
     i = 0 
     for j in range(0,ny): up[i,j] = bc(x[i], y[j], t) 
     j = 0 
     for i in range(0,nx): up[i,j] = bc(x[i], y[j], t) 
     i = nx 
     for j in range(0,ny): up[i,j] = bc(x[i], y[j], t) 
     j = ny 
     for i in range(0,nx): up[i,j] = bc(x[i], y[j], t) 

     if user_action is not None: 
      user_action(up, x, y, t) 

     um, u, up = u, up, um #update data structures 
    return up #dt might be computed in this function 

Lx = 10 
Ly = 10 
c = 1.0 
dt = 0 
nx = 40 
ny = 40 
tstop = 20 
x = linspace(0, Lx, nx+1) #grid points in x dir 
y = linspace(0, Ly, ny+1) #grid points in y dir 

canvas = scene.SceneCanvas(keys='interactive') 
view = canvas.central_widget.add_view() 
view.camera = scene.TurntableCamera(up='z') 

p1 = scene.visuals.SurfacePlot(z=solver0(I, f, c, bc, Lx, Ly, nx, ny, dt, 0, user_action=None),color=(0,0,1,1),shading='smooth') 
p1.transform = scene.transforms.AffineTransform() 
p1.transform.scale([1/29.,1/29.,1.0]) 
p1.transform.translate([-1.0,-0.5,0]) 

view.add(p1) 
t = 0.0 
def update(ev): 
    global p1 
    global t 
    t += 1.0 
    p1.set_data(z=solver0(I, f, c, bc, Lx, Ly, nx, ny, dt, t, user_action=None)) 

timer = app.Timer() 
timer.connect(update) 
timer.start(0) 
if __name__ == '__main__': 
    canvas.show() 
    if sys.flags.interactive == 0: 
     app.run() 

回答

1

尝试使用这样的:

up[1:nx-1,1:nx-1]=-um[1:nx-1,1:ny-1]+2*u[1:nx-1,1:ny-1]+ \ 
       Cx2*(u[0:nx-2,1:ny-1]-2*u[1:nx-1,1:ny-1]+ u[2:nx,1:ny-1]) + \ 
       Cy2*(u[1:nx-1,0:ny-2]-2*u[1:nx-1,1:ny-1]+ u[1:nx-1,2:ny]) + \ 
       [[dt2*f(x[i],y[j],t_old) for i in range(1,nx-1)] for j in range(1,ny-1)] 

,而不是

#update all inner points: 
    for i in range(1,nx-1): 
     for j in range(1,ny-1): 
      up[i,j] = -um[i,j] + 2*u[i,j] + \ 
         Cx2*(u[i-1,j] - 2*u[i,j] + u[i+1,j]) + \ 
         Cy2*(u[i,j-1] - 2*u[i,j] + u[i,j+1]) + \ 
         dt2*f(x[i], y[j], t_old) 

切片是非常快的在Python等都是列表理解。我基本上通过添加矩阵来计算单个表达式中的整个double循环。

告诉我,如果有任何错误(可能有intendation错误)

编辑:一些错误的指标,一些由矩阵测试这在100次;我的变体比double for循环快大约10倍。您可以用类似的方式更改程序中的其他double for循环。

EDITEDIT:

列表综合是几乎总是比for循环快,所以不是

for j in range(0,ny): up[i,j] = bc(x[i], y[j], t) 

up[i,0:ny] = [bc(x[i],y[j],t) for j in range(0,ny)] 

这应该是快一点 另外,如果你是使用Python 2.7使用xrange而不是范围来节省内存。

这里有一些更多的例子说明如何为环插入列表解析: http://www.u.arizona.edu/~erdmann/mse350/topics/list_comprehensions.html

这是一个有点复杂在第一,但definitly值得!

+0

改变double for loop后,它确实变得更加平滑。我是一条新鱼,我甚至不知道这一点。非常感谢您的建议。这是我能改进的唯一部分吗?我应该用我的单循环做些什么吗? –

+0

我只知道切片,因为我已经在细胞自动机中使用它;) 但是,您应该阅读列表解析,它们只是更快。我将编辑答案,为您的单个循环添加列表组件。给我一点时间。同时你可以看看这个:https://wiki.python.org/moin/PythonSpeed/PerformanceTips – JeD

+0

非常感谢你:) –