コロケート格子の一般直交座標のPISO法のコード（Julia）

お詫び

Jacobianの計算がおかしかったのと，速度の収束判定がおかしかったのを修正しました。
ついでに格子点の生成も修正しました。

----------


普通有限体積法のコードだと，一般直交座標とかは使わずに，四面体とかで複雑形状に適用すると思います。それはややこしいので，ここでは一般直交座標を適用します。コロケート格子のPISO法のコードのベースはここにあります。　

https://qiita.com/ishigaki/items/dec6d42f78b3f8a06bad

これに一般直交座標に適用します。カルマン渦の計算しました。下記の結果を出します。
![uv-vec143.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/43609/92a9468e-c4b9-72e4-4761-024ac5277a75.jpeg)


# 基礎式
説明ややこしいので，下の画像を見てください。基本的には一般直交座標でのSIMPLE法の説明です。RhieとChowの論文が参考になります。
Numerical study of the turbulent flow past an airfoil with trailing edge separation
C. M. Rhie and W. L. Chow
AIAA Journal 1983 21:11, 1525-1532

![1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/43609/5ed98f45-07c1-59b6-39c4-e9b1443ebc80.png)
![2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/43609/54bfaff8-4dd5-1b8e-d5e0-ab0b2ef7677a.png)
![3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/43609/b7581910-88e9-9093-e030-30a2d1a0aec8.png)
![4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/43609/98db68b7-1513-f43d-5f36-7e9f783554e4.png)
![5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/43609/18d083ef-4599-2f33-24ec-754302a95ed1.png)
![6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/43609/3188fadc-39e5-88e7-4ab2-b2dd58000bfc.png)
![7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/43609/551c8001-e2b2-e0fa-4301-21244827090a.png)
![8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/43609/83789e8f-6011-063b-4b71-a6e594c8eb43.png)

# グリッドの生成
まずはグリッドを作ります。グリッド生成のコード。周期境界条件を適用するため，j方向のグリッドが一部重なります。

高見，河村「偏微分方程式の差分解法」を参考にしました。

```julia
using Parameters
using HDF5
using PyPlot

@with_kw struct Param{T1,T2}
    R1::T1 = 0.5
    R2::T1 = 10.0
    xr::T1 = 3.
    ytop::T1 = 3.
    
    ε::T1 = 1e-6
    itermax::T2 = 100
    
    nx::T2 = 150
    ny::T2 = 146
    
    
    outputname::String = "Grid/grid-largecirc2"
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, R2, xr, ytop = para
    
    # initial grid
    for j = 1:ny
        r = R1
        θ = 2pi*(j-1)/(ny-2)
        x[1, j] = r*cos(θ) 
        y[1, j] = r*sin(θ)
    end 

    
    for j = 1:ny
        for i = 2:nx
            r = R1 + (R2-R1)*(i-1.5)/(nx-1)
            θ = 2pi*(j-1)/(ny-2)
            x[i, j] = r*cos(θ) 
            y[i, j] = r*sin(θ) 
        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)
        write(file, "y", y)
        write(file, "xT", xT)
        write(file, "yT", yT)
    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.05
    
    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
        
        for i = 1:nx
            x[i, 1] = x[i, ny-1]
            y[i, 1] = y[i, ny-1]
            x[i, ny] = x[i, 2]
            y[i, ny] = y[i, 2]
        end
       
        if error < ε
            println("iteration, error=", i, ", ", error)
            break
        elseif i==itermax
            println("error= ",error)
        end 
    end
end
```

実行と可視化

```julia
para = Param(nx=150, ny=146)
gridgen(para)

@unpack outputname = para
    
filename = outputname*".h5"

file = h5open(filename, "r") 
x = read(file, "x")
y = read(file, "y")
xT = read(file, "xT")
yT = read(file, "yT")
nx = read(file, "nx")
ny = read(file, "ny")
close(file)

fig = plt.figure(figsize = (13, 13))
ax1 = fig.add_subplot(111)

ax1.plot(xT[:,:], yT[:,:], "*" )
#ax1.plot(x[:,:], y[:,:], color="black")
#ax1.plot(xT[:,:], yT[:,:], color="black")

#ax1.set_xlim(-1,1)
#ax1.set_ylim(-1,1)
```

こんなグリッドです。
![grid.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/43609/85c9d163-faba-ffb8-e9c3-7df9cc31eb7d.png)

# 流体コード
流体の解析は以下です。グリッドの点数を表す，nx, nyと数字がずれてます。これは流体は1ーnx+2, 1-ny+2まで点数があるためです。

