2次元非線形長波方程式による津波伝播計算

draw_surface_all_images.py

# https://github.com/kaigan-osu/2D_Tsunami/blob/main/draw_surface_all_images.py
# を改造したもの
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm
from natsort import natsorted
import glob

mx = 750  # 東西方向の格子数
my = 495  # 南北方向の格子数
dx = 1620.0  # 格子間隔 [m]
dy = 1620.0  # 格子間隔 [m]

# 読込データのファイルリストを読み込む
# case = np.loadtxt('flist.csv',delimiter=',', dtype='str')
case = natsorted(glob.glob("./data/*.csv"))

for kt, _ in enumerate(case):
    print(case[kt], " in process")
    eta = np.loadtxt(case[kt], delimiter=",")
    etanew = np.flipud(eta[1 : my + 1, 1 : mx + 1])

    # 図形描画の準備
    fig = plt.figure(figsize=(11.2, 6))  # 図の枠の準備
    ax = fig.add_subplot(111)

    x = [dx * i / 1000.0 for i in range(mx)]  # x座標の値
    y = [dy * j / 1000.0 for j in range(my)]  # y座標の値
    X, Y = np.meshgrid(x, y)

    # 結果の出力
    ax.set_xlabel("X-axis [km]")
    ax.set_ylabel("Y-axis [km]")

    # カラーバーの範囲を設定
    vmin, vmax = -3.0, 3.0

    # 等高線の間隔を細かく設定
    levels = np.linspace(vmin, vmax, 100)

    # BoundaryNormを使用して、範囲外の色を指定
    cmap = plt.get_cmap("bwr")
    norm = BoundaryNorm(np.linspace(vmin, vmax, 101), cmap.N, clip=True)

    map = ax.contourf(
        X, Y, etanew, cmap=cmap, levels=levels, norm=norm, extend="both"
    )  # 水位分布をプロット
    fig.colorbar(map, ax=ax, label="elevation [m]", ticks=np.arange(-3, 3.1, 0.5))

    # fig内でのaxes座標を取得，戻り値はBbox
    ax_pos = ax.get_position()
    # fig内座標でテキストを表示 Bboxは Bbox.x0, Bbox.x1, Bbox.y0, Bbox.y1で座標を取得できる
    fig.text(ax_pos.x1 - 0.11, ax_pos.y1 + 0.01, "file= " + case[kt])

    h = np.loadtxt("depth_1620-01.csv", delimiter=",")  # 水深データの読み込み
    hnew = np.flipud(h)
    ax.contour(X, Y, hnew, levels=[0.01], colors="black", linewidths=0.5)  # 地形の形状をプロット
    print("images/" + str(kt).zfill(4) + ".jpg")
    plt.savefig("images/" + str(kt).zfill(4) + ".jpg")
    plt.cla()
    plt.close()

    # plt.show()

main.jl

using OffsetArrays: Origin
using DelimitedFiles

function setupdepth()
    # 水深データ
    depth = readdlm(joinpath("./depth_1620-01.csv"), ',')::Matrix{Float64}
    depth = vcat(depth[begin, :]', depth, depth[end, :]')
    depth = hcat(depth[:, begin], depth, depth[:, end])
    # 陸地の水深はゼロ
    depth .= max.(depth, 0.0)
    # 水深10mよりも浅い海はすべて水深10mにセット
    depth .= min.(depth, 10.0)
    return Origin(0)(Matrix{Float64}(depth))
end


function setupη()
    # 初期水位データ
    η = readdlm(joinpath("./surface_1620.csv"), ',')::Matrix{Float64}
    η = vcat(η[begin, :]', η, η[end, :]')
    η = hcat(η[:, begin], η, η[:, end])
    return Origin(0)(Matrix{Float64}(η))
end

