#!/usr/bin/env python # -*- coding: utf-8 -*- """ dclaw_debris_flow_demo_template.py D-Clawの土石流向けモデルを使うための、デモサイト掲載用テンプレートです。 - 1m DEM、発生源ポリゴン、河川・堤防・砂防施設データをD-Clawアプリの入力へ整理する想定です。 - D-Claw本体はFortran計算コードを含むため、このスクリプト単体では土石流計算を実行しません。 - 代わりに、D-Clawで必要になる地形・初期土砂・固相濃度・間隙圧の入力ファイル雛形を作ります。 - 実務版では rasterio / geopandas で実DEM・ベクトルデータを読み込み、この雛形の配列生成部分を置き換えます。 実行例: python scripts/dclaw_debris_flow_demo_template.py --outdir out_dclaw_debris_app """ from __future__ import annotations import argparse from pathlib import Path import numpy as np def create_synthetic_debris_case(nx: int, ny: int, dx: float) -> dict[str, np.ndarray]: """Webデモ用の渓流谷・発生源・砂防施設を作ります。""" x = np.arange(nx) * dx y = np.arange(ny) * dx xx, yy = np.meshgrid(x, y) centerline = y.mean() + 18.0 * np.sin(xx / 85.0) longitudinal = 180.0 - 0.105 * xx cross_valley = 0.0035 * (yy - centerline) ** 2 bank_roughness = 0.8 * np.sin(xx / 30.0) * np.sin(yy / 26.0) topo = longitudinal + cross_valley + bank_roughness # 河道を少し掘り下げる channel = np.exp(-((yy - centerline) / 6.5) ** 2) topo -= 1.2 * channel # 砂防堰堤風の横断構造物 dam_x = x[int(nx * 0.58)] dam = np.exp(-((xx - dam_x) / 4.0) ** 2) * np.exp(-((yy - centerline) / 48.0) ** 2) topo += 3.5 * dam # 発生源: 上流左岸寄りの楕円状土石流塊 sx = x[int(nx * 0.16)] sy = y[int(ny * 0.55)] source_shape = ((xx - sx) / 32.0) ** 2 + ((yy - sy) / 18.0) ** 2 h = 4.0 * np.exp(-source_shape) h[h < 0.04] = 0.0 # D-Clawの代表的な追加変数として、固相濃度と間隙圧を初期化する。 solid_fraction = np.where(h > 0.0, 0.60, 0.0) pore_pressure = np.where(h > 0.0, 0.8 * h, 0.0) return { "topo": topo, "initial_depth": h, "solid_fraction": solid_fraction, "pore_pressure": pore_pressure, "channel_mask": channel, "dam_mask": dam, } def write_esri_ascii(path: Path, arr: np.ndarray, dx: float, nodata: float = -9999.0) -> None: """ESRI ASCII grid形式で書き出します。""" ny, nx = arr.shape with path.open("w", encoding="utf-8") as f: f.write(f"ncols {nx}\n") f.write(f"nrows {ny}\n") f.write("xllcorner 0.0\n") f.write("yllcorner 0.0\n") f.write(f"cellsize {dx}\n") f.write(f"NODATA_value {nodata}\n") for row in arr[::-1, :]: f.write(" ".join(f"{v:.6f}" for v in row) + "\n") def write_setrun_template(path: Path, nx: int, ny: int, dx: float, duration: float) -> None: """D-Claw土石流アプリ用setrun.py雛形を書き出します。""" xupper = (nx - 1) * dx yupper = (ny - 1) * dx lines = [ "# -*- coding: utf-8 -*-", '"""D-Claw土石流アプリ用 setrun.py 雛形。"""', "", "", 'def setrun(claw_pkg="dclaw"):', " # 実際にはD-Claw公式ドキュメントのsetrun.pyテンプレートをベースにしてください。", f" xlower, xupper = 0.0, {xupper:.3f}", f" ylower, yupper = 0.0, {yupper:.3f}", f" mx, my = {nx}, {ny}", f" tfinal = {duration:.3f}", " nout = 30", "", ' topo_file = "topo_debris_valley.asc"', ' qinit_depth_file = "qinit_debris_depth.asc"', ' qinit_solid_fraction_file = "qinit_solid_fraction.asc"', ' qinit_pore_pressure_file = "qinit_pore_pressure.asc"', "", " # Example material parameters", " rho_f = 1000.0", " rho_s = 2700.0", " m0 = 0.60", " friction_angle = 35.0", " permeability = 1.0e-10", " compressibility = 1.0e-8", "", " # TODO: 公式D-Claw setrun.pyテンプレートに上記値を反映してください。", " return None", "", ] path.write_text("\n".join(lines), encoding="utf-8") def write_setplot_template(path: Path) -> None: """D-Claw結果可視化setplot.py雛形を書き出します。""" path.write_text( "# -*- coding: utf-8 -*-\n" '"""D-Claw結果可視化 setplot.py 雛形。"""\n\n' "def setplot(plotdata):\n" " # 実際のD-Claw/VisClaw環境で、深さ、速度、固相濃度、間隙圧などを描画します。\n" " plotdata.clearfigures()\n" " return plotdata\n", encoding="utf-8", ) def write_readme(path: Path) -> None: path.write_text( "# D-Claw debris-flow demo application template\n\n" "このフォルダは、D-Clawで土石流解析を行うためのデモ用アプリ雛形です。\n\n" "## 生成ファイル\n\n" "- topo_debris_valley.asc: 河道・堰堤を含む人工谷地形\n" "- qinit_debris_depth.asc: 初期土石流流動深\n" "- qinit_solid_fraction.asc: 固相体積濃度の初期分布\n" "- qinit_pore_pressure.asc: 間隙圧相当量の初期分布\n" "- setrun.py: D-Claw設定ファイルの雛形\n" "- setplot.py: 可視化設定の雛形\n" "- Makefile: 実行手順メモ\n", encoding="utf-8", ) def write_makefile(path: Path) -> None: path.write_text( "# D-Claw debris-flow application Makefile template\n\n" ".PHONY: help\n" "help:\n" "\t@echo \"D-Claw debris-flow demo template\"\n" "\t@echo \"Adjust setrun.py using the official D-Claw template, then run:\"\n" "\t@echo \" make .data\"\n" "\t@echo \" make .output\"\n" "\t@echo \" make .plots\"\n", encoding="utf-8", ) def main() -> None: parser = argparse.ArgumentParser(description="Create a D-Claw debris-flow demo app template") parser.add_argument("--outdir", default="out_dclaw_debris_app") parser.add_argument("--nx", type=int, default=220) parser.add_argument("--ny", type=int, default=120) parser.add_argument("--dx", type=float, default=1.0) parser.add_argument("--duration", type=float, default=120.0) args = parser.parse_args() outdir = Path(args.outdir) outdir.mkdir(parents=True, exist_ok=True) case = create_synthetic_debris_case(args.nx, args.ny, args.dx) write_esri_ascii(outdir / "topo_debris_valley.asc", case["topo"], args.dx) write_esri_ascii(outdir / "qinit_debris_depth.asc", case["initial_depth"], args.dx) write_esri_ascii(outdir / "qinit_solid_fraction.asc", case["solid_fraction"], args.dx) write_esri_ascii(outdir / "qinit_pore_pressure.asc", case["pore_pressure"], args.dx) write_setrun_template(outdir / "setrun.py", args.nx, args.ny, args.dx, args.duration) write_setplot_template(outdir / "setplot.py") write_makefile(outdir / "Makefile") write_readme(outdir / "README.md") volume = float(np.sum(case["initial_depth"]) * args.dx * args.dx) print(f"D-Claw土石流アプリ雛形を作成しました: {outdir}") print(f"初期土石流体積の概算: {volume:,.1f} m3") print("D-Claw本体の環境で setrun.py を公式テンプレートに合わせて調整してください。") if __name__ == "__main__": main()