できたばかりなので，間違っているところが含まれる可能性があるので，注意してください。

```julia


using Parameters

using HDF5

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

    nstep::T2 = 2000
    alphaua::T1 = 0.9
    alphava::T1 = 0.9
    alphapa::T1 = 0.5
    niter::T2 = 500
    ipter::T2 = 25
    ipter2::T2 = 100

    start::T2 =1

    # mesh size
    nx::T2 = 147
    ny::T2 = 144

    ndata::T2 = 50

    gridfile::String = "Grid/grid-rectcirc4.h5"
    outfile::String = "Data-piso-Kar/uv" 

end

function main(para)
    @unpack nx, ny, nstep, niter, ndata, Uext = 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 = ones(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)

    ξ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)

    rest = 0.0

    readgrid(para, x, y)
    jacobian(para, x, y, ξx, ξy, ηx, ηy, J)


    for l = 1:nstep
        BC(para, ua, va, pc, pa, ξx, ξy,  ηx, ηy, J, uold, vold, ua2, va2, pc2, l)
        
        for n = 1:niter
            BC(para, ua, va, pc, pa, ξx, ξy,  ηx, ηy, J, uold, vold, ua2, va2, pc2, l)
            makesource(bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J)
            
            solveUA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
            solveVA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, Vn, ξx, ξy, ηx, ηy, J, U, V)
            RhieChowInterpolation(para, utemp, vtemp, 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,)
            update(para, utemp, vtemp, pc, pa, apuva, ua2, va2, ξx, ξy, ηx, ηy, J)
            solvePC2(para, utemp, ua2, vtemp, 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)
            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)
            if rest <=1e-7 && n >=2
                if l%1 == 0
                    println(l, " th step, ", n, " th iteration, residual=", rest)
                end
                break
            end
        end

        for j=1:ny+2
            for i=1:nx+2
                uold[i,j] = ua[i,j]
            end
        end
    
        for j=2:ny+1
            for i=2:nx+1
                vold[i,j] = va[i,j]
            end
        end
        if l%ndata == 0
            println("residual=", rest)
            vor(para, omega, ua, va, ξx, ξy, ηx, ηy)
            output(para, ua, va, pa, uold, vold, x, y, l, omega)
        end
    end
    #streamfunction(para, u, v, psi)
    
end

function main(para, restartFile)
    @unpack nx, ny, nstep, niter, ndata = 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)

    ξ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)
    
    rest = 0.0

    readgrid(para, x, y)
    jacobian(para, x, y, ξx, ξy, ηx, ηy, J)

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

    for l = nstart+1:nstart+nstep
        BC(para, ua, va, pc, pa, ξx, ξy,  ηx, ηy, J, uold, vold, ua2, va2, pc2, l)
        
        for n = 1:niter
            BC(para, ua, va, pc, pa, ξx, ξy,  ηx, ηy, J, uold, vold, ua2, va2, pc2, l)
            makesource(bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J)
            
            solveUA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
            solveVA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, Vn, ξx, ξy, ηx, ηy, J, U, V)
            RhieChowInterpolation(para, utemp, vtemp, 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,)
            update(para, utemp, vtemp, pc, pa, apuva, ua2, va2, ξx, ξy, ηx, ηy, J)
            solvePC2(para, utemp, ua2, vtemp, 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)
            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)
 
            if rest <=1e-7 && n >=2
                if l%1 == 0
                    println(l, " th step, ", n, " th iteration, residual=", rest)
                end
                break
            end
        end
        for j=1:ny+2
            for i=1:nx+2
                uold[i,j] = ua[i,j]
            end
        end
    
        for j=2:ny+1
            for i=2:nx+1
                vold[i,j] = va[i,j]
            end
        end
        if l%ndata == 0
            println("residual=", rest)
            vor(para, omega, ua, va, ξx, ξy, ηx, ηy)
            output(para, ua, va, pa, uold, vold, x, y, 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)
    @unpack nx, ny, Uext, dt, start = para

    for j=1:ny+2
        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
    
    # outlet
    for j=1:2div(ny,8) +1
        ua[nx+2,j] = ua[nx+1,j]
        va[nx+2,j] = va[nx+1,j]
        ua2[nx+2,j] = ua2[nx+1,j]
        va2[nx+2,j] = va2[nx+1,j]
        pa[nx+2, j] = 0 #pa[nx+1, j]
        pc[nx+2, j] = 0 #pc[nx+1, j]
        pc2[nx+2, j] = 0
    end

    # inlet
    Ubound = min(Uext*l/start, Uext)
    
    for j=2div(ny,8)+2:ny -2div(ny,8)
        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

    # outlet
    for j=ny+1-2div(ny,8):ny+2
        ua[nx+2,j] = ua[nx+1,j]
        va[nx+2,j] = va[nx+1,j]
        ua2[nx+2,j] = ua2[nx+1,j]
        va2[nx+2,j] = va2[nx+1,j]
        pa[nx+2, j] = 0 #pa[nx+1, j]
        pc[nx+2, j] = 0 #pc[nx+1, j]
        pc2[nx+2, j] = 0
    end

    # periodic condition
    for i=1:nx+2
        ua[i, ny+2] = ua[i, 2]
        va[i, ny+2] = va[i, 2]
        pc[i, ny+2] = pc[i, 2]
        pa[i, ny+2] = pa[i, 2]
        ua2[i, ny+2] = ua2[i, 2]
        va2[i, ny+2] = va2[i, 2]
        pc2[i, ny+2] = pc2[i, 2]

        ua[i, 1] = ua[i, ny+1]
        va[i, 1] = va[i, ny+1]
        pc[i, 1] = pc[i, ny+1]        
        pa[i, 1] = pa[i, ny+1]
        ua2[i, 1] = ua2[i, ny+1]
        va2[i, 1] = va2[i, ny+1]
        pc2[i, 1] = pc2[i, ny+1]    
    end
end

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

    file = h5open(gridfile, "r") 
    x .= read(file, "x")
    y .= read(file, "y")
    close(file)
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
        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
        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] = 0.5(x[i,2] - x[i,ny])
        yη[i,1] = 0.5(y[i,2] - y[i,ny])
        xη[i, ny+2] = 0.5(x[i,3] - x[i,ny+1])
        yη[i, ny+2] = 0.5(y[i,3] - y[i,ny+1])
    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])
end

function makesource(bua, bva, uold, vold, pa, beta2, gam, ξx, ξy, ηx, ηy, J)
    @unpack nx, ny, dt, nu, = para

    for j=2:ny+1
        for i=2:nx+1
            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
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 solveUA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bua, uold, utemp, Ue, ξx, ξy, ηx, ηy, J, U, V)
    @unpack nx, ny, dt, nu, alphaua, alphava, alphapa, 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
                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])
            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]) + 0.25*(U[i,j]/J[i,j] + U[i-1,j]/J[i-1,j])
            end
            if i==nx+2
                aeuva[i,j] = 0
            elseif i==nx+1
                aeuva[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]) - 0.25*(U[i+1,j]/J[i+1,j] + U[i,j]/J[i,j])
            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]) - 0.25*(U[i+1,j]/J[i+1,j] + U[i,j]/J[i,j])
            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,ny+1]^2 + ηy[i,ny+1]^2)/J[i,ny+1]) + 0.25*(V[i,j]/J[i,j] + V[i,ny+1]/J[i,ny+1])
            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]) + 0.25*(V[i,j]/J[i,j] + V[i,j-1]/J[i,j-1])
            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,2]^2 + ηy[i,2]^2)/J[i,2]) - 0.25*(V[i,2]/J[i,2] + V[i,j]/J[i,j])
            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]) - 0.25*(V[i,j+1]/J[i,j+1] + V[i,j]/J[i,j])
            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=1:ny+2
        utemp[1,j] = ua[1,j]
        utemp[nx+2,j] = ua[nx+2,j]
    end
    for i=1:nx+2
        utemp[i,1] = ua[i,1]
        utemp[i,ny+2] = ua[i,ny+2]
    end

    #for j=2:ny+1
    #    for i=2:nx+1
    #        dξ[i,j] = (ξx[i,j]^2+ξy[i,j]^2)/apuva[i,j]/J[i,j]
    #        dη[i,j] = (ηx[i,j]^2+ηy[i,j]^2)/apuva[i,j]/J[i,j]
    #    end
    #end
end

function solveVA(para, ua, va, pc, pa, apuva, aeuva, awuva, anuva, asuva, bva, vold, vtemp, Vn, ξx, ξy, ηx, ηy, J, U, V)
    @unpack nx, ny, dt, nu, alphaua, alphava, alphapa, niter, ipter,  = para
    
    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=1:ny+2
        vtemp[1,j] = va[1,j]
        vtemp[nx+2,j] = va[nx+2,j]
    end
    for i=1:nx+2
        vtemp[i,1] = va[i,1]
        vtemp[i,ny+2] = va[i,ny+2]
    end
end

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

    makeCotravariant(para, utemp, vtemp, 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]
    end

    for j=2:ny+1
        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
        j=ny+1
        #Vn[i,j] = 0.5(V[i,2]+V[i,j]) + 0.5(dη[i,j]+dη[i,2])*(pa[i,j]-pa[i,2]) -0.25dη[i,j]*(pa[i,j-1]-pa[i,2]) -0.25dη[i,2]*(pa[i,j]-pa[i,3])    
        Vn[i,j] = 0.5(V[i,2]+V[i,j]) -0.5((ηx[i,j]*ηx[i,2]/J[i,2]+ηy[i,j]*ηy[i,2]/J[i,2])/apuva[i,j] + (ηx[i,2]*ηx[i,j]/J[i,j]+ηy[i,2]*ηy[i,j]/J[i,j])/apuva[i,2] )*(pa[i,2]-pa[i,j]) + 
        0.25((ηx[i,j]*ηx[i,2]/J[i,2]+ηy[i,j]*ηy[i,2]/J[i,2])*pa[i,2] - (η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,2]*ηx[i,3]/J[i,2]+ηy[i,2]*ηy[i,3]/J[i,3])*pa[i,3] - (ηx[i,2]*ηx[i,j]/J[i,j]+ηy[i,2]*ηy[i,j]/J[i,j])*pa[i,j])/apuva[i,2]
    end
    for i=2:nx+1
        Vn[i,1] = Vn[i,ny+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)
    @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.5((ηx[i,ny+1]*ηx[i,j]/J[i,j]+ηy[i,ny+1]*ηy[i,j]/J[i,j])/apuva[i,ny+1] + (ηx[i,j]*ηx[i,ny+1]/J[i,ny+1]+ηy[i,j]*ηy[i,ny+1]/J[i,ny+1])/apuva[i,j] )/(0.5(J[i,ny+1]+J[i,j]))
            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.5((ηx[i,j]*ηx[i,2]/J[i,2]+ηy[i,j]*ηy[i,2]/J[i,2])/apuva[i,j] + (ηx[i,2]*ηx[i,j]/J[i,j]+ηy[i,2]*ηy[i,j]/J[i,j])/apuva[i,2] )/(0.5(J[i,j]+J[i,2]))
            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=1:ny+2
            pc[1,j] = pc[2,j] 
        end

        for j=1:2div(ny,8) +1
            pc[nx+2,j] = 0 #pc[nx+1, j]
        end

        for j=2div(ny,8)+2:ny -2div(ny,8)
        #for i=1:nx+2
            pc[nx+2,j] = pc[nx+1, j] 
        end

        for j=ny+1-2div(ny,8):ny+2
            pc[nx+2,j] = 0 #pc[nx+1, j]
        end

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

        rest=convergP(para, pc, ptemp)
        if rest < 1e-5 && 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, utemp, vtemp, 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] = utemp[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] = vtemp[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, utemp, ua2, vtemp, 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)
    @unpack nx, ny, ipter2, alphapa, nu = para

    makeCotravariant(para, ua2 .- utemp, va2 .- vtemp, 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 i=1:nx+2
            pc2[i, 1] = pc2[i, ny+1]
            pc2[i, ny+2] = pc2[i, 2]
        end

        
        for j=1:ny+2
            pc2[1, j] = pc2[2, j] 
        end
        for j=1:2div(ny,8) +1
            pc2[nx+2, j] = 0
        end
        for j=2div(ny,8)+2:ny -2div(ny,8)
            pc2[nx+2,j] = pc2[nx+1, j]
        end
        for j=ny+1-2div(ny,8):ny+2
            pc2[nx+2, j] = 0
        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]-utemp[i+1,j]) + awuva[i,j]*(ua2[i-1,j]-utemp[i-1,j]) + anuva[i,j]*(ua2[i,j+1]-utemp[i,j+1]) + asuva[i,j]*(ua2[i,j-1]-utemp[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]-vtemp[i+1,j]) + awuva[i,j]*(va2[i-1,j]-vtemp[i-1,j]) + anuva[i,j]*(va2[i,j+1]-vtemp[i,j+1]) + asuva[i,j]*(va2[i,j-1]-vtemp[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 output(para, ua, va, pa, uold, vold, x, y, 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, "Psi", psi[2:nx+1, 2:ny+1])
        write(file, "X", x[1:nx+2, 1:ny+2])
        write(file, "Y", y[1:nx+2, 1:ny+2])
        write(file, "ua_old", uold)
        write(file, "va_old", vold)
        write(file, "n", l)
    end
    
end

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

function vor(para, omega, ua, va, ξx, ξy, ηx, ηy)
    @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 i=2:nx+1
        j=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,ny])) - 
        (ξy[i,j]*0.5*(ua[i+1,j] -ua[i-1,j]) + ηy[i,j]*0.5*(ua[i,j+1] -ua[i,ny]))
        j=ny+2
        omega[i,j] = (ξx[i,j]*0.5*(va[i+1,j] -va[i-1,j]) + ηx[i,j]*0.5*(va[i,3] -va[i,j-1])) - 
        (ξy[i,j]*0.5*(ua[i+1,j] -ua[i-1,j]) + ηy[i,j]*0.5*(ua[i,3] -ua[i,j-1]))
    end
end

#para = Param(nu=0.02, nstep=500, outfile="Data-simple/uv", nx=30, ny=30, ndata=50, alphapa=0.5)
#para = Param(nu=0.02, nstep=1, outfile="Data-simple/uv", nx=30, ny=30, ndata=1, niter=5000)
 
#@time main(para)

        

```

