Numerical Stability

#!/usr/bin/env python3
import argparse
import json
import math
import sys
from typing import Dict, List, Optional


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Compute CFL/Fourier numbers and suggest stable dt limits.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument("--dx", type=float, required=True, help="Grid spacing")
    parser.add_argument("--dt", type=float, required=True, help="Time step")
    parser.add_argument("--velocity", type=float, default=None, help="Advection velocity")
    parser.add_argument("--diffusivity", type=float, default=None, help="Diffusivity")
    parser.add_argument("--reaction-rate", type=float, default=None, help="Reaction rate")
    parser.add_argument("--dimensions", type=int, default=1, help="Spatial dimensions")
    parser.add_argument(
        "--scheme",
        choices=["explicit", "implicit"],
        default="explicit",
        help="Time integration scheme",
    )
    parser.add_argument("--advection-limit", type=float, default=None, help="CFL limit")
    parser.add_argument("--diffusion-limit", type=float, default=None, help="Fourier limit")
    parser.add_argument("--reaction-limit", type=float, default=None, help="Reaction limit")
    parser.add_argument("--safety", type=float, default=1.0, help="Safety factor for dt")
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def compute_cfl(
    dx: float,
    dt: float,
    velocity: Optional[float],
    diffusivity: Optional[float],
    reaction_rate: Optional[float],
    dimensions: int,
    scheme: str,
    advection_limit: Optional[float],
    diffusion_limit: Optional[float],
    reaction_limit: Optional[float],
    safety: float,
) -> Dict[str, object]:
    notes: List[str] = []
    if dx <= 0 or dt <= 0:
        raise ValueError("dx and dt must be positive")
    if dimensions <= 0:
        raise ValueError("dimensions must be positive")
    if safety <= 0:
        raise ValueError("safety must be positive")

    v = None
    if velocity is not None:
        if velocity < 0:
            notes.append("Negative velocity detected; using absolute value for stability analysis")
        v = abs(velocity)
    d = None
    if diffusivity is not None:
        if diffusivity < 0:
            notes.append("Negative diffusivity detected; may indicate spinodal decomposition; using absolute value")
        d = abs(diffusivity)
    k = None
    if reaction_rate is not None:
        if reaction_rate < 0:
            notes.append("Negative reaction rate detected; using absolute value for stability analysis")
        k = abs(reaction_rate)

    if advection_limit is None:
        advection_limit = 1.0 if scheme == "explicit" else math.inf
    if diffusion_limit is None:
        diffusion_limit = 1.0 / (2.0 * dimensions) if scheme == "explicit" else math.inf
    if reaction_limit is None:
        reaction_limit = 1.0 if scheme == "explicit" else math.inf

    cfl = None
    dt_max_adv = None
    if v is not None and v > 0:
        cfl = v * dt / dx
        if math.isfinite(advection_limit):
            dt_max_adv = advection_limit * dx / v
    elif v == 0:
        cfl = 0.0

    fo = None
    dt_max_diff = None
    if d is not None and d > 0:
        fo = d * dt / (dx ** 2)
        if math.isfinite(diffusion_limit):
            dt_max_diff = diffusion_limit * (dx ** 2) / d
    elif d == 0:
        fo = 0.0

    react = None
    dt_max_react = None
    if k is not None and k > 0:
        react = k * dt
        if math.isfinite(reaction_limit):
            dt_max_react = reaction_limit / k
    elif k == 0:
        react = 0.0

    criteria_applied: List[str] = []
    stable: Optional[bool] = True
    if scheme == "explicit":
        if cfl is not None and math.isfinite(advection_limit):
            criteria_applied.append("advection")
            stable = stable and (cfl <= advection_limit + 1e-12)
        if fo is not None and math.isfinite(diffusion_limit):
            criteria_applied.append("diffusion")
            stable = stable and (fo <= diffusion_limit + 1e-12)
        if react is not None and math.isfinite(reaction_limit):
            criteria_applied.append("reaction")
            stable = stable and (react <= reaction_limit + 1e-12)
    else:
        notes.append("Implicit scheme: stability limits are relaxed; check accuracy.")

    if not criteria_applied:
        stable = None
        notes.append("No stability criteria applied; provide velocity/diffusivity/reaction.")

    dt_candidates = [x for x in [dt_max_adv, dt_max_diff, dt_max_react] if x is not None]
    recommended_dt = min(dt_candidates) * safety if dt_candidates else None

    return {
        "inputs": {
            "dx": dx,
            "dt": dt,
            "velocity": velocity,
            "diffusivity": diffusivity,
            "reaction_rate": reaction_rate,
            "dimensions": dimensions,
            "scheme": scheme,
            "safety": safety,
        },
        "metrics": {
            "cfl": cfl,
            "fourier": fo,
            "reaction": react,
        },
        "limits": {
            "advection_limit": advection_limit,
            "diffusion_limit": diffusion_limit,
            "reaction_limit": reaction_limit,
        },
        "criteria_applied": criteria_applied,
        "recommended_dt": recommended_dt,
        "stable": stable,
        "notes": notes,
    }


