make1-param.py

#!/usr/bin/env python3
"""
CA DENDRITIC SOLIDIFICATION PARAMETER HELPER
D.R. Johnson

make1-param.py
Lean interactive parameter-file editor for the CA dendritic solidification
model.  Generated with AI assistance from an older F77 helper routine.

Purpose
-------
Create or edit parameter files read by the Fortran CA program.
The output file is intentionally plain:

PARAM    value

with no comments, so it remains easy for the Fortran reader to parse.

Conventions
 -----------
ATHERM
Effective thermal smoothing coefficient [m^2/s].  It controls the
explicit local spreading/smoothing of the latent-heat temperature
field generated by interface growth.

LHEAT
Effective latent-heat scale [J/kg].  In the CA program this sets the
temperature source strength associated with newly frozen subcell
capture, rather than representing a rigorously calibrated alloy
latent heat.

CPHEAT
Effective heat-capacity scale [J/(kg K)].  Together with LHEAT, it
sets the local temperature perturbation scale through roughly
LHEAT/CPHEAT.

 These terms are best viewed as optional/local latent-heat feedback controls
 for numerical experiments, not as validated thermal-property
 inputs.
"""

import sys
import re
import math
from pathlib import Path


# =====================================================================
# PARAMETER TABLE
# =====================================================================
# Each row is:
#   NAME, TYPE, DEFAULT, DESCRIPTION
#
# Keep this table as the main editable source of parameter definitions.
# The file writer preserves this order, then appends any extra parameters
# found in an existing file or added through the menu.
# =====================================================================

PARAM_TABLE = [
    # --- Alloy / phase diagram parameters ---
    ("TMP",     float, 327.502,       "Melting point (degC)"),
    ("MLIQ",    float, -2.326,        "Liquidus slope (degC/wt%)"),
    ("KALPHA",  float, 0.296,         "Partition coefficient k"),
    ("C0",      float, 8.0,           "Nominal composition (wt%)"),
    ("CEUT",    float, 61.9,          "Eutectic composition (wt%)"),

    # --- Interface / anisotropy parameters ---
    ("GIBT",    float, 2.0e-7,        "Gibbs-Thomson coefficient (K*m)"),
    ("DD",      float, 0.05,          "Anisotropy strength"),
    ("THETA1",  float, 1.22173,       "Preferred angle 1"),
    ("THETA2",  float, 1.91986,       "Preferred angle 2"),

    # --- Transport / grid parameters ---
    ("DCOEFF",  float, 1.0e-9,        "Solute diffusivity (m^2/s)"),
    ("DT",      float, 2.0e-4,        "Time step (s)"),
    ("DX",      float, 1.0e-6,        "Cell size (m)"),
    ("DMASK",   float, 0.60,          "Liquid mask threshold"),

    # --- Run-control parameters ---
    ("FIN",     float, 10.0,          "Simulation end time (s)"),
    ("DTSEED",  float, 1.0,           "Initial undercooling (K)"),
    ("SPARE",   float, 1.0,           "Spare multiplier"),

    # --- Thermal / latent-heat parameters ---
    ("GRAD",    float, 0.0150,        "Thermal gradient (K/cell)"),
    ("VEL",     float, 200.0,         "DS velocity (cells/s)"),
    ("ATHERM",  float, 2.5e-9,        "Thermal diffusivity"),
    ("LHEAT",   float, 1.0e4,         "Latent heat scaling"),
    ("CPHEAT",  float, 5.0e2,         "Heat capacity scale"),

    # --- Domain / option parameters ---
    ("XLEN",    int,   200,           "Domain width (cells)"),
    ("YLEN",    int,   400,           "Domain height (cells)"),
    ("IEMBED",  int,   0,             "Thermal embedding mode"),
    ("IBCTYPE", int,   1,             "Boundary condition type"),

    # --- File names ---
    ("PARFIL",  str,   "xparam.txt",   "Parameter file"),
    ("GRIDFIL", str,   "xgrid.txt",    "Starting grid"),
    ("THFIL",   str,   "xthermal.csv", "Thermal history file"),
]

PARAM_NAMES = [row[0] for row in PARAM_TABLE]
PARAM_TYPES = {row[0]: row[1] for row in PARAM_TABLE}
PARAM_DEFAULT = {row[0]: row[2] for row in PARAM_TABLE}
PARAM_DESC = {row[0]: row[3] for row in PARAM_TABLE}

SIGMA_STAR = 0.0253156


# =====================================================================
# BASIC HELPERS
# =====================================================================

def is_theta_param(name):
    return name.upper() in ("THETA1", "THETA2")


def parse_value(text):
    val = text.strip()
    try:
        if "." in val or "E" in val.upper():
            return float(val)
        return int(val)
    except ValueError:
        return val


