make3-grid.py

#!/usr/bin/env python3
"""
make3-grid.py
Starting-grid builder for the CA dendritic solidification workflow.
Generated with AI assistance from an older F77 helper routine.

Purpose
-------
This script creates the initial `xgrid.txt` file read by GETGRID in the
Fortran CA code.  It combines:

    1) a named CA parameter file, such as xparam.txt
    2) the first row of a thermal-history CSV, such as xthermal.csv

The result is a full DOTABLE/GETGRID-style grid containing the starting
liquid field, initial solid seed geometry, local temperature field, and
orientation information.

Grid column meanings
--------------------
The output columns are:

    X          1-based CA grid x-index
    Y          1-based CA grid y-index
    #liqCells  liquid subcell count in the CA cell
    Lcomp      liquid composition, usually initialized to C0
    Scomp      solid composition, initialized to KALPHA*C0 for seed cells
    Phase      phase flag: 0.0 = liquid, 1.0 = primary solid
    Temp       starting cell temperature in degC
    Orient     crystal/growth orientation in radians

The #liqCells column is the most important state flag at startup:
    #liqCells = NUNITS  -> fully liquid
    #liqCells = 0       -> fully solid seed

Temperature initialization
--------------------------
The script reads the first thermal CSV row and constructs the starting
macro-temperature field as:

    T(i,j) = Tref + Gx*(x-xref) + Gy*(y-yref)

where x and y are physical positions in meters.  The reference point is
controlled by xref_frac and yref_frac, usually 0.5 and 0.5 for the grid
center.

For uniform-cooling CSV files, Gx=Gy=0, so the initial temperature is
spatially uniform.

Programming approach
--------------------
This is a lean, self-contained console helper. Startup geometry is
kept in plain helper routines and command options are parsed with a
small local parser rather than a package-style interface.

Seed options
------------
The script supports several startup geometries:

    bottom       old SETUP-like bottom seed rows and bumps
    center       single  seed
    multicenter  random center seeded grains
    none         all-liquid field

# ----------------------------------------------------------------------
# Initial seed / grain-geometry options
# ----------------------------------------------------------------------
# This menu controls the starting solid/grain geometry written to xgrid.txt.
# It defines where initial solid/interface cells are placed and how their
# crystallographic growth angles are assigned.
#
# Bottom-seed options are mainly for directional-solidification/chill-growth
# cases, where grains start at the lower boundary and grow upward.
# Center-seed options are mainly for free-growth, equiaxed-growth, or
# symmetry/debugging tests, where grains start away from the boundary.
#
# Options:
#   1) single crystal bottom seed       uses theta1
#      One bottom grain with one growth orientation.
#
#   2) bicrystal bottom seed            uses theta1 and theta2
#      Two bottom grains with prescribed left/right orientations.
#
#   3) multigrain sequential bottom     spans theta1 to theta2
#      Multiple bottom grains placed side-by-side across the domain;
#      orientations are assigned in sequence from theta1 to theta2.
#
#   4) multigrain clustered bottom      clustered between theta1 and theta2
#      Multiple bottom grains placed side-by-side across the domain;
#      orientations are assigned from three representative angle groups
#      between theta1 and theta2.
#
#   5) multigrain fixed-angle bottom    fixed angle list in this script
#      Multiple bottom grains using a fixed repeating set of prescribed
#      angles, useful for repeatable orientation-competition tests:
#      62.4, 108.9, 79.3, 134.2, 48.7,
#      91.5, 122.1, 70.8, 100.4, 54.6 degrees.
#
#   6) multigrain random bottom         random between theta1 and theta2
#      Multiple bottom grains with randomized growth orientations.
#
#   7) single center seed               uses theta1
#      One seed placed near the center of the domain.
#
#   8) random centered grains           random from 0 to pi radians
#      Multiple grains placed near the center region with randomized
#      positions and growth orientations.
#
#   9) no seed
#      Initially fully liquid grid; growth requires later nucleation or
#      solid-cell creation elsewhere in the CA model.
# ----------------------------------------------------------------------
"""

import sys
import csv
import math
import random
from pathlib import Path


# CA startup-state conventions used by the grid writer.
NUNITS_DEFAULT = 1000.0
NUNITS = NUNITS_DEFAULT
PH_LIQUID = 0.0
PH_PRIMARY = 1.0


