代码优化问题,求大佬指点迷津

自己初学julia语言,终于把毕设的雏形做出来了,但是运行计算速度实在是很慢,又无从下手,等几个大佬来优化一下QwQ

# code
global cx = [0,1,0,-1,0,1,-1,-1,1] #水平方向速度分量
global cy = [0,0,1,0,-1,1,1,-1,-1] #锤子方向速度分量
global w = [4.0 / 9.0,1.0 / 9.0,1.0 / 9.0,1.0 / 9.0,1.0 / 9.0,1.0 / 36.0,1.0 / 36.0,1.0 / 36.0,1.0 / 36.0]
mstep=1000
nx=256
ny=256
Q=9
rho0=1.0
U=0.1
dx=1.0
dy=1.0
Lx=dx*nx
Ly=dy*ny
dt=dx
c=dx/dt
Re=1000
niu=U*Lx/Re
tau=3.0*niu+0.5
f=zeros(Float64,nx+1,ny+1,Q)
F=zeros(Float64,nx+1,ny+1,Q)
u = ones(Float64, nx+1, ny+1) #水平速度
u0 = ones(Float64, nx+1, ny+1)
v = ones(Float64, nx+1, ny+1) #垂直速度
v0 = ones(Float64, nx+1, ny+1)
rho = ones(Float64, nx+1, ny+1) #初试密度
function ph(k,rho,u,v)
    t1=cx[k]*u+cy[k]*v
    t2=u*u+v*v
    feq=w[k]*rho*(1.0+3.0*t1+4.5*t1*t1-1.5t2)
    return  feq
end
function evolution()
    for i=2:nx
        for j=2:nx
            for k=1:Q
                ip=i-cx[k]
                jp=j-cy[k]
                F[i,j,k]=f[ip,jp,k]+(ph(k,rho[ip,jp],u[ip,jp],v[ip,jp])-f[ip,jp,k])/tau
            end
        end
    end
    for i=2:nx
        for j=2:ny
            u0[i,j]=u[i,j]
            v0[i,j]=v[i,j]
            rho[i,j]=0
            u[i,j]=0
            v[i,j]=0
            for k=1:Q
                f[i,j,k]=F[i,j,k]
                rho[i,j]+=f[i,j,k]
                u[i,j]+=cx[k]*f[i,j,k]
                v[i,j]+=cy[k]*f[i,j,k]
            end
            u[i,j]/=rho[i,j]
            v[i,j]/=rho[i,j]
        end
    end
    for j=2:ny
        for k=1:Q
            rho[nx+1,j]=rho[nx,j]
            f[nx+1,j,k]=ph(k,rho[nx+1,j],u[nx+1,j],v[nx+1,j])+f[nx,j,k]-ph(k,rho[nx,j],u[nx,j],v[nx,j])
            rho[1,j]=rho[2,j]
            f[1,j,k]=ph(k,rho[1,j],u[1,j],v[1,j])+f[2,j,k]-ph(k,rho[2,j],u[2,j],v[2,j])
        end
    end
    for i=2:nx
        for k=1:Q
            rho[i,1]=rho[i,2]
            f[i,1,k]=ph(k,rho[i,1],u[i,1],v[i,1])+f[i,2,k]-ph(k,rho[i,2],u[i,2],v[i,2])
            rho[i,ny+1]=rho[i,ny]
            u[i,ny+1]=U
            f[i,ny+1,k]=ph(k,rho[i,ny+1],u[i,ny+1],v[i,ny+1])+f[i,ny,k]-ph(k,rho[i,ny],u[i,ny],v[i,ny])
        end
    end
end

for i=1:nx+1
    for j=1:ny+1
        u[i,j]=0
        v[i,j]=0
        rho[i,j]=rho0
        u[i,ny+1]=U
        for k=1:Q
            f[i,j,k]=ph(k,rho[i,j],u[i,j],v[i,j])
        end
    end
end
for kk=1:mstep
    evolution()
    if(kk%100==0)
        println("the u,v of point(nx/2,ny/2)is",u[129,129],",",v[129,129])

    end
end

