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

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

3 个赞

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

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

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

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

5 个赞

LBM计算量本身不小，不能要求秒算出结果

2 个赞