make2-thermal.py

#!/usr/bin/env python3
"""
make2_thermal_csv.py
Build a simple one-dimensional thermal-history CSV file for the CA dendrite code.
Generated with AI assistance from an older F77 helper routine.

Purpose
-------
The CA code can either impose a built-in moving thermal gradient or read an
external thermal history file through LOADTHERMCSV.  This helper creates the
external CSV form from the same compact parameter file used by the CA run.

It is mainly a convenience/consistency tool: the velocity, gradient, final time,
cell size, and starting temperature are taken from the CA parameter file so that
the thermal CSV and the solidification run remain synchronized.

Input parameter conventions
---------------------------
The parameter file is assumed to contain named value pairs, for example:

    TMP     327.502
    MLIQ    -2.326
    C0      8.0
    DX      1.0E-6
    FIN     10.0
    GRAD    0.015
    VEL     200.0
    DTSEED  1.0

The important CA-side meanings are:

    TMP, MLIQ, C0   define the nominal liquidus temperature
                    TL0 = TMP + MLIQ*C0

    DTSEED          initial undercooling used to set the default Tref0
                    Tref0 = TL0 - DTSEED

    DX              CA cell size, m/cell

    GRAD            thermal gradient stored in the CA parameter file as K/cell

    VEL             imposed thermal/front speed stored as cells/s

    FIN             final time for the thermal history, unless overridden

Output CSV columns and meaning
------------------------------
This helper writes the explicit six-column form accepted by LOADTHERMCSV:

    time_s,Tref_C,Gx_K_per_m,Gy_K_per_m,cooling_rate_K_per_s,Vthermal_m_per_s

Column meanings:

    time_s
        Simulation time in seconds.

    Tref_C
        Reference temperature in degrees C at the CA thermal reference row.
        The CA reconstructs the local macro-temperature field from this value
        plus the supplied gradient.

    Gx_K_per_m
        X-component of the imposed thermal gradient in K/m.  This simple
        one-dimensional helper sets Gx = 0.

    Gy_K_per_m
        Y-component of the imposed thermal gradient in K/m.  This is converted
        from the CA parameter GRAD by Gy = GRAD/DX.

    cooling_rate_K_per_s
        Cooling rate at the reference row.  For the moving-gradient case,
        cooling_rate = G*V = (GRAD/DX)*(VEL*DX) = GRAD*VEL.

    Vthermal_m_per_s
        Thermal/front speed in m/s, converted from the CA parameter VEL by
        Vthermal = VEL*DX.

Thermal field represented
-------------------------
The generated reference history is linear in time:

    Tref(t) = Tref0 - cooling_rate*t

and the CA-side macro field is equivalent to:

    T(j,t) = Tref(t) + G*(j - Jref)*DX

where Jref is the reference row used internally by the CA thermal logic.

IEMBED note
-----------
For IEMBED=0, the CA code can apply a startup TSHIFT to align the initial
seed/front row to TL0-DTSEED.  In that mode, the exact absolute Tref0 is less
critical.  For IEMBED=1, Tref0 is used directly at the CA reference row, so the
absolute starting temperature matters.
"""

from __future__ import annotations

import argparse
import csv
from pathlib import Path


def read_params(path: Path) -> dict[str, float]:
    """Read numeric named-value pairs from the CA parameter file."""
    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:
    """Return a required CA parameter or stop with a clear error."""
    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")


def print_input_summary(params: dict[str, float], fin: float, nrows: int, tref0: float) -> None:
    """Print the CA-to-thermal conversion summary before writing the CSV."""
    tmp = get_required(params, "TMP")
    mliq = get_required(params, "MLIQ")
    c0 = get_required(params, "C0")
    dx = get_required(params, "DX")
    vel_cells = get_required(params, "VEL")
    grad_cell = get_required(params, "GRAD")
    dtseed = params.get("DTSEED", 0.0)
    tl0 = tmp + mliq*c0
    v_m_s = vel_cells*dx
    g_k_m = grad_cell/dx
    cooling_rate = grad_cell*vel_cells

    print()
    print("Current 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"DX                    : {dx:.6e} m")
    print(f"VEL                   : {vel_cells:.6g} cells/s")
    print(f"Vthermal              : {v_m_s:.6e} m/s")
    print(f"GRAD                  : {grad_cell:.6g} K/cell")
    print(f"G                     : {g_k_m:.6e} K/m")
    print(f"Cooling rate          : {cooling_rate:.6e} K/s")
    print("----------------------------------------------")


def interactive_args(args: argparse.Namespace) -> argparse.Namespace:
    """Prompt for the few values most likely to change between runs."""
    print()
    print("==============================================")
    print("CA 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", "xthermal.csv")

    out_path = Path(args.output)

    tmp = get_required(params, "TMP")
    mliq = get_required(params, "MLIQ")
    c0 = get_required(params, "C0")
    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", "xthermal.csv")
    args.output = ask_string("Output 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 reference-row temperature at t=0.")
    print("Press Enter to use TL0 - DTSEED from the parameter file.")
    tref0_value = ask_optional_float("Initial Tref0 (degC)", args.tref0, current_tref0)
    args.tref0 = tref0_value

    actual_tref0 = args.tref0 if args.tref0 is not None else default_tref0
    print_input_summary(params, args.fin, args.nrows, actual_tref0)

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

    return args


def main() -> int:
    parser = argparse.ArgumentParser(
        description="Generate thermal CSV from CA parameter file."
    )
    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(
        "--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", "xthermal.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")
    vel_cells = get_required(params, "VEL")
    grad_cell = get_required(params, "GRAD")

    dtseed = params.get("DTSEED", 0.0)
    fin = args.fin if args.fin is not None else params.get("FIN", 1.0)

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

    # Nominal alloy liquidus and default starting reference temperature.
    tl0 = tmp + mliq*c0
    tref0 = args.tref0 if args.tref0 is not None else tl0 - dtseed

    # Convert CA units to physical units expected by LOADTHERMCSV.
    v_m_s = vel_cells*dx
    g_k_m = grad_cell/dx
    cooling_rate = g_k_m*v_m_s

    # Equivalent and often numerically clearer:
    # cooling_rate = grad_cell*vel_cells

    with out_path.open("w", newline="", encoding="utf-8") as f:
        writer = csv.writer(f)
        # Six-column explicit form: time, Tref, gradient vector, cooling rate, velocity.
        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

            # Reference-row temperature decreases linearly with imposed cooling.
            tref = tref0 - cooling_rate*t
            writer.writerow([
                f"{t:.8e}",
                f"{tref:.8e}",
                f"{0.0:.8e}",
                f"{g_k_m:.8e}",
                f"{cooling_rate:.8e}",
                f"{v_m_s:.8e}",
            ])

    if not args.quiet:
        print("Generated:", out_path)
        print("From parameter file:", param_path)
        print(f"TL0                  = {tl0:.6g} degC")
        print(f"Tref0                = {tref0:.6g} degC")
        print(f"VEL                  = {vel_cells:.6g} cells/s")
        print(f"DX                   = {dx:.6g} m")
        print(f"Vthermal_m_per_s     = {v_m_s:.6e} m/s")
        print(f"GRAD                 = {grad_cell:.6g} K/cell")
        print(f"G_K_per_m            = {g_k_m:.6e} K/m")
        print(f"cooling_rate         = {cooling_rate:.6e} K/s")
        print(f"final time           = {fin:.6g} s")
        print(f"rows                 = {args.nrows}")

    return 0


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