激光物理模拟
激光物理模拟
提问于 2022-02-08 14:28:33
这是用于激光物理学的计算。速度相当慢，这似乎被指数j、k和L的三个嵌套for循环所延迟(下面有注释)。计算时间似乎与nx的第四次方成正比，但我想把nx增加到100。
我也猜七大指数Numpy阵列S拖了很多计算。如果我能删除三个嵌套的for循环，那就更好了。然而，这似乎是不可能的。有没有提高速度的方法？
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np
import matplotlib as mpl
import math
from scipy import linalg, optimize
from scipy.integrate import quad
import scipy.integrate as integrate
import scipy.special as special
from matplotlib.colors import LogNorm
import time
from datetime import datetime
import copy
nz=70
nx=3
ny=nx
ntau=70
gamsp = 1/(160.*10.**-12)
Rincrease=10**0
R = Rincrease*400*10.**-9
nnu = 1.6*10.**25*Rincrease**-2
c = 3.*10.**8
n = nnu*math.pi*R**2
lambda0= 12.4/0.88*10**-10
delo = 4.*10.**-6
tau_max= 2./gamsp
z_max=3.*10**-4
L=z_max
dtau = tau_max/ntau
dx=2*R/nx
dy=dx
dz = z_max/nz
Omega=2*math.pi*c/lambda0
alpha=np.arctan(R/L)
S=np.empty((nx,ny,nz,nx,ny,nz,ntau),dtype='complex_')
S[:,:,:,:,:,:,0]=0
rho=np.empty((nx,ny,nz,ntau),dtype='complex_')
for i in range(nx):
for j in range(ny):
rho[i,j,:,0]=1 if math.sqrt((i-(nx-1)/2)**2+(j-(ny-1)/2)**2)nx-1:
g[i]=copy.copy(np.conjugate(g[2*nx-2-i]))
else:
for j in range (2*(ny-1)+1):
if j>ny-1:
g[i,j]=copy.copy(np.conjugate(g[i,2*ny-2-j]))
else:
for k in range (2*(nz-1)+1):
if k>nz-1:
g[i,j,k]=copy.copy(np.conjugate(g[i,j,2*nz-2-k]))
else:
xvalue=copy.copy(-2*R+4*R*i/(2*(nx-1)))
yvalue=copy.copy(-2*R+4*R*j/(2*(nx-1)))
radius=copy.copy(math.sqrt(xvalue**2+yvalue**2))
zvalue=copy.copy(-L+2*L*k/(2*(nz-1)))
Funcre = "2*math.pi*np.sin(x)*(np.cos(x))**2*np.cos(Omega/c*({})*(np.cos(x)-1))*special.jv(0,Omega/c*{}*math.sin(x))".format(zvalue,radius)
Funcim = "2*math.pi*np.sin(x)*(np.cos(x))**2*np.sin(Omega/c*({})*(np.cos(x)-1))*special.jv(0,Omega/c*{}*math.sin(x))".format(zvalue,radius)
def fre(x):
fre=eval(Funcre)
return fre
def fim(x):
fim=eval(Funcim)
return fim
g[i,j,k]=quad(fre,0,alpha)[0]+1j*quad(fim,0,alpha)[0]
for i in range(2*(nz-1)+1):
g2[:,:,i]=copy.copy(g[:,:,i])*(0 if i>nz else 1)
for i in range(ntau):
i1=np.zeros((nx,ny,nz),dtype='complex_')
i2=np.zeros((nx,ny,nz),dtype='complex_')
if ik else snoise2[:,:,j,:,:,k]
rho[:,:,:,i+1]=rho[:,:,:,i]+gamsp*(-2*(1+rho[:,:,:,i])/2-3*nnu/(8*math.pi)*2*np.real(s1))*dtau
S[:,:,:,:,:,:,i+1]=S[:,:,:,:,:,:,i]+gamsp*(-S[:,:,:,:,:,:,i]+3/(16*math.pi)*(nnu*(s2+s3)+snoise))*dtau
for i in range(ntau):
for jz in range (nz):
Intensity[jz,i]=np.sum(Intensitytrans[:,:,jz,i])/(nx*ny)
peak=np.empty((nz,2))
for i in range(nz):
photonn[i]=np.sum(Intensity[i,:])/ntau*tau_max*delo*math.pi*R**2
peak[i,0]=np.argmax(Intensity[i,:])/ntau
peak[i,1]=i/nz
if max(Intensity[i,:])!=0:
Intensityn[i,:]=Intensity[i,:]/max(Intensity[i,:])回答 1
Code Review用户
发布于 2022-02-08 18:20:08
速度相当慢，这似乎被指数j、k和L的三个嵌套for循环所延迟。
这是一个很好的开端，可以发现代码的哪一部分造成了瓶颈。猜测哪个部分是慢的可能有助于建立一个直觉，所以让我们测试代码，看看什么是减速。
我们可以用几种方法来衡量代码的性能。因为我对代码不了解，所以我会首先向它抛出一个通用分析器。Python附带了一个名为cProfile的分析器。我发现另一个有用的分析器是琵琶。
通过cProfile python -m cProfile -s time cr_laser_physics.py运行代码并在一分钟后将其杀死，这给了我一个概要文件。它的顶部看起来像这样
801051 function calls (786162 primitive calls) in 62.441 seconds
Ordered by: internal time
ncalls tottime percall cumtime percall filename:lineno(function)
1 45.036 45.036 62.804 62.804 cr_laser_physics.py:1()
13081/8182 13.425 0.001 15.052 0.002 {built-in method numpy.core._multiarray_umath.implement_array_function}
26476 1.239 0.000 1.612 0.000 {built-in method builtins.eval}
12644 1.195 0.000 1.195 0.000 {method 'reduce' of 'numpy.ufunc' objects}
4719 0.226 0.000 13.868 0.003 numeric.py:943(tensordot)
104/102 0.201 0.002 0.203 0.002 {built-in method _imp.create_dynamic}
448 0.073 0.000 0.073 0.000 {method 'read' of '_io.FileIO' objects}
448 0.072 0.000 0.072 0.000 {built-in method marshal.loads}
14217 0.057 0.000 0.057 0.000 {method 'reshape' of 'numpy.ndarray' objects}
1 0.034 0.034 0.034 0.034 {built-in method mkl._py_mkl_service.get_version}
...不幸的是，这看起来没有多大帮助。由于大部分时间都花在模块中(45秒)，剩下的是我们没有直接调用的numpy方法，所以我们没有学到多少。我们可以肯定地说，整个代码是缓慢的。
谷歌搜索把我带到了这个问题。他们遇到的问题看起来与我们的相似，cProfile并没有具体说明性能问题在哪里。建议的答案是将代码分解成更多的函数，这样cProfile就可以知道哪些函数比较慢。
因此，这似乎是一个很好的时间，通过代码，并作出任何改进，我们可以。一些美好的事情是
- 删除未使用的进口产品。他们给轮廓仪增加了一点噪音。皮兰斯建议，大多数进口产品都是不需要的。您可以随时在需要时将它们添加回。
special导入被错误标记为未使用，因为Pylance无法判断代码字符串将在这里运行。这是一个很好的指示，表明代码正在以一种非常奇怪的方式做一些事情。是否可以简化呢？ - 把事情做好。虽然每个函数调用都有开销，但我们需要它，以便cProfile能够很好地分析。
- 将长队分割成较小的块。这将使您更容易了解代码的哪些部分取决于什么。
我不会实现所有这些建议，但我将尝试概述文件内容将是什么样子。
# Imports
import copy
import math
import numpy as np
import scipy.special as special
from scipy.integrate import quad
# Global constants. (The usual convention is to make these all uppercase).
nz = 70
ny = nx = 3
ntau = 70
...
def main():
S=np.empty((nx,ny,nz,nx,ny,nz,ntau),dtype='complex_')
S[:,:,:,:,:,:,0]=0
rho=np.empty((nx,ny,nz,ntau),dtype='complex_')
...
if __name__ == "__main__":
main()从这一点开始，我将选择主方法的内容，直到它被分割成尽可能多的函数。这里最难的部分是跟踪哪些变量依赖于哪些其他变量。
这样做，我最终得到了一个主要的方法，它看起来像
def main():
rho = make_rho()
g = make_g()
g2 = make_g2(g)
Intensitytrans = big_loop(g, g2, rho)
Intensity = make_intensity(Intensitytrans)
Intensityn, photonn = smaller_loop(Intensity)在这里，我将变量保持不变，但尽可能多地将代码移到函数中。每个函数都遵循类似的模式：“创建变量->循环以填充它们的值，->返回稍后引用的变量”。举个例子，这里是make_g2。
def make_g2(g):
g2=np.zeros((2*(nx-1)+1,2*(ny-1)+1,2*(nz-1)+1),dtype ='complex_')
for i in range(2*(nz-1)+1):
g2[:,:,i]=copy.copy(g[:,:,i])*(0 if i>nz else 1)
return g2在代码上运行分析器现在给出：(只显示5个最昂贵的函数)
516776 function calls (504160 primitive calls) in 60.096 seconds
Ordered by: internal time
ncalls tottime percall cumtime percall filename:lineno(function)
1 44.271 44.271 59.010 59.010 cr_laser_physics.py:75(big_loop)
19107/11944 12.793 0.001 14.640 0.001 {built-in method numpy.core._multiarray_umath.implement_array_function}
19102 1.385 0.000 1.385 0.000 {method 'reduce' of 'numpy.ufunc' objects}
26460 0.803 0.000 1.060 0.000 {built-in method builtins.eval}
7163 0.236 0.000 13.265 0.002 numeric.py:943(tensordot)这告诉我们两条信息。您是对的，大多数情况下(代码运行60秒时为45秒)。花在重嵌套的循环中。与这些循环相比，其他函数的速度并不明显，因此我们现在不需要担心它们。
缓慢的函数如下所示
@profile
def big_loop(g, g2, rho):
# This has a lot going on, can it be simplified?
S=np.empty((nx,ny,nz,nx,ny,nz,ntau),dtype='complex_')
S[:,:,:,:,:,:,0]=0
Intensitytrans=np.zeros((nx,ny,nz,ntau),dtype='complex_')
for i in range(ntau):
i1=np.zeros((nx,ny,nz),dtype='complex_')
i2=np.zeros((nx,ny,nz),dtype='complex_')
if ik else snoise2[:,:,j,:,:,k]
rho[:,:,:,i+1]=rho[:,:,:,i]+gamsp*(-2*(1+rho[:,:,:,i])/2-3*nnu/(8*math.pi)*2*np.real(s1))*dtau
S[:,:,:,:,:,:,i+1]=S[:,:,:,:,:,:,i]+gamsp*(-S[:,:,:,:,:,:,i]+3/(16*math.pi)*(nnu*(s2+s3)+snoise))*dtau
return Intensitytrans因为拆分这个函数看起来并不容易，接下来要做的就是使用不同的分析器。在这里，逐行分析器应该更有帮助。
我通过pip安装了线路_轮廓仪，添加了@profile来指定要配置文件的函数，并使用kernprof -l cr_laser_physics.py运行它，60秒后将其终止。然后，我使用python -m line_profiler cr_laser_physics.py.lprof查看了配置文件。只保留超过5%的时间的线显示三行是真的很慢。
Timer unit: 1e-06 s
Total time: 60.4413 s
File: cr_laser_physics.py
Function: big_loop at line 74
Line # Hits Time Per Hit % Time Line Contents
==============================================================
74 @profile
75 def big_loop(g, g2, rho):
97 2433 22926489.0 9423.1 37.9 i1[j,k,l]=np.sum(S[:,:,:,:,:,:,i]*np.tensordot(np.conjugate(g2c),g2c,axes=0))*(dx*dy*dz)**2
101 2433 18182798.0 7473.4 30.1 s2+=np.tensordot(np.conjugate(g2c)*rho[:,:,:,i],S[j,k,l,:,:,:,i],axes = 0)*dx*dy*dz
102 2433 18384856.0 7556.5 30.4 s3+=np.tensordot(S[:,:,:,j,k,l,i],g2c*rho[:,:,:,i],axes = 0)*dx*dy*dz这些线路中的每一条都是“热”的，需要很长时间，而且运行了很多次。我对可以改变的东西不太自信，所以我会在下面留下一些建议。
i1[j,k,l] = np.sum(S[:,:,:,:,:,:,i]*np.tensordot(np.conjugate(g2c),g2c,axes=0))*(dx*dy*dz)**2np.sum(S[:,:,:,:,:,:,i]和(dx*dy*dz)**2都会给出相同的结果，即使j、k和l发生了变化。它们可以在三重嵌套循环之前完成一次。我怀疑这一总和是非常昂贵的计算太多。
每个慢行都有一个调用np.conjugate(g2c)。我建议这样做一次，并将其存储在一个变量中。
后两条慢线都是将张量乘积加到和中。您能为np.tensordot指定更多的轴以删除某些循环嵌套吗？通常，纯python中的循环和求和比对numpy函数的等效调用要慢。
页面原文内容由Code Review提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://codereview.stackexchange.com/questions/273874
复制相关文章
相似问题
腾讯云开发者