def cast_value(name, text):
    if name in PARAM_TYPES:
        return PARAM_TYPES[name](text)
    return parse_value(text)


def make_default_params():
    params = {}
    for name in PARAM_NAMES:
        params[name] = PARAM_DEFAULT[name]
    return params


def theta_display(rad_value):
    try:
        rad = float(rad_value)
        deg = math.degrees(rad)
        return f"{deg:.4f} deg ({rad:.6g} rad)"
    except Exception:
        return str(rad_value)

def theta_display(rad_value):
    try:
        rad = float(rad_value)
        deg = math.degrees(rad)
        return f"{deg:.4f} deg ({rad:.6g} rad)"
    except Exception:
        return str(rad_value) 
 
def grad_vel_display(name, val, params):
    """Show stored CA units plus physical units for GRAD and VEL."""
    try:
        dx = float(params.get("DX", PARAM_DEFAULT["DX"]))
        xval = float(val)

        if name.upper() == "GRAD" and dx > 0.0:
            grad_k_cm = xval / dx / 100.0
            return f"{xval:g}  ({grad_k_cm:.2f} K/cm)"

        if name.upper() == "VEL" and dx > 0.0:
            vel_um_s = xval * dx * 1.0e6
            return f"{xval:g}  ({vel_um_s:.1f} um/s)"

    except Exception:
        pass

    return str(val)


def format_value(val):
    if isinstance(val, (float, int)):
        return f"{val:14.8E}"
    return str(val)


def format_value(val):
    if isinstance(val, (float, int)):
        return f"{val:14.8E}"
    return str(val)


# =====================================================================
# FILE I/O
# =====================================================================

def load_param_file(filename):
    params = {}
    if not Path(filename).exists():
        return params

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

            parts = re.split(r"\s+", line, 1)
            if len(parts) != 2:
                continue

            name = parts[0].upper()
            params[name] = parse_value(parts[1])

    return params


def save_param_file(filename, params):
    with open(filename, "w", encoding="utf-8") as f:
        for name in PARAM_NAMES:
            if name in params:
                f.write(f"{name:<8} {format_value(params[name])}\n")

        for name, val in params.items():
            if name not in PARAM_NAMES:
                f.write(f"{name:<8} {format_value(val)}\n")


# =====================================================================
# DERIVED CHECKS
# =====================================================================

def estimate_kgt_radius(params):
    mliq = float(params.get("MLIQ", PARAM_DEFAULT["MLIQ"]))
    k = float(params.get("KALPHA", PARAM_DEFAULT["KALPHA"]))
    c0 = float(params.get("C0", PARAM_DEFAULT["C0"]))
    giebt = float(params.get("GIBT", PARAM_DEFAULT["GIBT"]))
    dcoeff = float(params.get("DCOEFF", PARAM_DEFAULT["DCOEFF"]))
    dx = float(params.get("DX", PARAM_DEFAULT["DX"]))
    grad = float(params.get("GRAD", PARAM_DEFAULT["GRAD"]))
    vel = float(params.get("VEL", PARAM_DEFAULT["VEL"]))

    if k <= 1.0e-8 or dx <= 0.0 or dcoeff <= 0.0 or giebt <= 0.0:
        return 0.0, 0.0

    delt0 = mliq * c0 - mliq * (c0 / k)
    if delt0 <= 1.0e-12:
        return 0.0, 0.0

    v_m_s = vel * dx
    g_k_m = abs(grad) / dx
    if v_m_s <= 0.0:
        return 0.0, 0.0

    lchar = giebt / (SIGMA_STAR * delt0)
    rtip0 = math.sqrt(dcoeff * lchar / v_m_s)
    bb = 1.0 + (g_k_m * giebt / (v_m_s * delt0 * delt0))
    if bb < 1.0:
        bb = 1.0

    rtip = rtip0 / bb
    return rtip, rtip / dx


def ivantsov_2d(peclet):
    # Same 2-D Ivantsov form used by the CA output/CSV diagnostics:
    #   OmegaKGT2D = sqrt(pi*P) * exp(P) * erfc(sqrt(P))
    if peclet <= 0.0:
        return 0.0
    rootp = math.sqrt(peclet)
    return math.sqrt(math.pi * peclet) * math.exp(peclet) * math.erfc(rootp)


