#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
dclaw_basic_granular_flow_template.py

D-Claw概要ページ用のPythonサンプルです。

D-Clawは、pip installだけで完結する純Pythonライブラリではなく、
Clawpack系のFortran計算コードとPython設定ファイルを組み合わせて使います。
このスクリプトは、D-Clawアプリケーションで必要になる
地形データ、初期土砂分布、setrun.py、Makefile、READMEを作るための
テンプレートです。

実際にD-Clawで計算するには、Linux/Unix上でD-Clawをインストールしたうえで、
生成されたアプリケーションフォルダ内で make .data / make .output などを実行します。

実行例:
    python scripts/dclaw_basic_granular_flow_template.py --outdir out_dclaw_basic_app
"""

from __future__ import annotations

import argparse
from pathlib import Path

import numpy as np


def make_valley_topography(nx: int, ny: int, dx: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """D-Claw入力に使う人工谷地形を作成します。"""
    x = np.arange(nx) * dx
    y = np.arange(ny) * dx
    xx, yy = np.meshgrid(x, y)

    # 下流方向へ低くなる地形 + 谷地形
    slope = -0.08 * xx
    valley = 0.0009 * (yy - y.mean()) ** 2
    roughness = 0.4 * np.sin(xx / 45.0) * np.cos(yy / 35.0)
    z = 120.0 + slope + valley + roughness
    return x, y, z


def make_initial_debris(nx: int, ny: int, dx: float) -> dict[str, np.ndarray]:
    """初期土石流塊の分布を作ります。D-Clawのqinit作成前処理を想定します。"""
    x = np.arange(nx) * dx
    y = np.arange(ny) * dx
    xx, yy = np.meshgrid(x, y)

    cx = x[int(nx * 0.18)]
    cy = y[int(ny * 0.52)]
    r2 = ((xx - cx) / 38.0) ** 2 + ((yy - cy) / 24.0) ** 2
    h = 3.0 * np.exp(-r2)
    h[h < 0.03] = 0.0

    # D-Clawで重要な考え方: 固相体積濃度と間隙圧を初期条件として用意する。
    solid_fraction = np.where(h > 0.0, 0.58, 0.0)
    pore_pressure = np.where(h > 0.0, 0.75 * h, 0.0)

    return {
        "h": h,
        "solid_fraction": solid_fraction,
        "pore_pressure": pore_pressure,
    }


def write_esri_ascii(path: Path, arr: np.ndarray, dx: float, nodata: float = -9999.0) -> None:
    """GeoClaw/D-Clawの前処理で使いやすいESRI ASCII grid形式で書き出します。"""
    ny, nx = arr.shape
    lines = [
        f"ncols {nx}",
        f"nrows {ny}",
        "xllcorner 0.0",
        "yllcorner 0.0",
        f"cellsize {dx}",
        f"NODATA_value {nodata}",
    ]
    for row in arr[::-1, :]:
        lines.append(" ".join(f"{v:.6f}" for v in row))
    path.write_text("\n".join(lines) + "\n", encoding="utf-8")


def write_setrun_template(path: Path, nx: int, ny: int, dx: float) -> None:
    """D-Claw用setrun.pyの雛形を書きます。値は実務データに合わせて調整します。"""
    xupper = (nx - 1) * dx
    yupper = (ny - 1) * dx

    # 実際のD-Claw環境では、公式setrun.pyテンプレートをベースにしてください。
    lines = [
        "# -*- coding: utf-8 -*-",
        '"""setrun.py - D-Claw application template."""',
        "",
        "",
        'def setrun(claw_pkg="dclaw"):',
        "    # 実際のD-Claw環境では、dclaw.data / clawpack.clawutil.data などを使って",
        "    # ClawRunDataを作成します。",
        "    lower = [0.0, 0.0]",
        f"    upper = [{xupper:.3f}, {yupper:.3f}]",
        f"    num_cells = [{nx}, {ny}]",
        "    tfinal = 120.0",
        "    num_output_times = 24",
        "",
        "    # D-Claw物性値の例。実務では再現計算で調整します。",
        "    rho_f = 1000.0      # fluid density kg/m3",
        "    rho_s = 2700.0      # solid density kg/m3",
        "    m_crit = 0.58       # representative solid volume fraction",
        "    delta = 35.0        # internal friction angle degree",
        "    kappa = 1.0e-10     # permeability",
        "",
        '    topography_file = "topo_valley.asc"',
        '    initial_depth_file = "qinit_depth.asc"',
        '    initial_solid_fraction_file = "qinit_solid_fraction.asc"',
        '    initial_pore_pressure_file = "qinit_pore_pressure.asc"',
        "",
        "    # TODO: 公式D-Claw setrun.pyテンプレートへ上記値を反映してください。",
        "    return None",
        "",
    ]
    path.write_text("\n".join(lines), encoding="utf-8")


def write_makefile(path: Path) -> None:
    path.write_text(
        "# D-Claw application Makefile template\n"
        "# 実際のD-Claw環境では、D-Claw/ClawpackのMakefile構成に合わせて調整してください。\n\n"
        ".PHONY: help\n"
        "help:\n"
        "\t@echo \"D-Claw demo application template\"\n"
        "\t@echo \"1) D-ClawをLinux/Unix環境へインストール\"\n"
        "\t@echo \"2) setrun.pyを公式テンプレートに合わせて調整\"\n"
        "\t@echo \"3) make .data\"\n"
        "\t@echo \"4) make .output\"\n"
        "\t@echo \"5) make .plots\"\n",
        encoding="utf-8",
    )


def write_readme(path: Path) -> None:
    path.write_text(
        "# D-Claw basic granular-flow application template\n\n"
        "このフォルダは、D-Clawの基本アプリケーションを作るための入力雛形です。\n\n"
        "## 含まれるファイル\n\n"
        "- topo_valley.asc: 人工谷地形\n"
        "- qinit_depth.asc: 初期流動深\n"
        "- qinit_solid_fraction.asc: 固相体積濃度の初期値\n"
        "- qinit_pore_pressure.asc: 間隙圧相当量の初期値\n"
        "- setrun.py: D-Claw設定ファイルの雛形\n"
        "- Makefile: 実行手順メモ\n\n"
        "D-Clawの実行にはLinux/Unix環境とD-Claw本体が必要です。\n",
        encoding="utf-8",
    )


def main() -> None:
    parser = argparse.ArgumentParser(description="Create a simple D-Claw application template")
    parser.add_argument("--outdir", default="out_dclaw_basic_app")
    parser.add_argument("--nx", type=int, default=160)
    parser.add_argument("--ny", type=int, default=90)
    parser.add_argument("--dx", type=float, default=2.0)
    args = parser.parse_args()

    outdir = Path(args.outdir)
    outdir.mkdir(parents=True, exist_ok=True)

    _, _, topo = make_valley_topography(args.nx, args.ny, args.dx)
    q = make_initial_debris(args.nx, args.ny, args.dx)

    write_esri_ascii(outdir / "topo_valley.asc", topo, args.dx)
    write_esri_ascii(outdir / "qinit_depth.asc", q["h"], args.dx)
    write_esri_ascii(outdir / "qinit_solid_fraction.asc", q["solid_fraction"], args.dx)
    write_esri_ascii(outdir / "qinit_pore_pressure.asc", q["pore_pressure"], args.dx)
    write_setrun_template(outdir / "setrun.py", args.nx, args.ny, args.dx)
    write_makefile(outdir / "Makefile")
    write_readme(outdir / "README.md")

    print(f"D-Clawアプリ雛形を作成しました: {outdir}")
    print("D-Claw実行環境で setrun.py を公式テンプレートに合わせて調整してください。")


if __name__ == "__main__":
    main()