function main()
    # 2次元非線形長波方程式による津波伝播計算（2024.6.17）
    h = setupdepth()
    eta = setupη()

    my = size(h, 1) - 2
    mx = size(h, 2) - 2
    # 初期条件（計算領域の外側にゼロの境界をすべての方向に設定する）
    etanew = Origin(0)(zeros(my + 2, mx + 2))
    Mnew = Origin(0)(zeros(my + 2, mx + 2))
    M = Origin(0)(zeros(my + 2, mx + 2))
    Nnew = Origin(0)(zeros(my + 2, mx + 2))
    N = Origin(0)(zeros(my + 2, mx + 2))

    Mp1 = Origin(0)(zeros(my + 2, mx + 2))
    Mp2 = Origin(0)(zeros(my + 2, mx + 2))
    Mp3 = Origin(0)(zeros(my + 2, mx + 2))
    Mp4 = Origin(0)(zeros(my + 2, mx + 2))
    dm = Origin(0)(zeros(my + 2, mx + 2))
    N4 = Origin(0)(zeros(my + 2, mx + 2))

    Np1 = Origin(0)(zeros(my + 2, mx + 2))
    Np2 = Origin(0)(zeros(my + 2, mx + 2))
    Np3 = Origin(0)(zeros(my + 2, mx + 2))
    Np4 = Origin(0)(zeros(my + 2, mx + 2))
    dn = Origin(0)(zeros(my + 2, mx + 2))
    M4 = Origin(0)(zeros(my + 2, mx + 2))
    C = Origin(0)(zeros(my + 2, mx + 2))

    main(
        h,
        eta,
        etanew,
        Mnew,
        M,
        Nnew,
        N,
        Mp1,
        Mp2,
        Mp3,
        Mp4,
        dm,
        N4,
        Np1,
        Np2,
        Np3,
        Np4,
        dn,
        M4,
        C,
    )
end