実行と可視化

```julia
include("Unsteady-PISO_colocation-gene-KarmanBC.jl")
para = Param(nu=0.002, nx=148, ny=144,  gridfile="Grid/grid-largecirc2.h5", outfile="Data-piso-Kar-largeCircle/uv", ndata=50, niter=500)

@time main(para)

using HDF5
using PyPlot

dir="Data-piso-Kar-largeCircle/"
name = "uv7000"
filename = dir*name*".h5"
file = h5open(filename, "r") 
u = read(file, "U")
v = read(file, "V")
x = read(file, "X")
y = read(file, "Y")
pa = read(file, "P")
omega = read(file, "omega")
close(file)

skipx = 4
skipy = 1

X = x[:,1]
Y = y[1,:]
fig = plt.figure(figsize = (15, 20))

ax1 = fig.add_subplot(321)
heatmap = ax1.pcolormesh(x, y, u, cmap="jet", )
fig.colorbar(heatmap, ax=ax1)
ax1.set_title("u")
#heatmap.set_clim(vmin=0, vmax=1)
ax1.set_xlim(-3,9)
ax1.set_ylim(-6,6)

ax2 = fig.add_subplot(322)
heatmap2 = ax2.pcolormesh(x, y, v, cmap="jet", )
fig.colorbar(heatmap2, ax=ax2)
ax2.set_title("v")
#heatmap2.set_clim(vmin=-1, vmax=1)
ax2.set_xlim(-3,9)
ax2.set_ylim(-6,6)

ax3 = fig.add_subplot(323)
umag = sqrt.(u.*u .+ v.*v)
#image = ax3.quiver(x[:,40:40], y[:,40:40], u[:,40:40], v[:,40:40], scale=4, cmap="jet")
image = ax3.quiver(x[1:skipx:size(x)[1], 1:skipy:size(x)[2]], y[1:skipx:size(x)[1], 1:skipy:size(x)[2]], u[1:skipx:size(x)[1], 1:skipy:size(x)[2]], v[1:skipx:size(x)[1], 1:skipy:size(x)[2]], 
    umag[1:skipx:size(x)[1], 1:skipy:size(x)[2]], scale=20, cmap="jet")
#image = ax3.streamplot(x', y', u, v)
cb=fig.colorbar(image, ax=ax3)
ax3.set_title("u,v vector")
ax3.set_xlim(-3,9)
ax3.set_ylim(-6,6)


ax4 = fig.add_subplot(324)
heatmap3 = ax4.pcolormesh(x, y, pa, cmap="viridis", )
fig.colorbar(heatmap3, ax=ax4)
ax4.set_title("p")
#heatmap3.set_clim(vmin=0.35, vmax=0.5)
ax4.set_xlim(-3,9)
ax4.set_ylim(-6,6)

ax5 = fig.add_subplot(325)
heatmap4 = ax5.pcolormesh(x, y, omega, cmap="bwr", )
fig.colorbar(heatmap4, ax=ax5)
ax5.set_title("vorticity")
heatmap4.set_clim(vmin=-1, vmax=1)
ax5.set_xlim(-3,9)
ax5.set_ylim(-6,6)

plt.savefig(dir*"Fig/"*name*".pdf")
```