make2-uniform.py

#!/usr/bin/env python3
"""
make2_thermal_uniform_csv.py
Spatially uniform thermal-history CSV helper for the CA dendrite model.
Generated with AI assistance from an older F77 helper routine.

Purpose
-------
This is the uniform/equiaxed-growth companion to make_thermal_csv.py.
The moving-gradient helper writes a thermal field with a nonzero spatial
gradient.  This helper writes a spatially uniform macro-temperature history:

    Tref(t) = Tref0 - cooling_rate*t
    Gx      = 0
    Gy      = 0

Thus every CA cell receives the same imposed macro temperature before any
local perturbation, latent-heat, or cell-scale thermal model is applied inside
CAden.

Parameter-file conventions
--------------------------
The input file uses named value pairs read from the CA parameter file.
The values used here are:

    TMP, MLIQ, C0   -> TL0 = TMP + MLIQ*C0
    DTSEED          -> default initial undercooling
    FIN             -> default final time
    DX              -> CA cell size in m
    GRAD, VEL       -> optional default cooling rate, GRAD*VEL

Here GRAD is the CA moving-gradient value in K/cell and VEL is in cells/s.
Their product is therefore K/s and can be reused as a convenient default
uniform cooling rate.

CSV column meaning
------------------
The output file uses the six-column LOADTHERMCSV handshake form:

    time_s
        Simulation time in seconds.

    Tref_C
        Imposed reference macro-temperature in degC.  For this uniform file,
        this is also the macro-temperature seen everywhere in the CA window.

    Gx_K_per_m, Gy_K_per_m
        Macro thermal-gradient components in K/m.  Both are intentionally
        zero for uniform cooling.  This is what makes the imposed field
        equiaxed/isotropic rather than directional.

    cooling_rate_K_per_s
        Positive cooling rate used to decrease Tref with time:
            Tref(t) = Tref0 - cooling_rate*t

    Vthermal_m_per_s
        Velocity scale passed to the CA thermal/KGT machinery.  In this
        uniform helper it does NOT create a moving gradient because Gx=Gy=0.
        It is retained because the present CA code can use Vthermal as the
        runtime KGT/staged-growth closure velocity scale.

Important distinction
---------------------
For uniform cooling, cooling_rate_K_per_s controls the imposed macro cooling.
Vthermal_m_per_s controls the analytical tip/freezing closure scale.  They are
allowed to be conceptually separate.  The default Vthermal is VEL*DX from the
parameter file; set --kgt-velocity-cells 0.0 only when intentionally disabling
that KGT/staged-growth velocity scale.

Cooling-rate choices
--------------------
Default:
    cooling_rate = GRAD*VEL, if both are present in the parameter file.

Manual alternatives:
    --cooling-rate R      use R directly in K/s
    --total-drop DTEMP    use cooling_rate = DTEMP/FIN

For equiaxed tests, --cooling-rate is often the clearest choice.
"""

from __future__ import annotations

import argparse
import csv
from pathlib import Path


# Read CA-style named-value pairs and ignore nonnumeric entries.
def read_params(path: Path) -> dict[str, float]:
    params: dict[str, float] = {}

    with path.open("r", encoding="utf-8", errors="replace") as f:
        for line in f:
            s = line.strip()
            if not s or s.startswith("#"):
                continue

            # Allow trailing comments.
            if "#" in s:
                s = s.split("#", 1)[0].strip()
            if not s:
                continue

            parts = s.split()
            if len(parts) < 2:
                continue

            name = parts[0].upper()
            try:
                value = float(parts[1])
            except ValueError:
                continue

            params[name] = value

    return params


def read_param_text(path: Path, name: str, default: str) -> str:
    """Read one named text value from the CA parameter file."""
    key = name.upper()
    if not path.exists():
        return default

    with path.open("r", encoding="utf-8", errors="replace") as f:
        for line in f:
            s = line.strip()
            if not s or s.startswith("#"):
                continue
            if "#" in s:
                s = s.split("#", 1)[0].strip()
            if not s:
                continue

            parts = s.split(None, 1)
            if len(parts) == 2 and parts[0].upper() == key:
                return parts[1].strip()

    return default


def get_required(params: dict[str, float], name: str) -> float:
    key = name.upper()
    if key not in params:
        raise KeyError(f"Required parameter {name!r} not found.")
    return params[key]


def ask_string(prompt: str, default: str) -> str:
    ans = input(f"{prompt} [{default}]: ").strip()
    return ans if ans else default


def ask_float(prompt: str, default: float) -> float:
    while True:
        ans = input(f"{prompt} [{default:g}]: ").strip()
        if ans == "":
            return default
        try:
            return float(ans)
        except ValueError:
            print("Please enter a numeric value.")


