#訂正
リスタートのためのデータ書き出しで，無駄に同じデータを別名で保存していたところを修正しました。
すいませんでした。

カルマン渦の計算

下記のようなカルマン渦の計算をするためのコードです。

背景

以前下記の記事を書きましたが，それにTVDスキームも実装しました。

TVD

Total Validation Diminishing だったかな。実装はVersteeg, Malalasekeraの「数値流体力学」を参考にしました。

下記のコードだと，makesourceの後半でTVDの実装してます。PISO法のa_P等にもいれる必要がありますが，
これはsolveUVAでやってます。

Paramのなかのパラメータで，TVDlimitterの値でスキームを選択できます。デフォルトでvan Albadaにしています。

グリッド生成

グリッドも前回から修正しました。境界面を指定するためのインデックスがコードに直接書き込んでありますが，
そこはご了承ください。

using Parameters
using HDF5

@with_kw struct Param{T1,T2}
    Rout::T1=10.0
    R1::T1 = 0.5
    xr::T1 = 20.
    ytop::T1 = 10.

    x1::T1 = 8.0
    x2::T1 = 10.0
    
    ε::T1 = 1e-5
    itermax::T2 = 100
    
    nx::T2 = 140
    ny::T2 = 561    
    
    outputname::String = "Grid/grid-cylinder-c"
end

function gridgen(para)
    @unpack nx, ny, = para
    
    x = zeros(nx, ny)
    y = zeros(nx, ny)
    xT = zeros(ny, nx)
    yT = zeros(ny, nx)
    # control function
    P = zeros(nx, ny)
    Q = zeros(nx, ny)
    
    initialize(para, x, y)
    
    xT .= x'
    yT .= y'
    output(para, x, y, xT, yT)
    
    # solve elliptic equation
    ellip(para, x, y, P, Q)
    

    xT .= x'
    yT .= y'
    
    output(para, x, y, xT, yT)
    
end

function initialize(para, x, y)
    @unpack nx, ny, R1, Rout, xr, ytop, x1, x2 = para
    
    dx = 0.05
    
    for j=1:div(ny,2)
        x[1, j] = xr - dx*(j-1)
    end
    
    x[1,div(ny,2)+1] = 7.5
    
    for j=div(ny,2)+2:ny
        x[1,j] = 7.5 + dx*(j-div(ny,2)-31)
    end
    
      for j=241:281
        for i=2:nx
            θ = 0.5*pi*((j-241)/(40)) + 0.5*pi
            R=Rout/1.001^(nx-i)
            x[i,j] = R*cos(θ) + x1
            y[i,j] = R*sin(θ)
            x[i,562-j] = R*cos(-θ) + x1
            y[i,562-j] = R*sin(-θ)
            x[i,281] = R*cos(pi) + x1
            y[i,281] = R*sin(pi)
        end
    end 
    
    for j=1:230
        y[1,j] = 0
    end
    
    for j=231:281
        θ = pi*(j-231)/50
        x[1,j] = R1*cos(θ) + x1
        y[1,j] = R1*sin(θ)
        x[1,562-j] = R1*cos(-θ) + x1
        y[1,562-j] = R1*sin(-θ)
    end
    x[1,281] = R1*cos(pi) + x1
    y[1,281] = R1*sin(pi)
    
    for j=332:ny
        y[1,j]=0
    end

    #for j=224:230
    #    x[1,j] = x1+R1 +0.005*1.8^(231-j)
    #end
    #for j=332:338
    #    x[1,j] = x1+R1 +0.005*1.8^(j-331)
    #end
    
    #外周
    for j=1:190
        for i=2:nx
            R = Rout
            x[i, j] = xr- dx*(j-1)
            y[i,j] = R/(nx-1)*(i-1)
        end
    end
    
    for j=190:241
        for i=2:nx-1
            R = Rout
            x[i, j] = xr- dx*(j-1)
            y[i,j] = 0.5R +0.5R/(nx-1)*(i-1)
        end
        i=nx
        R = Rout
        x[i,j] = xr- dx*(j-1)
        y[i,j] = R
    end
    
    for j= 321:372
        for i=2:nx-1
            R = Rout
            x[i,j] = 7.5 +dx*(j-div(ny,2)-31)
            y[i,j] = -0.5R/(nx-1)*(i-1) - 0.5R
        end
        i=nx
        R = Rout
        x[i,j] = 7.5 +dx*(j-div(ny,2)-31)
        y[i,j] = -R
    end
    
    for j=373:ny
       for i=2:nx
            R = Rout
            x[i,j] = 7.5 +dx*(j-div(ny,2)-31)
            y[i,j] = -R/(nx-1)*(i-1)
        end
    end
        
end

function output(para, x, y, xT, yT)
    @unpack nx, ny, outputname = para
    
    filename = outputname*".h5"
    h5open(filename, "w") do file
        write(file, "nx", nx)
        write(file, "ny", ny)
        write(file, "x", x[1:nx-1,:])
        write(file, "y", y[1:nx-1,:])
        write(file, "xall", x[1:nx,:])
        write(file, "yall", y[1:nx,:])
        write(file, "xT", xT[:,1:nx-1])
        write(file, "yT", yT[:,1:nx-1])
        write(file, "xTall", xT[:,1:nx])
        write(file, "yTall", yT[:,1:nx])
    end
end

function ellip(para, x, y, P, Q)
    @unpack nx, ny, itermax, ε = para
    
    itermax3 = div(itermax, 3)
    itermax6 = itermax3*2
    
    println(x[1:3,1:3])
    
    cons = 0.01
    
    for i = 1:itermax
        error = 0.0
        if i >= itermax3 
            cons = 0.25
        elseif i>=itermax6
            cons = 1.0
        end
        
        for l = 2:ny-1
            for k = 2:nx-1
                xξ = 0.5(x[k+1, l] - x[k-1, l])
                yξ = 0.5(y[k+1, l] - y[k-1, l])
                xη = 0.5(x[k, l+1] - x[k, l-1])
                yη = 0.5(y[k, l+1] - y[k, l-1])
                α = xη^2 + yη^2
                β = xξ*xη + yξ*yη
                γ = xξ^2 + yξ^2
                J = xξ*yη - xη*yξ
                
                PP = P[k, l]
                QQ = Q[k, l]
                ag = 0.5/(α + γ)
                
                sx = ( α*(x[k+1,l] + x[k-1,l]) -
                    0.5β*(x[k+1,l+1] - x[k-1,l+1] - x[k+1,l-1] + x[k-1,l-1]) + 
                    γ*(x[k,l+1] + x[k,l-1]) + J^2*(PP*xξ + QQ*xη) )*ag - x[k,l]
                sy = (α*(y[k+1,l] + y[k-1,l]) -
                    0.5β*(y[k+1,l+1] - y[k-1,l+1] - y[k+1,l-1] + y[k-1,l-1]) + 
                    γ*(y[k,l+1] + y[k,l-1]) + J^2*(PP*yξ + QQ*yη) )*ag - y[k,l]
                
                
                x[k,l] = x[k,l] + cons*sx
                y[k,l] = y[k,l] + cons*sy
                
                error += sx^2 + sy^2
            end
        end
        
        
        if error < ε
            println("iteration, error=", i, ", ", error)
            break
        elseif i==itermax
            println("error= ",error)
        end
        
        if mod(i,2000) == 0
            println("iteration=", i, ", error=", error)
        end
        
    end
