改写成loop以后，我稍微加了点输出到面板和画图的代码，如下
using ShiftedArrays, Plots
function mainFlow(iterStepMax::Int64)
Q = 9 #D2Q9
nx = 1200 #lattice number in horizontal direction
ny = 600 #lattice number in vertical direction
U = 0.1 #inlet velocity
cx = [1 0 -1 0 1 -1 -1 1 0] #velocity components in horizontal direction
cy = [0 1 0 -1 1 1 -1 -1 0] #velocity components in vertical direction
w = [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 4.0 / 9.0] #weight parameters
rel_err_temp = ones(Float64, nx, ny)
ρ = ones(Float64, nx, ny) #initial density
u = zeros(Float64, nx, ny) #horizontal velocity component
v = zeros(Float64, nx, ny) #vertical velocity component
u0 = zeros(Float64, nx, ny) #horizontal velocity component before iteration
v0 = zeros(Float64, nx, ny) #vertical velocity component befre iteration
uv = zeros(Float64, nx, ny) #resultant velocity
f = zeros(Float64, nx, ny, Q) #distribution function
F = zeros(Float64, nx, ny, Q)
feq = zeros(Float64, nx, ny, Q) #equilibrium distribution function
dx = 1.0 #horizontal lattice length
dy = 1.0 #vertical lattice length
lx = dx * nx #horizontal domain length
ly = dy * ny #vertical domain length
dt = dx #dt = 1.0
c = dx / dt #lattice sound speed
Re = 1000.0 #Reynoldz number
ν = U * ly / Re #kinematic viscousity
τ = 3.0ν + 0.5 #relaxation time
ω = 1.0 / τ #relaxation frequency
println("Re = ", Re, ", τ = ", τ) #instead of "Re = $Re, τ = $τ"
ε = 1.0e-6 #convergence requirement
@views u[1, 2:ny-1] .= U
for iterStep = 1:iterStepMax
collision(nx, ny, Q, u, v, ρ, cx, cy, w, ω, f, feq)
streaming(Q, cx, cy, f, F)
boundaryCondition(nx, ny, u, v, ρ, U, f)
macroQuantity(nx, ny, Q, u, v, uv, u0, v0, ρ, cx, cy, f)
@. rel_err_temp = (u - u0)^2 + (v - v0)^2
temp1 = sum(rel_err_temp)
@. rel_err_temp = u^2 + v^2
temp2 = sum(rel_err_temp)
err = sqrt(temp1 / (temp2 + 1e-30))
if fld(iterStep, 100) * 100 == iterStep #instead of "floor"
println(
"iterStep = $iterStep",
", residual = ",
round(err, digits = 8),
", top-left velocity = ",
round(uv[300, 450], digits = 8),
", bottom-left velocity = ",
round(uv[300, 150], digits = 8),
", top-right velocity = ",
round(uv[900, 450], digits = 8),
", bottom-right velocity = ",
round(uv[900, 150], digits = 8),
)
end
if err <= ε
break
end
end
clibrary(:colorcet)
contourf(uv, color = :rainbow, aspect_ratio = 1.0)
end
function collision(nx, ny, Q, u, v, ρ, cx, cy, w, ω, f, feq)
for j = 1:ny, i = 1:nx
t1 = u[i, j] * u[i, j] + v[i, j] * v[i, j]
for k = 1:Q
t2 = u[i, j] * cx[k] + v[i, j] * cy[k]
feq[i, j, k] = ρ[i, j] * w[k] * (1.0 + 3.0t2 + 4.5 * t2^2 - 1.5t1)
f[i, j, k] = (1.0 - ω) * f[i, j, k] + ω * feq[i, j, k]
end
end
end
function streaming(Q, cx, cy, f, F)
copy!(F, f)
for k = 1:Q
@views f[:, :, k] = ShiftedArrays.circshift(F[:, :, k], (cx[k], cy[k]))
end
end
function boundaryCondition(nx, ny, u, v, ρ, U, f)
#left and right
for j = 2:ny-1
f[1, j, 1] = f[1, j, 3] + ρ[1, j] * U * 2.0 / 3.0
f[1, j, 5] = f[1, j, 7] - 0.5 * (f[1, j, 2] - f[1, j, 4]) + ρ[1, j] * U / 6.0
f[1, j, 8] = f[1, j, 6] + 0.5 * (f[1, j, 2] - f[1, j, 4]) + ρ[1, j] * U / 6.0
f[nx, j, 3] = f[nx-1, j, 3]
f[nx, j, 6] = f[nx-1, j, 6]
f[nx, j, 7] = f[nx-1, j, 7]
end
u[1, 2:ny-1] .= U
v[1, 2:ny-1] .= 0.0
#bottom and top
for i = 1:nx
f[i, 1, 2] = f[i, 1, 4]
f[i, 1, 5] = f[i, 1, 7]
f[i, 1, 6] = f[i, 1, 8]
f[i, ny, 4] = f[i, ny, 2]
f[i, ny, 8] = f[i, ny, 6]
f[i, ny, 7] = f[i, ny, 5]
end
u[:, 1] .= 0.0
v[:, 1] .= 0.0
u[:, ny] .= 0.0
v[:, ny] .= 0.0
end
function macroQuantity(nx, ny, Q, u, v, uv, u0, v0, ρ, cx, cy, f)
sum!(ρ, f)
copy!(u0, u)
copy!(v0, v)
for j = 1:ny, i = 1:nx
uSum = 0.0
vSum = 0.0
for k = 1:Q
uSum += f[i, j, k] * cx[k]
vSum += f[i, j, k] * cy[k]
end
u[i, j] = uSum / ρ[i, j]
v[i, j] = vSum / ρ[i, j]
end
@. uv = sqrt(u^2 + v^2)
end
.mem文件
- using ShiftedArrays, Plots
-
- function mainFlow(iterStepMax::Int64)
211946304 Q = 9 #D2Q9
- nx = 1200 #lattice number in horizontal direction
- ny = 600 #lattice number in vertical direction
- U = 0.1 #inlet velocity
160 cx = [1 0 -1 0 1 -1 -1 1 0] #velocity components in horizontal direction
160 cy = [0 1 0 -1 1 1 -1 -1 0] #velocity components in vertical direction
160 w = [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 4.0 / 9.0] #weight parameters
5760080 rel_err_temp = ones(Float64, nx, ny)
5760080 ρ = ones(Float64, nx, ny) #initial density
5760080 u = zeros(Float64, nx, ny) #horizontal velocity component
5760080 v = zeros(Float64, nx, ny) #vertical velocity component
5760080 u0 = zeros(Float64, nx, ny) #horizontal velocity component before iteration
5760080 v0 = zeros(Float64, nx, ny) #vertical velocity component befre iteration
5760080 uv = zeros(Float64, nx, ny) #resultant velocity
51840080 f = zeros(Float64, nx, ny, Q) #distribution function
51840080 F = zeros(Float64, nx, ny, Q)
51840080 feq = zeros(Float64, nx, ny, Q) #equilibrium distribution function
- dx = 1.0 #horizontal lattice length
- dy = 1.0 #vertical lattice length
- lx = dx * nx #horizontal domain length
- ly = dy * ny #vertical domain length
- dt = dx #dt = 1.0
- c = dx / dt #lattice sound speed
- Re = 1000.0 #Reynoldz number
- ν = U * ly / Re #kinematic viscousity
- τ = 3.0ν + 0.5 #relaxation time
- ω = 1.0 / τ #relaxation frequency
4941680 println("Re = ", Re, ", τ = ", τ) #instead of "Re = $Re, τ = $τ"
- ε = 1.0e-6 #convergence requirement
-
0 @views u[1, 2:ny-1] .= U
0 for iterStep = 1:iterStepMax
0 collision(nx, ny, Q, u, v, ρ, cx, cy, w, ω, f, feq)
0 streaming(Q, cx, cy, f, F)
0 boundaryCondition(nx, ny, u, v, ρ, U, f)
0 macroQuantity(nx, ny, Q, u, v, uv, u0, v0, ρ, cx, cy, f)
0 @. rel_err_temp = (u - u0)^2 + (v - v0)^2
0 temp1 = sum(rel_err_temp)
0 @. rel_err_temp = u^2 + v^2
0 temp2 = sum(rel_err_temp)
0 err = sqrt(temp1 / (temp2 + 1e-30))
0 if fld(iterStep, 100) * 100 == iterStep #instead of "floor"
0 println(
- "iterStep = $iterStep",
- ", residual = ",
- round(err, digits = 8),
- ", top-left velocity = ",
- round(uv[300, 450], digits = 8),
- ", bottom-left velocity = ",
- round(uv[300, 150], digits = 8),
- ", top-right velocity = ",
- round(uv[900, 450], digits = 8),
- ", bottom-right velocity = ",
- round(uv[900, 150], digits = 8),
- )
- end
0 if err <= ε
0 break
- end
- end
0 clibrary(:colorcet)
64 contourf(uv, color = :rainbow, aspect_ratio = 1.0)
- end
-
- function collision(nx, ny, Q, u, v, ρ, cx, cy, w, ω, f, feq)
0 for j = 1:ny, i = 1:nx
0 t1 = u[i, j] * u[i, j] + v[i, j] * v[i, j]
0 for k = 1:Q
0 t2 = u[i, j] * cx[k] + v[i, j] * cy[k]
0 feq[i, j, k] = ρ[i, j] * w[k] * (1.0 + 3.0t2 + 4.5 * t2^2 - 1.5t1)
0 f[i, j, k] = (1.0 - ω) * f[i, j, k] + ω * feq[i, j, k]
- end
- end
- end
-
- function streaming(Q, cx, cy, f, F)
0 copy!(F, f)
0 for k = 1:Q
43200 @views f[:, :, k] = ShiftedArrays.circshift(F[:, :, k], (cx[k], cy[k]))
- end
- end
-
- function boundaryCondition(nx, ny, u, v, ρ, U, f)
- #left and right
0 for j = 2:ny-1
0 f[1, j, 1] = f[1, j, 3] + ρ[1, j] * U * 2.0 / 3.0
0 f[1, j, 5] = f[1, j, 7] - 0.5 * (f[1, j, 2] - f[1, j, 4]) + ρ[1, j] * U / 6.0
0 f[1, j, 8] = f[1, j, 6] + 0.5 * (f[1, j, 2] - f[1, j, 4]) + ρ[1, j] * U / 6.0
0 f[nx, j, 3] = f[nx-1, j, 3]
0 f[nx, j, 6] = f[nx-1, j, 6]
0 f[nx, j, 7] = f[nx-1, j, 7]
- end
0 u[1, 2:ny-1] .= U
0 v[1, 2:ny-1] .= 0.0
-
- #bottom and top
0 for i = 1:nx
0 f[i, 1, 2] = f[i, 1, 4]
0 f[i, 1, 5] = f[i, 1, 7]
0 f[i, 1, 6] = f[i, 1, 8]
0 f[i, ny, 4] = f[i, ny, 2]
0 f[i, ny, 8] = f[i, ny, 6]
0 f[i, ny, 7] = f[i, ny, 5]
- end
0 u[:, 1] .= 0.0
0 v[:, 1] .= 0.0
0 u[:, ny] .= 0.0
2400 v[:, ny] .= 0.0
- end
-
-
- function macroQuantity(nx, ny, Q, u, v, uv, u0, v0, ρ, cx, cy, f)
0 sum!(ρ, f)
0 copy!(u0, u)
0 copy!(v0, v)
0 for j = 1:ny, i = 1:nx
- uSum = 0.0
- vSum = 0.0
0 for k = 1:Q
0 uSum += f[i, j, k] * cx[k]
0 vSum += f[i, j, k] * cy[k]
- end
0 u[i, j] = uSum / ρ[i, j]
0 v[i, j] = vSum / ρ[i, j]
- end
0 @. uv = sqrt(u^2 + v^2)
- end
-
除了 2400 v[:, ny] .= 0.0 有点不一致以外，又改善了一些
性能测试如下
@benchmark mainFlow(50)
BenchmarkTools.Trial:
memory estimate: 197.99 MiB
allocs estimate: 3231
--------------
minimum time: 4.342 s (3.28% GC)
median time: 4.362 s (1.65% GC)
mean time: 4.362 s (1.65% GC)
maximum time: 4.382 s (0.03% GC)
--------------
samples: 2
evals/sample: 1
@benchmark mainFlow(100)
BenchmarkTools.Trial:
memory estimate: 198.04 MiB
allocs estimate: 4244
--------------
minimum time: 8.771 s (0.01% GC)
median time: 8.771 s (0.01% GC)
mean time: 8.771 s (0.01% GC)
maximum time: 8.771 s (0.01% GC)
--------------
samples: 1
evals/sample: 1