我试图用Python计算许多粒子位置的演化。最终我将有许多粒子(大约10000个)在100000个时间步上被计算出来。由于我正在不惜一切代价避免使用Fortran,我正在努力加快这个过程。你知道吗
我要解的方程是
d X_i/dt = u
d Y_i/dt = v
所以理想情况下我会使用二维数组:一个在粒子之间变化,另一个在x
和y
之间变化。如果我有100个粒子,我会有一个100x2
数组。你知道吗
问题的一部分是因为scipy.integrate.odeint
只接受一维数组,所以我必须将初始条件展平,在导数函数(RHS_im
)中拆分,然后在输出时再次展平,这很慢(这占RHS_im
调用的20%左右)。你知道吗
我可以用一种丑陋的方式来做。这是我贴的一张MWE
import numpy as np
from scipy.interpolate import RegularGridInterpolator
from scipy.integrate import odeint
x=np.arange(2.5, 500, 5)
y=np.arange(2.5, 500, 5)
X, Y = np.meshgrid(x, y, indexing='xy')
U=np.full_like(X, -0.1)
V=0.2*np.sin(2*np.pi*X/200)
tend=100
dt=1
tsteps=np.arange(0, tend, dt)
#------
# Create interpolation
U_int = RegularGridInterpolator((x, y), U.T)
V_int = RegularGridInterpolator((x, y), V.T)
#------
#------
# Initial conditions
x0=np.linspace(200,300,5)
y0=np.linspace(200,300,5)
#------
#------
# Calculation for many
def RHS_im(XY, t):
X, Y = np.split(XY, 2, axis=0)
pts=np.array([X,Y]).T
return np.concatenate([ U_int(pts), V_int(pts) ])
XY0 = np.concatenate([x0, y0])
XY, info = odeint(RHS_im, XY0, tsteps[:-1], args=(), hmax=dt, hmin=dt, atol=0.1, full_output=True)
X, Y = np.split(XY, 2, axis=1)
有没有办法避免拆分过程?你知道吗
此外,虽然我选择了odeint
来集成这个系统,但我并不致力于这个选择。如果有另一个函数做得更快,我会很容易改变(即使是用一个简单的Euler格式)。我只是没找到。(同样的事情也适用于插值方案,它大约占RHS_im
所需时间的80%)。你知道吗
编辑
稍微改进了分裂过程(通过使用np.split
),但是整个程序仍然需要改进。你知道吗
感谢您对MWE的澄清。通过将
U_int
和V_int
组合成UV_int
并改变XY的布局,我可以让MWE运行快50%左右np.重塑可以用来代替np.拆分. 这两种变化都有助于就地和连续地访问内存。你知道吗相关问题 更多 >
编程相关推荐