end

para = Param( itermax=100000,)
@time gridgen(para)

2分くらいでグリッドができます。

流体のほうのコード


using Parameters
using HDF5

@with_kw struct Param{T1,T2}
    dt::T1 = 0.01
    # viscosity
    nu::T1 = 0.01
    Uext::T1 = 1.0

    nstep::T2 = 400
    alphaua::T1 = 0.5
    alphava::T1 = 0.5
    alphapa::T1 = 0.25
    niter::T2 = 20
    ipter::T2 = 5
    ipter2::T2 = 5

    start::T2 =1

    TVDlimitter::T2 = 1 #van Albada

    # mesh size
    nx::T2 = 137
    ny::T2 = 559

    ndata::T2 = 50

    gridfile::String = "Grid/grid-cylinder-c.h5"
    outfile::String = "Data-cyl-Re100/uv" 

end

function main(para)
    @unpack nx, ny, nstep, niter, ndata, Uext, TVDlimitter = para

    # auxiliary velocity
    ua = ones(Float64, nx+2, ny+2)
    ua2 = ones(Float64, nx+2, ny+2)
    ua3 = ones(Float64, nx+2, ny+2)
    utemp = zeros(Float64, nx+2, ny+2)
    uold = ones(Float64, nx+2, ny+2)
    apuva = zeros(Float64, nx+2, ny+2)
    aeuva = zeros(Float64, nx+2, ny+2)
    awuva  = zeros(Float64, nx+2, ny+2)
    anuva = zeros(Float64, nx+2, ny+2)
    asuva = zeros(Float64, nx+2, ny+2)
    bua  = zeros(Float64, nx+2, ny+2)
    va = zeros(Float64, nx+2, ny+2)
    va2 = zeros(Float64, nx+2, ny+2)
    va3 = zeros(Float64, nx+2, ny+2)
    vtemp = zeros(Float64, nx+2, ny+2)
    vold = zeros(Float64, nx+2, ny+2)
    bva = zeros(Float64, nx+2, ny+2)
    Ue = zeros(Float64, nx+2, ny+2)
    Vn = zeros(Float64, nx+2, ny+2)

    UE = zeros(Float64, nx+2, ny+2)
    UW = zeros(Float64, nx+2, ny+2)
    VN = zeros(Float64, nx+2, ny+2)
    VS = zeros(Float64, nx+2, ny+2)
    
    omega = zeros(Float64, nx+2, ny+2)

    beta2 = zeros(Float64, nx+2, ny+2)
    gam = zeros(Float64, nx+2, ny+2)
    U = zeros(Float64, nx+2, ny+2)
    V = zeros(Float64, nx+2, ny+2)
    #dξ = zeros(Float64, nx+2, ny+2)
    #dη = zeros(Float64, nx+2, ny+2)
    
    # pressure, corrected pressure
    pa = zeros(Float64, nx+2, ny+2)
    pc = zeros(Float64, nx+2, ny+2)
    pc2 = zeros(Float64, nx+2, ny+2)
    ptemp = zeros(Float64, nx+2, ny+2)
    appc = zeros(Float64, nx+2, ny+2)
    aepc = zeros(Float64, nx+2, ny+2)
    awpc = zeros(Float64, nx+2, ny+2)
    anpc = zeros(Float64, nx+2, ny+2)
    aspc = zeros(Float64, nx+2, ny+2)
    bpc = zeros(Float64, nx+2, ny+2)
    # coordinate
    x = zeros(Float64, nx+2, ny+2)
    y = zeros(Float64, nx+2, ny+2)
    st = zeros(Float64, ny+2)

    ξx = zeros(Float64, nx+2, ny+2)
    ξy = zeros(Float64, nx+2, ny+2) 
    ηx = zeros(Float64, nx+2, ny+2) 
    ηy = zeros(Float64, nx+2, ny+2) 
    J = zeros(Float64, nx+2, ny+2)

    Su = zeros(Float64, nx+2, ny+2)
    Sv = zeros(Float64, nx+2, ny+2)
    re⁺ = zeros(Float64, nx+2, ny+2)
    re⁻ = zeros(Float64, nx+2, ny+2)
    rn⁺ = zeros(Float64, nx+2, ny+2)
    rn⁻ = zeros(Float64, nx+2, ny+2)
    αe = zeros(Int64, nx+2, ny+2)
    αn = zeros(Int64, nx+2, ny+2)

    rest = 0.0

    readgrid(para, x, y, st)
    jacobian(para, x, y, ξx, ξy, ηx, ηy, J)
    ψ = selectTVDlimitter(TVDlimitter)


    for l = 1:nstep
    
        for j=1:ny+2
            for i=1:nx+2
                uold[i,j] = ua[i,j]
                vold[i,j] = va[i,j]
            end
        end
    
        BC(para, ua, va, pc, pa, ξx, ξy,  ηx, ηy, J, uold, vold, ua2, va2, pc2, l, st)
        
        for n = 1:niter
            BC(para, ua, va, pc, pa, ξx, ξy,  ηx, ηy, J, uold, vold, ua2, va2, pc2, l, st)
            makesource(para, bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J, U, V, Su, Sv, re⁺, re⁻, rn⁺, rn⁻, αe, αn, Ue, Vn, ψ)
            
            for k = 1:5
                solveUVA(para, ua, va, apuva, aeuva, awuva, anuva, asuva, bua, utemp, bva, vtemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
            end
            RhieChowInterpolation(para, ua, va, pa, apuva, U, V,  Ue, Vn, ξx, ξy, ηx, ηy, J)
            #println("ua, va=", ua[2,2],", ", va[2,2])
            solvePC(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, ptemp, Ue, Vn, ξx, ξy, ηx, ηy, J, st)
            update(para, ua, va, pc, pa, apuva, ua2, va2, ξx, ξy, ηx, ηy, J)
            solvePC2(para, ua, ua2, va, va2, pc, pc2, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, UE, UW, VN, VS, ptemp, ξx, ξy, ηx, ηy, U, V, J , st)
            update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva, ξx, ξy, ηx, ηy, J, utemp, vtemp)

            rest = converg(para, ua, va, ua2, va2)
            divu = divergence(para, ua, va, U, V, ξx, ξy, ηx, ηy, J) 
            
            if rest <=1e-7 && n >=2
                if l%20 == 0
                    println(l, " th step, ", n, " th iteration, u residual=", rest, ", div=", divu)
                end
                break
            elseif n==niter
                println(l, " th step, ", n, " th iteration, u residual=", rest, ", div=", divu)
            end
        end

        if l%ndata == 0
            println("residual=", rest)
            vor(para, omega, ua, va, ξx, ξy, ηx, ηy, st)
            output(para, ua, va, pa, l, omega)
        end
    end
    #streamfunction(para, u, v, psi)
    
end