def estimate_kgt_values(params):
    mliq = float(params.get("MLIQ", PARAM_DEFAULT["MLIQ"]))
    k = float(params.get("KALPHA", PARAM_DEFAULT["KALPHA"]))
    c0 = float(params.get("C0", PARAM_DEFAULT["C0"]))
    giebt = float(params.get("GIBT", PARAM_DEFAULT["GIBT"]))
    dcoeff = float(params.get("DCOEFF", PARAM_DEFAULT["DCOEFF"]))
    dx = float(params.get("DX", PARAM_DEFAULT["DX"]))
    grad = float(params.get("GRAD", PARAM_DEFAULT["GRAD"]))
    vel = float(params.get("VEL", PARAM_DEFAULT["VEL"]))

    rtip, rhotip = estimate_kgt_radius(params)
    v_m_s = vel * dx
    g_k_m = abs(grad) / dx if dx > 0.0 else 0.0

    if rtip <= 0.0 or dcoeff <= 0.0 or v_m_s <= 0.0 or k <= 0.0:
        return None

    peclet = v_m_s * rtip / (2.0 * dcoeff)
    omega2d = ivantsov_2d(peclet)

    denom = 1.0 - (1.0 - k) * omega2d
    if denom <= 1.0e-12:
        return None

    cl_tip = c0 / denom
    dt_solute = abs(mliq) * (cl_tip - c0)
    dt_curve = 2.0 * giebt / rtip
    dt_thermal = g_k_m * dcoeff / v_m_s
    dt_kgt = dt_solute + dt_curve + dt_thermal

    bl_cells = 0.0
    if dx > 0.0 and v_m_s > 0.0:
        bl_cells = dcoeff / v_m_s / dx

    return {
        "rtip_m": rtip,
        "rhotip": rhotip,
        "peclet": peclet,
        "omega2d": omega2d,
        "cl_tip": cl_tip,
        "dt_solute": dt_solute,
        "dt_curve": dt_curve,
        "dt_thermal": dt_thermal,
        "dt_kgt": dt_kgt,
        "boundary_cells": bl_cells,
    }


def estimate_eutectic_values(params):
    tmp = float(params.get("TMP", PARAM_DEFAULT["TMP"]))
    mliq = float(params.get("MLIQ", PARAM_DEFAULT["MLIQ"]))
    k = float(params.get("KALPHA", PARAM_DEFAULT["KALPHA"]))
    c0 = float(params.get("C0", PARAM_DEFAULT["C0"]))
    ceut = float(params.get("CEUT", PARAM_DEFAULT["CEUT"]))

    teut = tmp + mliq * ceut
    fs_scheil = None
    if c0 > 0.0 and ceut > c0 and k != 1.0:
        try:
            fs_scheil = 1.0 - (ceut / c0) ** (1.0 / (k - 1.0))
        except Exception:
            fs_scheil = None

    return teut, fs_scheil


def derived_checks_text(params):
    lines = ["Derived checks:"]

    try:
        dcoeff = float(params.get("DCOEFF", PARAM_DEFAULT["DCOEFF"]))
        dt = float(params.get("DT", PARAM_DEFAULT["DT"]))
        dx = float(params.get("DX", PARAM_DEFAULT["DX"]))
        grad = float(params.get("GRAD", PARAM_DEFAULT["GRAD"]))
        vel = float(params.get("VEL", PARAM_DEFAULT["VEL"]))

        alpha = dcoeff * dt / (dx * dx)
        grad_k_cm = grad / dx / 100.0
        vel_um_s = vel * dx * 1.0e6
        vel_mm_h = vel * dx * 1.0e3 * 3600.0
        kgt = estimate_kgt_values(params)
        teut, fs_scheil = estimate_eutectic_values(params)

        lines.append(f"D*DT/DX^2={alpha:10.3E}   Grad={grad_k_cm:6.2f} K/cm")
        lines.append(f"Vel={vel_um_s:7.1f} um/s  ({vel_mm_h:6.1f} mm/h)")

        if kgt is not None:
            lines.append(
                f"RtipKGT={kgt['rtip_m']*1.0e6:7.3f} um  "
                f"Rtip/DX={kgt['rhotip']:6.2f}"
            )
            lines.append(
                f"KGTdelT2D={kgt['dt_kgt']:7.3f} K "
                f"Peclet={kgt['peclet']:9.3E}"
            )
#            lines.append(
#                f"OmegaKGT2D={kgt['omega2d']:7.4f} "
#                f"BL={kgt['boundary_cells']:6.2f} cells"
#            )
#            lines.append(
#                f"dTsol={kgt['dt_solute']:6.3f} "
#                f"dTcurv={kgt['dt_curve']:6.3f} "
#                f"dTtherm={kgt['dt_thermal']:6.3f}"
#            )
            if kgt['rhotip'] < 0.9:
                lines.append("(!!) RtipKGT<0.9*DX: limited resolution")
        else:
            lines.append("KGT estimates unavailable")

        if fs_scheil is not None:
            lines.append(f"Teut={teut:7.2f} C   Scheil fs @Teut={fs_scheil:6.3f}")
        else:
            lines.append(f"Teut={teut:7.2f} C")

    except Exception:
        lines.append("(checks unavailable)")

    return "\n".join(lines)


