Last active
June 19, 2024 14:35
-
-
Save terasakisatoshi/3e79334bae96ddf7901446223d5b4ea2 to your computer and use it in GitHub Desktop.
2次元非線形長波方程式による津波伝播計算
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
output.mp4
もう少し長い版
output.mp4
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
元実装は Python 実装.
https://github.com/kaigan-osu/2D_Tsunami
Julia の方が明らかに高速
Usage