function main(para, restartFile)
    @unpack nx, ny, nstep, niter, ndata, TVDlimitter = para

    # auxiliary velocity
    ua = zeros(Float64, nx+2, ny+2)
    ua2 = zeros(Float64, nx+2, ny+2)
    ua3 = zeros(Float64, nx+2, ny+2)
    utemp = zeros(Float64, nx+2, ny+2)
    uold = zeros(Float64, nx+2, ny+2)
    apuva = zeros(Float64, nx+2, ny+2)
    aeuva = zeros(Float64, nx+2, ny+2)
    awuva  = zeros(Float64, nx+2, ny+2)
    anuva = zeros(Float64, nx+2, ny+2)
    asuva = zeros(Float64, nx+2, ny+2)
    bua  = zeros(Float64, nx+2, ny+2)
    va = zeros(Float64, nx+2, ny+2)
    va2 = zeros(Float64, nx+2, ny+2)
    va3 = zeros(Float64, nx+2, ny+2)
    vtemp = zeros(Float64, nx+2, ny+2)
    vold = zeros(Float64, nx+2, ny+2)
    bva = zeros(Float64, nx+2, ny+2)
    Ue = zeros(Float64, nx+2, ny+2)
    Vn = zeros(Float64, nx+2, ny+2)

    omega = zeros(Float64, nx+2, ny+2)

    UE = zeros(Float64, nx+2, ny+2)
    UW = zeros(Float64, nx+2, ny+2)
    VN = zeros(Float64, nx+2, ny+2)
    VS = zeros(Float64, nx+2, ny+2)

    beta2 = zeros(Float64, nx+2, ny+2)
    gam = zeros(Float64, nx+2, ny+2)
    U = zeros(Float64, nx+2, ny+2)
    V = zeros(Float64, nx+2, ny+2)
    dξ = zeros(Float64, nx+2, ny+2)
    dη = zeros(Float64, nx+2, ny+2)
    
    # pressure, corrected pressure
    pa = zeros(Float64, nx+2, ny+2)
    pc = zeros(Float64, nx+2, ny+2)
    pc2 = zeros(Float64, nx+2, ny+2)
    ptemp = zeros(Float64, nx+2, ny+2)
    appc = zeros(Float64, nx+2, ny+2)
    aepc = zeros(Float64, nx+2, ny+2)
    awpc = zeros(Float64, nx+2, ny+2)
    anpc = zeros(Float64, nx+2, ny+2)
    aspc = zeros(Float64, nx+2, ny+2)
    bpc = zeros(Float64, nx+2, ny+2)
    # coordinate
    x = zeros(Float64, nx+2, ny+2)
    y = zeros(Float64, nx+2, ny+2)
    st = zeros(Float64, ny+2)

    ξx = zeros(Float64, nx+2, ny+2)
    ξy = zeros(Float64, nx+2, ny+2) 
    ηx = zeros(Float64, nx+2, ny+2) 
    ηy = zeros(Float64, nx+2, ny+2) 
    J = zeros(Float64, nx+2, ny+2)
    
    Su = zeros(Float64, nx+2, ny+2)
    Sv = zeros(Float64, nx+2, ny+2)
    re⁺ = zeros(Float64, nx+2, ny+2)
    re⁻ = zeros(Float64, nx+2, ny+2)
    rn⁺ = zeros(Float64, nx+2, ny+2)
    rn⁻ = zeros(Float64, nx+2, ny+2)
    αe = zeros(Int64, nx+2, ny+2)
    αn = zeros(Int64, nx+2, ny+2)

    rest = 0.0

    readgrid(para, x, y, st)
    jacobian(para, x, y, ξx, ξy, ηx, ηy, J)
    ψ = selectTVDlimitter(TVDlimitter)

    nstart = restart(para, restartFile, ua, va, pa,  )
    println("nstart=",nstart)

    for l = nstart+1:nstart+nstep
    
        for j=1:ny+2
            for i=1:nx+2
                uold[i,j] = ua[i,j]
                vold[i,j] = va[i,j]
            end
        end
        
        BC(para, ua, va, pc, pa, ξx, ξy,  ηx, ηy, J, uold, vold, ua2, va2, pc2, l, st)
        
        for n = 1:niter
            BC(para, ua, va, pc, pa, ξx, ξy,  ηx, ηy, J, uold, vold, ua2, va2, pc2, l, st)
            makesource(para, bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J, U, V, Su, Sv, re⁺, re⁻, rn⁺, rn⁻, αe, αn, Ue, Vn, ψ)
            
            for k = 1:5
                solveUVA(para, ua, va, apuva, aeuva, awuva, anuva, asuva, bua, utemp, bva, vtemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
            end
            RhieChowInterpolation(para, ua, va, pa, apuva, U, V,  Ue, Vn, ξx, ξy, ηx, ηy, J)
            #println("ua, va=", ua[2,2],", ", va[2,2])
            solvePC(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, ptemp, Ue, Vn, ξx, ξy, ηx, ηy, J, st)
            update(para, ua, va, pc, pa, apuva, ua2, va2, ξx, ξy, ηx, ηy, J)
            solvePC2(para, ua, ua2, va, va2, pc, pc2, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, UE, UW, VN, VS, ptemp, ξx, ξy, ηx, ηy, U, V, J , st)
            update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva, ξx, ξy, ηx, ηy, J, utemp, vtemp)

            rest = converg(para, ua, va, ua2, va2)
            divu = divergence(para, ua, va, U, V, ξx, ξy, ηx, ηy, J) 
 
            if rest <=1e-7 && n >=2
            #if divu <=1e-7 && n>=2
                if l%20 == 0
                    println(l, " th step, ", n, " th iteration, residual=", rest, ", div=", divu)
                end
                break
            elseif n==niter
                println(l, " th step, ", n, " th iteration, u residual=", rest, ", div=", divu)
            end
        end

        if l%ndata == 0
            println("residual=", rest)
            vor(para, omega, ua, va, ξx, ξy, ηx, ηy, st)
            output(para, ua, va, pa, l, omega)
        end
    end
    #streamfunction(para, u, v, psi)
    
end

function BC(para, ua, va, pc, pa, ξx, ξy,  ηx, ηy, J, uold, vold, ua2, va2, pc2, l, st)
    @unpack nx, ny, Uext, dt, start = para

    # cylinder 
    for j=232:330
        ua[1, j] = 0.0
        va[1, j] = 0.0
        ua2[1,j] = 0.0
        va2[1,j] = 0.0
        pc[1, j] = pc[2,j] 
        pc2[1, j] = pc2[2,j] 
        pa[1, j] = pa[2,j] 
    end

    # j=231は円柱表面だが，外挿で圧力を設定する。
    for j=1:231
        ua[1, j] = 0.5*(ua[2,j]-st[j]*(ua[3,j]-ua[2,j]) + ua[2,ny+3-j]-st[ny+3-j]*(ua[3,ny+3-j]-ua[2,ny+3-j]))
        ua[1, ny+3-j] = ua[1,j]
        va[1, j] = 0.5*(va[2,j]-st[j]*(va[3,j]-va[2,j]) + va[2,ny+3-j]-st[ny+3-j]*(va[3,ny+3-j]-va[2,ny+3-j]))
        va[1, ny+3-j] = va[1,j]
        ua2[1,j] = 0.5*(ua2[2,j]-st[j]*(ua2[3,j]-ua2[2,j]) + ua2[2,ny+3-j]-st[ny+3-j]*(ua2[3,ny+3-j]-ua2[2,ny+3-j]))
        ua2[1,ny+3-j] = ua2[1,j]
        va2[1,j] = 0.5*(va2[2,j]-st[j]*(va2[3,j]-va2[2,j]) + va2[2,ny+3-j]-st[ny+3-j]*(va2[3,ny+3-j]-va2[2,ny+3-j]))
        va2[1,ny+3-j] = va2[1,j]
        pc[1, j] = 0.5*(pc[2,j]-st[j]*(pc[3,j]-pc[2,j]) + pc[2,ny+3-j]-st[ny+3-j]*(pc[3,ny+3-j]-pc[2,ny+3-j]))
        pc[1, ny+3-j] = pc[1,j]
        pc2[1, j] = 0.5*(pc2[2,j]-st[j]*(pc2[3,j]-pc2[2,j]) + pc2[2,ny+3-j]-st[ny+3-j]*(pc2[3,ny+3-j]-pc2[2,ny+3-j]))
        pc2[1, ny+3-j] = pc2[1,j]
        pa[1, j] = 0.5*(pa[2,j]-st[j]*(pa[3,j]-pa[2,j]) + pa[2,ny+3-j]-st[ny+3-j]*(pa[3,ny+3-j]-pa[2,ny+3-j]))
        pa[1, ny+3-j] = pa[1,j]
    end