def print_dt_hint(params):
    try:
        dcoeff = float(params.get("DCOEFF", PARAM_DEFAULT["DCOEFF"]))
        dx = float(params.get("DX", PARAM_DEFAULT["DX"]))
        if dcoeff > 0.0 and dx > 0.0:
            dt_alpha = 0.20 * dx * dx / dcoeff
            print(f"DT for ALPHA=0.20: {dt_alpha:.6E} s")
    except Exception:
        pass


# =====================================================================
# MENU / EDITING
# =====================================================================

def print_menu(params, target_file):
    print("\n" + "=" * 80)
    print("CA PARAMETER EDITOR  |  Target: " + target_file)
    print("=" * 80)
    print("Parameters:")

    for i, name in enumerate(PARAM_NAMES, 1):
        val = params.get(name, PARAM_DEFAULT[name])
        if is_theta_param(name):
            shown = theta_display(val)
            line = f"{i:2d}) {name:<8} = {shown}"
        elif name.upper() in ("GRAD", "VEL"):
            shown = grad_vel_display(name, val, params)
            line = f"{i:2d}) {name:<8} = {shown:<14}"            
        else:
            line = f"{i:2d}) {name:<8} = {str(val):<14}"

        desc = PARAM_DESC[name]
        if len(desc) > 46:
            desc = desc[:46] + "..."
        print(f"{line:<46}# {desc}")

    print("-" * 80)
    print_command_and_checks(params)


def print_command_and_checks(params):
    cmd_lines = [
        "Commands:",
        "   a)  Add/edit arbitrary parameter",
        "   s)  Save and quit",
        "   sa) Save-As and continue",
        "   r)  Reload from file",
        "   q)  Quit without saving",
    ]

    check_lines = derived_checks_text(params).split("\n")
    width = max(len(line) for line in cmd_lines)

    for i in range(max(len(cmd_lines), len(check_lines))):
        cmd = cmd_lines[i] if i < len(cmd_lines) else ""
        chk = check_lines[i] if i < len(check_lines) else ""
        print(f"{cmd:<{width+2}}   {chk}")


def edit_value(name, params):
    meta_default = PARAM_DEFAULT[name] if name in PARAM_DEFAULT else ""
    default = params.get(name, meta_default)

    print(f"\nEditing: {name}")

    if is_theta_param(name):
        try:
            default_deg = math.degrees(float(default))
            print(f"Current: {default_deg:.6g} deg  ({float(default):.8g} rad)")
            new_val = input(f"New value in degrees [{default_deg:.6g}]: ").strip()
            if new_val == "":
                return default
            return math.radians(float(new_val))
        except Exception:
            pass

    print(f"Current: {default}")

    if name.upper() == "DT":
        print_dt_hint(params)

    new_val = input(f"New value [{default}]: ").strip()
    if new_val == "":
        return default

    try:
        return cast_value(name, new_val)
    except ValueError:
        return new_val


def save_as_prompt(target_file, force_new_name):
    if force_new_name:
        newfile = input("Save as: ").strip()
    else:
        newfile = input(f"Save as [{target_file}]: ").strip()

    if newfile:
        return newfile
    return target_file


# =====================================================================
# MAIN PROGRAM
# =====================================================================

def main():
    default_file = sys.argv[1] if len(sys.argv) > 1 else "xparam.txt"
    filename = input(f"Load parameter file [{default_file}]: ").strip()
    if filename == "":
        filename = default_file

    params = make_default_params()
    file_params = load_param_file(filename)
    params.update(file_params)

    if file_params:
        print(f"Loaded {len(file_params)} parameters from {filename}.")
    else:
        print(f"File not found/empty -> using defaults (save to {filename}).")

    target_file = filename

    while True:
        print_menu(params, target_file)
        choice = input("Choice: ").strip().lower()

        if choice == "q":
            print("Exited without saving.")
            break

        elif choice in ("s", "sa"):
            target_file = save_as_prompt(target_file, choice == "sa")
            save_param_file(target_file, params)
            print(f"Saved: {target_file}")
            if choice == "s":
                break

        elif choice == "r":
            file_params = load_param_file(filename)
            params.update(file_params)
            print(f"Reloaded from {filename}.")

        elif choice == "a":
            name = input("Parameter name: ").strip().upper()
            if name:
                params[name] = edit_value(name, params)

        elif choice.isdigit():
            idx = int(choice) - 1
            if 0 <= idx < len(PARAM_NAMES):
                name = PARAM_NAMES[idx]
                params[name] = edit_value(name, params)
            else:
                print("Invalid number.")

        else:
            print("Unknown command.")

    print("Done!")


if __name__ == "__main__":
    main()