概要
上記のコードにTVDスキームも実装しました。
前半は配列等のデータを直接引数渡しするコードで，同じ計算を配列等を構造体でまとめて，構造体を引数渡しするコードです。
構造体 struct を使えば，最初の配列の宣言の部分がすっきりするかとおもって作りましたが，計算が遅くなりました。配列を直接渡すほうが5.6秒くらいに対して，struct渡しのほうは8秒くらい掛かる感じです。
改善点があれば，教えて下さい。
せっかく書いたので，供養のためにおいておきます。
PISO法の説明は下記の後半にあります。
配列で全部データを渡すコード
using Parameters
using HDF5
@with_kw struct Param{T1,T2}
dt::T1 = 0.002
# Grashoff number
Ra::T1 = 1e8
Pr::T1 = 5.3
Gr::T1 = Ra/Pr
nu::T1 = 1.0/sqrt(Gr)
alpha::T1 = 1.0/Pr/sqrt(Gr)
# boundary temp
Th::T1 = 1.0
Tc::T1 = 0.0
nstep::T2 = 20000
niter::T2 = 15
ipter::T2 = 10
tol::T1 = 1e-5
alphapa::T1 = 0.5
TVDlimitter::T2 = 1 #van Albada
# mesh size
nx::T2 = 512
ny::T2 = 512
# length
Lx::T1 = 1.0
Ly::T1 = 1.0
dx::T1 = Lx/nx
dy::T1 = Ly/ny
ndata::T2 = 2000
outfile::String = "Piso-Ra1e8-TVD/uv"
end
function main(para)
@unpack nx, ny, nstep, niter, ndata, tol, 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)
uW = zeros(Float64, nx+2, ny+2)
vN = zeros(Float64, nx+2, ny+2)
vS = zeros(Float64, nx+2, ny+2)
uf = zeros(Float64, nx+2, ny+2)
vf = zeros(Float64, nx+2, ny+2)
# temperature
Ta = zeros(Float64, nx+2, ny+2)
Told = zeros(Float64, nx+2, ny+2)
apTa = zeros(Float64, nx+2, ny+2)
aeTa = zeros(Float64, nx+2, ny+2)
awTa = zeros(Float64, nx+2, ny+2)
anTa = zeros(Float64, nx+2, ny+2)
asTa = zeros(Float64, nx+2, ny+2)
bTa = zeros(Float64, nx+2, ny+2)
Ttemp = zeros(Float64, nx+2, ny+2)
# pressure, corrected pressure
pa = zeros(Float64, nx+2, ny+2)
pc = zeros(Float64, nx+2, ny+2)
ptemp = zeros(Float64, nx+2, ny+2)
pc2 = 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+1, ny+1)
y = zeros(Float64, nx+1, ny+1)
# stream function
psi = zeros(Float64, nx+2, ny+2)
Su = zeros(Float64, nx+2, ny+2)
Sv = zeros(Float64, nx+2, ny+2)
ST = 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)
ψ = selectTVDlimitter(TVDlimitter)
rest = 0.0
init(para, x, y)
for l = 1:nstep
BC(para, ua, va, ua2, va2, pa, pc, pc2, Ta)
for n = 1:niter
BC(para, ua, va, ua2, va2, pa, pc, pc2, Ta)
makesource(para, bua, bva, bTa, uold, vold, Told, pa, Ta, Su, Sv, ST, re⁺, re⁻, rn⁺, rn⁻, αe, αn, uf, vf, ψ)
solveUA(para, ua, va, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, uf)
solveVA(para, ua, va, Ta, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, vf)
solvePC(para, ua, va, pc, ptemp, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, uf, vf)
update(para, ua, va, pc, apuva, ua2, va2)
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)
update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva)
solveTA(para, ua, va, Ta, apTa, aeTa, awTa, anTa, asTa, bTa, Told, Ttemp)
rest = converg(para, ua, va, utemp, vtemp)
#println("iter, redisual ", n,", ",rest)
if rest <=tol && n>=2
if l%10 == 0
println(l, " th step, ", n, " th iteration, residual=", rest)
end
break
elseif n==niter
println(l, " th step, ", n, " th iteration, u residual=", rest,)
end
end
for j=1:ny+2
for i=1:nx+2
uold[i,j] = ua[i,j]
vold[i,j] = va[i,j]
Told[i,j] = Ta[i,j]
end
end
if l%ndata == 0
println("residual=", rest)
output(para, ua, va, pa, Ta, uold, vold, Told, x, y, l)
end
end
#streamfunction(para, u, v, psi)
end
function main(para, restartFile)
@unpack nx, ny, nstep, niter, ndata, tol, 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)
uW = zeros(Float64, nx+2, ny+2)
vN = zeros(Float64, nx+2, ny+2)
vS = zeros(Float64, nx+2, ny+2)
uf = zeros(Float64, nx+2, ny+2)
vf = zeros(Float64, nx+2, ny+2)
# temperature
Ta = zeros(Float64, nx+2, ny+2)
Told = zeros(Float64, nx+2, ny+2)
apTa = zeros(Float64, nx+2, ny+2)
aeTa = zeros(Float64, nx+2, ny+2)
awTa = zeros(Float64, nx+2, ny+2)
anTa = zeros(Float64, nx+2, ny+2)
asTa = zeros(Float64, nx+2, ny+2)
bTa = zeros(Float64, nx+2, ny+2)
Ttemp = zeros(Float64, nx+2, ny+2)
# pressure, corrected pressure
pa = zeros(Float64, nx+2, ny+2)
pc = zeros(Float64, nx+2, ny+2)
ptemp = zeros(Float64, nx+2, ny+2)
pc2 = 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+1, ny+1)
y = zeros(Float64, nx+1, ny+1)
# stream function
psi = zeros(Float64, nx+2, ny+2)
Su = zeros(Float64, nx+2, ny+2)
Sv = zeros(Float64, nx+2, ny+2)
ST = 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)
ψ = selectTVDlimitter(TVDlimitter)
rest = 0.0
nstart = restart(para, restartFile, x, y, ua, va, pa, Ta, uold, vold, Told)
println("nstart=",nstart)
for l = nstart+1:nstart+nstep
BC(para, ua, va, ua2, va2, pa, pc, pc2, Ta)
for n = 1:niter
BC(para, ua, va, ua2, va2, pa, pc, pc2, Ta)
makesource(para, bua, bva, bTa, uold, vold, Told, pa, Ta, Su, Sv, ST, re⁺, re⁻, rn⁺, rn⁻, αe, αn, uf, vf, ψ)
solveUA(para, ua, va, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, uf)
solveVA(para, ua, va, Ta, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, vf)
solvePC(para, ua, va, pc, ptemp, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, uf, vf)
update(para, ua, va, pc, apuva, ua2, va2)
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)
update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva)
solveTA(para, ua, va, Ta, apTa, aeTa, awTa, anTa, asTa, bTa, Told, Ttemp)
rest = converg(para, ua, va, utemp, vtemp)
#println("iter, redisual ", n,", ",rest)
if rest <=tol && n >=2
if l%10 == 0
println(l, " th step, ", n, " th iteration, residual=", rest)
end
break
elseif n==niter
println(l, " th step, ", n, " th iteration, u residual=", rest,)
end
#println(rest)
end
for j=2:ny+1
for i=2:nx+1
uold[i,j] = ua[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
vold[i,j] = va[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
Told[i,j] = Ta[i,j]
end
end
if l%ndata == 0
println("residual=", rest)
output(para, ua, va, pa, Ta, uold, vold, Told, x, y, l)
end
end
#streamfunction(para, u, v, psi)
end
function init(para, x, y)
@unpack nx, ny, dx, dy = para
for j in 2:ny+1
for i in 2:nx+1
x[i,j] = (i-2+0.5)*dx
end
end
for j in 2:ny+1
for i in 2:nx+1
y[i,j] = (j-2+0.5)*dy
end
end
end
function BC(para, ua, va, ua2, va2, pa, pc, pc2, Ta)
@unpack nx, ny, dy, Th, Tc = para
for i=1:nx+2
ua[i, ny+2] = 0.0
va[i, ny+2] = 0.0
ua[i, 1] = 0.0
va[i, 1] = 0.0
# ua2[i, ny+2]をupdate2でつかうが，updateではここの値は更新されないので，BCで与える。
ua2[i, ny+2] = 0.0
va2[i, ny+2] = 0.0
ua2[i, 1] = 0.0
va2[i, 1] = 0.0
pa[i, 1] = pa[i, 2]-Ta[i,1]*dy
pa[i, ny+2] = pa[i, ny+1]+Ta[i,ny+2]*dy
pc[i, 1] = pc[i, 2]
pc[i, ny+2] = pc[i, ny+1]
pc2[i, 1] = pc2[i, 2]
pc2[i, ny+2] = pc2[i, ny+1]
Ta[i, ny+2] = Tc
Ta[i, 1] = Th
end
for j=1:ny+2
ua[nx+2, j] = 0.0
va[nx+2, j] = 0.0
ua[1, j] = 0.0
va[1, j] = 0.0
ua2[nx+2, j] = 0.0
va2[nx+2, j] = 0.0
ua2[1, j] = 0.0
va2[1, j] = 0.0
pa[1, j] = pa[2, j]
pa[nx+2, j] = pa[nx+1, j]
pc[nx+2, j] = pc[nx+1, j]
pc[1, j] = pc[2, j]
pc2[nx+2, j] = pc2[nx+1, j]
pc2[1, j] = pc2[2, j]
Ta[nx+2, j] = Ta[nx+1, j]
Ta[1, j] = Ta[2,j]
end
end
function makesource(para, bua, bva, bTa, uold, vold, Told, pa, Ta, Su, Sv, ST, re⁺, re⁻, rn⁺, rn⁻, αe, αn, uf, vf, ψ)
@unpack nx, ny, dx, dy, dt = para
for j=2:ny+1
for i=2:nx+1
bua[i,j] = -0.5dy*(pa[i+1,j] - pa[i-1,j]) + dx*dy/dt*uold[i,j]
bva[i,j] = -0.5dx*(pa[i,j+1] - pa[i,j-1]) + dx*dy/dt*vold[i,j] + dx*dy*Told[i,j]
bTa[i,j] = dx*dy/dt*Told[i,j]
end
end
for j=2:ny
for i=2:nx
uf[i,j] = 0.5*(uold[i+1,j] + uold[i,j])
vf[i,j] = 0.5*(vold[i,j+1] + vold[i,j])
αe[i,j] = 1 - signbit(uf[i,j])
αn[i,j] = 1 - signbit(vf[i,j])
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] = dy*0.5uf[i,j]*((1-αe[i,j])*ψ(re⁻[i,j]) - αe[i,j]*ψ(re⁺[i,j]))*(uold[i+1,j] - uold[i,j]) -
dy*0.5uf[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]) +
dx*0.5vf[i,j]*((1-αn[i,j])*ψ(rn⁻[i,j]) - αn[i,j]*ψ(rn⁺[i,j]))*(uold[i,j+1] - uold[i,j]) -
dx*0.5vf[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] =
dy*0.5uf[i,j]*((1-αe[i,j])*ψ(re⁻[i,j]) - αe[i,j]*ψ(re⁺[i,j]))*(vold[i+1,j] - vold[i,j]) -
dy*0.5uf[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]) +
dx*0.5vf[i,j]*((1-αn[i,j])*ψ(rn⁻[i,j]) - αn[i,j]*ψ(rn⁺[i,j]))*(vold[i,j+1] - vold[i,j]) -
dx*0.5vf[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
# ST
for j=2:ny
for i=2:nx
re⁺[i,j] = (Told[i,j]-Told[i-1,j])/(Told[i+1,j]-Told[i,j]+ 1e-10)
re⁻[i,j] = (Told[i+2,j]-Told[i+1,j])/(Told[i+1,j]-Told[i,j]+ 1e-10)
rn⁺[i,j] = (Told[i,j]-Told[i,j-1])/(Told[i,j+1]-Told[i,j]+ 1e-10)
rn⁻[i,j] = (Told[i,j+2]-Told[i,j+1])/(Told[i,j+1]-Told[i,j]+ 1e-10)
end
end
for j=3:ny
for i=3:nx
ST[i,j] =
dy*0.5uf[i,j]*((1-αe[i,j])*ψ(re⁻[i,j]) - αe[i,j]*ψ(re⁺[i,j]))*(Told[i+1,j] - Told[i,j]) -
dy*0.5uf[i-1,j]*((1-αe[i-1,j])*ψ(re⁻[i-1,j]) - αe[i-1,j]*ψ(re⁺[i-1,j]))*(Told[i,j] - Told[i-1,j]) +
dx*0.5vf[i,j]*((1-αn[i,j])*ψ(rn⁻[i,j]) - αn[i,j]*ψ(rn⁺[i,j]))*(Told[i,j+1] - Told[i,j]) -
dx*0.5vf[i,j-1]*((1-αn[i,j-1])*ψ(rn⁻[i,j-1]) - αn[i,j-1]*ψ(rn⁺[i,j-1]))*(Told[i,j] - Told[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]
bTa[i,j] += ST[i,j]
end
end
end
function selectTVDlimitter(TVDlimitter)
# van Albada
function vanAlbada(r)
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 solveUA(para, ua, va, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, uf)
@unpack nx, ny, dx, dy, dt, nu, = para
for j=1:ny+2
for i=1:nx+2
if i==nx+2
aeuva[i,j] = 0.0
elseif i==nx+1
aeuva[i,j] = 2nu*dy/dx
else
aeuva[i,j] = nu*dy/dx + dy*max(-0.5*(ua[i+1,j] + ua[i,j]), 0.0)#- 0.25*dy*(ua[i+1, j] + ua[i,j])
end
if i==1
awuva[i,j] = 0.0
elseif i==2
awuva[i,j] = 2nu*dy/dx
else
awuva[i,j] = nu*dy/dx + dy*max(0.5*(ua[i,j] + ua[i-1,j]), 0.0)
end
if j==ny+2
anuva[i,j] = 0.0
elseif j == ny+1
anuva[i,j] = 2nu*dx/dy
else
anuva[i,j] = nu*dx/dy + dx*max(-0.5*(va[i,j+1] + va[i,j]), 0.0)#- 0.25dx*(va[i,j+1] + va[i,j])
end
if j==1
asuva[i,j] = 0.0
elseif j==2
asuva[i,j] = 2nu*dx/dy
else
asuva[i,j] = nu*dx/dy + dx*max(0.5*(va[i,j] + va[i,j-1]), 0.0)
end
apuva[i,j] = aeuva[i,j] + awuva[i,j] + anuva[i,j] + asuva[i,j] + dx*dy/dt
end
end
for j=2:ny+1
for i=2:nx+1
utemp[i,j] = (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=1:ny+1
for i=2:nx
uf[i,j] = 0.5(ua[i+1,j]+ua[i,j]) + 0.5(dy/apuva[i,j] + dy/apuva[i+1,j])*(pa[i,j]-pa[i+1,j]) - 0.25dy/apuva[i,j]*(pa[i-1,j]-pa[i+1,j]) -0.25dy/apuva[i+1,j]*(pa[i,j]-pa[i+2,j])
end
uf[1,j] = 0.0
uf[nx+1,j] = 0.0
end
end
function solveVA(para, ua, va, Ta, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, vf)
@unpack nx, ny, dx, dy, dt, nu, = para
for j=2:ny+1
for i=2:nx+1
vtemp[i,j] = (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
for j=2:ny
for i=1:nx+1
vf[i,j] = 0.5(va[i,j+1]+va[i,j]) + 0.5(dx/apuva[i,j] + dx/apuva[i,j+1])*(pa[i,j]-pa[i,j+1]) - 0.25dx/apuva[i,j]*(pa[i,j-1]-pa[i,j+1]) -0.25dx/apuva[i,j+1]*(pa[i,j]-pa[i,j+2])
end
end
for i=1:nx+1
vf[i,1] = 0.0
vf[i,ny+1] = 0.0
end
end
function solvePC(para, ua, va, pc, ptemp, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, uf, vf)
@unpack nx, ny, dx, dy, nu, ipter, alphapa = para
for j=2:ny+1
for i=2:nx+1
pc[i,j] = 0.0
bpc[i,j] = 0.0
end
end
for j=2:ny+1
for i=2:nx+1
if i == 2
awpc[i,j] = 0#0.5dy*(dy/apuva[i,j])
else
awpc[i,j] = 0.5dy*(dy/apuva[i-1,j]+dy/apuva[i,j])
end
if i == nx+1
aepc[i,j] = 0#0.5dy*(dy/apuva[i,j])
else
aepc[i,j] = 0.5dy*(dy/apuva[i,j]+dy/apuva[i+1,j])
end
if j == 2
aspc[i,j] = 0#0.5dx*(dx/apuva[i,j])
else
aspc[i,j] = 0.5dx*(dx/apuva[i,j-1]+dx/apuva[i,j])
end
if j == ny+1
anpc[i,j] = 0#0.5dx*(dx/apuva[i,j])
else
anpc[i,j] = 0.5dx*(dx/apuva[i,j]+dx/apuva[i,j+1])
end
appc[i,j] = aepc[i,j] + awpc[i,j] + anpc[i,j] + aspc[i,j]
bpc[i,j] = -dy*(uf[i,j] - uf[i-1,j]) - dx*(vf[i,j] - vf[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] = ptemp[i,j]*(1.0 - alphapa) + alphapa*(aepc[i,j]*ptemp[i+1,j] + awpc[i,j]*ptemp[i-1,j] + anpc[i,j]*ptemp[i,j+1] + aspc[i,j]*ptemp[i,j-1] + bpc[i,j])/appc[i,j]
end
end
for i=1:nx+2
pc[i, 1] = pc[i, 2]
pc[i, ny+2] = pc[i, ny+1]
end
for j=1:ny+2
pc[nx+2, j] = pc[nx+1, j]
pc[1, j] = pc[2, j]
end
rest=convergP(para, pc, ptemp)
if rest < 1e-5 && ip>=2
#println("rest PC ", rest, ", ip=", ip)
break
end
end
#println("pc=",pc[15,15])
#println("max pc=",findmax(abs.(pc)))
end
function update(para, ua, va, pc, apuva, ua2, va2)
@unpack nx, ny, dx, dy, = para
# u** = u* + u'
for j=2:ny+1
for i=2:nx+1
ua2[i,j] = ua[i,j] -0.5dy*(pc[i+1,j] - pc[i-1,j])/apuva[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
va2[i,j] = va[i,j] -0.5dx*(pc[i,j+1] - pc[i,j-1])/apuva[i,j]
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)
@unpack nx, ny, dx, dy, nu, ipter, alphapa = para
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]*(ua2[i,j]-ua[i,j])
else
uE[i,j] = aeuva[i+1,j]*(ua2[i+2,j]-ua[i+2,j]) +awuva[i+1,j]*(ua2[i,j]-ua[i,j]) +anuva[i+1,j]*(ua2[i+1,j+1]-ua[i+1,j+1]) +asuva[i+1,j]*(ua2[i+1,j-1]-ua[i+1,j-1])
end
if i==2
uW[i,j] = aeuva[i-1,j]*(ua2[i,j]-ua[i,j])
else
uW[i,j] = aeuva[i-1,j]*(ua2[i,j]-ua[i,j]) +awuva[i-1,j]*(ua2[i-2,j]-ua[i-2,j]) +anuva[i-1,j]*(ua2[i-1,j+1]-ua[i-1,j+1]) +asuva[i-1,j]*(ua2[i-1,j-1]-ua[i-1,j-1])
end
if j==ny+1
vN[i,j] = asuva[i,j+1]*(va2[i,j]-va[i,j])
else
vN[i,j] = aeuva[i,j+1]*(va2[i+1,j+1]-va[i+1,j+1]) +awuva[i,j+1]*(va2[i-1,j+1]-va[i-1,j+1]) +anuva[i,j+1]*(va2[i,j+2]-va[i,j+2]) +asuva[i,j+1]*(va2[i,j]-va[i,j])
end
if j==2
vS[i,j] = anuva[i,j-1]*(va2[i,j]-va[i,j])
else
vS[i,j] = aeuva[i,j-1]*(va2[i+1,j-1]-va[i+1,j-1]) +awuva[i,j-1]*(va2[i-1,j-1]-va[i-1,j-1]) +anuva[i,j-1]*(va2[i,j]-va[i,j]) +asuva[i,j-1]*(va2[i,j-2]-va[i,j-2])
end
end
end
for j=2:ny+1
for i=2:nx+1
bpc[i,j] = -0.5(dy/apuva[i+1,j]*uE[i,j] -dy/apuva[i-1,j]*uW[i,j]) -0.5(dx/apuva[i,j+1]*vN[i,j] -dx/apuva[i,j-1]*vS[i,j])
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:ipter
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] = ptemp[i,j]*(1.0 - alphapa) + alphapa*(aepc[i,j]*ptemp[i+1,j] + awpc[i,j]*ptemp[i-1,j] + anpc[i,j]*ptemp[i,j+1] + aspc[i,j]*ptemp[i,j-1] + bpc[i,j])/appc[i,j]
end
end
for i=1:nx+2
pc2[i, 1] = pc2[i, 2]
pc2[i, ny+2] = pc2[i, ny+1]
end
for j=1:ny+2
pc2[nx+2, j] = pc2[nx+1, j]
pc2[1, j] = pc2[2, j]
end
rest=convergP(para, pc2, ptemp)
if rest < 1e-5 && ip>=2
#println("rest PC2 ", rest, ", ip=", ip)
break
end
end
#println("max pc2=",findmax(abs.(pc2)))
end
function update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva)
@unpack nx, ny, dx, dy, = para
#println("max ua2-ua=",findmax(abs.(ua2[2:nx,2:ny+1]-ua[2:nx,2:ny+1])))
#println("max va2-va=",findmax(abs.(va2[2:nx+1,2:ny]-va[2:nx+1,2:ny])))
#println("ua2[2,ny+2], ua[2,ny+2]", ua2[2,ny+2], ",",ua[2,ny+2])
# 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.5dy*(pc2[i+1,j] - pc2[i-1,j]))/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.5dx*(pc2[i,j+1] - pc2[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
pref = pa[2,2]
for j=2:ny+1
for i=2:nx+1
pa[i,j] -= pref
end
end
end
function solveTA(para, ua, va, Ta, apTa, aeTa, awTa, anTa, asTa, bTa, Told, Ttemp)
@unpack nx, ny, dx, dy, dt, alpha, = para
for j=2:ny+1
for i=2:nx+1
if i==2
awTa[i,j] = 2alpha*dy/dx
else
awTa[i,j] = alpha*dy/dx + dy*max(0.5*(ua[i,j] + ua[i-1,j]), 0.0)
end
if i==nx+1
aeTa[i,j]=2alpha*dy/dx
else
aeTa[i,j] = alpha*dy/dx + dy*max(-0.5*(ua[i+1,j] + ua[i,j]), 0.0)
end
if j == ny+1
anTa[i,j] = 2alpha*dx/dy
else
anTa[i,j] = alpha*dx/dy + dx*max(-0.5*(va[i,j+1] + va[i,j]), 0.0)
end
if j == 2
asTa[i,j] = 2alpha*dx/dy
else
asTa[i,j] = alpha*dx/dy + dx*max(0.5*(va[i,j] + va[i,j-1]), 0.0)
end
apTa[i,j] = aeTa[i,j] + awTa[i,j] + anTa[i,j] + asTa[i,j] + dx*dy/dt
Ttemp[i,j] = (aeTa[i,j]*Ta[i+1,j] + awTa[i,j]*Ta[i-1,j] + anTa[i,j]*Ta[i,j+1] + asTa[i,j]*Ta[i,j-1] + bTa[i,j])/apTa[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
Ta[i,j] = Ttemp[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 = max(restua, restva)
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 output(para, ua, va, pa, Ta, uold, vold, Told, x, y, l)
@unpack outfile, dt, nx, ny = para
#println("time=", time)
filename = outfile*string(l)*".h5"
h5open(filename, "w") do file
write(file, "U", ua[2:nx+1, 2:ny+1])
write(file, "V", va[2:nx+1, 2:ny+1])
write(file, "T", Ta[2:nx+1, 2:ny+1])
write(file, "P", pa[2:nx+1, 2:ny+1])
#write(file, "Psi", psi[2:nx+1, 2:ny+1])
write(file, "X", x[2:nx+1, 2:ny+1])
write(file, "Y", y[2:nx+1, 2:ny+1])
end
h5open(outfile*"restart.h5", "w") do file
write(file, "ua_latest", ua)
write(file, "va_latest", va)
write(file, "pa_latest", pa)
write(file, "Ta_latest", Ta)
write(file, "ua_old", uold)
write(file, "va_old", vold)
write(file, "Ta_old", Told)
write(file, "n", l)
end
#open( outfile*time*".csv", "w") do f
# println(f, "x, y, u, v, pa")
# for j=2:ny+1
# for i=2:nx+1
# println(f, x[i,j], ",", y[i,j], ",",(u[i,j]+u[i-1,j])/2., ",", (v[i,j]+v[i,j-1])/2.0, ", ", pa[i,j] )
# end
# end
#end
end
function restart(para, restartFile, x, y, ua, va, pa, Ta, uold, vold, Told)
@unpack nx, ny, dx, dy = para
for j in 2:ny+1
for i in 2:nx+1
x[i,j] = (i-2+0.5)*dx
end
end
for j in 2:ny+1
for i in 2:nx+1
y[i,j] = (j-2+0.5)*dy
end
end
file = h5open(restartFile, "r")
ua .= read(file, "ua_latest")
va .= read(file, "va_latest")
pa .= read(file, "pa_latest")
Ta .= read(file, "Ta_latest")
uold .= read(file, "ua_old")
vold .= read(file, "va_old")
Told .= read(file, "Ta_old")
nstart = read(file, "n")
close(file)
return nstart
end
function streamfunction(para, u, v, psi)
@unpack dx, dy, nx, ny = para
for j=2:ny+1
for i=2:nx+1
psi[i,j] = 0.0
end
end
for j=2:ny+1
for i=2:nx+1
for k=2:j
psi[i,j] += u[i,k]*dy
end
end
end
end
#para = Param()
#@time main(para)
structでデータを渡すコード
using HDF5
using Parameters
@with_kw struct Param{T1,T2}
dt::T1 = 0.002
# Grashoff number
Ra::T1 = 1e8
Pr::T1 = 5.3
Gr::T1 = Ra/Pr
nu::T1 = 1.0/sqrt(Gr)
alpha::T1 = 1.0/Pr/sqrt(Gr)
# boundary temp
Th::T1 = 1.0
Tc::T1 = 0.0
nstep::T2 = 2000
niter::T2 = 15
ipter::T2 = 10
tol::T1 = 1e-5
alphapa::T1 = 0.5
TVDlimitter::T2 = 1 #van Albada
# mesh size
nx::T2 = 512
ny::T2 = 512
# length
Lx::T1 = 1.0
Ly::T1 = 1.0
dx::T1 = Lx/nx
dy::T1 = Ly/ny
ndata::T2 = 200
outfile::String = "Piso-Ra1e8-TVD/uv"
end
mutable struct uvTps
ua::Array{Float64, 2}
ua2::Array{Float64, 2}
ua3::Array{Float64, 2}
utemp::Array{Float64, 2}
uold::Array{Float64, 2}
va::Array{Float64, 2}
va2::Array{Float64, 2}
va3::Array{Float64, 2}
vtemp::Array{Float64, 2}
vold::Array{Float64, 2}
uE::Array{Float64, 2}
uW::Array{Float64, 2}
vN::Array{Float64, 2}
vS::Array{Float64, 2}
uf::Array{Float64, 2}
vf::Array{Float64, 2}
Ta::Array{Float64, 2}
Told::Array{Float64, 2}
Ttemp::Array{Float64, 2}
pa::Array{Float64, 2}
pc::Array{Float64, 2}
ptemp::Array{Float64, 2}
pc2::Array{Float64, 2}
bua::Array{Float64, 2}
bva::Array{Float64, 2}
bTa::Array{Float64, 2}
bpc::Array{Float64, 2}
end
mutable struct coeffs
apuva::Array{Float64, 2}
aeuva::Array{Float64, 2}
awuva::Array{Float64, 2}
anuva::Array{Float64, 2}
asuva::Array{Float64, 2}
apTa::Array{Float64, 2}
aeTa::Array{Float64, 2}
awTa::Array{Float64, 2}
anTa::Array{Float64, 2}
asTa::Array{Float64, 2}
appc::Array{Float64, 2}
aepc::Array{Float64, 2}
awpc::Array{Float64, 2}
anpc::Array{Float64, 2}
aspc::Array{Float64, 2}
end
mutable struct TVDcoeffs
Su::Array{Float64, 2}
Sv::Array{Float64, 2}
ST::Array{Float64, 2}
re⁺::Array{Float64, 2}
re⁻::Array{Float64, 2}
rn⁺::Array{Float64, 2}
rn⁻::Array{Float64, 2}
αe::Array{Int64, 2}
αn::Array{Int64, 2}
end
uvTps(para::Param) = uvTps(
zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2))
coeffs(para::Param) = coeffs(
zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2))
TVDcoeffs(para::Param) = TVDcoeffs(
zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2), zeros(Float64, para.nx+2, para.ny+2),
zeros(Int64, para.nx+2, para.ny+2), zeros(Int64, para.nx+2, para.ny+2))
function main(para)
@unpack nx, ny, nstep, niter, ndata, tol, TVDlimitter = para
# coordinate
x = zeros(Float64, nx+1, ny+1)
y = zeros(Float64, nx+1, ny+1)
uvTp = uvTps(para)
coeff = coeffs(para)
TVDcoeff = TVDcoeffs(para)
ψ = selectTVDlimitter(TVDlimitter)
rest = 0.0
init(para, x, y)
for l = 1:nstep
BC(para, uvTp)
for n = 1:niter
BC(para, uvTp)
makesource(para, uvTp, TVDcoeff, ψ)
solveUA(para, uvTp, coeff)
solveVA(para, uvTp, coeff)
solvePC(para, uvTp, coeff)
update(para, uvTp, coeff)
solvePC2(para, uvTp, coeff)
update2(para, uvTp, coeff)
solveTA(para, uvTp, coeff)
rest = converg(para, uvTp.ua, uvTp.va, uvTp.utemp, uvTp.vtemp)
#println("iter, redisual ", n,", ",rest)
if rest <=tol && n>=2
if l%10 == 0
println(l, " th step, ", n, " th iteration, residual=", rest)
end
break
elseif n==niter
println(l, " th step, ", n, " th iteration, u residual=", rest,)
end
end
for j=1:ny+2
for i=1:nx+2
uvTp.uold[i,j] = uvTp.ua[i,j]
uvTp.vold[i,j] = uvTp.va[i,j]
uvTp.Told[i,j] = uvTp.Ta[i,j]
end
end
if l%ndata == 0
println("residual=", rest)
output(para, uvTp, x, y, l)
end
end
#streamfunction(para, u, v, psi)
end
function init(para, x, y)
@unpack nx, ny, dx, dy = para
for j in 2:ny+1
for i in 2:nx+1
x[i,j] = (i-2+0.5)*dx
end
end
for j in 2:ny+1
for i in 2:nx+1
y[i,j] = (j-2+0.5)*dy
end
end
end
function BC(para, uvTp)
@unpack nx, ny, dy, Th, Tc = para
for i=1:nx+2
uvTp.ua[i, ny+2] = 0.0
uvTp.va[i, ny+2] = 0.0
uvTp.ua[i, 1] = 0.0
uvTp.va[i, 1] = 0.0
# ua2[i, ny+2]をupdate2でつかうが，updateではここの値は更新されないので，BCで与える。
uvTp.ua2[i, ny+2] = 0.0
uvTp.va2[i, ny+2] = 0.0
uvTp.ua2[i, 1] = 0.0
uvTp.va2[i, 1] = 0.0
uvTp.pa[i, 1] = uvTp.pa[i, 2] - uvTp.Ta[i,1]*dy
uvTp.pa[i, ny+2] = uvTp.pa[i, ny+1] + uvTp.Ta[i,ny+2]*dy
uvTp.pc[i, 1] = uvTp.pc[i, 2]
uvTp.pc[i, ny+2] = uvTp.pc[i, ny+1]
uvTp.pc2[i, 1] = uvTp.pc2[i, 2]
uvTp.pc2[i, ny+2] = uvTp.pc2[i, ny+1]
uvTp.Ta[i, ny+2] = Tc
uvTp.Ta[i, 1] = Th
end
for j=1:ny+2
uvTp.ua[nx+2, j] = 0.0
uvTp.va[nx+2, j] = 0.0
uvTp.ua[1, j] = 0.0
uvTp.va[1, j] = 0.0
uvTp.ua2[nx+2, j] = 0.0
uvTp.va2[nx+2, j] = 0.0
uvTp.ua2[1, j] = 0.0
uvTp.va2[1, j] = 0.0
uvTp.pa[1, j] = uvTp.pa[2, j]
uvTp.pa[nx+2, j] = uvTp.pa[nx+1, j]
uvTp.pc[nx+2, j] = uvTp.pc[nx+1, j]
uvTp.pc[1, j] = uvTp.pc[2, j]
uvTp.pc2[nx+2, j] = uvTp.pc2[nx+1, j]
uvTp.pc2[1, j] = uvTp.pc2[2, j]
uvTp.Ta[nx+2, j] = uvTp.Ta[nx+1, j]
uvTp.Ta[1, j] = uvTp.Ta[2,j]
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 makesource(para, uvTp, TVDcoeff, ψ)
@unpack nx, ny, dx, dy, dt = para
for j=2:ny+1
for i=2:nx+1
uvTp.bua[i,j] = -0.5dy*(uvTp.pa[i+1,j] - uvTp.pa[i-1,j]) + dx*dy/dt*uvTp.uold[i,j]
uvTp.bva[i,j] = -0.5dx*(uvTp.pa[i,j+1] - uvTp.pa[i,j-1]) + dx*dy/dt*uvTp.vold[i,j] + dx*dy*uvTp.Told[i,j]
uvTp.bTa[i,j] = dx*dy/dt*uvTp.Told[i,j]
end
end
for j=2:ny
for i=2:nx
uvTp.uf[i,j] = 0.5*(uvTp.uold[i+1,j] + uvTp.uold[i,j])
uvTp.vf[i,j] = 0.5*(uvTp.vold[i,j+1] + uvTp.vold[i,j])
TVDcoeff.αe[i,j] = 1 - signbit(uvTp.uf[i,j])
TVDcoeff.αn[i,j] = 1 - signbit(uvTp.vf[i,j])
end
end
# Su
for j=2:ny
for i=2:nx
TVDcoeff.re⁺[i,j] = (uvTp.uold[i,j] - uvTp.uold[i-1,j])/(uvTp.uold[i+1,j] - uvTp.uold[i,j] + 1e-10)
TVDcoeff.re⁻[i,j] = (uvTp.uold[i+2,j] - uvTp.uold[i+1,j])/(uvTp.uold[i+1,j] - uvTp.uold[i,j]+ 1e-10)
TVDcoeff.rn⁺[i,j] = (uvTp.uold[i,j] - uvTp.uold[i,j-1])/(uvTp.uold[i,j+1] - uvTp.uold[i,j]+ 1e-10)
TVDcoeff.rn⁻[i,j] = (uvTp.uold[i,j+2] - uvTp.uold[i,j+1])/(uvTp.uold[i,j+1] - uvTp.uold[i,j]+ 1e-10)
end
end
for j=3:ny
for i=3:nx
TVDcoeff.Su[i,j] =
dy*0.5uvTp.uf[i,j]*((1-TVDcoeff.αe[i,j])*ψ(TVDcoeff.re⁻[i,j]) - TVDcoeff.αe[i,j]*ψ(TVDcoeff.re⁺[i,j]))*(uvTp.uold[i+1,j] - uvTp.uold[i,j]) -
dy*0.5uvTp.uf[i-1,j]*((1-TVDcoeff.αe[i-1,j])*ψ(TVDcoeff.re⁻[i-1,j]) - TVDcoeff.αe[i-1,j]*ψ(TVDcoeff.re⁺[i-1,j]))*(uvTp.uold[i,j] - uvTp.uold[i-1,j]) +
dx*0.5uvTp.vf[i,j]*((1-TVDcoeff.αn[i,j])*ψ(TVDcoeff.rn⁻[i,j]) - TVDcoeff.αn[i,j]*ψ(TVDcoeff.rn⁺[i,j]))*(uvTp.uold[i,j+1] - uvTp.uold[i,j]) -
dx*0.5uvTp.vf[i,j-1]*((1-TVDcoeff.αn[i,j-1])*ψ(TVDcoeff.rn⁻[i,j-1]) - TVDcoeff.αn[i,j-1]*ψ(TVDcoeff.rn⁺[i,j-1]))*(uvTp.uold[i,j] - uvTp.uold[i,j-1])
end
end
# Sv
for j=2:ny
for i=2:nx
TVDcoeff.re⁺[i,j] = (uvTp.vold[i,j] - uvTp.vold[i-1,j])/(uvTp.vold[i+1,j]-uvTp.vold[i,j]+ 1e-10)
TVDcoeff.re⁻[i,j] = (uvTp.vold[i+2,j] - uvTp.vold[i+1,j])/(uvTp.vold[i+1,j]-uvTp.vold[i,j]+ 1e-10)
TVDcoeff.rn⁺[i,j] = (uvTp.vold[i,j] - uvTp.vold[i,j-1])/(uvTp.vold[i,j+1]-uvTp.vold[i,j]+ 1e-10)
TVDcoeff.rn⁻[i,j] = (uvTp.vold[i,j+2] - uvTp.vold[i,j+1])/(uvTp.vold[i,j+1]-uvTp.vold[i,j]+ 1e-10)
end
end
for j=3:ny
for i=3:nx
TVDcoeff.Sv[i,j] =
dy*0.5uvTp.uf[i,j]*((1-TVDcoeff.αe[i,j])*ψ(TVDcoeff.re⁻[i,j]) - TVDcoeff.αe[i,j]*ψ(TVDcoeff.re⁺[i,j]))*(uvTp.vold[i+1,j] - uvTp.vold[i,j]) -
dy*0.5uvTp.uf[i-1,j]*((1-TVDcoeff.αe[i-1,j])*ψ(TVDcoeff.re⁻[i-1,j]) - TVDcoeff.αe[i-1,j]*ψ(TVDcoeff.re⁺[i-1,j]))*(uvTp.vold[i,j] - uvTp.vold[i-1,j]) +
dx*0.5uvTp.vf[i,j]*((1-TVDcoeff.αn[i,j])*ψ(TVDcoeff.rn⁻[i,j]) - TVDcoeff.αn[i,j]*ψ(TVDcoeff.rn⁺[i,j]))*(uvTp.vold[i,j+1] - uvTp.vold[i,j]) -
dx*0.5uvTp.vf[i,j-1]*((1-TVDcoeff.αn[i,j-1])*ψ(TVDcoeff.rn⁻[i,j-1]) - TVDcoeff.αn[i,j-1]*ψ(TVDcoeff.rn⁺[i,j-1]))*(uvTp.vold[i,j] - uvTp.vold[i,j-1])
end
end
# ST
for j=2:ny
for i=2:nx
TVDcoeff.re⁺[i,j] = (uvTp.Told[i,j]-uvTp.Told[i-1,j])/(uvTp.Told[i+1,j]-uvTp.Told[i,j]+ 1e-10)
TVDcoeff.re⁻[i,j] = (uvTp.Told[i+2,j]-uvTp.Told[i+1,j])/(uvTp.Told[i+1,j]-uvTp.Told[i,j]+ 1e-10)
TVDcoeff.rn⁺[i,j] = (uvTp.Told[i,j]-uvTp.Told[i,j-1])/(uvTp.Told[i,j+1]-uvTp.Told[i,j]+ 1e-10)
TVDcoeff.rn⁻[i,j] = (uvTp.Told[i,j+2]-uvTp.Told[i,j+1])/(uvTp.Told[i,j+1]-uvTp.Told[i,j]+ 1e-10)
end
end
for j=3:ny
for i=3:nx
TVDcoeff.ST[i,j] =
dy*0.5uvTp.uf[i,j]*((1-TVDcoeff.αe[i,j])*ψ(TVDcoeff.re⁻[i,j]) - TVDcoeff.αe[i,j]*ψ(TVDcoeff.re⁺[i,j]))*(uvTp.Told[i+1,j] - uvTp.Told[i,j]) -
dy*0.5uvTp.uf[i-1,j]*((1-TVDcoeff.αe[i-1,j])*ψ(TVDcoeff.re⁻[i-1,j]) - TVDcoeff.αe[i-1,j]*ψ(TVDcoeff.re⁺[i-1,j]))*(uvTp.Told[i,j] - uvTp.Told[i-1,j]) +
dx*0.5uvTp.vf[i,j]*((1-TVDcoeff.αn[i,j])*ψ(TVDcoeff.rn⁻[i,j]) - TVDcoeff.αn[i,j]*ψ(TVDcoeff.rn⁺[i,j]))*(uvTp.Told[i,j+1] - uvTp.Told[i,j]) -
dx*0.5uvTp.vf[i,j-1]*((1-TVDcoeff.αn[i,j-1])*ψ(TVDcoeff.rn⁻[i,j-1]) - TVDcoeff.αn[i,j-1]*ψ(TVDcoeff.rn⁺[i,j-1]))*(uvTp.Told[i,j] - uvTp.Told[i,j-1])
end
end
for j=3:ny
for i=3:nx
uvTp.bua[i,j] += TVDcoeff.Su[i,j]
uvTp.bva[i,j] += TVDcoeff.Sv[i,j]
uvTp.bTa[i,j] += TVDcoeff.ST[i,j]
end
end
end
function solveUA(para, uvTp, coeff)
@unpack nx, ny, dx, dy, dt, nu, = para
for j=1:ny+2
for i=1:nx+2
if i==nx+2
coeff.aeuva[i,j] = 0.0
elseif i==nx+1
coeff.aeuva[i,j] = 2nu*dy/dx
else
coeff.aeuva[i,j] = nu*dy/dx + dy*max(-0.5*(uvTp.ua[i+1,j] + uvTp.ua[i,j]), 0.0)#- 0.25*dy*(ua[i+1, j] + ua[i,j])
end
if i==1
coeff.awuva[i,j] = 0.0
elseif i==2
coeff.awuva[i,j] = 2nu*dy/dx
else
coeff.awuva[i,j] = nu*dy/dx + dy*max(0.5*(uvTp.ua[i,j] + uvTp.ua[i-1,j]), 0.0)
end
if j==ny+2
coeff.anuva[i,j] = 0.0
elseif j == ny+1
coeff.anuva[i,j] = 2nu*dx/dy
else
coeff.anuva[i,j] = nu*dx/dy + dx*max(-0.5*(uvTp.va[i,j+1] + uvTp.va[i,j]), 0.0)#- 0.25dx*(va[i,j+1] + va[i,j])
end
if j==1
coeff.asuva[i,j] = 0.0
elseif j==2
coeff.asuva[i,j] = 2nu*dx/dy
else
coeff.asuva[i,j] = nu*dx/dy + dx*max(0.5*(uvTp.va[i,j] + uvTp.va[i,j-1]), 0.0)
end
coeff.apuva[i,j] = coeff.aeuva[i,j] + coeff.awuva[i,j] + coeff.anuva[i,j] + coeff.asuva[i,j] + dx*dy/dt
end
end
for j=2:ny+1
for i=2:nx+1
uvTp.utemp[i,j] = (coeff.aeuva[i,j]*uvTp.ua[i+1,j] + coeff.awuva[i,j]*uvTp.ua[i-1,j] + coeff.anuva[i,j]*uvTp.ua[i,j+1] + coeff.asuva[i,j]*uvTp.ua[i,j-1] + uvTp.bua[i,j])/coeff.apuva[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
uvTp.ua[i,j] = uvTp.utemp[i,j]
end
end
for j=1:ny+1
for i=2:nx
uvTp.uf[i,j] = 0.5(uvTp.ua[i+1,j] + uvTp.ua[i,j]) + 0.5(dy/coeff.apuva[i,j] + dy/coeff.apuva[i+1,j])*(uvTp.pa[i,j] - uvTp.pa[i+1,j]) - 0.25dy/coeff.apuva[i,j]*(uvTp.pa[i-1,j] - uvTp.pa[i+1,j]) -0.25dy/coeff.apuva[i+1,j]*(uvTp.pa[i,j] - uvTp.pa[i+2,j])
end
uvTp.uf[1,j] = 0.0
uvTp.uf[nx+1,j] = 0.0
end
end
function solveVA(para, uvTp, coeff)
@unpack nx, ny, dx, dy, dt, nu, = para
for j=2:ny+1
for i=2:nx+1
uvTp.vtemp[i,j] = (coeff.aeuva[i,j]*uvTp.va[i+1,j] + coeff.awuva[i,j]*uvTp.va[i-1,j] + coeff.anuva[i,j]*uvTp.va[i,j+1] + coeff.asuva[i,j]*uvTp.va[i,j-1] + uvTp.bva[i,j])/coeff.apuva[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
uvTp.va[i,j] = uvTp.vtemp[i,j]
end
end
for j=2:ny
for i=1:nx+1
uvTp.vf[i,j] = 0.5(uvTp.va[i,j+1] + uvTp.va[i,j]) + 0.5(dx/coeff.apuva[i,j] + dx/coeff.apuva[i,j+1])*(uvTp.pa[i,j] - uvTp.pa[i,j+1]) - 0.25dx/coeff.apuva[i,j]*(uvTp.pa[i,j-1] - uvTp.pa[i,j+1]) -0.25dx/coeff.apuva[i,j+1]*(uvTp.pa[i,j] - uvTp.pa[i,j+2])
end
end
for i=1:nx+1
uvTp.vf[i,1] = 0.0
uvTp.vf[i,ny+1] = 0.0
end
end
function solvePC(para, uvTp, coeff)
@unpack nx, ny, dx, dy, nu, ipter, alphapa = para
for j=2:ny+1
for i=2:nx+1
uvTp.pc[i,j] = 0.0
uvTp.bpc[i,j] = 0.0
end
end
for j=2:ny+1
for i=2:nx+1
if i == 2
coeff.awpc[i,j] = 0#0.5dy*(dy/apuva[i,j])
else
coeff.awpc[i,j] = 0.5dy*(dy/coeff.apuva[i-1,j] + dy/coeff.apuva[i,j])
end
if i == nx+1
coeff.aepc[i,j] = 0#0.5dy*(dy/apuva[i,j])
else
coeff.aepc[i,j] = 0.5dy*(dy/coeff.apuva[i,j] + dy/coeff.apuva[i+1,j])
end
if j == 2
coeff.aspc[i,j] = 0#0.5dx*(dx/apuva[i,j])
else
coeff.aspc[i,j] = 0.5dx*(dx/coeff.apuva[i,j-1] + dx/coeff.apuva[i,j])
end
if j == ny+1
coeff.anpc[i,j] = 0#0.5dx*(dx/apuva[i,j])
else
coeff.anpc[i,j] = 0.5dx*(dx/coeff.apuva[i,j] + dx/coeff.apuva[i,j+1])
end
coeff.appc[i,j] = coeff.aepc[i,j] + coeff.awpc[i,j] + coeff.anpc[i,j] + coeff.aspc[i,j]
uvTp.bpc[i,j] = -dy*(uvTp.uf[i,j] - uvTp.uf[i-1,j]) - dx*(uvTp.vf[i,j] - uvTp.vf[i,j-1])
end
end
# SOR法で圧力の解を求めてる
for ip=1:ipter
for j=2:ny+1
for i=2:nx+1
uvTp.ptemp[i,j] = uvTp.pc[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
uvTp.pc[i,j] = uvTp.ptemp[i,j]*(1.0 - alphapa) + alphapa*(coeff.aepc[i,j]*uvTp.ptemp[i+1,j] + coeff.awpc[i,j]*uvTp.ptemp[i-1,j] + coeff.anpc[i,j]*uvTp.ptemp[i,j+1] + coeff.aspc[i,j]*uvTp.ptemp[i,j-1] + uvTp.bpc[i,j])/coeff.appc[i,j]
end
end
for i=1:nx+2
uvTp.pc[i, 1] = uvTp.pc[i, 2]
uvTp.pc[i, ny+2] = uvTp.pc[i, ny+1]
end
for j=1:ny+2
uvTp.pc[nx+2, j] = uvTp.pc[nx+1, j]
uvTp.pc[1, j] = uvTp.pc[2, j]
end
rest=convergP(para, uvTp.pc, uvTp.ptemp)
if rest < 1e-5 && ip>=2
#println("rest PC ", rest, ", ip=", ip)
break
end
end
#println("pc=",pc[15,15])
#println("max pc=",findmax(abs.(pc)))
end
function update(para, uvTp, coeff)
@unpack nx, ny, dx, dy, = para
# u** = u* + u'
for j=2:ny+1
for i=2:nx+1
uvTp.ua2[i,j] = uvTp.ua[i,j] -0.5dy*(uvTp.pc[i+1,j] - uvTp.pc[i-1,j])/coeff.apuva[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
uvTp.va2[i,j] = uvTp.va[i,j] -0.5dx*(uvTp.pc[i,j+1] - uvTp.pc[i,j-1])/coeff.apuva[i,j]
end
end
end
function solvePC2(para, uvTp, coeff)
@unpack nx, ny, dx, dy, nu, ipter, alphapa = para
for j=2:ny+1
for i=2:nx+1
uvTp.pc2[i,j] = 0.0
uvTp.bpc[i,j] = 0.0
end
end
for j=2:ny+1
for i=2:nx+1
if i==nx+1
uvTp.uE[i,j] = coeff.awuva[i+1,j]*(uvTp.ua2[i,j] - uvTp.ua[i,j])
else
uvTp.uE[i,j] = coeff.aeuva[i+1,j]*(uvTp.ua2[i+2,j] - uvTp.ua[i+2,j]) + coeff.awuva[i+1,j]*(uvTp.ua2[i,j] - uvTp.ua[i,j]) + coeff.anuva[i+1,j]*(uvTp.ua2[i+1,j+1] - uvTp.ua[i+1,j+1]) + coeff.asuva[i+1,j]*(uvTp.ua2[i+1,j-1] - uvTp.ua[i+1,j-1])
end
if i==2
uvTp.uW[i,j] = coeff.aeuva[i-1,j]*(uvTp.ua2[i,j] - uvTp.ua[i,j])
else
uvTp.uW[i,j] = coeff.aeuva[i-1,j]*(uvTp.ua2[i,j] - uvTp.ua[i,j]) + coeff.awuva[i-1,j]*(uvTp.ua2[i-2,j] - uvTp.ua[i-2,j]) + coeff.anuva[i-1,j]*(uvTp.ua2[i-1,j+1] - uvTp.ua[i-1,j+1]) + coeff.asuva[i-1,j]*(uvTp.ua2[i-1,j-1] - uvTp.ua[i-1,j-1])
end
if j==ny+1
uvTp.vN[i,j] = coeff.asuva[i,j+1]*(uvTp.va2[i,j] - uvTp.va[i,j])
else
uvTp.vN[i,j] = coeff.aeuva[i,j+1]*(uvTp.va2[i+1,j+1] - uvTp.va[i+1,j+1]) + coeff.awuva[i,j+1]*(uvTp.va2[i-1,j+1] - uvTp.va[i-1,j+1]) + coeff.anuva[i,j+1]*(uvTp.va2[i,j+2] - uvTp.va[i,j+2]) + coeff.asuva[i,j+1]*(uvTp.va2[i,j] - uvTp.va[i,j])
end
if j==2
uvTp.vS[i,j] = coeff.anuva[i,j-1]*(uvTp.va2[i,j] - uvTp.va[i,j])
else
uvTp.vS[i,j] = coeff.aeuva[i,j-1]*(uvTp.va2[i+1,j-1] - uvTp.va[i+1,j-1]) + coeff.awuva[i,j-1]*(uvTp.va2[i-1,j-1] - uvTp.va[i-1,j-1]) + coeff.anuva[i,j-1]*(uvTp.va2[i,j] - uvTp.va[i,j]) + coeff.asuva[i,j-1]*(uvTp.va2[i,j-2] - uvTp.va[i,j-2])
end
end
end
for j=2:ny+1
for i=2:nx+1
uvTp.bpc[i,j] = -0.5(dy/coeff.apuva[i+1,j]*uvTp.uE[i,j] -dy/coeff.apuva[i-1,j]*uvTp.uW[i,j]) -0.5(dx/coeff.apuva[i,j+1]*uvTp.vN[i,j] -dx/coeff.apuva[i,j-1]*uvTp.vS[i,j])
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:ipter
for j=2:ny+1
for i=2:nx+1
uvTp.ptemp[i,j] = uvTp.pc2[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
uvTp.pc2[i,j] = uvTp.ptemp[i,j]*(1.0 - alphapa) + alphapa*(coeff.aepc[i,j]*uvTp.ptemp[i+1,j] + coeff.awpc[i,j]*uvTp.ptemp[i-1,j] + coeff.anpc[i,j]*uvTp.ptemp[i,j+1] + coeff.aspc[i,j]*uvTp.ptemp[i,j-1] + uvTp.bpc[i,j])/coeff.appc[i,j]
end
end
for i=1:nx+2
uvTp.pc2[i, 1] = uvTp.pc2[i, 2]
uvTp.pc2[i, ny+2] = uvTp.pc2[i, ny+1]
end
for j=1:ny+2
uvTp.pc2[nx+2, j] = uvTp.pc2[nx+1, j]
uvTp.pc2[1, j] = uvTp.pc2[2, j]
end
rest=convergP(para, uvTp.pc2, uvTp.ptemp)
if rest < 1e-5 && ip>=2
#println("rest PC2 ", rest, ", ip=", ip)
break
end
end
#println("max pc2=",findmax(abs.(pc2)))
end
function update2(para, uvTp, coeff)
@unpack nx, ny, dx, dy, = para
#println("max ua2-ua=",findmax(abs.(ua2[2:nx,2:ny+1]-ua[2:nx,2:ny+1])))
#println("max va2-va=",findmax(abs.(va2[2:nx+1,2:ny]-va[2:nx+1,2:ny])))
#println("ua2[2,ny+2], ua[2,ny+2]", ua2[2,ny+2], ",",ua[2,ny+2])
# u** = u* + u'
for j=2:ny+1
for i=2:nx+1
uvTp.ua3[i,j] = uvTp.ua2[i,j] + (coeff.aeuva[i,j]*(uvTp.ua2[i+1,j] - uvTp.ua[i+1,j]) + coeff.awuva[i,j]*(uvTp.ua2[i-1,j] - uvTp.ua[i-1,j]) + coeff.anuva[i,j]*(uvTp.ua2[i,j+1] - uvTp.ua[i,j+1]) + coeff.asuva[i,j]*(uvTp.ua2[i,j-1] - uvTp.ua[i,j-1]) -
0.5dy*(uvTp.pc2[i+1,j] - uvTp.pc2[i-1,j]))/coeff.apuva[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
uvTp.ua[i,j] = uvTp.ua3[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
uvTp.va3[i,j] = uvTp.va2[i,j] + (coeff.aeuva[i,j]*(uvTp.va2[i+1,j] - uvTp.va[i+1,j]) + coeff.awuva[i,j]*(uvTp.va2[i-1,j] - uvTp.va[i-1,j]) + coeff.anuva[i,j]*(uvTp.va2[i,j+1] - uvTp.va[i,j+1]) + coeff.asuva[i,j]*(uvTp.va2[i,j-1] - uvTp.va[i,j-1]) -
0.5dx*(uvTp.pc2[i,j+1] - uvTp.pc2[i,j-1]))/coeff.apuva[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
uvTp.va[i,j] = uvTp.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
uvTp.pa[i,j] = uvTp.pa[i,j] + uvTp.pc[i,j] + uvTp.pc2[i,j]
end
end
pref = uvTp.pa[2,2]
for j=2:ny+1
for i=2:nx+1
uvTp.pa[i,j] -= pref
end
end
end
function solveTA(para, uvTp, coeff)
@unpack nx, ny, dx, dy, dt, alpha, = para
for j=2:ny+1
for i=2:nx+1
if i==2
coeff.awTa[i,j] = 2alpha*dy/dx
else
coeff.awTa[i,j] = alpha*dy/dx + dy*max(0.5*(uvTp.ua[i,j] + uvTp.ua[i-1,j]), 0.0)
end
if i==nx+1
coeff.aeTa[i,j] = 2alpha*dy/dx
else
coeff.aeTa[i,j] = alpha*dy/dx + dy*max(-0.5*(uvTp.ua[i+1,j] + uvTp.ua[i,j]), 0.0)
end
if j == ny+1
coeff.anTa[i,j] = 2alpha*dx/dy
else
coeff.anTa[i,j] = alpha*dx/dy + dx*max(-0.5*(uvTp.va[i,j+1] + uvTp.va[i,j]), 0.0)
end
if j == 2
coeff.asTa[i,j] = 2alpha*dx/dy
else
coeff.asTa[i,j] = alpha*dx/dy + dx*max(0.5*(uvTp.va[i,j] + uvTp.va[i,j-1]), 0.0)
end
coeff.apTa[i,j] = coeff.aeTa[i,j] + coeff.awTa[i,j] + coeff.anTa[i,j] + coeff.asTa[i,j] + dx*dy/dt
uvTp.Ttemp[i,j] = (coeff.aeTa[i,j]*uvTp.Ta[i+1,j] + coeff.awTa[i,j]*uvTp.Ta[i-1,j] + coeff.anTa[i,j]*uvTp.Ta[i,j+1] + coeff.asTa[i,j]*uvTp.Ta[i,j-1] + uvTp.bTa[i,j])/coeff.apTa[i,j]
end
end
for j=2:ny+1
for i=2:nx+1
uvTp.Ta[i,j] = uvTp.Ttemp[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 = max(restua, restva)
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 output(para, uvTp, x, y, l)
@unpack outfile, dt, nx, ny = para
#println("time=", time)
filename = outfile*string(l)*".h5"
h5open(filename, "w") do file
write(file, "U", uvTp.ua[2:nx+1, 2:ny+1])
write(file, "V", uvTp.va[2:nx+1, 2:ny+1])
write(file, "T", uvTp.Ta[2:nx+1, 2:ny+1])
write(file, "P", uvTp.pa[2:nx+1, 2:ny+1])
#write(file, "Psi", psi[2:nx+1, 2:ny+1])
write(file, "X", x[2:nx+1, 2:ny+1])
write(file, "Y", y[2:nx+1, 2:ny+1])
end
h5open(outfile*"restart.h5", "w") do file
write(file, "ua_latest", uvTp.ua)
write(file, "va_latest", uvTp.va)
write(file, "pa_latest", uvTp.pa)
write(file, "Ta_latest", uvTp.Ta)
write(file, "ua_old", uvTp.uold)
write(file, "va_old", uvTp.vold)
write(file, "Ta_old", uvTp.Told)
write(file, "n", l)
end
#open( outfile*time*".csv", "w") do f
# println(f, "x, y, u, v, pa")
# for j=2:ny+1
# for i=2:nx+1
# println(f, x[i,j], ",", y[i,j], ",",(u[i,j]+u[i-1,j])/2., ",", (v[i,j]+v[i,j-1])/2.0, ", ", pa[i,j] )
# end
# end
#end
end
#para = Param()
#@time main(para)