"""
TorsoCAE Validation -- thermal_only case 2: cylindrical wall radial conduction.

Classic textbook steady-state radial conduction through a hollow cylinder
(pipe) wall, inner radius a at T_in, outer radius b at T_out. Exact
solution (1D radial Fourier conduction in cylindrical coordinates):

    T(r) = T_in + (T_out - T_in) * ln(r/a) / ln(b/a)
    q(r) = k*(T_in - T_out) / (r * ln(b/a))          [W/m^2, radial flux]
    Q    = 2*pi*k*L*(T_in - T_out) / ln(b/a)          [W, total radial heat rate]

Geometry: a 90-degree WEDGE of the tube (same CSG construction as
structural/elastoplastic/case_02_thick_cylinder_pressure -- cylinder minus
cylinder, two csg_slice cuts), which for a purely radial-conduction problem
needs no explicit symmetry BC: the default adiabatic (zero-flux) condition
on the wedge's flat radial-cut faces IS the correct physical condition
(no heat flows across a radial plane in an axisymmetric problem), and
insulating the axial (z) end faces approximates the long/2D-radial
idealization the exact solution assumes.

Run: HWLOC_COMPONENTS=-gl python3 validation/thermal/thermal_only/case_02_cylindrical_wall_conduction/run_validation.py
"""
import json
import math
import sys
import time
from pathlib import Path

import numpy as np

CASE_DIR = Path(__file__).resolve().parent
sys.path.insert(0, str(CASE_DIR.parents[3] / "server"))
from pytorsocae import TorsoCAESession

A, B, LZ = 0.05, 0.10, 0.02     # m -- inner r, outer r, axial slice length
K = 50.0                        # W/m-K
T_IN, T_OUT = 400.0, 300.0      # K


def t_exact(r: np.ndarray) -> np.ndarray:
    return T_IN + (T_OUT - T_IN) * np.log(r / A) / math.log(B / A)


def q_exact(r: np.ndarray) -> np.ndarray:
    return K * (T_IN - T_OUT) / (r * math.log(B / A))


def build_wedge(session):
    session.csg_make_shape("cylinder", {"r": B, "h": LZ, "cx": 0, "cy": 0, "cz": 0}, name="outer", body_id="body_0")
    session.csg_make_shape("cylinder", {"r": A, "h": LZ, "cx": 0, "cy": 0, "cz": 0}, name="inner", body_id="body_1")
    session.csg_boolean("cut", "body_0", "body_1", name="tube", body_id="body_2")
    session.csg_delete("body_0")
    session.csg_delete("body_1")
    session.csg_slice("body_2", p=[0, 0, 0], n=[1, 0, 0], name="q1", result_body_id="body_2")
    session.csg_slice("body_2", p=[0, 0, 0], n=[0, 1, 0], name="q2", result_body_id="body_2")
    session.csg_select("body_2")
    session.mesh(mesh_id="mesh_0", algo_id="hex_transfinite", size_factor=0.4, size_min=0, size_max=0, order=1, dim=3)
    # Face tags verified empirically (same construction as the proven structural
    # thick-cylinder case): 1 = outer cyl (r=b), 6 = inner cyl (r=a); 2,3,4,5
    # (z-ends, radial-cut symmetry planes) left at the default adiabatic BC.
    session.surface(1, scope_id="body_2").bc("fixed_temp", values=[T_OUT])
    session.surface(6, scope_id="body_2").bc("fixed_temp", values=[T_IN])
    return "body_2"


def run_case() -> dict:
    session = TorsoCAESession()
    build_wedge(session)
    session.solid("tube").material(k=K)
    session.set_physics("thermal", submodel="thermal_only", backend="dolfinx")
    session.set_solver_options(algo="cg", precond="gamg", tol=1e-12, max_iter=1000, n_cores=2, device="cpu")
    t0 = time.time()
    result = session.compute(mesh_ids=["mesh_0"])
    elapsed = time.time() - t0

    npz = np.load(result["npz_path"])
    coords = npz["coordinates"]
    temp = npz["temperature"]
    r = np.hypot(coords[:, 0], coords[:, 1])
    t_ref = t_exact(r)
    temp_err_pct = np.abs(temp - t_ref) / (T_IN - T_OUT) * 100.0
    max_temp_err_pct = float(temp_err_pct.max())

    flux = npz["heat_flux"]
    flux_coords = npz["heat_flux_coords"]
    r_flux = np.hypot(flux_coords[:, 0], flux_coords[:, 1])
    q_ref = q_exact(r_flux)
    flux_err_pct = np.abs(flux - q_ref) / q_ref * 100.0
    max_flux_err_pct = float(np.median(flux_err_pct))   # median: flux singular-ish near sharp mesh corners, see README

    return {
        "case": "thermal_only / cylindrical wall radial conduction",
        "citation": "Standard 1D radial Fourier conduction in a hollow cylinder (textbook exact solution)",
        "geometry": {"a_m": A, "b_m": B, "axial_slice_m": LZ, "model": "90-degree wedge"},
        "material": {"k_W_mK": K},
        "bc": {"T_in_K": T_IN, "T_out_K": T_OUT},
        "method": "FEM temperature at every mesh node compared to the exact logarithmic "
                  "profile T(r) = T_in + (T_out-T_in)*ln(r/a)/ln(b/a); FEM heat flux "
                  "magnitude (per-cell) compared to the exact q(r) = k*(T_in-T_out)/(r*ln(b/a)) "
                  "(median relative error reported for flux, since DG0 cell-center flux "
                  "samples very close to the curved boundaries are mesh-resolution-limited, "
                  "not a solver defect -- see README).",
        "max_temp_error_pct": max_temp_err_pct,
        "median_flux_error_pct": max_flux_err_pct,
        "n_nodes": int(coords.shape[0]),
        "pass": bool(max_temp_err_pct < 2.0 and max_flux_err_pct < 5.0),
        "elapsed_s": elapsed,
        "profile": {
            "r_m": r.tolist(),
            "t_fem_K": temp.tolist(),
        },
    }


def plot_profile(summary: dict) -> None:
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    prof = summary["profile"]
    r = np.array(prof["r_m"])
    t_fem = np.array(prof["t_fem_K"])
    r_curve = np.linspace(A, B, 200)

    fig, ax = plt.subplots(figsize=(6.4, 4.4), dpi=150)
    ax.plot(r_curve, t_exact(r_curve), "-", color="#dc2626", linewidth=2, label="Exact log profile")
    ax.plot(r, t_fem, "o", color="#2563eb", markersize=3, alpha=0.4, label="TorsoCAE FEM (all nodes)")
    ax.set_xlabel("r (m)")
    ax.set_ylabel("Temperature (K)")
    ax.set_title(f"Cylindrical wall conduction: FEM vs. exact (max err {summary['max_temp_error_pct']:.2e}%)")
    ax.legend(loc="best")
    ax.grid(alpha=0.3)
    fig.tight_layout()
    fig.savefig(CASE_DIR / "plots" / "temperature_profile.png")


def main():
    summary = run_case()
    profile = summary.pop("profile")
    (CASE_DIR / "results.json").write_text(json.dumps(summary, indent=2))
    print(json.dumps(summary, indent=2), flush=True)
    summary["profile"] = profile
    plot_profile(summary)


if __name__ == "__main__":
    main()