#!/usr/bin/env python # -*- coding: utf-8 -*- """ anuga_debris_flow_demo.py ANUGAを使った「土石流風」平面2次元デモプログラムです。 重要: ANUGAは本来、2次元浅水流方程式を解く洪水・津波・氾濫解析用の ライブラリです。このサンプルでは、土石流を厳密な2相流として扱う のではなく、以下のように「高粗度・高濃度の混合流体」とみなして 流下・氾濫範囲のデモを行います。 - 地形: 渓流谷 + 河道 + 堤防/盛土 + 簡易砂防堰堤 - 流体: 発生源に初期流動深を与えるダムブレイク型 - 抵抗: Manning粗度を大きめに設定して土石流らしい減衰を表現 - 出力: ANUGA標準のSWWファイル + 最大流動深PNG/NPZ 実データへ拡張する場合: 1. GeoTIFF DEMを読み込んで elevation に補間する 2. 河川ライン・堤防ライン・砂防施設を地形補正として反映する 3. 発生源ポリゴンから初期体積またはハイドログラフを与える 4. 実績範囲と比較しながら friction, 初期流動深, 供給時間を調整する インストール例: python -m venv .venv source .venv/bin/activate # Windowsは .venv\\Scripts\\activate pip install anuga numpy matplotlib 実行例: python scripts/anuga_debris_flow_demo.py --outdir out_anuga_debris --duration 80 """ from __future__ import annotations import argparse from pathlib import Path import numpy as np try: import anuga except ImportError as exc: raise SystemExit( "ANUGAがインストールされていません。\n" "例: pip install anuga numpy matplotlib\n" "ANUGAはC/Cython拡張を含むため、環境によってはビルドツールが必要です。" ) from exc def make_topography(length: float, width: float): """ANUGAの set_quantity('elevation', func) に渡す人工地形関数を作る。""" center_y = width / 2.0 def topography(x: np.ndarray, y: np.ndarray) -> np.ndarray: # 上流 x=0 が高く、下流 x=length が低い縦断勾配 longitudinal = 80.0 - 0.11 * x # 谷地形: 中央河道から離れるほど高くする valley = 0.0045 * (y - center_y) ** 2 # 河道掘り込み: 中央部をガウス型に低くする channel = -2.0 * np.exp(-((y - center_y) / 10.0) ** 2) # 下流部の左右に堤防・道路盛土を置く levee_left = 2.5 * np.exp(-((y - (center_y - 32.0)) / 3.2) ** 2) levee_right = 2.5 * np.exp(-((y - (center_y + 32.0)) / 3.2) ** 2) levee_zone = ((x > 280.0) & (x < 560.0)).astype(float) levee = (levee_left + levee_right) * levee_zone # 簡易砂防堰堤: 谷を横断する高まり。ただし中央に越流部を残す。 dam_x = np.exp(-((x - 210.0) / 4.0) ** 2) dam_y = 1.0 - 0.75 * np.exp(-((y - center_y) / 11.0) ** 2) check_dam = 4.0 * dam_x * dam_y return longitudinal + valley + channel + levee + check_dam return topography def make_initial_stage(topography, length: float, width: float, source_depth: float): """発生源に初期流動深を与え、その他はドライベッドとする。""" center_y = width / 2.0 def stage(x: np.ndarray, y: np.ndarray) -> np.ndarray: z = topography(x, y) # 上流の楕円形発生源。 # 実務では崩壊源ポリゴンを読み込み、この条件を置き換える。 sx = (x - 58.0) / 46.0 sy = (y - center_y) / 24.0 source = sx * sx + sy * sy <= 1.0 # 上流側ほど厚く、端部でなめらかに薄くする。 depth = np.zeros_like(x, dtype=float) depth[source] = source_depth * (1.0 - 0.35 * np.clip(sx[source], -1.0, 1.0)) # ANUGAの stage は「水面標高」なので、stage = elevation + depth。 return z + depth return stage def run_demo(args: argparse.Namespace) -> None: outdir = Path(args.outdir) outdir.mkdir(parents=True, exist_ok=True) # ANUGAの矩形三角形メッシュ。 # 実務では create_domain_from_regions() で解析範囲ポリゴン・細分化領域を与える。 domain = anuga.rectangular_cross_domain( args.nx, args.ny, len1=args.length, len2=args.width, ) domain.set_name("anuga_debris_flow_demo") domain.set_datadir(str(outdir)) domain.set_store_vertices_smoothly(False) topography = make_topography(args.length, args.width) initial_stage = make_initial_stage(topography, args.length, args.width, args.source_depth) # 地形、粗度、初期水位を設定。 domain.set_quantity("elevation", topography, location="centroids") domain.set_quantity("friction", args.friction, location="centroids") domain.set_quantity("stage", initial_stage, location="centroids") # 境界条件。左・上下は反射境界、右は可能なら透過境界とする。 Br = anuga.Reflective_boundary(domain) try: Bout = anuga.Transmissive_boundary(domain) except AttributeError: # 古い環境・API差異への保険。下流端を十分低いドライ境界にする。 Bout = anuga.Dirichlet_boundary([0.0, 0.0, 0.0]) domain.set_boundary({"left": Br, "right": Bout, "top": Br, "bottom": Br}) # 計算中に最大流動深・最大流速を記録する。 max_depth = None max_speed = None for t in domain.evolve(yieldstep=args.yieldstep, finaltime=args.duration): stage = domain.quantities["stage"].centroid_values elev = domain.quantities["elevation"].centroid_values xmom = domain.quantities["xmomentum"].centroid_values ymom = domain.quantities["ymomentum"].centroid_values depth = np.maximum(stage - elev, 0.0) speed = np.zeros_like(depth) wet = depth > 1.0e-6 speed[wet] = np.sqrt(xmom[wet] ** 2 + ymom[wet] ** 2) / depth[wet] max_depth = depth.copy() if max_depth is None else np.maximum(max_depth, depth) max_speed = speed.copy() if max_speed is None else np.maximum(max_speed, speed) print(domain.timestepping_statistics()) # 後処理用に重心座標と最大値を保存。 xy = domain.get_centroid_coordinates() np.savez( outdir / "anuga_debris_result_centroids.npz", x=xy[:, 0], y=xy[:, 1], max_depth=max_depth, max_speed=max_speed, ) # 簡易PNGを作成。SWWはQGIS/CrayfishやANUGA Viewerで確認可能。 try: import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(9, 3.8), constrained_layout=True) tri = ax.tricontourf(xy[:, 0], xy[:, 1], max_depth, levels=20) ax.set_aspect("equal") ax.set_title("ANUGA debris-flow style demo: maximum flow depth") ax.set_xlabel("x [m]") ax.set_ylabel("y [m]") fig.colorbar(tri, ax=ax, label="maximum depth [m]") fig.savefig(outdir / "anuga_debris_max_depth.png", dpi=180) plt.close(fig) except Exception as exc: # matplotlibなしでも計算結果は保存する print(f"PNG作成をスキップしました: {exc}") print("\n完了しました。") print(f" SWW: {outdir / 'anuga_debris_flow_demo.sww'}") print(f" NPZ: {outdir / 'anuga_debris_result_centroids.npz'}") print(f" PNG: {outdir / 'anuga_debris_max_depth.png'}") def build_parser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser(description="ANUGAを使った土石流風2Dデモ") parser.add_argument("--outdir", default="out_anuga_debris", help="出力フォルダ") parser.add_argument("--length", type=float, default=650.0, help="解析領域の長さ m") parser.add_argument("--width", type=float, default=220.0, help="解析領域の幅 m") parser.add_argument("--nx", type=int, default=130, help="x方向分割数") parser.add_argument("--ny", type=int, default=44, help="y方向分割数") parser.add_argument("--duration", type=float, default=80.0, help="計算時間 秒") parser.add_argument("--yieldstep", type=float, default=2.0, help="出力・表示間隔 秒") parser.add_argument("--source-depth", type=float, default=5.0, help="発生源の代表初期流動深 m") parser.add_argument("--friction", type=float, default=0.08, help="Manning粗度係数。土石流風に大きめを設定") return parser if __name__ == "__main__": run_demo(build_parser().parse_args())