def read_params(path):
    """Read named-value CA parameters, allowing trailing comments."""
    params = {}
    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:
                continue
            key = parts[0].upper()
            val = parts[1].strip()
            try:
                params[key] = float(val)
            except ValueError:
                params[key] = val
    return params


def get_float(params, key, default=None):
    k = key.upper()
    if k in params:
        return float(params[k])
    if default is not None:
        return float(default)
    raise KeyError(f"Required parameter {key!r} not found in parameter file.")


def get_int(params, key, default=None):
    k = key.upper()
    if k in params:
        return int(round(float(params[k])))
    if default is not None:
        return int(default)
    raise KeyError(f"Required parameter {key!r} not found in parameter file.")


def norm_name(name):
    """Normalize CSV column names for forgiving header matching."""
    return name.strip().lower().replace(" ", "").replace("_", "")


def find_col(fieldnames, choices):
    lookup = {norm_name(c): c for c in fieldnames}
    for c in choices:
        got = lookup.get(norm_name(c))
        if got is not None:
            return got
    return None


def read_thermal_first_row(path):
    """Return time, Tref_C, Gx_K_per_m, Gy_K_per_m from first CSV data row."""
    with path.open("r", newline="", encoding="utf-8", errors="replace") as f:
        reader = csv.DictReader(f)
        if reader.fieldnames is None:
            raise ValueError(f"No CSV header found in {path}")
        names = reader.fieldnames
        c_time = find_col(names, ["time_s", "time(s)", "time"])
        c_tref = find_col(names, ["Tref_C", "Tref(C)", "Tref"])
        c_gx = find_col(names, ["Gx_K_per_m", "Gx", "Gx(K/m)"])
        c_gy = find_col(names, ["Gy_K_per_m", "Gy", "Gy(K/m)"])
        c_g = find_col(names, ["G_K_per_m", "G", "grad", "gradient"])

        row = next(reader, None)
        if row is None:
            raise ValueError(f"No data rows found in {path}")

        if c_tref is None:
            raise ValueError("Thermal CSV needs Tref_C/Tref(C)/Tref column.")

        time = float(row[c_time]) if c_time is not None and row[c_time] != "" else 0.0
        tref = float(row[c_tref])

        if c_gx is not None or c_gy is not None:
            gx = float(row[c_gx]) if c_gx is not None and row[c_gx] != "" else 0.0
            gy = float(row[c_gy]) if c_gy is not None and row[c_gy] != "" else 0.0
        elif c_g is not None:
            gx = 0.0
            gy = float(row[c_g])
        else:
            gx = 0.0
            gy = 0.0

        return time, tref, gx, gy


# -------------------------------------------------------------------
# Bottom-seed orientation helpers
#
# These reproduce the old IVAR-style SETUP behavior in Python so the
# Fortran CA code no longer has to own startup geometry generation.
# -------------------------------------------------------------------
def theta_pick(
    i,
    imid,
    mode,
    theta1,
    theta2,
    thmin,
    thmax,
    w,
    nstr,
    random_angles = None,
):
    # Single crystal and bicrystal are the two simplest IVAR modes.
    if mode == -1:
        return theta1
    if mode == 0:
        return theta1 if i <= imid else theta2

    # For polycrystal-style bottom seeds, assign orientation by stripe.
    s = 1 + (i - 1) // max(1, w)
    s = max(1, min(nstr, s))

    if random_angles is not None:
        return random_angles[s - 1]

    if mode == 1:
        if nstr <= 1:
            return thmin
        return thmin + float(s - 1) * (thmax - thmin) / float(nstr - 1)

    if mode == 2:
        r = s % 3
        if r == 1:
            return thmin + 0.15 * (thmax - thmin)
        if r == 2:
            return thmin + 0.50 * (thmax - thmin)
        return thmin + 0.85 * (thmax - thmin)

    if mode == 3:
        deg = [62.4, 108.9, 79.3, 134.2, 48.7,
               91.5, 122.1, 70.8, 100.4, 54.6]
        return math.radians(deg[(s - 1) % len(deg)])

    return theta1


def centered_span(center, width, lo, hi):
    """Return a clipped centered integer span with at least one point."""
    w = max(1, int(width))
    start = int(center) - w // 2
    end = start + w - 1
    if start < lo:
        end += lo - start
        start = lo
    if end > hi:
        start -= end - hi
        end = hi
    start = max(lo, start)
    end = min(hi, end)
    return range(start, end + 1)


