#訂正
リスタートのためのデータ書き出しで，無駄に同じデータを別名で保存していたところを修正しました。
すいませんでした。
カルマン渦の計算
下記のようなカルマン渦の計算をするためのコードです。
背景
以前下記の記事を書きましたが，それにTVDスキームも実装しました。
TVD
Total Validation Diminishing だったかな。実装はVersteeg, Malalasekeraの「数値流体力学」を参考にしました。
下記のコードだと，makesourceの後半でTVDの実装してます。PISO法のa_P等にもいれる必要がありますが，
これはsolveUVAでやってます。
Paramのなかのパラメータで，TVDlimitterの値でスキームを選択できます。デフォルトでvan Albadaにしています。
グリッド生成
グリッドも前回から修正しました。境界面を指定するためのインデックスがコードに直接書き込んでありますが，
そこはご了承ください。
using Parameters
using HDF5
@with_kw struct Param{T1,T2}
Rout::T1=10.0
R1::T1 = 0.5
xr::T1 = 20.
ytop::T1 = 10.
x1::T1 = 8.0
x2::T1 = 10.0
ε::T1 = 1e-5
itermax::T2 = 100
nx::T2 = 140
ny::T2 = 561
outputname::String = "Grid/grid-cylinder-c"
end
function gridgen(para)
@unpack nx, ny, = para
x = zeros(nx, ny)
y = zeros(nx, ny)
xT = zeros(ny, nx)
yT = zeros(ny, nx)
# control function
P = zeros(nx, ny)
Q = zeros(nx, ny)
initialize(para, x, y)
xT .= x'
yT .= y'
output(para, x, y, xT, yT)
# solve elliptic equation
ellip(para, x, y, P, Q)
xT .= x'
yT .= y'
output(para, x, y, xT, yT)
end
function initialize(para, x, y)
@unpack nx, ny, R1, Rout, xr, ytop, x1, x2 = para
dx = 0.05
for j=1:div(ny,2)
x[1, j] = xr - dx*(j-1)
end
x[1,div(ny,2)+1] = 7.5
for j=div(ny,2)+2:ny
x[1,j] = 7.5 + dx*(j-div(ny,2)-31)
end
for j=241:281
for i=2:nx
θ = 0.5*pi*((j-241)/(40)) + 0.5*pi
R=Rout/1.001^(nx-i)
x[i,j] = R*cos(θ) + x1
y[i,j] = R*sin(θ)
x[i,562-j] = R*cos(-θ) + x1
y[i,562-j] = R*sin(-θ)
x[i,281] = R*cos(pi) + x1
y[i,281] = R*sin(pi)
end
end
for j=1:230
y[1,j] = 0
end
for j=231:281
θ = pi*(j-231)/50
x[1,j] = R1*cos(θ) + x1
y[1,j] = R1*sin(θ)
x[1,562-j] = R1*cos(-θ) + x1
y[1,562-j] = R1*sin(-θ)
end
x[1,281] = R1*cos(pi) + x1
y[1,281] = R1*sin(pi)
for j=332:ny
y[1,j]=0
end
#for j=224:230
# x[1,j] = x1+R1 +0.005*1.8^(231-j)
#end
#for j=332:338
# x[1,j] = x1+R1 +0.005*1.8^(j-331)
#end
#外周
for j=1:190
for i=2:nx
R = Rout
x[i, j] = xr- dx*(j-1)
y[i,j] = R/(nx-1)*(i-1)
end
end
for j=190:241
for i=2:nx-1
R = Rout
x[i, j] = xr- dx*(j-1)
y[i,j] = 0.5R +0.5R/(nx-1)*(i-1)
end
i=nx
R = Rout
x[i,j] = xr- dx*(j-1)
y[i,j] = R
end
for j= 321:372
for i=2:nx-1
R = Rout
x[i,j] = 7.5 +dx*(j-div(ny,2)-31)
y[i,j] = -0.5R/(nx-1)*(i-1) - 0.5R
end
i=nx
R = Rout
x[i,j] = 7.5 +dx*(j-div(ny,2)-31)
y[i,j] = -R
end
for j=373:ny
for i=2:nx
R = Rout
x[i,j] = 7.5 +dx*(j-div(ny,2)-31)
y[i,j] = -R/(nx-1)*(i-1)
end
end
end
function output(para, x, y, xT, yT)
@unpack nx, ny, outputname = para
filename = outputname*".h5"
h5open(filename, "w") do file
write(file, "nx", nx)
write(file, "ny", ny)
write(file, "x", x[1:nx-1,:])
write(file, "y", y[1:nx-1,:])
write(file, "xall", x[1:nx,:])
write(file, "yall", y[1:nx,:])
write(file, "xT", xT[:,1:nx-1])
write(file, "yT", yT[:,1:nx-1])
write(file, "xTall", xT[:,1:nx])
write(file, "yTall", yT[:,1:nx])
end
end
function ellip(para, x, y, P, Q)
@unpack nx, ny, itermax, ε = para
itermax3 = div(itermax, 3)
itermax6 = itermax3*2
println(x[1:3,1:3])
cons = 0.01
for i = 1:itermax
error = 0.0
if i >= itermax3
cons = 0.25
elseif i>=itermax6
cons = 1.0
end
for l = 2:ny-1
for k = 2:nx-1
xξ = 0.5(x[k+1, l] - x[k-1, l])
yξ = 0.5(y[k+1, l] - y[k-1, l])
xη = 0.5(x[k, l+1] - x[k, l-1])
yη = 0.5(y[k, l+1] - y[k, l-1])
α = xη^2 + yη^2
β = xξ*xη + yξ*yη
γ = xξ^2 + yξ^2
J = xξ*yη - xη*yξ
PP = P[k, l]
QQ = Q[k, l]
ag = 0.5/(α + γ)
sx = ( α*(x[k+1,l] + x[k-1,l]) -
0.5β*(x[k+1,l+1] - x[k-1,l+1] - x[k+1,l-1] + x[k-1,l-1]) +
γ*(x[k,l+1] + x[k,l-1]) + J^2*(PP*xξ + QQ*xη) )*ag - x[k,l]
sy = (α*(y[k+1,l] + y[k-1,l]) -
0.5β*(y[k+1,l+1] - y[k-1,l+1] - y[k+1,l-1] + y[k-1,l-1]) +
γ*(y[k,l+1] + y[k,l-1]) + J^2*(PP*yξ + QQ*yη) )*ag - y[k,l]
x[k,l] = x[k,l] + cons*sx
y[k,l] = y[k,l] + cons*sy
error += sx^2 + sy^2
end
end
if error < ε
println("iteration, error=", i, ", ", error)
break
elseif i==itermax
println("error= ",error)
end
if mod(i,2000) == 0
println("iteration=", i, ", error=", error)
end
end
end
para = Param( itermax=100000,)
@time gridgen(para)
2分くらいでグリッドができます。
流体のほうのコード
using Parameters
using HDF5
@with_kw struct Param{T1,T2}
dt::T1 = 0.01
# viscosity
nu::T1 = 0.01
Uext::T1 = 1.0
nstep::T2 = 400
alphaua::T1 = 0.5
alphava::T1 = 0.5
alphapa::T1 = 0.25
niter::T2 = 20
ipter::T2 = 5
ipter2::T2 = 5
start::T2 =1
TVDlimitter::T2 = 1 #van Albada
# mesh size
nx::T2 = 137
ny::T2 = 559
ndata::T2 = 50
gridfile::String = "Grid/grid-cylinder-c.h5"
outfile::String = "Data-cyl-Re100/uv"
end
function main(para)
@unpack nx, ny, nstep, niter, ndata, Uext, TVDlimitter = para
# auxiliary velocity
ua = ones(Float64, nx+2, ny+2)
ua2 = ones(Float64, nx+2, ny+2)
ua3 = ones(Float64, nx+2, ny+2)
utemp = zeros(Float64, nx+2, ny+2)
uold = ones(Float64, nx+2, ny+2)
apuva = zeros(Float64, nx+2, ny+2)
aeuva = zeros(Float64, nx+2, ny+2)
awuva = zeros(Float64, nx+2, ny+2)
anuva = zeros(Float64, nx+2, ny+2)
asuva = zeros(Float64, nx+2, ny+2)
bua = zeros(Float64, nx+2, ny+2)
va = zeros(Float64, nx+2, ny+2)
va2 = zeros(Float64, nx+2, ny+2)
va3 = zeros(Float64, nx+2, ny+2)
vtemp = zeros(Float64, nx+2, ny+2)
vold = zeros(Float64, nx+2, ny+2)
bva = zeros(Float64, nx+2, ny+2)
Ue = zeros(Float64, nx+2, ny+2)
Vn = zeros(Float64, nx+2, ny+2)
UE = zeros(Float64, nx+2, ny+2)
UW = zeros(Float64, nx+2, ny+2)
VN = zeros(Float64, nx+2, ny+2)
VS = zeros(Float64, nx+2, ny+2)
omega = zeros(Float64, nx+2, ny+2)
beta2 = zeros(Float64, nx+2, ny+2)
gam = zeros(Float64, nx+2, ny+2)
U = zeros(Float64, nx+2, ny+2)
V = zeros(Float64, nx+2, ny+2)
#dξ = zeros(Float64, nx+2, ny+2)
#dη = zeros(Float64, nx+2, ny+2)
# pressure, corrected pressure
pa = zeros(Float64, nx+2, ny+2)
pc = zeros(Float64, nx+2, ny+2)
pc2 = zeros(Float64, nx+2, ny+2)
ptemp = zeros(Float64, nx+2, ny+2)
appc = zeros(Float64, nx+2, ny+2)
aepc = zeros(Float64, nx+2, ny+2)
awpc = zeros(Float64, nx+2, ny+2)
anpc = zeros(Float64, nx+2, ny+2)
aspc = zeros(Float64, nx+2, ny+2)
bpc = zeros(Float64, nx+2, ny+2)
# coordinate
x = zeros(Float64, nx+2, ny+2)
y = zeros(Float64, nx+2, ny+2)
st = zeros(Float64, ny+2)
ξx = zeros(Float64, nx+2, ny+2)
ξy = zeros(Float64, nx+2, ny+2)
ηx = zeros(Float64, nx+2, ny+2)
ηy = zeros(Float64, nx+2, ny+2)
J = zeros(Float64, nx+2, ny+2)
Su = zeros(Float64, nx+2, ny+2)
Sv = zeros(Float64, nx+2, ny+2)
re⁺ = zeros(Float64, nx+2, ny+2)
re⁻ = zeros(Float64, nx+2, ny+2)
rn⁺ = zeros(Float64, nx+2, ny+2)
rn⁻ = zeros(Float64, nx+2, ny+2)
αe = zeros(Int64, nx+2, ny+2)
αn = zeros(Int64, nx+2, ny+2)
rest = 0.0
readgrid(para, x, y, st)
jacobian(para, x, y, ξx, ξy, ηx, ηy, J)
ψ = selectTVDlimitter(TVDlimitter)
for l = 1:nstep
for j=1:ny+2
for i=1:nx+2
uold[i,j] = ua[i,j]
vold[i,j] = va[i,j]
end
end
BC(para, ua, va, pc, pa, ξx, ξy, ηx, ηy, J, uold, vold, ua2, va2, pc2, l, st)
for n = 1:niter
BC(para, ua, va, pc, pa, ξx, ξy, ηx, ηy, J, uold, vold, ua2, va2, pc2, l, st)
makesource(para, bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J, U, V, Su, Sv, re⁺, re⁻, rn⁺, rn⁻, αe, αn, Ue, Vn, ψ)
for k = 1:5
solveUVA(para, ua, va, apuva, aeuva, awuva, anuva, asuva, bua, utemp, bva, vtemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
end
RhieChowInterpolation(para, ua, va, pa, apuva, U, V, Ue, Vn, ξx, ξy, ηx, ηy, J)
#println("ua, va=", ua[2,2],", ", va[2,2])
solvePC(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, ptemp, Ue, Vn, ξx, ξy, ηx, ηy, J, st)
update(para, ua, va, pc, pa, apuva, ua2, va2, ξx, ξy, ηx, ηy, J)
solvePC2(para, ua, ua2, va, va2, pc, pc2, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, UE, UW, VN, VS, ptemp, ξx, ξy, ηx, ηy, U, V, J , st)
update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva, ξx, ξy, ηx, ηy, J, utemp, vtemp)
rest = converg(para, ua, va, ua2, va2)
divu = divergence(para, ua, va, U, V, ξx, ξy, ηx, ηy, J)
if rest <=1e-7 && n >=2
if l%20 == 0
println(l, " th step, ", n, " th iteration, u residual=", rest, ", div=", divu)
end
break
elseif n==niter
println(l, " th step, ", n, " th iteration, u residual=", rest, ", div=", divu)
end
end
if l%ndata == 0
println("residual=", rest)
vor(para, omega, ua, va, ξx, ξy, ηx, ηy, st)
output(para, ua, va, pa, l, omega)
end
end
#streamfunction(para, u, v, psi)
end
function main(para, restartFile)
@unpack nx, ny, nstep, niter, ndata, TVDlimitter = para
# auxiliary velocity
ua = zeros(Float64, nx+2, ny+2)
ua2 = zeros(Float64, nx+2, ny+2)
ua3 = zeros(Float64, nx+2, ny+2)
utemp = zeros(Float64, nx+2, ny+2)
uold = zeros(Float64, nx+2, ny+2)
apuva = zeros(Float64, nx+2, ny+2)
aeuva = zeros(Float64, nx+2, ny+2)
awuva = zeros(Float64, nx+2, ny+2)
anuva = zeros(Float64, nx+2, ny+2)
asuva = zeros(Float64, nx+2, ny+2)
bua = zeros(Float64, nx+2, ny+2)
va = zeros(Float64, nx+2, ny+2)
va2 = zeros(Float64, nx+2, ny+2)
va3 = zeros(Float64, nx+2, ny+2)
vtemp = zeros(Float64, nx+2, ny+2)
vold = zeros(Float64, nx+2, ny+2)
bva = zeros(Float64, nx+2, ny+2)
Ue = zeros(Float64, nx+2, ny+2)
Vn = zeros(Float64, nx+2, ny+2)
omega = zeros(Float64, nx+2, ny+2)
UE = zeros(Float64, nx+2, ny+2)
UW = zeros(Float64, nx+2, ny+2)
VN = zeros(Float64, nx+2, ny+2)
VS = zeros(Float64, nx+2, ny+2)
beta2 = zeros(Float64, nx+2, ny+2)
gam = zeros(Float64, nx+2, ny+2)
U = zeros(Float64, nx+2, ny+2)
V = zeros(Float64, nx+2, ny+2)
dξ = zeros(Float64, nx+2, ny+2)
dη = zeros(Float64, nx+2, ny+2)
# pressure, corrected pressure
pa = zeros(Float64, nx+2, ny+2)
pc = zeros(Float64, nx+2, ny+2)
pc2 = zeros(Float64, nx+2, ny+2)
ptemp = zeros(Float64, nx+2, ny+2)
appc = zeros(Float64, nx+2, ny+2)
aepc = zeros(Float64, nx+2, ny+2)
awpc = zeros(Float64, nx+2, ny+2)
anpc = zeros(Float64, nx+2, ny+2)
aspc = zeros(Float64, nx+2, ny+2)
bpc = zeros(Float64, nx+2, ny+2)
# coordinate
x = zeros(Float64, nx+2, ny+2)
y = zeros(Float64, nx+2, ny+2)
st = zeros(Float64, ny+2)
ξx = zeros(Float64, nx+2, ny+2)
ξy = zeros(Float64, nx+2, ny+2)
ηx = zeros(Float64, nx+2, ny+2)
ηy = zeros(Float64, nx+2, ny+2)
J = zeros(Float64, nx+2, ny+2)
Su = zeros(Float64, nx+2, ny+2)
Sv = zeros(Float64, nx+2, ny+2)
re⁺ = zeros(Float64, nx+2, ny+2)
re⁻ = zeros(Float64, nx+2, ny+2)
rn⁺ = zeros(Float64, nx+2, ny+2)
rn⁻ = zeros(Float64, nx+2, ny+2)
αe = zeros(Int64, nx+2, ny+2)
αn = zeros(Int64, nx+2, ny+2)
rest = 0.0
readgrid(para, x, y, st)
jacobian(para, x, y, ξx, ξy, ηx, ηy, J)
ψ = selectTVDlimitter(TVDlimitter)
nstart = restart(para, restartFile, ua, va, pa, )
println("nstart=",nstart)
for l = nstart+1:nstart+nstep
for j=1:ny+2
for i=1:nx+2
uold[i,j] = ua[i,j]
vold[i,j] = va[i,j]
end
end
BC(para, ua, va, pc, pa, ξx, ξy, ηx, ηy, J, uold, vold, ua2, va2, pc2, l, st)
for n = 1:niter
BC(para, ua, va, pc, pa, ξx, ξy, ηx, ηy, J, uold, vold, ua2, va2, pc2, l, st)
makesource(para, bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J, U, V, Su, Sv, re⁺, re⁻, rn⁺, rn⁻, αe, αn, Ue, Vn, ψ)
for k = 1:5
solveUVA(para, ua, va, apuva, aeuva, awuva, anuva, asuva, bua, utemp, bva, vtemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
end
RhieChowInterpolation(para, ua, va, pa, apuva, U, V, Ue, Vn, ξx, ξy, ηx, ηy, J)
#println("ua, va=", ua[2,2],", ", va[2,2])
solvePC(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, ptemp, Ue, Vn, ξx, ξy, ηx, ηy, J, st)
update(para, ua, va, pc, pa, apuva, ua2, va2, ξx, ξy, ηx, ηy, J)
solvePC2(para, ua, ua2, va, va2, pc, pc2, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, UE, UW, VN, VS, ptemp, ξx, ξy, ηx, ηy, U, V, J , st)
update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva, ξx, ξy, ηx, ηy, J, utemp, vtemp)
rest = converg(para, ua, va, ua2, va2)
divu = divergence(para, ua, va, U, V, ξx, ξy, ηx, ηy, J)
if rest <=1e-7 && n >=2
#if divu <=1e-7 && n>=2
if l%20 == 0
println(l, " th step, ", n, " th iteration, residual=", rest, ", div=", divu)
end
break
elseif n==niter
println(l, " th step, ", n, " th iteration, u residual=", rest, ", div=", divu)
end
end
if l%ndata == 0
println("residual=", rest)
vor(para, omega, ua, va, ξx, ξy, ηx, ηy, st)
output(para, ua, va, pa, l, omega)
end
end
#streamfunction(para, u, v, psi)
end
function BC(para, ua, va, pc, pa, ξx, ξy, ηx, ηy, J, uold, vold, ua2, va2, pc2, l, st)
@unpack nx, ny, Uext, dt, start = para
# cylinder
for j=232:330
ua[1, j] = 0.0
va[1, j] = 0.0
ua2[1,j] = 0.0
va2[1,j] = 0.0
pc[1, j] = pc[2,j]
pc2[1, j] = pc2[2,j]
pa[1, j] = pa[2,j]
end
# j=231は円柱表面だが，外挿で圧力を設定する。
for j=1:231
ua[1, j] = 0.5*(ua[2,j]-st[j]*(ua[3,j]-ua[2,j]) + ua[2,ny+3-j]-st[ny+3-j]*(ua[3,ny+3-j]-ua[2,ny+3-j]))
ua[1, ny+3-j] = ua[1,j]
va[1, j] = 0.5*(va[2,j]-st[j]*(va[3,j]-va[2,j]) + va[2,ny+3-j]-st[ny+3-j]*(va[3,ny+3-j]-va[2,ny+3-j]))
va[1, ny+3-j] = va[1,j]
ua2[1,j] = 0.5*(ua2[2,j]-st[j]*(ua2[3,j]-ua2[2,j]) + ua2[2,ny+3-j]-st[ny+3-j]*(ua2[3,ny+3-j]-ua2[2,ny+3-j]))
ua2[1,ny+3-j] = ua2[1,j]
va2[1,j] = 0.5*(va2[2,j]-st[j]*(va2[3,j]-va2[2,j]) + va2[2,ny+3-j]-st[ny+3-j]*(va2[3,ny+3-j]-va2[2,ny+3-j]))
va2[1,ny+3-j] = va2[1,j]
pc[1, j] = 0.5*(pc[2,j]-st[j]*(pc[3,j]-pc[2,j]) + pc[2,ny+3-j]-st[ny+3-j]*(pc[3,ny+3-j]-pc[2,ny+3-j]))
pc[1, ny+3-j] = pc[1,j]
pc2[1, j] = 0.5*(pc2[2,j]-st[j]*(pc2[3,j]-pc2[2,j]) + pc2[2,ny+3-j]-st[ny+3-j]*(pc2[3,ny+3-j]-pc2[2,ny+3-j]))
pc2[1, ny+3-j] = pc2[1,j]
pa[1, j] = 0.5*(pa[2,j]-st[j]*(pa[3,j]-pa[2,j]) + pa[2,ny+3-j]-st[ny+3-j]*(pa[3,ny+3-j]-pa[2,ny+3-j]))
pa[1, ny+3-j] = pa[1,j]
end
#
# outlet
for i=1:nx+2
ua[i, 1] = ua[i, 2]
va[i, 1] = va[i, 2]
ua[i, ny+2] = ua[i, ny+1]
va[i, ny+2] = va[i, ny+1]
ua2[i, 1] = ua2[i, 2]
va2[i, 1] = va2[i, 2]
ua2[i, ny+2] = ua2[i, ny+1]
va2[i, ny+2] = va2[i, ny+1]
pa[i, 1] = 0
pc[i, 1] = 0
pc2[i, 1] = 0
pa[i, ny+2] = 0
pc[i, ny+2] = 0
pc2[i, ny+2] = 0
end
# inlet
Ubound = min(Uext*l/start, Uext)
for j=1:ny+2
ua[nx+2,j] = Ubound
va[nx+2,j] = 0.0
ua2[nx+2,j] = Ubound
va2[nx+2,j] = 0.0
pa[nx+2,j] = pa[nx+1, j]
pc[nx+2,j] = pc[nx+1, j]
pc2[nx+2,j] = pc2[nx+1, j]
end
end
function readgrid(para, x, y, st)
@unpack gridfile, nx, ny = para
file = h5open(gridfile, "r")
x .= read(file, "x")
y .= read(file, "y")
close(file)
for j=1:ny+2
st[j] = sqrt(((x[1,j]-x[2,j])^2 + (y[1,j]-y[2,j])^2)/((x[2,j]-x[3,j])^2 + (y[2,j]-y[3,j])^2))
end
end
function jacobian(para, x, y, ξx, ξy, ηx, ηy, J)
@unpack nx, ny = para
xξ = zeros(Float64, nx+2, ny+2)
yξ = zeros(Float64, nx+2, ny+2)
xη = zeros(Float64, nx+2, ny+2)
yη = zeros(Float64, nx+2, ny+2)
for j = 1:ny+2
for i = 2:nx+1
xξ[i,j] = 0.5(x[i+1,j] - x[i-1,j])
yξ[i,j] = 0.5(y[i+1,j] - y[i-1,j])
end
end
for j=1:ny+2
if (231<=j<=331)
xξ[1,j] = -(3*x[1,j] - 4*x[2,j] + x[3,j])*0.5
yξ[1,j] = -(3*y[1,j] - 4*y[2,j] + y[3,j])*0.5
else
xξ[1,j] = 0.5(x[2,j] - x[2,ny+3-j])
yξ[1,j] = 0.5(y[2,j] - y[2,ny+3-j])
end
xξ[nx+2,j] = (3*x[nx+2,j] - 4*x[nx+1,j] + x[nx,j])*0.5
yξ[nx+2,j] = (3*y[nx+2,j] - 4*y[nx+1,j] + y[nx,j])*0.5
end
for j = 2:ny+1
for i = 1:nx+2
xη[i,j] = 0.5(x[i,j+1] - x[i,j-1])
yη[i,j] = 0.5(y[i,j+1] - y[i,j-1])
end
end
for i=1:nx+2
xη[i,1] = -(3*x[i,1] - 4*x[i,2] + x[i,3])*0.5
yη[i,1] = -(3*y[i,1] - 4*y[i,2] + y[i,3])*0.5
xη[i, ny+2] = (3*x[i,ny+2] - 4*x[i,ny+1] + x[i,ny])*0.5
yη[i, ny+2] = (3*y[i,ny+2] - 4*y[i,ny+1] + y[i,ny])*0.5
end
for j = 1:ny+2
for i = 1:nx+2
J[i,j] = 1.0/(xξ[i,j]*yη[i,j] - xη[i,j]*yξ[i,j])
ξx[i,j] = J[i,j]*yη[i,j]
ηx[i,j] = -J[i,j]*yξ[i,j]
ξy[i,j] = -J[i,j]*xη[i,j]
ηy[i,j] = J[i,j]*xξ[i,j]
end
end
println("metrics ", ξx[1:3,1:2], ξy[1:3,1:2], ηx[1:3,1:2], ηy[1:3,1:2], J[1:3,1:2])
println("negative volume, ", findall(x->(x<=0), J), ", ",J[findall(x->(x<=0), J)])
end
function makesource(para, bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J, U, V, Su, Sv, re⁺, re⁻, rn⁺, rn⁻,αe, αn, Ue, Vn, ψ)
@unpack nx, ny, dt, nu, = para
for j=1:ny+2
for i=1:nx+2
beta2[i,j] = nu/J[i,j]*(ξx[i,j]*ηx[i,j] + ξy[i,j]*ηy[i,j])
end
end
# bua
for j=2:ny+1
for i=2:nx+1
bua[i,j] = 0.125( (uold[i,j+1] - uold[i,j-1])*(beta2[i+1,j] - beta2[i-1,j]) + (uold[i+1,j] - uold[i-1,j])*(beta2[i,j+1] - beta2[i,j-1]) + uold[i+1,j+1]*(2beta2[i,j] + beta2[i+1,j] +beta2[i,j+1]) - uold[i-1,j+1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j+1]) - uold[i+1,j-1]*(2beta2[i,j] + beta2[i+1,j] + beta2[i,j-1]) + uold[i-1,j-1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j-1]))
bua[i,j] += -0.5(ξx[i+1,j]*pa[i+1,j]/J[i+1,j] - ξx[i-1,j]*pa[i-1,j]/J[i-1,j] + ηx[i,j+1]*pa[i,j+1]/J[i,j+1] - ηx[i,j-1]*pa[i,j-1]/J[i,j-1]) + uold[i,j]/J[i,j]/dt
end
end
#bva
for j=2:ny+1
for i=2:nx+1
bva[i,j] = 0.125( (vold[i,j+1] - vold[i,j-1])*(beta2[i+1,j] - beta2[i-1,j]) + (vold[i+1,j] - vold[i-1,j])*(beta2[i,j+1] - beta2[i,j-1]) + vold[i+1,j+1]*(2beta2[i,j] + beta2[i+1,j] +beta2[i,j+1]) - vold[i-1,j+1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j+1]) - vold[i+1,j-1]*(2beta2[i,j] + beta2[i+1,j] + beta2[i,j-1]) + vold[i-1,j-1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j-1]))
bva[i,j] += -0.5(ξy[i+1,j]*pa[i+1,j]/J[i+1,j] - ξy[i-1,j]*pa[i-1,j]/J[i-1,j] + ηy[i,j+1]*pa[i,j+1]/J[i,j+1] - ηy[i,j-1]*pa[i,j-1]/J[i,j-1]) + vold[i,j]/J[i,j]/dt
end
end
makeCotravariant(para, uold, vold, U, V, ξx, ξy, ηx, ηy)
for j=2:ny
for i=2:nx
Ue[i,j] = 0.5*(U[i+1,j]/J[i+1,j] + U[i,j]/J[i,j])
Vn[i,j] = 0.5*(V[i,j+1]/J[i,j+1] + V[i,j]/J[i,j])
αe[i,j] = 1 - signbit(Ue[i,j]) #0.5*(sign(Ue[i,j]) + 1.)
αn[i,j] = 1 - signbit(Vn[i,j]) #0.5*(sign(Vn[i,j]) + 1.)
end
end
# Su
for j=2:ny
for i=2:nx
re⁺[i,j] = (uold[i,j]-uold[i-1,j])/(uold[i+1,j]-uold[i,j] + 1e-10)
re⁻[i,j] = (uold[i+2,j]-uold[i+1,j])/(uold[i+1,j]-uold[i,j]+ 1e-10)
rn⁺[i,j] = (uold[i,j]-uold[i,j-1])/(uold[i,j+1]-uold[i,j]+ 1e-10)
rn⁻[i,j] = (uold[i,j+2]-uold[i,j+1])/(uold[i,j+1]-uold[i,j]+ 1e-10)
end
end
for j=3:ny
for i=3:nx
Su[i,j] = 0.5Ue[i,j]*((1-αe[i,j])*ψ(re⁻[i,j]) - αe[i,j]*ψ(re⁺[i,j]))*(uold[i+1,j] - uold[i,j]) -
0.5Ue[i-1,j]*((1-αe[i-1,j])*ψ(re⁻[i-1,j]) - αe[i-1,j]*ψ(re⁺[i-1,j]))*(uold[i,j] - uold[i-1,j]) +
0.5Vn[i,j]*((1-αn[i,j])*ψ(rn⁻[i,j]) - αn[i,j]*ψ(rn⁺[i,j]))*(uold[i,j+1] - uold[i,j]) -
0.5Vn[i,j-1]*((1-αn[i,j-1])*ψ(rn⁻[i,j-1]) - αn[i,j-1]*ψ(rn⁺[i,j-1]))*(uold[i,j] - uold[i,j-1])
end
end
# Sv
for j=2:ny
for i=2:nx
re⁺[i,j] = (vold[i,j]-vold[i-1,j])/(vold[i+1,j]-vold[i,j]+ 1e-10)
re⁻[i,j] = (vold[i+2,j]-vold[i+1,j])/(vold[i+1,j]-vold[i,j]+ 1e-10)
rn⁺[i,j] = (vold[i,j]-vold[i,j-1])/(vold[i,j+1]-vold[i,j]+ 1e-10)
rn⁻[i,j] = (vold[i,j+2]-vold[i,j+1])/(vold[i,j+1]-vold[i,j]+ 1e-10)
end
end
for j=3:ny
for i=3:nx
Sv[i,j] = 0.5Ue[i,j]*((1-αe[i,j])*ψ(re⁻[i,j]) - αe[i,j]*ψ(re⁺[i,j]))*(vold[i+1,j] - vold[i,j]) -
0.5Ue[i-1,j]*((1-αe[i-1,j])*ψ(re⁻[i-1,j]) - αe[i-1,j]*ψ(re⁺[i-1,j]))*(vold[i,j] - vold[i-1,j]) +
0.5Vn[i,j]*((1-αn[i,j])*ψ(rn⁻[i,j]) - αn[i,j]*ψ(rn⁺[i,j]))*(vold[i,j+1] - vold[i,j]) -
0.5Vn[i,j-1]*((1-αn[i,j-1])*ψ(rn⁻[i,j-1]) - αn[i,j-1]*ψ(rn⁺[i,j-1]))*(vold[i,j] - vold[i,j-1])
end
end
for j=3:ny
for i=3:nx
bua[i,j] += Su[i,j]
bva[i,j] += Sv[i,j]
end
end
end
function selectTVDlimitter(TVDlimitter)
# van Albada
function vanAlbada(r)::Float64
return (r+r^2)/(1.0+r^2)
end
function vanLeer(r)
return (r + abs(r))/(1.0 + r)
end
function QUICK(r)
return max(0.0, min(2*r, (3+r)*0.25, 2.0))
end
function superbee(r)
return max(0.0, min(2*r, 1.0), min(r,2.0))
end
if TVDlimitter == 1
return vanAlbada
elseif TVDlimitter == 2
return vanLeer
elseif TVDlimitter == 3
return QUICK
elseif TVDlimitter == 4
return superbee
end
end
function makeCotravariant(para, ua, va, U, V, ξx, ξy, ηx, ηy)
@unpack nx, ny, = para
for j=1:ny+2
for i=1:nx+2
U[i,j] = ξx[i,j]*ua[i,j] + ξy[i,j]*va[i,j]
V[i,j] = ηx[i,j]*ua[i,j] + ηy[i,j]*va[i,j]
end
end
end
function solveUVA(para, ua, va, apuva, aeuva, awuva, anuva, asuva, bua, utemp, bva, vtemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
@unpack nx, ny, dt, nu, alphaua, alphava, niter, ipter, Uext = para
makeCotravariant(para, ua, va, U, V, ξx, ξy, ηx, ηy)
for j=1:ny+2
for i=1:nx+2
if i==1
awuva[i,j] = 0
elseif i==2 && (231<=j<=331)
awuva[i,j] = nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i-1,j]^2 + ξy[i-1,j]^2)/J[i-1,j])
elseif i==2 && (j<=230 || 332<=j)
awuva[i,j] = 0.5nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i-1,j]^2 + ξy[i-1,j]^2)/J[i-1,j]) + max(0.5*(U[i,j]/J[i,j] + U[i-1,j]/J[i-1,j]), 0.0)
else
awuva[i,j] = 0.5nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i-1,j]^2 + ξy[i-1,j]^2)/J[i-1,j]) + max(0.5*(U[i,j]/J[i,j] + U[i-1,j]/J[i-1,j]), 0.0)
end
if i==nx+2
aeuva[i,j] = 0
elseif i==nx+1
aeuva[i,j] = 0.5nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i+1,j]^2 + ξy[i+1,j]^2)/J[i+1,j]) + max(-0.5*(U[i+1,j]/J[i+1,j] + U[i,j]/J[i,j]), 0.0)
else
aeuva[i,j] = 0.5nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i+1,j]^2 + ξy[i+1,j]^2)/J[i+1,j]) + max(-0.5*(U[i+1,j]/J[i+1,j] + U[i,j]/J[i,j]), 0.0)
end
if j==1
asuva[i,j] = 0
elseif j==2
asuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,j-1]^2 + ηy[i,j-1]^2)/J[i,j-1]) + max(0.5*(V[i,j]/J[i,j] + V[i,j-1]/J[i,j-1]), 0.0)
else
asuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,j-1]^2 + ηy[i,j-1]^2)/J[i,j-1]) + max(0.5*(V[i,j]/J[i,j] + V[i,j-1]/J[i,j-1]), 0.0)
end
if j==ny+2
anuva[i,j] = 0
elseif j == ny+1
anuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,j+1]^2 + ηy[i,j+1]^2)/J[i,j+1]) + max(-0.5*(V[i,j+1]/J[i,j+1] + V[i,j]/J[i,j]), 0.0)
else
anuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,j+1]^2 + ηy[i,j+1]^2)/J[i,j+1]) + max(-0.5*(V[i,j+1]/J[i,j+1] + V[i,j]/J[i,j]), 0.0)
end
apuva[i,j] = aeuva[i,j] + awuva[i,j] + anuva[i,j] + asuva[i,j] + 1/dt/J[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
utemp[i,j] = ua[i,j]*(1.0-alphaua) + alphaua*(
aeuva[i,j]*ua[i+1,j] + awuva[i,j]*ua[i-1,j] + anuva[i,j]*ua[i,j+1] + asuva[i,j]*ua[i,j-1] + bua[i,j])/apuva[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
ua[i,j] = utemp[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
vtemp[i,j] = va[i,j] * (1.0 - alphava) + alphava * (
aeuva[i,j]*va[i+1,j] + awuva[i,j]*va[i-1,j] + anuva[i,j]*va[i,j+1] + asuva[i,j]*va[i,j-1] + bva[i,j])/apuva[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
va[i,j] = vtemp[i,j]
end
end
end
function RhieChowInterpolation(para, ua, va, pa, apuva, U, V, Ue, Vn, ξx, ξy, ηx, ηy, J)
@unpack nx, ny, = para
makeCotravariant(para, ua, va, U, V, ξx, ξy, ηx, ηy)
# ue: interface velocity
for j=2:ny+1
for i=2:nx
#Ue[i,j] = 0.5(U[i+1,j]+U[i,j]) + 0.5(dξ[i,j]+dξ[i+1,j])*(pa[i,j]-pa[i+1,j]) -0.25dξ[i,j]*(pa[i-1,j]-pa[i+1,j]) -0.25*dξ[i+1,j]*(pa[i,j]-pa[i+2,j])
Ue[i,j] = 0.5(U[i+1,j]+U[i,j]) -0.5((ξx[i,j]*ξx[i+1,j]/J[i+1,j]+ξy[i,j]*ξy[i+1,j]/J[i+1,j])/apuva[i,j] + (ξx[i+1,j]*ξx[i,j]/J[i,j]+ξy[i+1,j]*ξy[i,j]/J[i,j])/apuva[i+1,j])*(pa[i+1,j]-pa[i,j]) +
0.25((ξx[i,j]*ξx[i+1,j]/J[i+1,j]+ξy[i,j]*ξy[i+1,j]/J[i+1,j])*pa[i+1,j] -(ξx[i,j]*ξx[i-1,j]/J[i-1,j]+ξy[i,j]*ξy[i-1,j]/J[i-1,j])*pa[i-1,j] )/apuva[i,j] +
0.25((ξx[i+1,j]*ξx[i+2,j]/J[i+2,j]+ξy[i+1,j]*ξy[i+2,j]/J[i+2,j])*pa[i+2,j] -(ξx[i+1,j]*ξx[i,j]/J[i,j]+ξy[i+1,j]*ξy[i,j]/J[i,j])*pa[i,j] )/apuva[i+1,j]
end
end
for j=2:ny+1
Ue[1,j] = U[1,j]
Ue[nx+1,j] = U[nx+2,j]
end
for j=2:ny
for i=2:nx+1
#Vn[i,j] = 0.5(V[i,j+1]+V[i,j]) + 0.5(dη[i,j]+dη[i,j+1])*(pa[i,j]-pa[i,j+1]) -0.25dη[i,j]*(pa[i,j-1]-pa[i,j+1]) -0.25dη[i,j+1]*(pa[i,j]-pa[i,j+2])
Vn[i,j] = 0.5(V[i,j+1]+V[i,j]) -0.5((ηx[i,j]*ηx[i,j+1]/J[i,j+1]+ηy[i,j]*ηy[i,j+1]/J[i,j+1])/apuva[i,j] + (ηx[i,j+1]*ηx[i,j]/J[i,j]+ηy[i,j+1]*ηy[i,j]/J[i,j])/apuva[i,j+1] )*(pa[i,j+1]-pa[i,j]) +
0.25((ηx[i,j]*ηx[i,j+1]/J[i,j+1]+ηy[i,j]*ηy[i,j+1]/J[i,j+1])*pa[i,j+1] - (ηx[i,j]*ηx[i,j-1]/J[i,j-1]+ηy[i,j]*ηy[i,j-1]/J[i,j-1])*pa[i,j-1])/apuva[i,j] +
0.25((ηx[i,j+1]*ηx[i,j+2]/J[i,j+2]+ηy[i,j+1]*ηy[i,j+2]/J[i,j+2])*pa[i,j+2] - (ηx[i,j+1]*ηx[i,j]/J[i,j]+ηy[i,j+1]*ηy[i,j]/J[i,j])*pa[i,j])/apuva[i,j+1]
end
end
for i=2:nx+1
Vn[i, ny+1] = V[i, ny+2]
Vn[i, 1] = V[i, 1]
end
end
function solvePC(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua,
bva, appc, aepc, awpc, anpc, aspc, bpc, ptemp, Ue, Vn, ξx, ξy, ηx, ηy, J, st)
@unpack nx, ny, nu, alphaua, alphava, alphapa, niter, ipter, nu = para
for j=2:ny+1
for i=2:nx+1
pc[i,j] = 0.0
end
end
for j=2:ny+1
for i=2:nx+1
if i == 2
awpc[i,j] =0
else
awpc[i,j] = 0.5((ξx[i-1,j]*ξx[i,j]/J[i,j]+ξy[i-1,j]*ξy[i,j]/J[i,j])/apuva[i-1,j] + (ξx[i,j]*ξx[i-1,j]/J[i-1,j]+ξy[i,j]*ξy[i-1,j]/J[i-1,j])/apuva[i,j])/(0.5(J[i-1,j]+J[i,j]))
end
if i == nx+1
aepc[i,j] = 0
else
aepc[i,j] = 0.5((ξx[i,j]*ξx[i+1,j]/J[i+1,j]+ξy[i,j]*ξy[i+1,j]/J[i+1,j])/apuva[i,j] + (ξx[i+1,j]*ξx[i,j]/J[i,j]+ξy[i+1,j]*ξy[i,j]/J[i,j])/apuva[i+1,j])/(0.5(J[i,j]+J[i+1,j]))
end
if j == 2
aspc[i,j] = 0
else
aspc[i,j] = 0.5((ηx[i,j-1]*ηx[i,j]/J[i,j]+ηy[i,j-1]*ηy[i,j]/J[i,j])/apuva[i,j-1] + (ηx[i,j]*ηx[i,j-1]/J[i,j-1]+ηy[i,j]*ηy[i,j-1]/J[i,j-1])/apuva[i,j] )/(0.5(J[i,j-1]+J[i,j]))
end
if j == ny+1
anpc[i,j] = 0
else
anpc[i,j] = 0.5((ηx[i,j]*ηx[i,j+1]/J[i,j+1]+ηy[i,j]*ηy[i,j+1]/J[i,j+1])/apuva[i,j] + (ηx[i,j+1]*ηx[i,j]/J[i,j]+ηy[i,j+1]*ηy[i,j]/J[i,j])/apuva[i,j+1] )/(0.5(J[i,j]+J[i,j+1]))
end
appc[i,j] = aepc[i,j] + awpc[i,j] + anpc[i,j] + aspc[i,j]
bpc[i,j] = -(Ue[i,j]/(0.5(J[i,j]+J[i+1,j])) - Ue[i-1,j]/(0.5(J[i,j]+J[i-1,j]))) - (Vn[i,j]/(0.5(J[i,j]+J[i,j+1])) - Vn[i,j-1]/(0.5(J[i,j]+J[i,j-1])))
end
end
# SOR法で圧力の解を求めてる
for ip=1:ipter
for j=2:ny+1
for i=2:nx+1
ptemp[i,j] = pc[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
pc[i,j] = pc[i,j]*(1.0 - alphapa) + alphapa*(
aepc[i,j]*pc[i+1,j] + awpc[i,j]*pc[i-1,j] + anpc[i,j]*pc[i,j+1] + aspc[i,j]*pc[i,j-1] + bpc[i,j])/appc[i,j]
end
end
for j=232:330
pc[1, j] = pc[2,j]
end
for j=1:231
pc[1, j] = 0.5*(pc[2,j]-st[j]*(pc[3,j]-pc[2,j]) + pc[2,ny+3-j]-st[ny+3-j]*(pc[3,ny+3-j]-pc[2,ny+3-j]))
pc[1, ny+3-j] = pc[1,j]
end
# outlet
for i=1:nx+2
pc[i, 1] = 0
pc[i, ny+2] = 0
end
for j=1:231
pc[nx+2,j] = 0
end
for j=331:ny+2
pc[nx+2,j] = 0
end
# inlet
for j=232:330
pc[nx+2,j] = pc[nx+1, j]
end
rest=convergP(para, pc, ptemp)
if rest < 1e-6 && ip >=2
#println("converged rest PC ", rest, ", ip=", ip)
break
end
if ip ==ipter
println("rest PC ", rest, ", ip=", ip)
end
end
end
function update(para, ua, va, pc, pa, apuva, ua2, va2, ξx, ξy, ηx, ηy, J)
@unpack nx, ny, alphapa, = para
# u = u* + u'
for j=2:ny+1
for i=2:nx+1
ua2[i,j] = ua[i,j] -0.5(ξx[i+1,j]*pc[i+1,j]/J[i+1,j] - ξx[i-1,j]*pc[i-1,j]/J[i-1,j] + ηx[i,j+1]*pc[i,j+1]/J[i,j+1] - ηx[i,j-1]*pc[i,j-1]/J[i,j-1])/apuva[i,j]
va2[i,j] = va[i,j] -0.5(ξy[i+1,j]*pc[i+1,j]/J[i+1,j] - ξy[i-1,j]*pc[i-1,j]/J[i-1,j] + ηy[i,j+1]*pc[i,j+1]/J[i,j+1] - ηy[i,j-1]*pc[i,j-1]/J[i,j-1])/apuva[i,j]
end
end
#pref = pa[2,2]
#for j=1:ny+2
# for i=1:nx+2
# pa[i,j] -= pref
# end
#end
end
function solvePC2(para, ua, ua2, va, va2, pc, pc2, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, UE, UW, VN, VS, ptemp, ξx, ξy, ηx, ηy, U, V, J, st)
@unpack nx, ny, ipter2, alphapa, = para
makeCotravariant(para, ua2 .- ua, va2 .- va, U, V, ξx, ξy, ηx, ηy)
for j=2:ny+1
for i=2:nx+1
pc2[i,j] = 0.0
bpc[i,j] = 0.0
end
end
for j=2:ny+1
for i=2:nx+1
if i==nx+1
UE[i,j] = awuva[i+1,j]*(U[i,j])
else
UE[i,j] = aeuva[i+1,j]*(U[i+2,j]) +awuva[i+1,j]*(U[i,j]) +anuva[i+1,j]*(U[i+1,j+1]) +asuva[i+1,j]*(U[i+1,j-1])
end
if i==2
UW[i,j] = aeuva[i-1,j]*(U[i,j])
else
UW[i,j] = aeuva[i-1,j]*(U[i,j]) +awuva[i-1,j]*(U[i-2,j]) +anuva[i-1,j]*(U[i-1,j+1]) +asuva[i-1,j]*(U[i-1,j-1])
end
if j==ny+1
VN[i,j] = asuva[i,j+1]*(V[i,j])
else
VN[i,j] = aeuva[i,j+1]*(V[i+1,j+1]) +awuva[i,j+1]*(V[i-1,j+1]) +anuva[i,j+1]*(V[i,j+2]) +asuva[i,j+1]*(V[i,j])
end
if j==2
VS[i,j] = anuva[i,j-1]*(V[i,j])
else
VS[i,j] = aeuva[i,j-1]*(V[i+1,j-1]) +awuva[i,j-1]*(V[i-1,j-1]) +anuva[i,j-1]*(V[i,j]) +asuva[i,j-1]*(V[i,j-2])
end
end
end
for j=2:ny+1
for i=2:nx+1
bpc[i,j] = -0.5(UE[i,j]/apuva[i+1,j]/(0.5(J[i+1,j]+J[i,j])) - UW[i,j]/apuva[i-1,j]/(0.5(J[i,j]+J[i-1,j]))) -0.5(VN[i,j]/apuva[i,j+1]/(0.5(J[i,j]+J[i,j+1])) -VS[i,j]/apuva[i,j-1]/(0.5(J[i,j]+J[i,j-1])))
end
end
#println("ue, uw, vn, vs ", findmax(abs.(UE)), ",", findmax(abs.(UW)), ", ", findmax(abs.(VN)), ",", findmax(abs.(VS)))
#println("apuva ", findmax(apuva))
#println("bpc ",findmax(abs.(bpc)))
# SOR法で圧力の解を求めてる
for ip=1:ipter2
for j=2:ny+1
for i=2:nx+1
ptemp[i,j] = pc2[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
pc2[i,j] = pc2[i,j]*(1.0 - alphapa) + alphapa*(aepc[i,j]*pc2[i+1,j] + awpc[i,j]*pc2[i-1,j] + anpc[i,j]*pc2[i,j+1] + aspc[i,j]*pc2[i,j-1] + bpc[i,j])/appc[i,j]
end
end
for j=232:330
pc2[1, j] = pc2[2,j]
end
for j=1:231
pc2[1, j] = 0.5*(pc2[2,j]-st[j]*(pc2[3,j]-pc2[2,j]) + pc2[2,ny+3-j]-st[ny+3-j]*(pc2[3,ny+3-j]-pc2[2,ny+3-j]))
pc2[1, ny+3-j] = pc2[1, j]
end
# outlet
for i=1:nx+2
pc2[i, 1] = 0
pc2[i, ny+2] = 0
end
for j=1:231
pc2[nx+2,j] = 0
end
for j=331:ny+2
pc2[nx+2,j] = 0
end
# inlet
for j=232:330
pc2[nx+2,j] = pc2[nx+1, j]
end
rest=convergP(para, pc2, ptemp)
if rest < 1e-7 && ip>=2
#println("converged rest PC2 ", rest, ", ip=", ip)
break
end
if ip == ipter2
println("rest PC2 ", rest, ", ip=", ip)
end
end
end
function update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva, ξx, ξy, ηx, ηy, J, utemp, vtemp)
@unpack nx, ny, = para
# u** = u* + u'
for j=2:ny+1
for i=2:nx+1
ua3[i,j] = ua2[i,j] + (aeuva[i,j]*(ua2[i+1,j]-ua[i+1,j]) + awuva[i,j]*(ua2[i-1,j]-ua[i-1,j]) + anuva[i,j]*(ua2[i,j+1]-ua[i,j+1]) + asuva[i,j]*(ua2[i,j-1]-ua[i,j-1]) -
0.5(ξx[i+1,j]*pc2[i+1,j]/J[i+1,j] - ξx[i-1,j]*pc2[i-1,j]/J[i-1,j] + ηx[i,j+1]*pc2[i,j+1]/J[i,j+1] - ηx[i,j-1]*pc2[i,j-1]/J[i,j-1]))/apuva[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
ua[i,j] = ua3[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
va3[i,j] = va2[i,j] + (aeuva[i,j]*(va2[i+1,j]-va[i+1,j]) + awuva[i,j]*(va2[i-1,j]-va[i-1,j]) + anuva[i,j]*(va2[i,j+1]-va[i,j+1]) + asuva[i,j]*(va2[i,j-1]-va[i,j-1]) -
0.5(ξy[i+1,j]*pc2[i+1,j]/J[i+1,j] - ξy[i-1,j]*pc2[i-1,j]/J[i-1,j] + ηy[i,j+1]*pc2[i,j+1]/J[i,j+1] - ηy[i,j-1]*pc2[i,j-1]/J[i,j-1]))/apuva[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
va[i,j] = va3[i,j]
end
end
#println("max ua3-ua2=",findmax(abs.(ua3[2:nx,2:ny+1]-ua2[2:nx,2:ny+1])))
#println("max va3-va2=",findmax(abs.(va3[2:nx+1,2:ny]-va2[2:nx+1,2:ny])))
for j=2:ny+1
for i=2:nx+1
pa[i,j] = pa[i,j] + pc[i,j] + pc2[i,j]
end
end
end
function converg(para, ua, va, utemp, vtemp)
@unpack nx, ny = para
restua = 0.0
restva = 0.0
for j=2:ny+1
for i=2:nx+1
restua = restua + abs(ua[i,j] - utemp[i,j])
end
end
for j=2:ny+1
for i=2:nx+1
restva = restva + abs(va[i,j] - vtemp[i,j])
end
end
rest = (restua + restva)/2.0
return rest/(nx*ny)
end
function convergP(para, pc, ptemp)
@unpack nx, ny = para
rest = 0.0
for j=2:ny+1
for i=2:nx+1
rest = rest + abs(pc[i,j] - ptemp[i,j])
end
end
return rest/(nx*ny)
end
function divergence(para, ua, va, U, V, ξx, ξy, ηx, ηy, J)
@unpack nx, ny = para
makeCotravariant(para, ua, va, U, V, ξx, ξy, ηx, ηy)
divu = 0.0
for j=2:ny+1
for i=2:nx+1
divu += abs((0.5*(U[i+1,j] + U[i,j])/(0.5(J[i,j]+J[i+1,j])) - 0.5*(U[i,j] + U[i-1,j])/(0.5(J[i-1,j]+J[i,j]))) + (0.5(V[i,j+1]+V[i,j])/(0.5(J[i,j]+J[i,j+1])) - 0.5(V[i,j]+V[i,j-1])/(0.5(J[i,j]+J[i,j-1]))))
end
end
return divu/(nx*ny)
end
function output(para, ua, va, pa, l, omega)
@unpack outfile, dt, nx, ny = para
#println("time=", time)
filename = outfile*string(l)*".h5"
h5open(filename, "w") do file
write(file, "U", ua)
write(file, "V", va)
write(file, "P", pa)
write(file, "omega", omega)
write(file, "n", l)
end
end
function restart(para, restartFile, ua, va, pa, )
@unpack nx, ny, = para
file = h5open(restartFile, "r")
ua .= read(file, "U")
va .= read(file, "V")
pa .= read(file, "P")
nstart = read(file, "n")
close(file)
return nstart
end
function vor(para, omega, ua, va, ξx, ξy, ηx, ηy, st)
@unpack nx, ny, = para
for j=2:ny+1
for i=2:nx+1
omega[i,j] = (ξx[i,j]*0.5*(va[i+1,j] -va[i-1,j]) + ηx[i,j]*0.5*(va[i,j+1] -va[i,j-1])) -
(ξy[i,j]*0.5*(ua[i+1,j] -ua[i-1,j]) + ηy[i,j]*0.5*(ua[i,j+1] -ua[i,j-1]))
end
end
for j=1:231
omega[1, j] = 0.5*(omega[2,j]-st[j]*(omega[3,j]-omega[2,j]) + omega[2,ny+3-j]-st[ny+3-j]*(omega[3,ny+3-j]-omega[2,ny+3-j]))
omega[1, ny+3-j] = omega[1,j]
end
for j=232:330
i=1
omega[i,j] = (ξx[i,j]*(va[i+1,j] -va[i,j]) + ηx[i,j]*0.5*(va[i,j+1] -va[i,j-1])) -
(ξy[i,j]*(ua[i+1,j] -ua[i,j]) + ηy[i,j]*0.5*(ua[i,j+1] -ua[i,j-1]))
end
end
実行は
para = Param()
@time main(para)
Mac miniで10分くらいで終わります。