#
    # outlet
    for i=1:nx+2
        ua[i, 1] = ua[i, 2]
        va[i, 1] = va[i, 2]
        ua[i, ny+2] = ua[i, ny+1]
        va[i, ny+2] = va[i, ny+1]
        ua2[i, 1] = ua2[i, 2]
        va2[i, 1] = va2[i, 2]
        ua2[i, ny+2] = ua2[i, ny+1]
        va2[i, ny+2] = va2[i, ny+1]

        pa[i, 1] = 0 
        pc[i, 1] = 0 
        pc2[i, 1] = 0
        pa[i, ny+2] = 0 
        pc[i, ny+2] = 0 
        pc2[i, ny+2] = 0
    end

    # inlet
    Ubound = min(Uext*l/start, Uext)
    
    for j=1:ny+2
        ua[nx+2,j] = Ubound
        va[nx+2,j] = 0.0
        ua2[nx+2,j] = Ubound
        va2[nx+2,j] = 0.0
        pa[nx+2,j] = pa[nx+1, j] 
        pc[nx+2,j] = pc[nx+1, j]
        pc2[nx+2,j] = pc2[nx+1, j]
    end

end

function readgrid(para, x, y, st)
    @unpack gridfile, nx, ny = para

    file = h5open(gridfile, "r") 
    x .= read(file, "x")
    y .= read(file, "y")
    close(file)

    for j=1:ny+2
        st[j] = sqrt(((x[1,j]-x[2,j])^2 + (y[1,j]-y[2,j])^2)/((x[2,j]-x[3,j])^2 + (y[2,j]-y[3,j])^2))
    end
end

function jacobian(para, x, y, ξx, ξy, ηx, ηy, J)
    @unpack nx, ny = para

    xξ = zeros(Float64, nx+2, ny+2)
    yξ = zeros(Float64, nx+2, ny+2) 
    xη = zeros(Float64, nx+2, ny+2) 
    yη = zeros(Float64, nx+2, ny+2) 

    for j = 1:ny+2
        for i = 2:nx+1
            xξ[i,j] = 0.5(x[i+1,j] - x[i-1,j])
            yξ[i,j] = 0.5(y[i+1,j] - y[i-1,j])
        end
    end
    for j=1:ny+2
        if (231<=j<=331)
            xξ[1,j] = -(3*x[1,j] - 4*x[2,j] + x[3,j])*0.5
            yξ[1,j] = -(3*y[1,j] - 4*y[2,j] + y[3,j])*0.5
        else
            xξ[1,j] = 0.5(x[2,j] - x[2,ny+3-j])
            yξ[1,j] = 0.5(y[2,j] - y[2,ny+3-j])
        end
        xξ[nx+2,j] = (3*x[nx+2,j] - 4*x[nx+1,j] + x[nx,j])*0.5
        yξ[nx+2,j] = (3*y[nx+2,j] - 4*y[nx+1,j] + y[nx,j])*0.5
    end

    for j = 2:ny+1
        for i = 1:nx+2
            xη[i,j] = 0.5(x[i,j+1] - x[i,j-1])
            yη[i,j] = 0.5(y[i,j+1] - y[i,j-1])
        end
    end
    
    for i=1:nx+2
        xη[i,1] = -(3*x[i,1] - 4*x[i,2] + x[i,3])*0.5
        yη[i,1] = -(3*y[i,1] - 4*y[i,2] + y[i,3])*0.5
        xη[i, ny+2] = (3*x[i,ny+2] - 4*x[i,ny+1] + x[i,ny])*0.5
        yη[i, ny+2] = (3*y[i,ny+2] - 4*y[i,ny+1] + y[i,ny])*0.5
    end

    for j = 1:ny+2
        for i = 1:nx+2
            J[i,j] =  1.0/(xξ[i,j]*yη[i,j] - xη[i,j]*yξ[i,j])
            ξx[i,j] = J[i,j]*yη[i,j]
            ηx[i,j] = -J[i,j]*yξ[i,j]
            ξy[i,j] = -J[i,j]*xη[i,j]
            ηy[i,j] = J[i,j]*xξ[i,j]
        end
    end

    println("metrics ", ξx[1:3,1:2], ξy[1:3,1:2], ηx[1:3,1:2], ηy[1:3,1:2], J[1:3,1:2])
    println("negative volume, ", findall(x->(x<=0), J), ", ",J[findall(x->(x<=0), J)])
end

