浅水方程式の先端の取り扱いに関する研究
https://www.jstage.jst.go.jp/article/jscejhe/72/4/72_I_325/_pdf
// シミュレーションパラメータ
float g = 9.81; // 重力加速度
float dt = 0.01; // 時間ステップ
float dx = 0.1; // x方向の空間ステップ
float dy = 0.1; // y方向の空間ステップ
float nu = 0.001; // 数値的粘性（拡散）係数
int rows = chi("rows"); // グリッドの行数
int cols = chi("cols"); // グリッドの列数
// 現在のポイントの水の高さと流速
float h = @h;
float u = @u;
float v = @v;
// 境界条件の確認
int i = @ptnum % cols;
int j = @ptnum / cols;
if (i == 0 || j == 0 || i == cols - 1 || j == rows - 1) {
// 境界での処理
h = 0.0;
u = 0.0;
v = 0.0;
}
else {
// 隣接ポイントの水の高さと流速
float h_right = point(0, "h", @ptnum + 1);
float u_right = point(0, "u", @ptnum + 1);
float h_left = point(0, "h", @ptnum - 1);
float u_left = point(0, "u", @ptnum - 1);
float h_up = point(0, "h", @ptnum + cols);
float v_up = point(0, "v", @ptnum + cols);
float h_down = point(0, "h", @ptnum - cols);
float v_down = point(0, "v", @ptnum - cols);
// 二次元浅水方程式の計算
float h_new = h - dt * ( (u_right * h_right - u_left * h_left) / (2 * dx) +
(v_up * h_up - v_down * h_down) / (2 * dy));
float u_new = u - dt * (u * (u_right - u_left) / (2 * dx) + g * (h_right - h_left) / (2 * dx));
float v_new = v - dt * (v * (v_up - v_down) / (2 * dy) + g * (h_up - h_down) / (2 * dy));
// 数値的拡散の追加
h_new += nu * ( (h_right - 2 * h + h_left) / (dx * dx) + (h_up - 2 * h + h_down) / (dy * dy));
u_new += nu * ( (u_right - 2 * u + u_left) / (dx * dx) + (u - 2 * u + u) / (dy * dy));
v_new += nu * ( (v - 2 * v + v) / (dx * dx) + (v_up - 2 * v + v_down) / (dy * dy));
// 負の高さを防ぐ
if(h_new < 0) h_new = 0;
// 更新
h = h_new;
u = u_new;
v = v_new;
}
// ポイント属性の更新
setpointattrib(0, "h", @ptnum, h, "set");
setpointattrib(0, "u", @ptnum, u, "set");
setpointattrib(0, "v", @ptnum, v, "set");
// 水の高さでY座標を更新
@P.y = @h;