def main() -> None:
    args = parse_args()
    try:
        payload = compute_cfl(
            dx=args.dx,
            dt=args.dt,
            velocity=args.velocity,
            diffusivity=args.diffusivity,
            reaction_rate=args.reaction_rate,
            dimensions=args.dimensions,
            scheme=args.scheme,
            advection_limit=args.advection_limit,
            diffusion_limit=args.diffusion_limit,
            reaction_limit=args.reaction_limit,
            safety=args.safety,
        )
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("CFL check")
    print(f"  scheme: {args.scheme}")
    cfl = payload["metrics"]["cfl"]
    if cfl is not None:
        print(f"  CFL: {cfl:.6g} (limit {payload['limits']['advection_limit']:.6g})")
    else:
        print("  CFL: n/a")
    fo = payload["metrics"]["fourier"]
    if fo is not None:
        print(f"  Fourier: {fo:.6g} (limit {payload['limits']['diffusion_limit']:.6g})")
    else:
        print("  Fourier: n/a")
    react = payload["metrics"]["reaction"]
    if react is not None:
        print(f"  Reaction: {react:.6g} (limit {payload['limits']['reaction_limit']:.6g})")
    else:
        print("  Reaction: n/a")
    stable = payload["stable"]
    stable_label = "unknown" if stable is None else str(stable)
    print(f"  stable: {stable_label}")
    recommended_dt = payload["recommended_dt"]
    if recommended_dt is not None:
        print(f"  recommended_dt: {recommended_dt:.6g}")
    for note in payload["notes"]:
        print(f"  note: {note}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
import argparse
import json
import os
import sys
from typing import Dict, Optional, Union

import numpy as np


def load_matrix(path: str, delimiter: Optional[str]) -> np.ndarray:
    _, ext = os.path.splitext(path)
    if ext == ".npy":
        return np.load(path)
    return np.loadtxt(path, delimiter=delimiter)


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Compute condition number and eigenvalue spread.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument("--matrix", required=True, help="Path to matrix file (.npy or text)")
    parser.add_argument(
        "--delimiter",
        default=None,
        help="Delimiter for text matrices (default: any whitespace)",
    )
    parser.add_argument(
        "--norm",
        default="2",
        help="Condition number norm: 2, fro, inf, -inf, 1, -1",
    )
    parser.add_argument(
        "--symmetry-tol",
        type=float,
        default=1e-8,
        help="Tolerance for symmetry check",
    )
    parser.add_argument(
        "--skip-eigs",
        action="store_true",
        help="Skip eigenvalue spread for large matrices",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def parse_norm(value: str) -> Union[float, str]:
    value = value.strip().lower()
    if value in {"2", "1", "-1"}:
        return float(value)
    if value in {"inf", "-inf", "fro"}:
        return value
    raise ValueError("norm must be one of: 2, 1, -1, inf, -inf, fro")


def compute_condition(
    matrix: np.ndarray,
    norm: Union[float, str],
    symmetry_tol: float,
    skip_eigs: bool,
) -> Dict[str, object]:
    if matrix.ndim != 2 or matrix.shape[0] != matrix.shape[1]:
        raise ValueError("matrix must be square")
    if not np.all(np.isfinite(matrix)):
        raise ValueError("matrix contains non-finite values")

    cond = float(np.linalg.cond(matrix, p=norm))
    is_symmetric = bool(np.allclose(matrix, matrix.T, atol=symmetry_tol, rtol=0.0))

    spread = None
    eig_min = None
    eig_max = None
    if not skip_eigs:
        eigvals = np.linalg.eigvals(matrix)
        abs_eigs = np.abs(eigvals)
        nonzero = abs_eigs[abs_eigs > 0]
        if nonzero.size:
            eig_min = float(np.min(nonzero))
            eig_max = float(np.max(nonzero))
            spread = float(eig_max / eig_min)
        else:
            spread = float("inf")

    status = "ok"
    note = None
    if cond > 1e10:
        status = "ill-conditioned"
        note = "Consider preconditioning or scaling."
    elif cond > 1e8:
        status = "poorly-conditioned"
        note = "Preconditioning likely needed."

    return {
        "condition_number": cond,
        "eigenvalue_spread": spread,
        "eigenvalue_min_abs": eig_min,
        "eigenvalue_max_abs": eig_max,
        "is_symmetric": is_symmetric,
        "status": status,
        "note": note,
    }


def main() -> None:
    args = parse_args()
    if not os.path.exists(args.matrix):
        print(f"Matrix not found: {args.matrix}", file=sys.stderr)
        sys.exit(2)

    try:
        norm = parse_norm(args.norm)
        matrix = load_matrix(args.matrix, args.delimiter)
        results = compute_condition(matrix, norm, args.symmetry_tol, args.skip_eigs)
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    payload = {
        "inputs": {
            "matrix": args.matrix,
            "shape": list(matrix.shape),
            "norm": args.norm,
            "symmetry_tol": args.symmetry_tol,
            "skip_eigs": args.skip_eigs,
        },
        "results": results,
    }

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Matrix conditioning")
    print(f"  shape: {matrix.shape}")
    print(f"  condition number: {results['condition_number']:.6g}")
    if results["eigenvalue_spread"] is not None:
        print(f"  eigenvalue spread: {results['eigenvalue_spread']:.6g}")
    else:
        print("  eigenvalue spread: skipped")
    print(f"  symmetric: {results['is_symmetric']}")
    print(f"  status: {results['status']}")
    if results["note"]:
        print(f"  note: {results['note']}")


if __name__ == "__main__":
    main()

numpy>=1.21

#!/usr/bin/env python3
import argparse
import json
import os
import sys
from typing import Dict, Optional

import numpy as np


def parse_eigs(raw: str) -> np.ndarray:
    parts = [p.strip() for p in raw.split(",") if p.strip()]
    if not parts:
        raise ValueError("eigs must be a comma-separated list")
    return np.array([complex(p) for p in parts], dtype=complex)


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Detect stiffness from eigenvalues or a Jacobian matrix.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    group = parser.add_mutually_exclusive_group(required=True)
    group.add_argument("--eigs", help="Comma-separated eigenvalues")
    group.add_argument("--jacobian", help="Path to Jacobian matrix (.npy or text)")
    parser.add_argument(
        "--delimiter",
        default=None,
        help="Delimiter for text Jacobians (default: any whitespace)",
    )
    parser.add_argument(
        "--threshold",
        type=float,
        default=1e3,
        help="Stiffness ratio threshold",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def load_matrix(path: str, delimiter: Optional[str]) -> np.ndarray:
    _, ext = os.path.splitext(path)
    if ext == ".npy":
        return np.load(path)
    return np.loadtxt(path, delimiter=delimiter)


def compute_stiffness(eigs: np.ndarray, threshold: float) -> Dict[str, object]:
    if threshold <= 0:
        raise ValueError("threshold must be positive")
    if eigs.size == 0:
        raise ValueError("eigs must be non-empty")
    if not np.all(np.isfinite(eigs)):
        raise ValueError("eigs contain non-finite values")

    abs_eigs = np.abs(eigs)
    nonzero = abs_eigs[abs_eigs > 0]
    ratio = float(np.max(nonzero) / np.min(nonzero)) if nonzero.size else float("inf")
    stiff = ratio >= threshold
    recommendation = "implicit (BDF/Radau)" if stiff else "explicit (RK/Adams)"

    return {
        "stiffness_ratio": ratio,
        "stiff": stiff,
        "recommendation": recommendation,
        "nonzero_count": int(nonzero.size),
        "total_count": int(eigs.size),
    }


def main() -> None:
    args = parse_args()
    try:
        if args.eigs is not None:
            eigs = parse_eigs(args.eigs)
            source = "eigs"
        else:
            if not os.path.exists(args.jacobian):
                print(f"Jacobian not found: {args.jacobian}", file=sys.stderr)
                sys.exit(2)
            jacobian = load_matrix(args.jacobian, args.delimiter)
            if jacobian.ndim != 2 or jacobian.shape[0] != jacobian.shape[1]:
                print("Jacobian must be square.", file=sys.stderr)
                sys.exit(2)
            eigs = np.linalg.eigvals(jacobian)
            source = "jacobian"
        results = compute_stiffness(eigs, args.threshold)
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    payload = {
        "inputs": {
            "source": source,
            "threshold": args.threshold,
        },
        "results": results,
    }

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Stiffness detection")
    print(f"  stiffness ratio: {results['stiffness_ratio']:.6g}")
    print(f"  stiff: {results['stiff']}")
    print(f"  recommendation: {results['recommendation']}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
import argparse
import json
import math
import sys
from typing import Dict, Optional

import numpy as np


def parse_coeffs(raw: str) -> np.ndarray:
    parts = [p.strip() for p in raw.split(",") if p.strip()]
    if not parts:
        raise ValueError("coeffs must be a comma-separated list")
    return np.array([float(p) for p in parts], dtype=float)


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Compute amplification factor for a linear update stencil.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--coeffs",
        required=True,
        help="Comma-separated stencil coefficients from negative to positive index",
    )
    parser.add_argument("--dx", type=float, default=1.0, help="Grid spacing")
    parser.add_argument(
        "--offset",
        type=int,
        default=None,
        help="Index of coefficient corresponding to j=0 (default: center index)",
    )
    parser.add_argument("--kmin", type=float, default=None, help="Minimum wavenumber")
    parser.add_argument("--kmax", type=float, default=None, help="Maximum wavenumber")
    parser.add_argument("--nk", type=int, default=256, help="Number of k samples")
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def compute_amplification(
    coeffs: np.ndarray,
    dx: float,
    nk: int,
    offset: Optional[int],
    kmin: Optional[float],
    kmax: Optional[float],
) -> Dict[str, object]:
    if dx <= 0:
        raise ValueError("dx must be positive")
    if nk <= 1:
        raise ValueError("nk must be > 1")

    n = len(coeffs)
    if n == 0:
        raise ValueError("coeffs must be non-empty")
    if not np.all(np.isfinite(coeffs)):
        raise ValueError("coeffs must be finite (no NaN or Inf values)")

    if offset is None:
        offset = n // 2
    if offset < 0 or offset >= n:
        raise ValueError("offset must be within coefficient indices")

    if kmin is None:
        kmin = -math.pi / dx
    if kmax is None:
        kmax = math.pi / dx
    if kmin >= kmax:
        raise ValueError("kmin must be < kmax")

    j = np.arange(-offset, n - offset)
    ks = np.linspace(kmin, kmax, nk)
    phase = np.exp(1j * np.outer(ks, j) * dx)
    amplification = phase @ coeffs
    amp_mag = np.abs(amplification)

    max_idx = int(np.argmax(amp_mag))
    max_amp = float(amp_mag[max_idx])
    k_at_max = float(ks[max_idx])
    stable = max_amp <= 1.0 + 1e-12
    warning = None
    if n % 2 == 0 and offset == n // 2:
        warning = "Even-length stencil: confirm offset aligns with j=0."

    return {
        "inputs": {
            "coeffs": coeffs.tolist(),
            "dx": dx,
            "kmin": kmin,
            "kmax": kmax,
            "nk": nk,
            "offset": offset,
        },
        "results": {
            "max_amplification": max_amp,
            "k_at_max": k_at_max,
            "stable": stable,
            "warning": warning,
        },
    }


def main() -> None:
    args = parse_args()
    try:
        coeffs = parse_coeffs(args.coeffs)
        payload = compute_amplification(
            coeffs=coeffs,
            dx=args.dx,
            nk=args.nk,
            offset=args.offset,
            kmin=args.kmin,
            kmax=args.kmax,
        )
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Von Neumann analysis")
    print(f"  coeffs: {payload['inputs']['coeffs']}")
    print(
        "  k range: [{:.6g}, {:.6g}] with {} samples".format(
            payload["inputs"]["kmin"],
            payload["inputs"]["kmax"],
            payload["inputs"]["nk"],
        )
    )
    print(
        "  max amplification: {:.6g} at k={:.6g}".format(
            payload["results"]["max_amplification"],
            payload["results"]["k_at_max"],
        )
    )
    print(f"  stable: {payload['results']['stable']}")
    if payload["results"]["warning"]:
        print(f"  warning: {payload['results']['warning']}")


if __name__ == "__main__":
    main()