function makesource(para, bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J, U, V, Su, Sv, re⁺, re⁻, rn⁺, rn⁻,αe, αn, Ue, Vn, ψ)
    @unpack nx, ny, dt, nu, = para

    for j=1:ny+2
        for i=1:nx+2
            beta2[i,j] = nu/J[i,j]*(ξx[i,j]*ηx[i,j] + ξy[i,j]*ηy[i,j])
        end
    end

    # bua
    for j=2:ny+1
        for i=2:nx+1
            bua[i,j] = 0.125( (uold[i,j+1] - uold[i,j-1])*(beta2[i+1,j] - beta2[i-1,j]) + (uold[i+1,j] - uold[i-1,j])*(beta2[i,j+1] - beta2[i,j-1]) + uold[i+1,j+1]*(2beta2[i,j] + beta2[i+1,j] +beta2[i,j+1]) - uold[i-1,j+1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j+1]) - uold[i+1,j-1]*(2beta2[i,j] + beta2[i+1,j] + beta2[i,j-1]) + uold[i-1,j-1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j-1]))

            bua[i,j] += -0.5(ξx[i+1,j]*pa[i+1,j]/J[i+1,j] - ξx[i-1,j]*pa[i-1,j]/J[i-1,j] + ηx[i,j+1]*pa[i,j+1]/J[i,j+1] - ηx[i,j-1]*pa[i,j-1]/J[i,j-1]) + uold[i,j]/J[i,j]/dt
        end
    end

    #bva
    for j=2:ny+1
        for i=2:nx+1
            bva[i,j] = 0.125( (vold[i,j+1] - vold[i,j-1])*(beta2[i+1,j] - beta2[i-1,j]) + (vold[i+1,j] - vold[i-1,j])*(beta2[i,j+1] - beta2[i,j-1]) + vold[i+1,j+1]*(2beta2[i,j] + beta2[i+1,j] +beta2[i,j+1]) - vold[i-1,j+1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j+1]) - vold[i+1,j-1]*(2beta2[i,j] + beta2[i+1,j] + beta2[i,j-1]) + vold[i-1,j-1]*(2beta2[i,j] + beta2[i-1,j] + beta2[i,j-1]))

            bva[i,j] += -0.5(ξy[i+1,j]*pa[i+1,j]/J[i+1,j] - ξy[i-1,j]*pa[i-1,j]/J[i-1,j] + ηy[i,j+1]*pa[i,j+1]/J[i,j+1] - ηy[i,j-1]*pa[i,j-1]/J[i,j-1]) + vold[i,j]/J[i,j]/dt
        end
    end

    makeCotravariant(para, uold, vold, U, V, ξx, ξy, ηx, ηy)
    for j=2:ny
        for i=2:nx
            Ue[i,j] = 0.5*(U[i+1,j]/J[i+1,j] + U[i,j]/J[i,j])
            Vn[i,j] = 0.5*(V[i,j+1]/J[i,j+1] + V[i,j]/J[i,j])
            αe[i,j] = 1 - signbit(Ue[i,j]) #0.5*(sign(Ue[i,j]) + 1.)
            αn[i,j] = 1 - signbit(Vn[i,j]) #0.5*(sign(Vn[i,j]) + 1.)
        end
    end
    # Su
    for j=2:ny
        for i=2:nx
            re⁺[i,j] = (uold[i,j]-uold[i-1,j])/(uold[i+1,j]-uold[i,j] + 1e-10)
            re⁻[i,j] = (uold[i+2,j]-uold[i+1,j])/(uold[i+1,j]-uold[i,j]+ 1e-10)
            rn⁺[i,j] = (uold[i,j]-uold[i,j-1])/(uold[i,j+1]-uold[i,j]+ 1e-10)
            rn⁻[i,j] = (uold[i,j+2]-uold[i,j+1])/(uold[i,j+1]-uold[i,j]+ 1e-10)
        end
    end
    for j=3:ny
        for i=3:nx
            Su[i,j] = 0.5Ue[i,j]*((1-αe[i,j])*ψ(re⁻[i,j]) - αe[i,j]*ψ(re⁺[i,j]))*(uold[i+1,j] - uold[i,j]) - 
            0.5Ue[i-1,j]*((1-αe[i-1,j])*ψ(re⁻[i-1,j]) - αe[i-1,j]*ψ(re⁺[i-1,j]))*(uold[i,j] - uold[i-1,j]) +
            0.5Vn[i,j]*((1-αn[i,j])*ψ(rn⁻[i,j]) - αn[i,j]*ψ(rn⁺[i,j]))*(uold[i,j+1] - uold[i,j]) -
            0.5Vn[i,j-1]*((1-αn[i,j-1])*ψ(rn⁻[i,j-1]) - αn[i,j-1]*ψ(rn⁺[i,j-1]))*(uold[i,j] - uold[i,j-1])
        end
    end
    # Sv
    for j=2:ny
        for i=2:nx
            re⁺[i,j] = (vold[i,j]-vold[i-1,j])/(vold[i+1,j]-vold[i,j]+ 1e-10)
            re⁻[i,j] = (vold[i+2,j]-vold[i+1,j])/(vold[i+1,j]-vold[i,j]+ 1e-10)
            rn⁺[i,j] = (vold[i,j]-vold[i,j-1])/(vold[i,j+1]-vold[i,j]+ 1e-10)
            rn⁻[i,j] = (vold[i,j+2]-vold[i,j+1])/(vold[i,j+1]-vold[i,j]+ 1e-10)
        end
    end
    for j=3:ny
        for i=3:nx
            Sv[i,j] = 0.5Ue[i,j]*((1-αe[i,j])*ψ(re⁻[i,j]) - αe[i,j]*ψ(re⁺[i,j]))*(vold[i+1,j] - vold[i,j]) - 
            0.5Ue[i-1,j]*((1-αe[i-1,j])*ψ(re⁻[i-1,j]) - αe[i-1,j]*ψ(re⁺[i-1,j]))*(vold[i,j] - vold[i-1,j]) +
            0.5Vn[i,j]*((1-αn[i,j])*ψ(rn⁻[i,j]) - αn[i,j]*ψ(rn⁺[i,j]))*(vold[i,j+1] - vold[i,j]) -
            0.5Vn[i,j-1]*((1-αn[i,j-1])*ψ(rn⁻[i,j-1]) - αn[i,j-1]*ψ(rn⁺[i,j-1]))*(vold[i,j] - vold[i,j-1])
        end
    end
    for j=3:ny
        for i=3:nx
            bua[i,j] += Su[i,j]
            bva[i,j] += Sv[i,j]
        end
    end
end

function selectTVDlimitter(TVDlimitter)
    # van Albada
    function vanAlbada(r)::Float64
        return (r+r^2)/(1.0+r^2)
    end

    function vanLeer(r)
        return (r + abs(r))/(1.0 + r)
    end

    function QUICK(r)
        return max(0.0, min(2*r, (3+r)*0.25, 2.0))
    end

    function superbee(r)
        return max(0.0, min(2*r, 1.0), min(r,2.0))
    end

    if TVDlimitter == 1
        return vanAlbada
    elseif TVDlimitter == 2
        return vanLeer
    elseif TVDlimitter == 3
        return QUICK
    elseif TVDlimitter == 4
        return superbee
    end
end

function makeCotravariant(para, ua, va, U, V, ξx, ξy, ηx, ηy)
    @unpack nx, ny, = para

    for j=1:ny+2
        for i=1:nx+2
            U[i,j] = ξx[i,j]*ua[i,j] + ξy[i,j]*va[i,j]
            V[i,j] = ηx[i,j]*ua[i,j] + ηy[i,j]*va[i,j]
        end
    end
end