def ask_int(prompt: str, default: int, min_value: int | None = None) -> int:
    while True:
        ans = input(f"{prompt} [{default}]: ").strip()
        if ans == "":
            value = default
        else:
            try:
                value = int(ans)
            except ValueError:
                print("Please enter an integer value.")
                continue
        if min_value is not None and value < min_value:
            print(f"Please enter a value >= {min_value}.")
            continue
        return value


def ask_optional_float(prompt: str, default: float | None, shown_default: float) -> float | None:
    label = f"{shown_default:g}" if default is None else f"{default:g}"
    ans = input(f"{prompt} [{label}]: ").strip()
    if ans == "":
        return default
    try:
        return float(ans)
    except ValueError:
        print("Please enter a numeric value; using default.")
        return default


def ask_yes_no(prompt: str, default: bool = False) -> bool:
    suffix = "[Y/n]" if default else "[y/N]"
    ans = input(f"{prompt} {suffix}: ").strip().lower()
    if ans == "":
        return default
    return ans in ("y", "yes")


# Interactive choice for the imposed uniform macro cooling rate.
def ask_cooling_choice(params: dict[str, float], fin: float) -> tuple[float | None, float | None]:
    default_rate = params.get("GRAD", 0.0) * params.get("VEL", 0.0)
    default_drop = default_rate * fin

    print()
    print("Uniform cooling-rate choice")
    print("  1) use GRAD*VEL from parameter file")
    print("  2) enter cooling rate directly (K/s)")
    print("  3) enter total temperature drop over FIN (K)")

    while True:
        choice = input("Choice [1]: ").strip().lower()
        if choice in ("", "1"):
            return None, None
        if choice == "2":
            rate = ask_float("Cooling rate (K/s)", default_rate)
            return rate, None
        if choice == "3":
            drop = ask_float("Total temperature drop over FIN (K)", default_drop)
            return None, drop
        print("Please choose 1, 2, or 3.")


# Interactive choice for Vthermal: the KGT/staged-growth velocity scale.
def ask_kgt_velocity_choice(params: dict[str, float], dx: float) -> tuple[float | None, float | None]:
    default_cells = params.get("VEL", 0.0)
    default_um_s = default_cells * dx * 1.0e6

    print()
    print("KGT/staged-growth velocity scale written to Vthermal")
    print("  1) use VEL from parameter file (cells/s)")
    print("  2) enter velocity in cells/s")
    print("  3) enter velocity in um/s")
    print("  4) write zero velocity")

    while True:
        choice = input("Choice [1]: ").strip().lower()
        if choice in ("", "1"):
            return None, None
        if choice == "2":
            return ask_float("KGT velocity scale (cells/s)", default_cells), None
        if choice == "3":
            return None, ask_float("KGT velocity scale (um/s)", default_um_s)
        if choice == "4":
            return 0.0, None
        print("Please choose 1, 2, 3, or 4.")


def print_input_summary(
    params: dict[str, float],
    fin: float,
    nrows: int,
    tref0: float,
    cooling_rate: float,
    rate_source: str,
    kgt_velocity_cells: float,
    velocity_source: str,
) -> None:
    tmp = get_required(params, "TMP")
    mliq = get_required(params, "MLIQ")
    c0 = get_required(params, "C0")
    dx = get_required(params, "DX")
    dtseed = params.get("DTSEED", 0.0)
    tl0 = tmp + mliq*c0
    v_m_s = kgt_velocity_cells*dx

    print()
    print("Current uniform thermal CSV settings")
    print("----------------------------------------------")
    print(f"TL0 = TMP + MLIQ*C0   : {tl0:.6g} degC")
    print(f"DTSEED                : {dtseed:.6g} K")
    print(f"Tref0                 : {tref0:.6g} degC")
    print(f"FIN                   : {fin:.6g} s")
    print(f"Rows                  : {nrows}")
    print(f"Gx, Gy                : 0, 0 K/m")
    print(f"Cooling rate          : {cooling_rate:.6e} K/s")
    print(f"Cooling-rate source   : {rate_source}")
    print(f"Total temperature drop: {cooling_rate*fin:.6g} K")
    print(f"DX                    : {dx:.6e} m")
    print(f"KGT velocity scale    : {kgt_velocity_cells:.6g} cells/s")
    print(f"Vthermal              : {v_m_s:.6e} m/s")
    print(f"Velocity source       : {velocity_source}")
    print("----------------------------------------------")


