コロケート格子上でのPISO法によるベナール対流の解析コード（Julia）

はじめに

基礎方程式は下記と同様です。問題設定もほぼ同様で，底面を高温，天井を低温にします。

コロケート格子を使うので，Rhie-Chow補間を行います。そこの説明は下記。

圧力の境界条件が必要になります。以前，浮力流れの論文を書いたくせ¹ に，圧力の境界条件が浮力があるときは少し変わるのを忘れててひどい目にあいました。（境界上の速度が変に振動してこれを取るのにくろうしました）

浮力がなければ壁面上で，下記の圧力条件になります。nは壁面に対する法線ベクトル。

\frac{\partial p}{\partial n}=0

浮力があるときの壁面上の境界条件は，基礎方程式中の速度に関する項を無視して残った項で，下記のようになります。

\frac{\partial p}{\partial n}=T

浮力に相当するTの項が残ってきます。下記コードではy方向に重力があるので，j=1, ny+2での圧力の境界条件にTが入ってきます。下記の部分。サイドの壁にはTの項が入ってきません。

pa[i, 1] = pa[i, 2]-Ta[i,1]*dy
pa[i, ny+2] = pa[i, ny+1]+Ta[i,ny+2]*dy

コード

using LinearAlgebra
using Parameters
using HDF5

@with_kw struct Param{T1,T2}
    dt::T1 = 0.01
    # Grashoff number
    Ra::T1 = 1e4
    Pr::T1 = 1.0
    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 = 200
    niter::T2 = 100
    ipter::T2 = 100
    tol::T1 = 1e-7

    alphapa::T1 = 0.5
    
    # mesh size
    nx::T2 = 40
    ny::T2 = 40
    
    # length
    Lx::T1 = 1.0
    Ly::T1 = 1.0
    dx::T1 = Lx/nx
    dy::T1 = Ly/ny

    ndata::T2 = 50

    outfile::String = "Piso-Ra1e4-co/uv" 

end

function main(para)
    @unpack nx, ny, nstep, niter, ndata, tol = 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)
    
    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)
            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
            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 = 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)
    
    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)
            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
            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 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 - 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 + 0.25*dy*(ua[i, j] + ua[i-1,j])
            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 - 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 + 0.25dx*(va[i,j] + va[i,j-1])
            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
            bua[i,j] = -0.5dy*(pa[i+1,j] - pa[i-1,j]) +dx*dy/dt*uold[i,j]
            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 
            #bva[i,j] = -0.5dx*(pa[i,j+1] - pa[i,j-1]) + dx*dy/dt*vold[i,j] + dx*dy*Ta[i, j]
            bva[i,j] = -0.5dx*(pa[i,j+1] - pa[i,j-1]) + dx*dy/dt*vold[i,j] + dx*dy*Ta[i,j]
            
            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
            else
                awpc[i,j] = 0.5dy*(dy/apuva[i-1,j]+dy/apuva[i,j])
            end
            
            if i == nx+1
                aepc[i,j] = 0
            else          
                aepc[i,j] = 0.5dy*(dy/apuva[i,j]+dy/apuva[i+1,j])
            end
            
            if j == 2
                aspc[i,j] = 0
            else
                aspc[i,j] = 0.5dx*(dx/apuva[i,j-1]+dx/apuva[i,j])
            end
                    
            if j == ny+1
                anpc[i,j] = 0
            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 + 0.25*dy*(ua[i, j] + ua[i-1,j])
            end
            if i==nx+1
                aeTa[i,j]=2alpha*dy/dx
            else
                aeTa[i,j] = alpha*dy/dx - 0.25*dy*(ua[i+1, j] + ua[i,j])
            end

            if j == ny+1
                anTa[i,j] = 2alpha*dx/dy
            else
                anTa[i,j] = alpha*dx/dy - 0.25dx*(va[i,j+1] + va[i,j])
            end
            if j == 2
                asTa[i,j] = 2alpha*dx/dy
            else
                asTa[i,j] = alpha*dx/dy + 0.25dx*(va[i,j] + va[i,j-1])
            end
            
            apTa[i,j] = aeTa[i,j] + awTa[i,j] + anTa[i,j] + asTa[i,j] + dx*dy/dt
            bTa[i,j] = dx*dy/dt*Told[i,j]
            
            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])
        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")
    
    return nstart
end

1. Influence of mesh non-orthogonality on numerical simulation of buoyant jet flowsっていう論文なので，興味があったら見て下さい。 ↩