function solveUVA(para, ua, va, apuva, aeuva, awuva, anuva, asuva, bua, utemp, bva, vtemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
    @unpack nx, ny, dt, nu, alphaua, alphava,  niter, ipter, Uext = para
    
    makeCotravariant(para, ua, va, U, V, ξx, ξy, ηx, ηy)

    for j=1:ny+2
        for i=1:nx+2
            if i==1
                awuva[i,j] = 0
            elseif i==2 && (231<=j<=331)
                awuva[i,j] = nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i-1,j]^2 + ξy[i-1,j]^2)/J[i-1,j])
            elseif i==2 && (j<=230 || 332<=j)
                awuva[i,j] = 0.5nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i-1,j]^2 + ξy[i-1,j]^2)/J[i-1,j]) + max(0.5*(U[i,j]/J[i,j] + U[i-1,j]/J[i-1,j]), 0.0)
            else
                awuva[i,j] = 0.5nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i-1,j]^2 + ξy[i-1,j]^2)/J[i-1,j]) + max(0.5*(U[i,j]/J[i,j] + U[i-1,j]/J[i-1,j]), 0.0)
            end

            if i==nx+2
                aeuva[i,j] = 0
            elseif i==nx+1
                aeuva[i,j] = 0.5nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i+1,j]^2 + ξy[i+1,j]^2)/J[i+1,j]) + max(-0.5*(U[i+1,j]/J[i+1,j] + U[i,j]/J[i,j]), 0.0)
            else
                aeuva[i,j] = 0.5nu*((ξx[i,j]^2 + ξy[i,j]^2)/J[i,j] + (ξx[i+1,j]^2 + ξy[i+1,j]^2)/J[i+1,j]) + max(-0.5*(U[i+1,j]/J[i+1,j] + U[i,j]/J[i,j]), 0.0)
            end
            
            if j==1
                asuva[i,j] = 0
            elseif j==2
                asuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,j-1]^2 + ηy[i,j-1]^2)/J[i,j-1]) + max(0.5*(V[i,j]/J[i,j] + V[i,j-1]/J[i,j-1]), 0.0)
            else
                asuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,j-1]^2 + ηy[i,j-1]^2)/J[i,j-1]) + max(0.5*(V[i,j]/J[i,j] + V[i,j-1]/J[i,j-1]), 0.0)
            end
            if j==ny+2
                anuva[i,j] = 0
            elseif j == ny+1
                anuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,j+1]^2 + ηy[i,j+1]^2)/J[i,j+1]) + max(-0.5*(V[i,j+1]/J[i,j+1] + V[i,j]/J[i,j]), 0.0)
            else
                anuva[i,j] = 0.5nu*((ηx[i,j]^2 + ηy[i,j]^2)/J[i,j] + (ηx[i,j+1]^2 + ηy[i,j+1]^2)/J[i,j+1]) + max(-0.5*(V[i,j+1]/J[i,j+1] + V[i,j]/J[i,j]), 0.0)
            end
            
            apuva[i,j] = aeuva[i,j] + awuva[i,j] + anuva[i,j] + asuva[i,j] + 1/dt/J[i,j]
            
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            utemp[i,j] = ua[i,j]*(1.0-alphaua) + alphaua*(
                aeuva[i,j]*ua[i+1,j] + awuva[i,j]*ua[i-1,j] + anuva[i,j]*ua[i,j+1] + asuva[i,j]*ua[i,j-1] + bua[i,j])/apuva[i,j]
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            ua[i,j] = utemp[i,j]
        end
    end

    for j=2:ny+1
        for i=2:nx+1 
            vtemp[i,j] = va[i,j] * (1.0 - alphava) + alphava * (
                aeuva[i,j]*va[i+1,j] + awuva[i,j]*va[i-1,j] + anuva[i,j]*va[i,j+1] + asuva[i,j]*va[i,j-1] + bva[i,j])/apuva[i,j]
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            va[i,j] = vtemp[i,j]
        end
    end
end

function RhieChowInterpolation(para, ua, va, pa, apuva, U, V,  Ue, Vn, ξx, ξy, ηx, ηy, J)
    @unpack nx, ny, = para

    makeCotravariant(para, ua, va, U, V, ξx, ξy, ηx, ηy)

    # ue: interface velocity
    for j=2:ny+1
        for i=2:nx
            #Ue[i,j] = 0.5(U[i+1,j]+U[i,j]) + 0.5(dξ[i,j]+dξ[i+1,j])*(pa[i,j]-pa[i+1,j]) -0.25dξ[i,j]*(pa[i-1,j]-pa[i+1,j]) -0.25*dξ[i+1,j]*(pa[i,j]-pa[i+2,j])  
            Ue[i,j] =  0.5(U[i+1,j]+U[i,j]) -0.5((ξx[i,j]*ξx[i+1,j]/J[i+1,j]+ξy[i,j]*ξy[i+1,j]/J[i+1,j])/apuva[i,j] + (ξx[i+1,j]*ξx[i,j]/J[i,j]+ξy[i+1,j]*ξy[i,j]/J[i,j])/apuva[i+1,j])*(pa[i+1,j]-pa[i,j]) +
            0.25((ξx[i,j]*ξx[i+1,j]/J[i+1,j]+ξy[i,j]*ξy[i+1,j]/J[i+1,j])*pa[i+1,j] -(ξx[i,j]*ξx[i-1,j]/J[i-1,j]+ξy[i,j]*ξy[i-1,j]/J[i-1,j])*pa[i-1,j] )/apuva[i,j] +
            0.25((ξx[i+1,j]*ξx[i+2,j]/J[i+2,j]+ξy[i+1,j]*ξy[i+2,j]/J[i+2,j])*pa[i+2,j] -(ξx[i+1,j]*ξx[i,j]/J[i,j]+ξy[i+1,j]*ξy[i,j]/J[i,j])*pa[i,j] )/apuva[i+1,j]
        end
    end

    for j=2:ny+1
        Ue[1,j] = U[1,j]
        
        Ue[nx+1,j] = U[nx+2,j] 
    end

    for j=2:ny
        for i=2:nx+1
            #Vn[i,j] = 0.5(V[i,j+1]+V[i,j]) + 0.5(dη[i,j]+dη[i,j+1])*(pa[i,j]-pa[i,j+1]) -0.25dη[i,j]*(pa[i,j-1]-pa[i,j+1]) -0.25dη[i,j+1]*(pa[i,j]-pa[i,j+2])  
            Vn[i,j] = 0.5(V[i,j+1]+V[i,j]) -0.5((ηx[i,j]*ηx[i,j+1]/J[i,j+1]+ηy[i,j]*ηy[i,j+1]/J[i,j+1])/apuva[i,j] + (ηx[i,j+1]*ηx[i,j]/J[i,j]+ηy[i,j+1]*ηy[i,j]/J[i,j])/apuva[i,j+1] )*(pa[i,j+1]-pa[i,j]) + 
            0.25((ηx[i,j]*ηx[i,j+1]/J[i,j+1]+ηy[i,j]*ηy[i,j+1]/J[i,j+1])*pa[i,j+1] - (ηx[i,j]*ηx[i,j-1]/J[i,j-1]+ηy[i,j]*ηy[i,j-1]/J[i,j-1])*pa[i,j-1])/apuva[i,j] +
            0.25((ηx[i,j+1]*ηx[i,j+2]/J[i,j+2]+ηy[i,j+1]*ηy[i,j+2]/J[i,j+2])*pa[i,j+2] - (ηx[i,j+1]*ηx[i,j]/J[i,j]+ηy[i,j+1]*ηy[i,j]/J[i,j])*pa[i,j])/apuva[i,j+1]
        end
    end

    for i=2:nx+1
        Vn[i, ny+1] = V[i, ny+2]
        Vn[i, 1] = V[i, 1]
    end

end