def add_vertical_bump(
    solid,
    xcenter,
    j0,
    j1,
    xlen,
    ylen,
    width,
    theta,
):
    """Add a finite-width vertical startup bump to the solid seed map."""
    for k in range(max(1, j0), min(ylen, j1) + 1):
        for ii in centered_span(xcenter, width, 1, xlen):
            solid[(ii, k)] = theta


def seed_cells_for_bottom(
    xlen,
    ylen,
    ivar,
    theta1,
    theta2,
    nstr,
    random_poly,
    seed_random,
    bicrystal_left_bumps = 1,
    bicrystal_right_bumps = 1,
    seed_thickness = 5,
    seed_length = 20,
):
    """Return solid seed cells for old SETUP-like bottom seeding."""
    # Seed features need finite thickness so early CA growth has enough
    # active interface area.  A value of 5 gives >=4-cell-thick starts
    # while staying close to the old narrow-bump geometry.
    thick = max(4, int(seed_thickness))
    box2 = max(4, thick)
    imid = xlen // 2
    w = max(1, xlen // max(1, nstr))
    thmin = min(theta1, theta2)
    thmax = max(theta1, theta2)

    # IVAR controls the old bottom-seed style.  The protrusions are now
    # deliberately longer than the bottom base so they remain active and
    # visually distinct during low-undercooling startup.
    if ivar < 0:
        mode = -1
    elif ivar == 0:
        mode = 0
    else:
        mode = ivar

    bump_len = max(thick, int(seed_length))
    bump_top = min(ylen, box2 + bump_len)

    random_angles = None
    if random_poly:
        rng = random.Random(seed_random)
        lo, hi = (thmin, thmax) if thmax > thmin else (0.0, math.pi)
        random_angles = [rng.uniform(lo, hi) for _ in range(nstr)]

    solid = {}

    def pick(i):
        return theta_pick(i, imid, mode, theta1, theta2, thmin, thmax,
                          w, nstr, random_angles)

    # The bottom rows provide the initial solid base.
    for j in range(1, min(ylen, box2) + 1):
        for i in range(1, xlen + 1):
            solid[(i, j)] = pick(i)

    # One bump at center of each seed stripe/grain.
    # For bicrystal starts, the left/right bump counts are handled below.
    if mode != 0:
        for sctr in range(1, nstr + 1):
            ic = (sctr - 1) * w + w // 2
            ic = max(1, min(xlen, ic))
            add_vertical_bump(solid, ic, box2 + 1, bump_top, xlen, ylen,
                              thick, pick(ic))

    # Extra bumps for single crystal and configurable bicrystal starts.
    if mode == -1:
        isp = max(1, min(xlen, xlen // 2))
        add_vertical_bump(solid, isp, box2 + 1, bump_top, xlen, ylen,
                          thick, pick(isp))
    elif mode == 0:
        nl = max(1, int(bicrystal_left_bumps))
        nr = max(1, int(bicrystal_right_bumps))

        def spaced_positions(x0, x1, n) -> list:
            x0 = max(1, min(xlen, x0))
            x1 = max(1, min(xlen, x1))
            if x1 < x0:
                x0, x1 = x1, x0
            if n <= 1:
                return [(x0 + x1) // 2]
            width = float(x1 - x0 + 1)
            return [
                max(x0, min(x1, int(round(x0 + (m + 0.5) * width / n))))
                for m in range(n)
            ]

        left_bumps = spaced_positions(1, imid, nl)
        right_bumps = spaced_positions(imid + 1, xlen, nr)

        for isp in left_bumps:
            add_vertical_bump(solid, isp, box2 + 1, bump_top, xlen, ylen,
                              thick, theta1)
        for isp in right_bumps:
            add_vertical_bump(solid, isp, box2 + 1, bump_top, xlen, ylen,
                              thick, theta2)

    return solid


def seed_cells_center(
    xlen,
    ylen,
    theta,
    radius,
    seed_thickness = 5,
):
    """Return a single circular seed at the grid center."""
    cx = max(1, xlen // 2)
    cy = max(1, ylen // 2)
    r = max(int(radius), int(math.ceil(0.5 * max(4, seed_thickness))))
    r2 = r ** 2
    solid = {}
    for j in range(max(1, cy - r), min(ylen, cy + r) + 1):
        for i in range(max(1, cx - r), min(xlen, cx + r) + 1):
            if (i - cx) ** 2 + (j - cy) ** 2 <= r2:
                solid[(i, j)] = theta
    return solid


def seed_cells_multicenter(
    xlen,
    ylen,
    ngrains,
    radius,
    seed_random,
    theta_min = 0.0,
    theta_max = math.pi,
    region_frac = 0.60,
    max_attempts = 10000,
    seed_thickness = 5,
):
    """Random finite nuclei in a centered rectangular region.

    Each grain is a circular solid seed with its own random orientation.
    Centers are kept at least 2*radius+2 cells apart where possible so
    the starting grains do not overlap.
    """
    rng = random.Random(seed_random)
    r = max(int(radius), int(math.ceil(0.5 * max(4, seed_thickness))))
    n = max(1, int(ngrains))
    frac = max(0.05, min(1.0, float(region_frac)))

    cx = 0.5 * float(xlen + 1)
    cy = 0.5 * float(ylen + 1)
    half_w = 0.5 * frac * float(xlen)
    half_h = 0.5 * frac * float(ylen)

    xmin = max(1 + r, int(round(cx - half_w)))
    xmax = min(xlen - r, int(round(cx + half_w)))
    ymin = max(1 + r, int(round(cy - half_h)))
    ymax = min(ylen - r, int(round(cy + half_h)))

    if xmin > xmax or ymin > ymax:
        raise ValueError("Centered grain region is too small for the requested radius.")

    centers = []
    min_sep2 = float((2*r + 2) ** 2)

    attempts = 0
    while len(centers) < n and attempts < max_attempts:
        attempts += 1
        ix = rng.randint(xmin, xmax)
        jy = rng.randint(ymin, ymax)
        ok = True
        for ox, oy, _ in centers:
            if float((ix - ox)**2 + (jy - oy)**2) < min_sep2:
                ok = False
                break
        if not ok:
            continue
        theta = rng.uniform(0.0, math.pi)
        centers.append((ix, jy, theta))

    # If the requested density is high, relax the spacing and place the rest.
    while len(centers) < n:
        ix = rng.randint(xmin, xmax)
        jy = rng.randint(ymin, ymax)
        theta = rng.uniform(0.0, math.pi)
        centers.append((ix, jy, theta))


    solid = {}
    r2 = r*r
    for cx_i, cy_i, theta in centers:
        for j in range(max(1, cy_i - r), min(ylen, cy_i + r) + 1):
            for i in range(max(1, cx_i - r), min(xlen, cx_i + r) + 1):
                if (i - cx_i)**2 + (j - cy_i)**2 <= r2:
                    solid[(i, j)] = theta
    return solid


def grain_orientation_summary(solid_cells, max_items = 12):
    """Return compact orientation summary for the generated seed cells."""
    if not solid_cells:
        return []

    counts = {}
    for theta in solid_cells.values():
        # Round only for grouping/printing; the grid still keeps full values.
        key = round(float(theta), 10)
        counts[key] = counts.get(key, 0) + 1

    items = sorted(counts.items(), key=lambda x: (-x[1], x[0]))
    out = []
    for theta, count in items[:max_items]:
        out.append((theta, math.degrees(theta), count))
    if len(items) > max_items:
        out.append((None, None, len(items) - max_items))
    return out


def print_seed_orientation_summary(seed, ivar, seed_name, solid_cells,
                                   theta1, theta2, args):
    """Print the actual seed orientations used in xgrid.txt."""
    print(f"Seed            : {seed_name}, solid cells={len(solid_cells)}")

    if seed == "none" or not solid_cells:
        print("Seed geometry   : no initial solid seed")
        return

    if seed == "bottom":
        if ivar < 0:
            desc = "single crystal bottom seed"
        elif ivar == 0:
            desc = "bicrystal bottom seed"
        elif args.random_poly:
            desc = "multigrain random bottom seed"
        elif ivar == 1:
            desc = "multigrain sequential bottom seed"
        elif ivar == 2:
            desc = "multigrain clustered bottom seed"
        elif ivar == 3:
            desc = "multigrain fixed-angle bottom seed"
        else:
            desc = f"bottom seed, IVAR={ivar}"
        print(f"Seed geometry   : {desc}")
        print(f"Seed stripes    : {args.nstr}")
    elif seed == "center":
        print("Seed geometry   : single center seed")
    elif seed == "multicenter":
        print("Seed geometry   : random centered grains")

    print("Seed angles     : radians / degrees / solid cells")
    for theta, deg, count in grain_orientation_summary(solid_cells):
        if theta is None:
            print(f"                  ... {count} more unique angles")
        else:
            print(f"                  {theta:10.6f} rad  {deg:8.3f} deg  n={count}")

    if seed == "bottom" and ivar == 0:
        print(f"Theta1/Theta2   : {theta1:.6f} rad / {theta2:.6f} rad")
    elif seed in ("bottom", "center") and ivar < 0:
        print(f"Theta used      : {theta1:.6f} rad  {math.degrees(theta1):.3f} deg")


# -------------------------------------------------------------------
# Grid writer
#
# This is the final CA handoff: every cell receives its liquid/solid
# state, composition, phase flag, temperature, and orientation.
# -------------------------------------------------------------------
def write_grid(
    out_path,
    xlen,
    ylen,
    dx_m,
    c0,
    kalpha,
    nunit,
    tref,
    gx,
    gy,
    xref_frac,
    yref_frac,
    liquid_theta,
    solid_cells,
):
    # Local-equilibrium seed composition.
    cs = kalpha * c0
    xref = xref_frac * float(xlen - 1) * dx_m
    yref = yref_frac * float(ylen - 1) * dx_m

    with out_path.open("w", encoding="utf-8") as f:
        f.write(" X Y #liqCells Lcomp Scomp Phase Temp Orient\n")
        for j in range(1, ylen + 1):
            y = float(j - 1) * dx_m
            for i in range(1, xlen + 1):
                x = float(i - 1) * dx_m
                temp = tref + gx * (x - xref) + gy * (y - yref)

                # Default state: fully liquid cell.
                theta = liquid_theta
                lcells = nunit
                phase = PH_LIQUID
                scomp = 0.0

                # Seed override: fully solid primary phase cell.
                if (i, j) in solid_cells:
                    theta = solid_cells[(i, j)]
                    lcells = 0.0
                    phase = PH_PRIMARY
                    scomp = cs

                f.write(
                    f"{i:4d}{j:4d}"
                    f"{lcells:10.3f}{c0:10.3f}{scomp:10.3f}"
                    f"{phase:10.2f}{temp:10.2f}{theta:10.4f}\n"
                )


def yes_no(prompt, default=True):
    tag = "Y/n" if default else "y/N"
    ans = input(f"{prompt} [{tag}]: ").strip().lower()
    if ans == "":
        return default
    return ans in ("y", "yes")


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


def ask_int(prompt, default):
    ans = input(f"{prompt} [{default}]: ").strip()
    return int(ans) if ans else int(default)


def ask_float(prompt, default):
    ans = input(f"{prompt} [{default}]: ").strip()
    return float(ans) if ans else float(default)


# =====================================================================
# SIMPLE ARGUMENT HOLDER
# =====================================================================
class Args:
    def __init__(self, **kw):
        self.__dict__.update(kw)


def default_args():
    return Args(
        param_file=None,
        thermal_csv=None,
        output=None,
        seed="single",
        ivar=None,
        nstr=10,
        random_poly=False,
        random_seed=1,
        center_radius=4,
        ngrains=8,
        center_region_frac=0.60,
        seed_thickness=5,
        seed_length=20,
        bicrystal_left_bumps=1,
        bicrystal_right_bumps=1,
        xref_frac=0.5,
        yref_frac=0.5,
        nunits=NUNITS_DEFAULT,
        quiet=False,
        interactive=False,
    )


def print_usage():
    print("Usage:")
    print("  python3 make3-grid.py")
    print("  python3 make3-grid.py xparam.txt [xthermal.csv] [options]")
    print()
    print("Options:")
    print("  -o, --output FILE")
    print("  --seed single|bicrystal|sequential|clustered|fixed|random|center|multicenter|none")
    print("  --ivar N")
    print("  --nstr N")
    print("  --random-poly")
    print("  --random-seed N")
    print("  --center-radius N")
    print("  --ngrains N")
    print("  --center-region-frac X")
    print("  --seed-thickness N")
    print("  --seed-length N")
    print("  --bicrystal-left-bumps N")
    print("  --bicrystal-right-bumps N")
    print("  --xref-frac X")
    print("  --yref-frac Y")
    print("  --quiet")


# =====================================================================
# INTERACTIVE FRONT END
# =====================================================================
def interactive_args(default_param="xparam.txt"):
    print("\n" + "="*72)
    print("CA STARTING GRID CREATOR")
    print("="*72)

    param_file = ask_str("Parameter file", default_param)
    params = read_params(Path(param_file))

    thermal_default = str(params.get("THFIL", "thermal.csv"))
    thermal_csv = ask_str("Thermal CSV", thermal_default)

    grid_default = str(params.get("GRIDFIL", "xgrid.txt"))
    output = ask_str("Output grid file", grid_default)

    #   1) single crystal bottom seed       uses theta1
    #      One bottom grain with one growth orientation.
    #
    #   2) bicrystal bottom seed            uses theta1 and theta2
    #      Two bottom grains with prescribed left/right orientations.
    #
    #   3) multigrain sequential bottom     spans theta1 to theta2
    #      Multiple bottom grains placed side-by-side across the domain;
    #      orientations are assigned in sequence from theta1 to theta2.
    #
    #   4) multigrain clustered bottom      clustered between theta1 and theta2
    #      Multiple bottom grains placed side-by-side across the domain;
    #      orientations are assigned from three representative angle groups
    #      between theta1 and theta2.
    #
    #   5) multigrain fixed-angle bottom    fixed angle list in this script
    #      Multiple bottom grains using a fixed repeating set of prescribed
    #      angles, useful for repeatable orientation-competition tests:
    #      62.4, 108.9, 79.3, 134.2, 48.7,
    #      91.5, 122.1, 70.8, 100.4, 54.6 degrees.
    #
    #   6) multigrain random bottom         random between theta1 and theta2
    #      Multiple bottom grains with randomized growth orientations.
    #
    #   7) single center seed               uses theta1
    #      One seed placed near the center of the domain.
    #
    #   8) random centered grains           random from 0 to pi radians
    #      Multiple grains placed near the center region with randomized
    #      positions and growth orientations.
    #
    #   9) no seed
    #      Initially fully liquid grid; growth requires later nucleation or
    #      solid-cell creation elsewhere in the CA model.
    # ----------------------------------------------------------------------


    print("\nSeed option:")
    print("  1) single crystal bottom seed    (theta1)")
    print("  2) bicrystal bottom seed         (theta1 & theta2)")
    print("  3) multigrain sequential bottom  (spans theta1 to theta2)")
    print("  4) multigrain clustered bottom   (within theta1 & theta2)")
    print("  5) multigrain fixed-angle bottom (fixed set of angles)")
    print("  6) multigrain random bottom      (random theta1 to theta2)")
    print("  7) center seed                   (theta1)")
    print("  8) random centered grains        (random; 0 to pi)")
    print("  9) no seed")

    seed_choice = ask_str("Choice", "1").lower()

    seed = "bottom"
    ivar = -1
    random_poly = False
    random_seed = 1
    nstr = 10
    center_radius = 4
    bicrystal_left_bumps = 1
    bicrystal_right_bumps = 1
    ngrains = 8
    center_region_frac = 0.60
    seed_thickness = 5
    seed_length = 20

    if seed_choice in ("1", "single", "s"):
        seed = "bottom"
        ivar = -1
    elif seed_choice in ("2", "bi", "bicrystal", "b"):
        seed = "bottom"
        ivar = 0
        bicrystal_left_bumps = ask_int("Number of left bicrystal seed bumps", 1)
        bicrystal_right_bumps = ask_int("Number of right bicrystal seed bumps", 1)
    elif seed_choice in ("3", "seq", "sequential"):
        seed = "bottom"
        ivar = 1
    elif seed_choice in ("4", "cluster", "clustered"):
        seed = "bottom"
        ivar = 2
    elif seed_choice in ("5", "fixed", "fixed-angle", "fixedangle"):
        seed = "bottom"
        ivar = 3
    elif seed_choice in ("6", "random", "rand"):
        seed = "bottom"
        ivar = 1
        random_poly = True
        random_seed = ask_int("Random seed", 1)
    elif seed_choice in ("7", "center", "c"):
        seed = "center"
        ivar = -1
        center_radius = ask_int("Center seed radius (cells)", 4)
    elif seed_choice in ("8", "multi", "multicenter", "random-center", "randomcenter"):
        seed = "multicenter"
        ivar = -1
        ngrains = ask_int("Number of centered grains", 8)
        center_radius = ask_int("Grain seed radius (cells)", 4)
        center_region_frac = ask_float("Centered region fraction", 0.60)
        random_seed = ask_int("Random seed", 1)
    elif seed_choice in ("9", "none", "n"):
        seed = "none"
        ivar = -1
    else:
        print("Unrecognized seed choice; using single crystal bottom seed.")
        seed = "bottom"
        ivar = -1

    if seed == "bottom" and ivar > 0:
        nstr = ask_int("Number of seed stripes/grains", 10)

    seed_thickness = ask_int("Minimum seed thickness (cells)", 5)
    if seed == "bottom":
        seed_length = ask_int("Bottom seed protrusion length (cells above base)", 20)

    xref_frac = ask_float("Thermal reference x fraction", 0.5)
    yref_frac = ask_float("Thermal reference y fraction", 0.5)
    quiet = not yes_no("Print summary after writing", True)

    return Args(
        param_file=param_file,
        thermal_csv=thermal_csv,
        output=output,
        seed=seed,
        ivar=ivar,
        nstr=nstr,
        random_poly=random_poly,
        random_seed=random_seed,
        center_radius=center_radius,
        ngrains=ngrains,
        center_region_frac=center_region_frac,
        seed_thickness=seed_thickness,
        seed_length=seed_length,
        bicrystal_left_bumps=bicrystal_left_bumps,
        bicrystal_right_bumps=bicrystal_right_bumps,
        xref_frac=xref_frac,
        yref_frac=yref_frac,
        nunits=NUNITS_DEFAULT,
        quiet=quiet,
        interactive=True,
    )


def parse_args(argv):
    if len(argv) == 1:
        return interactive_args()

    args = default_args()
    i = 1

    if argv[i] in ("-h", "--help"):
        print_usage()
        raise SystemExit(0)

    if argv[i] in ("-i", "--interactive"):
        return interactive_args()

    args.param_file = argv[i]
    i += 1

    if i < len(argv) and not argv[i].startswith("-"):
        args.thermal_csv = argv[i]
        i += 1

    while i < len(argv):
        opt = argv[i]

        if opt in ("-h", "--help"):
            print_usage()
            raise SystemExit(0)
        elif opt in ("-i", "--interactive"):
            return interactive_args(args.param_file or "xparam.txt")
        elif opt in ("-o", "--output"):
            i += 1; args.output = argv[i]
        elif opt == "--seed":
            i += 1; args.seed = argv[i].lower()
        elif opt == "--ivar":
            i += 1; args.ivar = int(argv[i])
        elif opt == "--nstr":
            i += 1; args.nstr = int(argv[i])
        elif opt == "--random-poly":
            args.random_poly = True
        elif opt == "--random-seed":
            i += 1; args.random_seed = int(argv[i])
        elif opt == "--center-radius":
            i += 1; args.center_radius = int(argv[i])
        elif opt == "--ngrains":
            i += 1; args.ngrains = int(argv[i])
        elif opt == "--center-region-frac":
            i += 1; args.center_region_frac = float(argv[i])
        elif opt == "--seed-thickness":
            i += 1; args.seed_thickness = int(argv[i])
        elif opt == "--seed-length":
            i += 1; args.seed_length = int(argv[i])
        elif opt == "--bicrystal-left-bumps":
            i += 1; args.bicrystal_left_bumps = int(argv[i])
        elif opt == "--bicrystal-right-bumps":
            i += 1; args.bicrystal_right_bumps = int(argv[i])
        elif opt == "--xref-frac":
            i += 1; args.xref_frac = float(argv[i])
        elif opt == "--yref-frac":
            i += 1; args.yref_frac = float(argv[i])
        elif opt == "--quiet":
            args.quiet = True
        else:
            raise ValueError("Unknown option: " + opt)

        i += 1

    valid = [
        "bottom", "single", "bicrystal", "sequential", "clustered",
        "fixed", "random", "center", "multicenter", "none",
    ]
    if args.seed not in valid:
        raise ValueError("Unknown seed option: " + args.seed)

    return args


# =====================================================================
# MAIN PROGRAM
# =====================================================================
def main():
    args = parse_args(sys.argv)
    params = read_params(Path(args.param_file))

    if args.thermal_csv is None:
        args.thermal_csv = str(params.get("THFIL", "thermal.csv"))

    if args.output is None:
        args.output = str(params.get("GRIDFIL", "xgrid.txt"))

    xlen = get_int(params, "XLEN", 200)
    ylen = get_int(params, "YLEN", 400)
    dx = get_float(params, "DX")
    c0 = get_float(params, "C0")
    kalpha = get_float(params, "KALPHA")
    theta1 = get_float(params, "THETA1", 0.0)
    theta2 = get_float(params, "THETA2", theta1)

    seed_name = args.seed
    random_poly = args.random_poly

    if seed_name == "single":
        seed = "bottom"
        ivar = -1
    elif seed_name == "bicrystal":
        seed = "bottom"
        ivar = 0
    elif seed_name == "sequential":
        seed = "bottom"
        ivar = 1
    elif seed_name == "clustered":
        seed = "bottom"
        ivar = 2
    elif seed_name == "fixed":
        seed = "bottom"
        ivar = 3
    elif seed_name == "random":
        seed = "bottom"
        ivar = 1
        random_poly = True
    elif seed_name == "bottom":
        seed = "bottom"
        ivar = args.ivar if args.ivar is not None else get_int(params, "IVAR", -1)
    else:
        seed = seed_name
        ivar = args.ivar if args.ivar is not None else get_int(params, "IVAR", -1)

    time, tref, gx, gy = read_thermal_first_row(Path(args.thermal_csv))

    if seed == "bottom":
        solid_cells = seed_cells_for_bottom(
            xlen, ylen, ivar, theta1, theta2, args.nstr,
            random_poly, args.random_seed,
            args.bicrystal_left_bumps, args.bicrystal_right_bumps,
            args.seed_thickness, args.seed_length,
        )
    elif seed == "center":
        solid_cells = seed_cells_center(
            xlen, ylen, theta1, args.center_radius, args.seed_thickness,
        )
    elif seed == "multicenter":
        solid_cells = seed_cells_multicenter(
            xlen, ylen, args.ngrains, args.center_radius,
            args.random_seed, 0.0, math.pi, args.center_region_frac,
            seed_thickness=args.seed_thickness,
        )
    else:
        solid_cells = {}

    write_grid(
        Path(args.output), xlen, ylen, dx, c0, kalpha, NUNITS_DEFAULT,
        tref, gx, gy, args.xref_frac, args.yref_frac,
        0.0, solid_cells,
    )

    if not args.quiet:
        print("Generated grid  :", args.output)
        print("Parameter file  :", args.param_file)
        print("Thermal CSV     :", args.thermal_csv)
        text = f"Thermal row     : time={time:.6g} s, Tref={tref:.6g} C, "
        text = text + f"Gx={gx:.6g} K/m, Gy={gy:.6g} K/m"
        print(text)
        print(f"Grid            : {xlen} x {ylen}, DX={dx:.6g} m")
        print_seed_orientation_summary(
            seed, ivar, seed_name, solid_cells, theta1, theta2, args
        )
        print(f"Seed thickness  : {max(4, int(args.seed_thickness))} cells minimum")
        if seed == "bottom":
            thick = max(4, int(args.seed_thickness))
            print(f"Seed length     : {max(thick, int(args.seed_length))} cells above base")
        if seed == "center":
            print(f"Center radius   : {args.center_radius} cells")
        if seed == "multicenter":
            print(f"Centered grains : {args.ngrains}")
            print(f"Grain radius    : {args.center_radius} cells")
            print(f"Region fraction : {args.center_region_frac}")
            print(f"Random seed     : {args.random_seed}")
        if seed == "bottom" and ivar == 0:
            text = f"Bicrystal bumps : left={args.bicrystal_left_bumps}, "
            text = text + f"right={args.bicrystal_right_bumps}"
            print(text)
        print(f"Compositions    : C0={c0:.6g}, KALPHA={kalpha:.6g}, CS(seed)={kalpha*c0:.6g}")

    return 0


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