建议你先读读文档的性能分析与性能建议部分吧,另外,请提供一个 benchmark.再附上你的 versioninfo() 也许会更好。
另外,你的代码我跑了一下,跑不出来。建议改一下例子(缩小规模),最好能够10s以内出结果。

3赞

emmm 我自己是跑的出来的,然后这是格子玻尔兹曼算法的一种应用

这是我迭代100次的benchmark
BenchmarkTools.Trial:
memory estimate: 19.11 GiB
allocs estimate: 1273520823

minimum time: 39.068 s (1.21% GC)
median time: 39.068 s (1.21% GC)
mean time: 39.068 s (1.21% GC)
maximum time: 39.068 s (1.21% GC)

samples: 1
evals/sample: 1

这是我更新后的代码 请大佬过目

# code
cx = [0,1,0,-1,0,1,-1,-1,1]::Array{Int64,1} #水平方向速度分量
cy = [0,0,1,0,-1,1,1,-1,-1]::Array{Int64,1}  #锤子方向速度分量
w = [4.0 / 9.0,1.0 / 9.0,1.0 / 9.0,1.0 / 9.0,1.0 / 9.0,1.0 / 36.0,1.0 / 36.0,1.0 / 36.0,1.0 / 36.0]::Array{Float64,1}
dx=1.0::Float64
dy=1.0::Float64
nx=256::Int64
ny=256::Int64
xx=129::Int64
Lx=dx*nx
Ly=Lx
dt=dx
c=dx/dt
Re=1000
U=0.1::Float64
niu=U*Lx/Re
tau=3.0*niu+0.5
function main(mstep::Int64)
    Q=9::Int64
    f=zeros(Float64,nx+1,ny+1,Q)
    F=zeros(Float64,nx+1,ny+1,Q)
    u = zeros(Float64, nx+1, ny+1) #水平速度
    u0= zeros(Float64, nx+1, ny+1)
    v = zeros(Float64, nx+1, ny+1) #垂直速度
    v0 = zeros(Float64, nx+1, ny+1)
    rho = ones(Float64, nx+1, ny+1) #初试密度
    rho0=1.0::Float64
    for i=1:nx+1
        for j=1:ny+1
            u[i,ny+1]=U
            for k=1:Q
                f[i,j,k]=ph(k,rho[i,j],u[i,j],v[i,j])
            end
        end
    end
    for kk=1:mstep
        evolution(nx, ny, Q, u, v, u0, v0, rho, cx, cy, f, F,U)
        hongguan(nx, ny, Q, u, v, u0, v0, rho, cx, cy, f, F,U)
        bianjie(nx, ny, Q, u, v, rho, cx, cy, f, F,U)
        if(kk%100==0)
            a=cuowu(nx,ny,u,v,u0,v0)
            println("the u,v of point(nx/2,ny/2)is:",kk,",",u[xx,xx],",",v[xx,xx])
            println("the max relative error of uv is:",a,"\n")
        end
    end
    o=open("dat.dat","w")
    write(o,"""Title="LMB Lid Driven Flow"\n""")
    write(o,"""Variables="x","y","u","v","rho"\n""")
    write(o,"""ZONE T= "BOX", I= 257, J=257, F=POINT\n""")
    for j=1:ny+1
        for i=1:nx+1
            print(o,i/(nx+1)," ",j/(ny+1)," ",u[i,j]," ",v[i,j]," ",rho[i,j],"\n")
        end
    end
    close(o)
end

function ph(k,rho,u,v)
    t1=cx[k]*u+cy[k]*v
    t2=u*u+v*v
    feq=w[k]*rho*(1.0+3.0*t1+4.5*t1*t1-1.5t2)
    return  feq
end

function evolution(nx, ny, Q, u, v, u0, v0, rho, cx, cy, f, F,U)
    for i=2:nx
        for j=2:nx
            for k=1:Q
                ip=i-cx[k]
                jp=j-cy[k]
                F[i,j,k]=f[ip,jp,k]+(ph(k,rho[ip,jp],u[ip,jp],v[ip,jp])-f[ip,jp,k])/tau
            end
        end
    end
