#!/usr/bin/env python # -*- coding: utf-8 -*- """ anuga_basic_shallow_water_sample.py ANUGA 入門サンプル: 三角形メッシュ上の2次元浅水流計算 目的: - ANUGAの基本的な使い方を確認する - 計算領域を三角形メッシュで作成する - 人工地形、初期水位、粗度、境界条件を設定する - evolve() で時間発展させる - 最大水深と最大流速をPNGとして保存する 実行例: pip install anuga numpy matplotlib python scripts/anuga_basic_shallow_water_sample.py --outdir out_anuga_basic 注意: - このサンプルは、ANUGAの基本構成を説明するための入門コードです。 - 土石流を厳密な2相流として扱うものではありません。 - 土石流風の概略デモは scripts/anuga_debris_flow_demo.py を参照してください。 """ 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', function) に渡す人工地形関数。""" center_y = width / 2.0 def topography(x: np.ndarray, y: np.ndarray) -> np.ndarray: # 全体として上流 x=0 から下流 x=length へ低くなる地形 longitudinal = 18.0 - 0.018 * x # 中央に浅い谷を作る valley = 0.0015 * (y - center_y) ** 2 channel = -0.8 * np.exp(-((y - center_y) / 18.0) ** 2) # 下流側に小さな盛土を作り、流れの回り込みを確認しやすくする mound = 1.2 * np.exp(-((x - 430.0) / 35.0) ** 2 - ((y - center_y) / 28.0) ** 2) return longitudinal + valley + channel + mound return topography def make_initial_stage(length: float, width: float, topography_func, initial_depth: float): """上流側の一部だけに初期水深を与える。""" center_y = width / 2.0 def stage(x: np.ndarray, y: np.ndarray) -> np.ndarray: z = topography_func(x, y) source = (x < 90.0) & (np.abs(y - center_y) < 24.0) h0 = np.where(source, initial_depth, 0.0) return z + h0 return stage def run(args: argparse.Namespace) -> None: outdir = Path(args.outdir) outdir.mkdir(parents=True, exist_ok=True) length = args.length width = args.width # ------------------------------------------------------------ # 1. 三角形メッシュ領域の作成 # ------------------------------------------------------------ # rectangular_cross_domain は、矩形領域を三角形メッシュに分割します。 # nx, ny を大きくすると解像度は上がりますが、計算時間も増えます。 domain = anuga.rectangular_cross_domain( args.nx, args.ny, len1=length, len2=width, ) domain.set_name("anuga_basic_shallow_water_sample") domain.set_datadir(str(outdir)) # ------------------------------------------------------------ # 2. 地形、粗度、初期水位の設定 # ------------------------------------------------------------ topography = make_topography(length, width) initial_stage = make_initial_stage(length, width, topography, args.initial_depth) domain.set_quantity("elevation", topography) domain.set_quantity("friction", args.friction) domain.set_quantity("stage", initial_stage) # ------------------------------------------------------------ # 3. 境界条件の設定 # ------------------------------------------------------------ # 上流・側方は反射境界、下流は水位ゼロ相当の流出境界にしています。 # 実務モデルでは、河川流量ハイドログラフや潮位・水位時系列を境界条件にします。 reflective = anuga.Reflective_boundary(domain) downstream = anuga.Dirichlet_boundary([0.0, 0.0, 0.0]) domain.set_boundary({ "left": reflective, "right": downstream, "top": reflective, "bottom": reflective, }) # ------------------------------------------------------------ # 4. 時間発展計算 # ------------------------------------------------------------ max_depth = np.zeros(domain.number_of_elements, dtype=float) max_speed = np.zeros(domain.number_of_elements, dtype=float) 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 = np.maximum(max_depth, depth) max_speed = np.maximum(max_speed, speed) if args.verbose: print(f"t={t:8.2f} sec, max_depth={max_depth.max():.3f} m") # ------------------------------------------------------------ # 5. 結果保存 # ------------------------------------------------------------ centroid_x = domain.get_centroid_coordinates()[:, 0] centroid_y = domain.get_centroid_coordinates()[:, 1] np.savez( outdir / "anuga_basic_result_centroids.npz", x=centroid_x, y=centroid_y, max_depth=max_depth, max_speed=max_speed, ) # matplotlibは、ANUGA本体の計算が終わった後だけ使います。 import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(9, 4.8)) sc = ax.tricontourf(centroid_x, centroid_y, max_depth, levels=24) ax.set_aspect("equal") ax.set_title("ANUGA sample: maximum water depth") ax.set_xlabel("x [m]") ax.set_ylabel("y [m]") fig.colorbar(sc, ax=ax, label="max depth [m]") fig.tight_layout() fig.savefig(outdir / "anuga_basic_max_depth.png", dpi=180) plt.close(fig) fig, ax = plt.subplots(figsize=(9, 4.8)) sc = ax.tricontourf(centroid_x, centroid_y, max_speed, levels=24) ax.set_aspect("equal") ax.set_title("ANUGA sample: maximum velocity") ax.set_xlabel("x [m]") ax.set_ylabel("y [m]") fig.colorbar(sc, ax=ax, label="max speed [m/s]") fig.tight_layout() fig.savefig(outdir / "anuga_basic_max_speed.png", dpi=180) plt.close(fig) print("完了しました。") print(f" {outdir / 'anuga_basic_shallow_water_sample.sww'}") print(f" {outdir / 'anuga_basic_result_centroids.npz'}") print(f" {outdir / 'anuga_basic_max_depth.png'}") print(f" {outdir / 'anuga_basic_max_speed.png'}") def build_parser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser(description="ANUGA入門 2次元浅水流サンプル") parser.add_argument("--outdir", default="out_anuga_basic", help="出力フォルダ") parser.add_argument("--length", type=float, default=600.0, help="領域長さ m") parser.add_argument("--width", type=float, default=180.0, help="領域幅 m") parser.add_argument("--nx", type=int, default=80, help="x方向分割数") parser.add_argument("--ny", type=int, default=28, help="y方向分割数") parser.add_argument("--duration", type=float, default=80.0, help="計算時間 秒") parser.add_argument("--yieldstep", type=float, default=1.0, help="出力・確認間隔 秒") parser.add_argument("--initial-depth", type=float, default=1.4, help="上流初期水深 m") parser.add_argument("--friction", type=float, default=0.045, help="Manning粗度係数") parser.add_argument("--verbose", action="store_true", help="各ステップの情報を表示") return parser if __name__ == "__main__": run(build_parser().parse_args())