function main(
    h,
    eta,
    etanew,
    Mnew,
    M,
    Nnew,
    N,
    Mp1,
    Mp2,
    Mp3,
    Mp4,
    dm,
    N4,
    Np1,
    Np2,
    Np3,
    Np4,
    dn,
    M4,
    C,
)
    # 計算条件
    ktend = 7200  # 総計算回数
    output_step = 30  # 計算結果の出力間隔
    mx = 750  # 東西方向の計算格子数
    my = 495  # 南北方向の計算格子数
    g = 9.80665  # 重力加速度 [m/s2]
    dt = 2.0  # 1ステップの時間 [sec]
    dx = 1620.0  # 東西方向の計算格子サイズ [m]
    dy = 1620.0  # 南北方向の計算格子サイズ [m]
    manning = 0.025  # マニングの粗度係数（摩擦の計算）
    omega = 7.29 * 10^(-5)  # 地球の角速度 [rad/s]
    f = 2.0 * omega * sin(pi / 6.0)  # コリオリ係数（北緯30°で代表）

    for j in range(0, my + 1)
        for i in range(0, mx + 1)
            if h[j, i] < 0.0  # 陸地の水深はゼロ
                h[j, i] = 0.0
            elseif h[j, i] < 10.0 && h[j, i] > 0.0  # 水深10mよりも浅い海はすべて水深10mにセット
                h[j, i] = 10.0
            end

            C[j, i] = √(g * h[j, i])  # 波速の計算
        end
    end

    # 最初の境界条件
    M[:, 1] .= 0.0
    M[:, mx+1] .= 0.0
    M[0, :] .= M[1, :]
    M[my+1, :] .= M[my, :]

    N[1, :] .= 0.0
    N[my+1, :] .= 0.0
    N[:, 0] .= N[:, 1]
    N[:, mx+1] .= N[:, mx]



    # 時間の計算--------------------------------------------------------------------------
    for kt in range(0, ktend - 1)
        # 水位の計算
        emax = -100.0
        emin = +100.0
        jemin = 1
        jemax = 1
        iemin = 1
        iemax = 1
        for i in range(1, mx)
            for j in range(1, my)
                if h[j, i] > 0.0
                    etanew[j, i] =
                        eta[j, i] -
                        dt * ((M[j, i+1] - M[j, i]) / dx + (N[j+1, i] - N[j, i]) / dy)
                    if etanew[j, i] < -h[j, i]  # 水位が水深よりも下回らないようにする
                        etanew[j, i] = -h[j, i]
                    end
                else
                    etanew[j, i] = 0.0
                end

                if etanew[j, i] > emax  # 水位の最大値を探す
                    emax = etanew[j, i]
                    jemax, iemax = j, i
                end
                if etanew[j, i] < emin  # 水位の最小値を探す
                    emin = etanew[j, i]
                    jemin, iemin = j, i
                end
            end
        end

        fn = lpad(kt, 4, "0") * ".csv"
        @info "info" kt fn emax jemax iemax emin jemin iemin

        # 水位の境界条件の更新（ゾンマーフェルトの境界条件）
        @views begin
            etanew[:, 0] .= eta[:, 0] .- dt .* C[:, 1] .* (eta[:, 0] .- eta[:, 1]) ./ dx
            etanew[:, mx+1] .=
                (eta[:, mx+1] .- dt .* C[:, mx] .* (eta[:, mx+1] .- eta[:, mx]) ./ dx)
            etanew[0, :] .= eta[0, :] .- dt .* C[1, :] .* (eta[0, :] .- eta[1, :]) ./ dy
            etanew[my+1, :] .=
                (eta[my+1, :] .- dt .* C[my, :] .* (eta[my+1, :] .- eta[my, :]) ./ dy)
        end
        # ========================================================================================東西の流れの計算
        for i in range(2, mx)
            for j in range(1, my)
                dm[j, i] = max(etanew[j, i], etanew[j, i-1]) + max(h[j, i], h[j, i-1])  # 計算格子間の水位
                N4[j, i] = (N[j, i] + N[j+1, i] + N[j, i-1] + N[j+1, i-1]) / 4.0  # Mの計算点におけるNの値
            end
        end

        for i in range(2, mx)
            for j in range(1, my)
                if h[j, i] > 0.0 && h[j, i-1] > 0.0  # 計算点の両側が海であること
                    if (M[j, i] > 0.0 && dm[j, i] > 0.0 && dm[j, i-1] > 0.0)  # 流れの方向の確認と計算格子間に水があること
                        Mp1[j, i] = (M[j, i]^2 / dm[j, i] - M[j, i-1]^2 / dm[j, i-1]) / dx

                    elseif M[j, i] < 0.0 && dm[j, i+1] > 0.0 && dm[j, i] > 0.0
                        Mp1[j, i] = (M[j, i+1]^2 / dm[j, i+1] - M[j, i]^2 / dm[j, i]) / dx

                    else
                        Mp1[j, i] = 0.0  # 流れの方向がない場合，計算格子間に水がない場合
                    end

                    if (N4[j, i] > 0.0 && dm[j, i] > 0.0 && dm[j-1, i] > 0.0)  # 流れの方向の確認と計算格子間に水があること
                        Mp2[j, i] =
                            (
                                M[j, i] * N4[j, i] / dm[j, i] -
                                M[j-1, i] * N4[j-1, i] / dm[j-1, i]
                            ) / dy

                    elseif N4[j, i] < 0.0 && dm[j+1, i] > 0.0 && dm[j, i] > 0.0
                        Mp2[j, i] =
                            (
                                M[j+1, i] * N4[j+1, i] / dm[j+1, i] -
                                M[j, i] * N4[j, i] / dm[j, i]
                            ) / dy

                    else
                        Mp2[j, i] = 0.0  # 流れの方向がない場合，計算格子間に水がない場合
                    end
                    # ------------------------------------------- 摩擦項の計算(浅すぎると摩擦項が発散する)
                    if dm[j, i] >= 1.0
                        Mp3[j, i] = (
                            g * manning^2 / dm[j, i]^(7.0 / 3.0) *
                            M[j, i] *
                            √(M[j, i]^2 + N4[j, i]^2)
                        )
                    elseif dm[j, i] < 1.0
                        Mp3_tmp = (
                            g * manning^2 / 1.0^(7.0 / 3.0) *
                            M[j, i] *
                            √(M[j, i]^2 + N4[j, i]^2)
                        )
                        Mp3[j, i] = Mp3_tmp * dm[j, i]^2
                    end
                    # --------------------------------------------------------------------------ここまで

                    Mp4[j, i] = -f * N4[j, i]  # コリオリ力の項

                else  # 計算点の少なくともどちらかが陸の場合は岸を抜ける流れはない
                    Mp1[j, i] = 0.0
                    Mp2[j, i] = 0.0
                    Mp3[j, i] = 0.0
                    Mp4[j, i] = 0.0
                end
            end
        end

        # ----------------------------------------------------------------------- 東西方向の流れを更新
        for i in range(2, mx)
            for j in range(1, my)
                if h[j, i] > 0.0 && h[j, i-1] > 0.0
                    if dm[j, i] > 0.01
                        Mnew[j, i] =
                            M[j, i] -
                            dt * (
                                Mp1[j, i] +
                                Mp2[j, i] +
                                Mp3[j, i] +
                                Mp4[j, i] +
                                g * h[j, i] * (etanew[j, i] - etanew[j, i-1]) / dx
                            )
                    else
                        Mnew[j, i] = 0.0
                    end
                else
                    Mnew[j, i] = 0.0
                end
            end
        end
        # -----------------------------------------------------------------------------------ここまで
        @views begin
            Mnew[0, :] .= Mnew[1, :]  # 境界に沿った方向の流れは境界の摩擦はない境界
            Mnew[my+1, :] .= Mnew[my, :]  # 境界に沿った方向の流れは境界の摩擦はない境界
            Mnew[:, 1] .= M[:, 1] .- dt .* C[:, 1] .* (M[:, 1] .- M[:, 2]) ./ dx  # 境界に垂直な流れ方向はゾンマーフェルト境界
            Mnew[:, mx+1] .=
                (M[:, mx+1] .- dt .* C[:, mx] .* (M[:, mx+1] .- M[:, mx]) ./ dx)  # 境界に垂直な流れ方向はゾンマーフェルト境界
        end
        # ========================================================================================南北の流れの計算
        for i in range(1, mx)
            for j in range(2, my)
                dn[j, i] = max(etanew[j, i], etanew[j-1, i]) + max(h[j, i], h[j-1, i])
                M4[j, i] = (M[j, i] + M[j-1, i] + M[j, i+1] + M[j-1, i+1]) / 4.0
            end
        end
        for i in range(1, mx)
            for j in range(2, my)
                if h[j, i] > 0.0 && h[j-1, i] > 0.0
                    if M4[j, i] > 0.0 && dn[j, i] > 0.0 && dn[j, i-1] > 0.0
                        Np1[j, i] =
                            (
                                N[j, i] * M4[j, i] / dn[j, i] -
                                N[j, i-1] * M4[j, i-1] / dn[j, i-1]
                            ) / dx

                    elseif M4[j, i] < 0.0 && dn[j, i+1] > 0.0 && dn[j, i] > 0.0
                        Np1[j, i] =
                            (
                                N[j, i+1] * M4[j, i+1] / dn[j, i+1] -
                                N[j, i] * M4[j, i] / dn[j, i]
                            ) / dx

                    else
                        Np1[j, i] = 0.0
                    end

                    if N[j, i] > 0.0 && dn[j, i] > 0.0 && dn[j-1, i] > 0.0
                        Np2[j, i] = (N[j, i]^2 / dn[j, i] - N[j-1, i]^2 / dn[j-1, i]) / dy

                    elseif N[j, i] < 0.0 && dn[j+1, i] > 0.0 && dn[j, i] > 0.0
                        Np2[j, i] = (N[j+1, i]^2 / dn[j+1, i] - N[j, i]^2 / dn[j, i]) / dy

                    else
                        Np2[j, i] = 0.0
                    end

                    # --------------------------------------------------------------------- 摩擦項の計算
                    if dn[j, i] >= 1.0
                        Np3[j, i] = (
                            g * manning^2 / dn[j, i]^(7.0 / 3.0) *
                            N[j, i] *
                            √(M4[j, i]^2 + N[j, i]^2)
                        )
                    elseif dn[j, i] < 1.0
                        Np3_tmp = (
                            g * manning^2 / 1.0^(7.0 / 3.0) *
                            N[j, i] *
                            √(M4[j, i]^2 + N[j, i]^2)
                        )
                        Np3[j, i] = Np3_tmp * dn[j, i]^2
                    end
                    # --------------------------------------------------------------------------ここまで

                    Np4[j, i] = f * M4[j, i]  # コリオリ力の項

                else
                    Np1[j, i] = 0.0
                    Np2[j, i] = 0.0
                    Np3[j, i] = 0.0
                    Np4[j, i] = 0.0
                end
            end
        end
        # --------------------------------------------------------------------------- 南北方向の流れを更新
        for i in range(1, mx)
            for j in range(2, my)
                if h[j, i] > 0.0 && h[j-1, i] > 0.0
                    if dn[j, i] > 0.01
                        Nnew[j, i] =
                            N[j, i] -
                            dt * (
                                Np1[j, i] +
                                Np2[j, i] +
                                Np3[j, i] +
                                Np4[j, i] +
                                g * dn[j, i] * (etanew[j, i] - etanew[j-1, i]) / dy
                            )
                    else
                        Nnew[j, i] = 0.0
                    end
                else
                    Nnew[j, i] = 0.0
                end
            end
        end
        # -----------------------------------------------------------------------------------ここまで
        @views begin
            Nnew[:, 0] .= Nnew[:, 1]
            Nnew[:, mx+1] .= Nnew[:, mx]
            Nnew[1, :] .= N[1, :] .- dt .* C[1, :] .* (N[1, :] .- N[2, :]) ./ dy
            Nnew[my+1, :] .= N[my+1, :] .- dt .* C[my, :] .* (N[my+1, :] .- N[my, :]) ./ dy
        end
        # ======================================================================================= 結果の出力
        if mod(kt, output_step) == 0
            writedlm("data/" * fn, etanew, ',')
        end
        # ====================================================================================結果の入れ替え
        @. eta = etanew
        @. M = Mnew
        @. N = Nnew
    end

end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

元実装は Python 実装．

https://github.com/kaigan-osu/2D_Tsunami

Julia の方が明らかに高速

Usage

$ mkdir data
$ julia main.jl
$ mkdir images
$ python draw_surface_all_images.py
$ cd images
$ ffmpeg -framerate 30 -i %04d.jpg -c:v libx264 -r 30 -pix_fmt yuv420p output.mp4

output.mp4

もう少し長い版

output.mp4