end
function hongguan(nx, ny, Q, u, v, u0, v0, rho, cx, cy, f, F,U)
    for i=2:nx
        for j=2:ny
            u0[i,j]=u[i,j]
            v0[i,j]=v[i,j]
            rho[i,j]=0
            u[i,j]=rho[i,j]
            v[i,j]=rho[i,j]
            for k=1:Q
                f[i,j,k]=F[i,j,k]
                rho[i,j]+=f[i,j,k]
                u[i,j]+=cx[k]*f[i,j,k]
                v[i,j]+=cy[k]*f[i,j,k]
            end
            u[i,j]/=rho[i,j]
            v[i,j]/=rho[i,j]
        end
    end
end
function bianjie(nx, ny, Q, u, v, rho, cx, cy, f, F,U)
    for j=2:ny
        for k=1:Q
            rho[nx+1,j]=rho[nx,j]
            f[nx+1,j,k]=ph(k,rho[nx+1,j],u[nx+1,j],v[nx+1,j])+f[nx,j,k]-ph(k,rho[nx,j],u[nx,j],v[nx,j])
            rho[1,j]=rho[2,j]
            f[1,j,k]=ph(k,rho[1,j],u[1,j],v[1,j])+f[2,j,k]-ph(k,rho[2,j],u[2,j],v[2,j])
        end
    end
    for i=2:nx
        for k=1:Q
            rho[i,1]=rho[i,2]
            f[i,1,k]=ph(k,rho[i,1],u[i,1],v[i,1])+f[i,2,k]-ph(k,rho[i,2],u[i,2],v[i,2])
            rho[i,ny+1]=rho[i,ny]
            u[i,ny+1]=U
            f[i,ny+1,k]=ph(k,rho[i,ny+1],u[i,ny+1],v[i,ny+1])+f[i,ny,k]-ph(k,rho[i,ny],u[i,ny],v[i,ny])
        end
    end
end
function cuowu(nx,ny,u,v,u0,v0)
    temp1=0.0::Float64
    temp2=0.0::Float64
    for i=2:nx
        for j=2:ny
            temp1+=((u[i,j]-u0[i,j])*(u[i,j]-u0[i,j])+(v[i,j]-v0[i,j])*(v[i,j]-v0[i,j]))
            temp2+=(u[i,j]*u[i,j])+(v[i,j]*v[i,j])
        end
    end
    temp1=sqrt(temp1)
    temp2=sqrt(temp2)
    error=temp1/(temp2+1e-30)

end

其中包含了一个输出文件 dat.dat,以及mstep需要自己去设置,因为时间太久,所以设置为main(100)
BenchmarkTools.Trial:
memory estimate: 19.11 GiB
allocs estimate: 1273520823

minimum time: 39.068 s (1.21% GC)
median time: 39.068 s (1.21% GC)
mean time: 39.068 s (1.21% GC)
maximum time: 39.068 s (1.21% GC)

samples: 1
evals/sample: 1

建议再写下数学表达式什么的,介绍下背景。我感觉很难理解你的代码。因为我没有这方面的知识。另外,在我的电脑上运行还是很慢。4800H,一次要跑27s

2333 比我的快~实在抱歉让你难以理解了,我明天来补一哈相关的东西

请介绍下那个 ph 函数,似乎是它非常慢。evolution 函数也非常慢。还有你写的时候最好可以用矩阵表达式。我现在完全没有思绪从哪里开始优化,但是我觉得你那个 ph 函数可以做成矩阵或者向量的形式。

这个代码还是写的有点可怕的。。。
最严重的问题是全局变量。请参考文章中的性能建议部分的第一条 避免使用全局变量

处理全局变量前(原代码)
original

处理后 (将全局变量放入main中,修改下个别函数)
new

如果我没理解错的话…这个大概是一个关于全局变量使用的非常好的反面教材…

4赞

LBM计算量本身不小,不能要求秒算出结果 :grinning:
不考虑GPU的话,Julia的循环和矢量化写法是一样快的。
对这个程序的优化而言,除了避免写global变量,还可以注意数组循环的顺序。保证循环的方向从数组内层开始,按照内存的顺序读写。
流体力学evolution和update都是并行粒度很高,可以用并行化加速循环,简单的做法是使用Threads.@threads或@distributed宏命令。

1赞

备案号:京ICP备17009874号-2