function solvePC(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, 
        bva, appc, aepc, awpc, anpc, aspc, bpc, ptemp, Ue, Vn, ξx, ξy, ηx, ηy, J, st)
    @unpack nx, ny, nu, alphaua, alphava, alphapa, niter, ipter, nu = para

    for j=2:ny+1
        for i=2:nx+1
            pc[i,j] = 0.0
        end
    end
    
    for j=2:ny+1
        for i=2:nx+1
            if i == 2
                awpc[i,j] =0
            else
                awpc[i,j] = 0.5((ξx[i-1,j]*ξx[i,j]/J[i,j]+ξy[i-1,j]*ξy[i,j]/J[i,j])/apuva[i-1,j] + (ξx[i,j]*ξx[i-1,j]/J[i-1,j]+ξy[i,j]*ξy[i-1,j]/J[i-1,j])/apuva[i,j])/(0.5(J[i-1,j]+J[i,j]))
            end
            
            if i == nx+1
                aepc[i,j] = 0
            else          
                aepc[i,j] = 0.5((ξx[i,j]*ξx[i+1,j]/J[i+1,j]+ξy[i,j]*ξy[i+1,j]/J[i+1,j])/apuva[i,j] + (ξx[i+1,j]*ξx[i,j]/J[i,j]+ξy[i+1,j]*ξy[i,j]/J[i,j])/apuva[i+1,j])/(0.5(J[i,j]+J[i+1,j]))
            end
            
            if j == 2
                aspc[i,j] = 0
            else
                aspc[i,j] = 0.5((ηx[i,j-1]*ηx[i,j]/J[i,j]+ηy[i,j-1]*ηy[i,j]/J[i,j])/apuva[i,j-1] + (ηx[i,j]*ηx[i,j-1]/J[i,j-1]+ηy[i,j]*ηy[i,j-1]/J[i,j-1])/apuva[i,j] )/(0.5(J[i,j-1]+J[i,j]))
            end
                    
            if j == ny+1
                anpc[i,j] = 0
            else
                anpc[i,j] = 0.5((ηx[i,j]*ηx[i,j+1]/J[i,j+1]+ηy[i,j]*ηy[i,j+1]/J[i,j+1])/apuva[i,j] + (ηx[i,j+1]*ηx[i,j]/J[i,j]+ηy[i,j+1]*ηy[i,j]/J[i,j])/apuva[i,j+1] )/(0.5(J[i,j]+J[i,j+1]))
            end
                        
            appc[i,j] = aepc[i,j] + awpc[i,j] + anpc[i,j] + aspc[i,j]
            bpc[i,j] = -(Ue[i,j]/(0.5(J[i,j]+J[i+1,j])) - Ue[i-1,j]/(0.5(J[i,j]+J[i-1,j]))) - (Vn[i,j]/(0.5(J[i,j]+J[i,j+1])) - Vn[i,j-1]/(0.5(J[i,j]+J[i,j-1])))
        end
    end

    
    # SOR法で圧力の解を求めてる
    for ip=1:ipter
        for j=2:ny+1
            for i=2:nx+1
                ptemp[i,j] = pc[i,j]
            end
        end
        for j=2:ny+1
            for i=2:nx+1
                pc[i,j] = pc[i,j]*(1.0 - alphapa) + alphapa*(
                    aepc[i,j]*pc[i+1,j] + awpc[i,j]*pc[i-1,j] + anpc[i,j]*pc[i,j+1] + aspc[i,j]*pc[i,j-1] + bpc[i,j])/appc[i,j]
            end
        end

        for j=232:330
            pc[1, j] = pc[2,j] 
        end
    
        for j=1:231
            pc[1, j] = 0.5*(pc[2,j]-st[j]*(pc[3,j]-pc[2,j]) + pc[2,ny+3-j]-st[ny+3-j]*(pc[3,ny+3-j]-pc[2,ny+3-j]))
            pc[1, ny+3-j] = pc[1,j]           
        end
        
        # outlet
        for i=1:nx+2
            pc[i, 1] = 0 
            pc[i, ny+2] = 0 
        end
        for j=1:231
            pc[nx+2,j] = 0
        end
        for j=331:ny+2
            pc[nx+2,j] = 0
        end
    
        # inlet
        for j=232:330
            pc[nx+2,j] = pc[nx+1, j]
        end

        rest=convergP(para, pc, ptemp)
        if rest < 1e-6 && ip >=2
            #println("converged rest PC ", rest, ", ip=", ip)
            break
        end
        if ip ==ipter
            println("rest PC ", rest, ", ip=", ip)
        end
    end
end    

function update(para, ua, va, pc, pa, apuva, ua2, va2, ξx, ξy, ηx, ηy, J)
    @unpack nx, ny, alphapa,  = para
    
    # u = u* + u'
    for j=2:ny+1
        for i=2:nx+1
            ua2[i,j] = ua[i,j] -0.5(ξx[i+1,j]*pc[i+1,j]/J[i+1,j] - ξx[i-1,j]*pc[i-1,j]/J[i-1,j] + ηx[i,j+1]*pc[i,j+1]/J[i,j+1] - ηx[i,j-1]*pc[i,j-1]/J[i,j-1])/apuva[i,j]
            
            va2[i,j] = va[i,j] -0.5(ξy[i+1,j]*pc[i+1,j]/J[i+1,j] - ξy[i-1,j]*pc[i-1,j]/J[i-1,j] + ηy[i,j+1]*pc[i,j+1]/J[i,j+1] - ηy[i,j-1]*pc[i,j-1]/J[i,j-1])/apuva[i,j]
        end
    end
    
    #pref = pa[2,2]
    #for j=1:ny+2
    #    for i=1:nx+2
    #        pa[i,j] -= pref
    #    end
    #end
end

function solvePC2(para, ua, ua2, va, va2, pc, pc2, apuva, aeuva, awuva, anuva, asuva, bua, bva, appc, aepc, awpc, anpc, aspc, bpc, UE, UW, VN, VS, ptemp, ξx, ξy, ηx, ηy, U, V, J, st)
    @unpack nx, ny, ipter2, alphapa, = para

    makeCotravariant(para, ua2 .- ua, va2 .- va, U, V, ξx, ξy, ηx, ηy)

    for j=2:ny+1
        for i=2:nx+1
            pc2[i,j] = 0.0
            bpc[i,j] = 0.0
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            if i==nx+1
                UE[i,j] = awuva[i+1,j]*(U[i,j])
            else
                UE[i,j] = aeuva[i+1,j]*(U[i+2,j]) +awuva[i+1,j]*(U[i,j]) +anuva[i+1,j]*(U[i+1,j+1]) +asuva[i+1,j]*(U[i+1,j-1])
            end

            if i==2
                UW[i,j] = aeuva[i-1,j]*(U[i,j]) 
            else
                UW[i,j] = aeuva[i-1,j]*(U[i,j]) +awuva[i-1,j]*(U[i-2,j]) +anuva[i-1,j]*(U[i-1,j+1]) +asuva[i-1,j]*(U[i-1,j-1])
            end

            if j==ny+1
                VN[i,j] = asuva[i,j+1]*(V[i,j])
            else
                VN[i,j] = aeuva[i,j+1]*(V[i+1,j+1]) +awuva[i,j+1]*(V[i-1,j+1]) +anuva[i,j+1]*(V[i,j+2]) +asuva[i,j+1]*(V[i,j])
            end

            if j==2
                VS[i,j] = anuva[i,j-1]*(V[i,j]) 
            else
                VS[i,j] = aeuva[i,j-1]*(V[i+1,j-1]) +awuva[i,j-1]*(V[i-1,j-1]) +anuva[i,j-1]*(V[i,j]) +asuva[i,j-1]*(V[i,j-2])
            end
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            bpc[i,j] = -0.5(UE[i,j]/apuva[i+1,j]/(0.5(J[i+1,j]+J[i,j])) - UW[i,j]/apuva[i-1,j]/(0.5(J[i,j]+J[i-1,j]))) -0.5(VN[i,j]/apuva[i,j+1]/(0.5(J[i,j]+J[i,j+1])) -VS[i,j]/apuva[i,j-1]/(0.5(J[i,j]+J[i,j-1])))           
        end
    end

    #println("ue, uw, vn, vs ", findmax(abs.(UE)), ",", findmax(abs.(UW)), ", ", findmax(abs.(VN)), ",", findmax(abs.(VS)))
    #println("apuva ", findmax(apuva))
    #println("bpc ",findmax(abs.(bpc)))

    # SOR法で圧力の解を求めてる
    for ip=1:ipter2
        for j=2:ny+1
            for i=2:nx+1
                ptemp[i,j] = pc2[i,j]
            end
        end
        for j=2:ny+1
            for i=2:nx+1
                pc2[i,j] = pc2[i,j]*(1.0 - alphapa) + alphapa*(aepc[i,j]*pc2[i+1,j] + awpc[i,j]*pc2[i-1,j] + anpc[i,j]*pc2[i,j+1] + aspc[i,j]*pc2[i,j-1] + bpc[i,j])/appc[i,j]
            end
        end
    
        for j=232:330
            pc2[1, j] = pc2[2,j] 
        end
    
        for j=1:231
            pc2[1, j] = 0.5*(pc2[2,j]-st[j]*(pc2[3,j]-pc2[2,j]) + pc2[2,ny+3-j]-st[ny+3-j]*(pc2[3,ny+3-j]-pc2[2,ny+3-j]))
            pc2[1, ny+3-j] = pc2[1, j]          
        end
        
        # outlet
        for i=1:nx+2
            pc2[i, 1] = 0 
            pc2[i, ny+2] = 0 
        end
        for j=1:231
            pc2[nx+2,j] = 0
        end
        for j=331:ny+2
            pc2[nx+2,j] = 0
        end
    
        # inlet
        for j=232:330
            pc2[nx+2,j] = pc2[nx+1, j]
        end

        rest=convergP(para, pc2, ptemp)
        if rest < 1e-7 && ip>=2
            #println("converged rest PC2 ", rest, ", ip=", ip)
            break
        end
        if ip == ipter2
            println("rest PC2 ", rest, ", ip=", ip)
        end
    end