def interactive_args(args: argparse.Namespace) -> argparse.Namespace:
    print()
    print("==============================================")
    print("CA UNIFORM THERMAL CSV HELPER")
    print("==============================================")

    default_param = args.param_file or "xparam.txt"
    args.param_file = ask_string("Parameter file", default_param)

    param_path = Path(args.param_file)
    params = read_params(param_path)

    if args.output is None:
        args.output = read_param_text(param_path, "THFIL", "thermal_uniform.csv")

    out_path = Path(args.output)

    tmp = get_required(params, "TMP")
    mliq = get_required(params, "MLIQ")
    c0 = get_required(params, "C0")
    dx = get_required(params, "DX")
    dtseed = params.get("DTSEED", 0.0)
    default_fin = args.fin if args.fin is not None else params.get("FIN", 1.0)
    default_tref0 = tmp + mliq*c0 - dtseed
    current_tref0 = args.tref0 if args.tref0 is not None else default_tref0

    default_output = args.output or read_param_text(param_path, "THFIL", "thermal_uniform.csv")
    args.output = ask_string("Output uniform thermal CSV", default_output)
    args.nrows = ask_int("Number of time rows", args.nrows, min_value=2)
    args.fin = ask_float("Final time FIN (s)", float(default_fin))

    print()
    print("Tref0 controls the uniform macro temperature at t=0.")
    print("Press Enter to use TL0 - DTSEED from the parameter file.")
    args.tref0 = ask_optional_float("Initial Tref0 (degC)", args.tref0, current_tref0)

    args.cooling_rate, args.total_drop = ask_cooling_choice(params, args.fin)
    args.kgt_velocity_cells, args.kgt_velocity_um_s = ask_kgt_velocity_choice(params, dx)

    actual_tref0 = args.tref0 if args.tref0 is not None else default_tref0
    cooling_rate, rate_source = choose_cooling_rate(
        params, args.fin, args.cooling_rate, args.total_drop
    )

    if args.kgt_velocity_um_s is not None:
        kgt_velocity_cells = args.kgt_velocity_um_s * 1.0e-6 / dx
        velocity_source = "interactive um/s input"
    elif args.kgt_velocity_cells is not None:
        kgt_velocity_cells = args.kgt_velocity_cells
        velocity_source = "interactive cells/s input"
    else:
        kgt_velocity_cells = params.get("VEL", 0.0)
        velocity_source = "parameter-file VEL"

    print_input_summary(
        params, args.fin, args.nrows, actual_tref0,
        cooling_rate, rate_source, kgt_velocity_cells, velocity_source
    )

    write_now = ask_yes_no("Write this uniform thermal CSV", default=True)
    if not write_now:
        print("No file written.")
        raise SystemExit(0)

    return args

# Resolve precedence for cooling-rate inputs: direct rate, total drop, or GRAD*VEL.
def choose_cooling_rate(
    params: dict[str, float],
    fin: float,
    cooling_rate_arg: float | None,
    total_drop_arg: float | None,
) -> tuple[float, str]:
    if cooling_rate_arg is not None and total_drop_arg is not None:
        raise ValueError("Use either --cooling-rate or --total-drop, not both.")

    if cooling_rate_arg is not None:
        return cooling_rate_arg, "command line --cooling-rate"

    if total_drop_arg is not None:
        return total_drop_arg / fin, "command line --total-drop / final time"

    if "GRAD" in params and "VEL" in params:
        return params["GRAD"] * params["VEL"], "parameter-file GRAD*VEL"

    raise KeyError(
        "No cooling rate available. Provide --cooling-rate, --total-drop, "
        "or include both GRAD and VEL in the parameter file."
    )


