概要

上記のコードに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)