end 

function update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apuva, aeuva, awuva, anuva, asuva, ξx, ξy, ηx, ηy, J, utemp, vtemp)
    @unpack nx, ny,  = para

    # u** = u* + u'
    for j=2:ny+1
        for i=2:nx+1
            ua3[i,j] = ua2[i,j] + (aeuva[i,j]*(ua2[i+1,j]-ua[i+1,j]) + awuva[i,j]*(ua2[i-1,j]-ua[i-1,j]) + anuva[i,j]*(ua2[i,j+1]-ua[i,j+1]) + asuva[i,j]*(ua2[i,j-1]-ua[i,j-1]) -
            0.5(ξx[i+1,j]*pc2[i+1,j]/J[i+1,j] - ξx[i-1,j]*pc2[i-1,j]/J[i-1,j] + ηx[i,j+1]*pc2[i,j+1]/J[i,j+1] - ηx[i,j-1]*pc2[i,j-1]/J[i,j-1]))/apuva[i,j]
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            ua[i,j] = ua3[i,j]
        end
    end
    
    for j=2:ny+1
        for i=2:nx+1
            va3[i,j] = va2[i,j] + (aeuva[i,j]*(va2[i+1,j]-va[i+1,j]) + awuva[i,j]*(va2[i-1,j]-va[i-1,j]) + anuva[i,j]*(va2[i,j+1]-va[i,j+1]) + asuva[i,j]*(va2[i,j-1]-va[i,j-1]) -
            0.5(ξy[i+1,j]*pc2[i+1,j]/J[i+1,j] - ξy[i-1,j]*pc2[i-1,j]/J[i-1,j] + ηy[i,j+1]*pc2[i,j+1]/J[i,j+1] - ηy[i,j-1]*pc2[i,j-1]/J[i,j-1]))/apuva[i,j]
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            va[i,j] = va3[i,j]
        end
    end

    #println("max ua3-ua2=",findmax(abs.(ua3[2:nx,2:ny+1]-ua2[2:nx,2:ny+1])))
    #println("max va3-va2=",findmax(abs.(va3[2:nx+1,2:ny]-va2[2:nx+1,2:ny])))
    for j=2:ny+1
        for i=2:nx+1
            pa[i,j] = pa[i,j] + pc[i,j] + pc2[i,j]
        end
    end
end

function converg(para, ua, va, utemp, vtemp)
    @unpack nx, ny = para
    
    restua = 0.0
    restva = 0.0
    
    for j=2:ny+1
        for i=2:nx+1
            restua = restua + abs(ua[i,j] - utemp[i,j])
            
        end
    end
    for j=2:ny+1
        for i=2:nx+1
            restva = restva + abs(va[i,j] - vtemp[i,j])
        end
    end
    rest = (restua + restva)/2.0
    return rest/(nx*ny)
end

function convergP(para, pc, ptemp)
    @unpack nx, ny = para
    
    rest = 0.0
    
    for j=2:ny+1
        for i=2:nx+1
            rest = rest + abs(pc[i,j] - ptemp[i,j])
        end
    end

    return rest/(nx*ny)
end

function divergence(para, ua, va, U, V, ξx, ξy, ηx, ηy, J) 
    @unpack nx, ny = para

    makeCotravariant(para, ua, va, U, V, ξx, ξy, ηx, ηy)
    divu = 0.0

    for j=2:ny+1
        for i=2:nx+1
            divu += abs((0.5*(U[i+1,j] + U[i,j])/(0.5(J[i,j]+J[i+1,j])) - 0.5*(U[i,j] + U[i-1,j])/(0.5(J[i-1,j]+J[i,j]))) + (0.5(V[i,j+1]+V[i,j])/(0.5(J[i,j]+J[i,j+1])) - 0.5(V[i,j]+V[i,j-1])/(0.5(J[i,j]+J[i,j-1]))))
        end
    end

    return divu/(nx*ny)

end 

function output(para, ua, va, pa,  l, omega)
    @unpack outfile, dt, nx, ny = para
    #println("time=", time)
    filename = outfile*string(l)*".h5"
    h5open(filename, "w") do file
        write(file, "U", ua)
        write(file, "V", va)
        write(file, "P", pa)
        write(file, "omega", omega)
        write(file, "n", l)
    end

end

function restart(para, restartFile, ua, va, pa, )
    @unpack nx, ny,  = para
    
    file = h5open(restartFile, "r") 
    ua .= read(file, "U")
    va .= read(file, "V")
    pa .= read(file, "P")
    nstart = read(file, "n")
    close(file)
    
    return nstart
end

function vor(para, omega, ua, va, ξx, ξy, ηx, ηy, st)
    @unpack nx, ny,  = para

    for j=2:ny+1
        for i=2:nx+1
            omega[i,j] = (ξx[i,j]*0.5*(va[i+1,j] -va[i-1,j]) + ηx[i,j]*0.5*(va[i,j+1] -va[i,j-1])) - 
            (ξy[i,j]*0.5*(ua[i+1,j] -ua[i-1,j]) + ηy[i,j]*0.5*(ua[i,j+1] -ua[i,j-1]))
        end    
    end
    for j=1:231
        omega[1, j] = 0.5*(omega[2,j]-st[j]*(omega[3,j]-omega[2,j]) + omega[2,ny+3-j]-st[ny+3-j]*(omega[3,ny+3-j]-omega[2,ny+3-j]))
        omega[1, ny+3-j] = omega[1,j]
    end
    for j=232:330
        i=1
        omega[i,j] = (ξx[i,j]*(va[i+1,j] -va[i,j]) + ηx[i,j]*0.5*(va[i,j+1] -va[i,j-1])) - 
            (ξy[i,j]*(ua[i+1,j] -ua[i,j]) + ηy[i,j]*0.5*(ua[i,j+1] -ua[i,j-1]))
    end
    
end

実行は

para = Param()
@time main(para)

Mac miniで10分くらいで終わります。