def main() -> int:
    parser = argparse.ArgumentParser(
        description="Generate uniform-cooling thermal CSV for equiaxed CA growth."
    )
    parser.add_argument(
        "param_file",
        nargs="?",
        help="CA parameter file, e.g. xparam.txt or CAparam.txt",
    )
    parser.add_argument(
        "-i", "--interactive",
        action="store_true",
        help="Use prompt-driven interactive mode.",
    )
    parser.add_argument(
        "-o", "--output",
        default=None,
        help="Output CSV filename. Default uses THFIL from the parameter file.",
    )
    parser.add_argument(
        "--nrows",
        type=int,
        default=501,
        help="Number of time rows to write [501]",
    )
    parser.add_argument(
        "--tref0",
        type=float,
        default=None,
        help="Initial Tref in degC. Default is TL0-DTSEED.",
    )
    parser.add_argument(
        "--fin",
        type=float,
        default=None,
        help="Override final time in seconds. Default uses FIN from parameter file.",
    )
    parser.add_argument(
        "--cooling-rate",
        type=float,
        default=None,
        help="Uniform cooling rate in K/s. Default uses GRAD*VEL from parameter file.",
    )
    parser.add_argument(
        "--total-drop",
        type=float,
        default=None,
        help="Total temperature drop over FIN in K. Sets cooling_rate=drop/FIN.",
    )
    parser.add_argument(
        "--kgt-velocity-cells",
        type=float,
        default=None,
        help=(
            "Velocity scale in cells/s written to Vthermal for KGT/staged "
            "growth closure. Default uses VEL from parameter file; use 0.0 "
            "to intentionally disable KGT-staged growth."
        ),
    )
    parser.add_argument(
        "--kgt-velocity-um-s",
        type=float,
        default=None,
        help=(
            "Velocity scale in um/s written to Vthermal for KGT/staged "
            "growth closure. Overrides --kgt-velocity-cells if supplied."
        ),
    )
    parser.add_argument(
        "--quiet",
        action="store_true",
        help="Suppress printed summary.",
    )

    args = parser.parse_args()

    if args.interactive or args.param_file is None:
        args = interactive_args(args)

    param_path = Path(args.param_file)

    params = read_params(param_path)

    if args.output is None:
        args.output = read_param_text(param_path, "THFIL", "thermal_uniform.csv")

    out_path = Path(args.output)

    tmp = get_required(params, "TMP")
    mliq = get_required(params, "MLIQ")
    c0 = get_required(params, "C0")
    dx = get_required(params, "DX")
    fin = args.fin if args.fin is not None else params.get("FIN", 1.0)
    dtseed = params.get("DTSEED", 0.0)

    if args.nrows < 2:
        raise ValueError("nrows must be at least 2.")
    if fin <= 0.0:
        raise ValueError("Final time must be positive.")
    if dx <= 0.0:
        raise ValueError("DX must be positive.")

    tl0 = tmp + mliq*c0
    tref0 = args.tref0 if args.tref0 is not None else tl0 - dtseed

    cooling_rate, rate_source = choose_cooling_rate(
        params, fin, args.cooling_rate, args.total_drop
    )
    if cooling_rate < 0.0:
        raise ValueError("Cooling rate should be non-negative for temperature decrease.")

    # Uniform/equiaxed thermal field: no imposed spatial gradient.
    gx_k_m = 0.0
    gy_k_m = 0.0

    if args.kgt_velocity_um_s is not None:
        kgt_velocity_cells = args.kgt_velocity_um_s * 1.0e-6 / dx
        velocity_source = "command line --kgt-velocity-um-s"
    elif args.kgt_velocity_cells is not None:
        kgt_velocity_cells = args.kgt_velocity_cells
        velocity_source = "command line --kgt-velocity-cells"
    else:
        kgt_velocity_cells = params.get("VEL", 0.0)
        velocity_source = "parameter-file VEL"

    if kgt_velocity_cells < 0.0:
        raise ValueError("KGT velocity scale should be non-negative.")

    # Vthermal is retained as the analytical KGT/staged-growth closure scale.
    vthermal_m_s = kgt_velocity_cells * dx

    with out_path.open("w", newline="", encoding="utf-8") as f:
        writer = csv.writer(f)
        writer.writerow([
            "time_s",
            "Tref_C",
            "Gx_K_per_m",
            "Gy_K_per_m",
            "cooling_rate_K_per_s",
            "Vthermal_m_per_s",
        ])

        for i in range(args.nrows):
            a = i / float(args.nrows - 1)
            t = a*fin
            tref = tref0 - cooling_rate*t
            writer.writerow([
                f"{t:.8e}",
                f"{tref:.8e}",
                f"{gx_k_m:.8e}",
                f"{gy_k_m:.8e}",
                f"{cooling_rate:.8e}",
                f"{vthermal_m_s:.8e}",
            ])

    if not args.quiet:
        print("Generated:", out_path)
        print("From parameter file:", param_path)
        print(f"Mode                 = uniform isotropic cooling")
        print(f"TL0                  = {tl0:.6g} degC")
        print(f"Tref0                = {tref0:.6g} degC")
        print(f"Gx_K_per_m            = {gx_k_m:.6e} K/m")
        print(f"Gy_K_per_m            = {gy_k_m:.6e} K/m")
        print(f"DX                    = {dx:.6e} m")
        print(f"KGT velocity scale    = {kgt_velocity_cells:.6g} cells/s")
        print(f"Vthermal_m_per_s      = {vthermal_m_s:.6e} m/s")
        print(f"velocity source       = {velocity_source}")
        print(f"cooling_rate          = {cooling_rate:.6e} K/s")
        print(f"cooling-rate source   = {rate_source}")
        print(f"total temperature drop= {cooling_rate*fin:.6g} K")
        print(f"final time            = {fin:.6g} s")
        print(f"rows                  = {args.nrows}")

    return 0


if __name__ == "__main__":
    raise SystemExit(main())