Convergence Study

#!/usr/bin/env python3
"""Temporal convergence study via dt-refinement analysis."""
import argparse
import json
import math
import sys


def parse_args():
    parser = argparse.ArgumentParser(
        description="Analyze temporal convergence by computing observed order from timestep refinement data.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--timesteps", type=str, required=True,
        help="Comma-separated timestep sizes (e.g. 0.04,0.02,0.01)",
    )
    parser.add_argument(
        "--values", type=str, required=True,
        help="Comma-separated solution values at each timestep level",
    )
    parser.add_argument(
        "--expected-order", type=float, default=None,
        help="Expected convergence order for assessment",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def compute_dt_refinement(timesteps, values, expected_order=None):
    """Compute temporal convergence order from timestep refinement data.

    Parameters
    ----------
    timesteps : list of float
        Timestep sizes (must be positive).
    values : list of float
        Solution values corresponding to each timestep.
    expected_order : float or None
        Expected convergence order for assessment.

    Returns
    -------
    dict
        Results with observed_orders, mean_order, richardson_extrapolated_value,
        in_asymptotic_range, convergence_assessment, and notes.
    """
    if len(timesteps) != len(values):
        raise ValueError("timesteps and values must have the same length")
    if len(timesteps) < 2:
        raise ValueError("At least 2 refinement levels required")
    for t in timesteps:
        if not math.isfinite(t) or t <= 0:
            raise ValueError("All timesteps must be positive finite numbers")
    for v in values:
        if not math.isfinite(v):
            raise ValueError("All values must be finite numbers")

    # Sort by timestep descending (coarsest first)
    paired = sorted(zip(timesteps, values), key=lambda x: -x[0])
    timesteps_sorted = [p[0] for p in paired]
    values_sorted = [p[1] for p in paired]

    notes = []
    observed_orders = []

    # Compute observed orders from consecutive triplets using log-ratio
    if len(timesteps_sorted) >= 3:
        for i in range(len(timesteps_sorted) - 2):
            dt_coarse = timesteps_sorted[i]
            dt_mid = timesteps_sorted[i + 1]
            dt_fine = timesteps_sorted[i + 2]
            f_coarse = values_sorted[i]
            f_mid = values_sorted[i + 1]
            f_fine = values_sorted[i + 2]

            e_coarse = abs(f_coarse - f_mid)
            e_fine = abs(f_mid - f_fine)

            if e_fine == 0 or e_coarse == 0:
                notes.append(
                    "Zero error difference at levels %d-%d; "
                    "cannot compute order for this pair" % (i, i + 2)
                )
                continue

            r_coarse = dt_coarse / dt_mid

            if r_coarse <= 0:
                continue

            p = math.log(e_coarse / e_fine) / math.log(r_coarse)
            observed_orders.append(p)

            if p < 0:
                notes.append(
                    "Negative observed order (%.2f) at levels %d-%d: "
                    "possible divergence" % (p, i, i + 2)
                )

    # Mean order
    mean_order = None
    if observed_orders:
        mean_order = sum(observed_orders) / len(observed_orders)

    # Check for pre-asymptotic behavior
    in_asymptotic_range = True
    if len(observed_orders) >= 2:
        for i in range(len(observed_orders) - 1):
            p1 = observed_orders[i]
            p2 = observed_orders[i + 1]
            avg = (abs(p1) + abs(p2)) / 2.0
            if avg > 0 and abs(p1 - p2) / avg > 0.5:
                in_asymptotic_range = False
                notes.append("Pre-asymptotic behavior detected: observed order varies >50%% between consecutive pairs")
                break
    elif len(observed_orders) == 1:
        in_asymptotic_range = True
    else:
        in_asymptotic_range = None

    # Richardson extrapolation from finest two levels
    richardson_extrapolated_value = None
    dt_fine = timesteps_sorted[-1]
    dt_next = timesteps_sorted[-2]
    f_fine = values_sorted[-1]
    f_next = values_sorted[-2]
    r = dt_next / dt_fine

    if mean_order is not None and mean_order > 0:
        richardson_extrapolated_value = f_fine + (f_fine - f_next) / (r ** mean_order - 1)
    elif expected_order is not None and expected_order > 0:
        richardson_extrapolated_value = f_fine + (f_fine - f_next) / (r ** expected_order - 1)
        notes.append("Richardson extrapolation used expected order (no observed order available)")

    # Convergence assessment
    convergence_assessment = "unknown"
    if expected_order is not None and mean_order is not None:
        if abs(mean_order - expected_order) / expected_order <= 0.1:
            convergence_assessment = "PASS: observed order (%.2f) within 10%% of expected (%.2f)" % (
                mean_order, expected_order,
            )
        else:
            convergence_assessment = "FAIL: observed order (%.2f) differs from expected (%.2f) by >10%%" % (
                mean_order, expected_order,
            )
    elif mean_order is not None:
        convergence_assessment = "Observed order: %.2f (no expected order given for comparison)" % mean_order
    elif len(timesteps_sorted) < 3:
        convergence_assessment = "Insufficient levels to compute observed order (need >= 3)"

    return {
        "inputs": {
            "timesteps": timesteps_sorted,
            "values": values_sorted,
            "expected_order": expected_order,
        },
        "results": {
            "observed_orders": observed_orders,
            "mean_order": mean_order,
            "richardson_extrapolated_value": richardson_extrapolated_value,
            "in_asymptotic_range": in_asymptotic_range,
            "convergence_assessment": convergence_assessment,
            "notes": notes,
        },
    }


def main():
    args = parse_args()
    try:
        timesteps = [float(x) for x in args.timesteps.split(",")]
        values = [float(x) for x in args.values.split(",")]
    except ValueError:
        print("Error: timesteps and values must be comma-separated numbers", file=sys.stderr)
        sys.exit(2)

    try:
        result = compute_dt_refinement(timesteps, values, args.expected_order)
    except ValueError as exc:
        print("Error: %s" % exc, file=sys.stderr)
        sys.exit(2)

    if args.json:
        print(json.dumps(result, indent=2))
    else:
        r = result["results"]
        print("Temporal Convergence Study (dt-refinement)")
        print("  levels: %d" % len(timesteps))
        if r["observed_orders"]:
            print("  observed orders: %s" % ", ".join("%.4f" % o for o in r["observed_orders"]))
        if r["mean_order"] is not None:
            print("  mean order: %.4f" % r["mean_order"])
        if r["richardson_extrapolated_value"] is not None:
            print("  Richardson extrapolated value: %.6g" % r["richardson_extrapolated_value"])
        if r["in_asymptotic_range"] is not None:
            print("  in asymptotic range: %s" % r["in_asymptotic_range"])
        print("  assessment: %s" % r["convergence_assessment"])
        for note in r["notes"]:
            print("  note: %s" % note)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Grid Convergence Index (GCI) calculator following Roache's method."""
import argparse
import json
import math
import sys


def parse_args():
    parser = argparse.ArgumentParser(
        description="Compute Grid Convergence Index (GCI) for solution verification.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--spacings", type=str, required=True,
        help="Comma-separated grid spacings for exactly 3 levels (e.g. 0.04,0.02,0.01)",
    )
    parser.add_argument(
        "--values", type=str, required=True,
        help="Comma-separated solution values at each spacing level",
    )
    parser.add_argument(
        "--safety-factor", type=float, default=1.25,
        help="Safety factor for GCI (default 1.25 for 3+ grids)",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def compute_gci(spacings, values, safety_factor=1.25):
    """Compute Grid Convergence Index.

    Parameters
    ----------
    spacings : list of float
        Exactly 3 grid spacings (positive).
    values : list of float
        Solution values at each spacing level.
    safety_factor : float
        Safety factor (default 1.25).

    Returns
    -------
    dict
        Results with observed_order, gci_fine, gci_coarse, asymptotic_ratio,
        in_asymptotic_range, refinement_ratio_21, refinement_ratio_32,
        and extrapolated_value.
    """
    if len(spacings) != 3 or len(values) != 3:
        raise ValueError("Exactly 3 refinement levels required for GCI calculation")
    if safety_factor <= 0:
        raise ValueError("safety_factor must be positive")
    for s in spacings:
        if not math.isfinite(s) or s <= 0:
            raise ValueError("All spacings must be positive finite numbers")
    for v in values:
        if not math.isfinite(v):
            raise ValueError("All values must be finite numbers")

    # Sort by spacing ascending: h1 (fine) < h2 (medium) < h3 (coarse)
    paired = sorted(zip(spacings, values), key=lambda x: x[0])
    h1, f1 = paired[0]  # fine
    h2, f2 = paired[1]  # medium
    h3, f3 = paired[2]  # coarse

    r21 = h2 / h1
    r32 = h3 / h2

    e21 = f2 - f1
    e32 = f3 - f2

    # Check for oscillatory convergence
    if e21 != 0 and e32 != 0 and (e32 * e21) < 0:
        raise ValueError("Oscillatory convergence detected")

    # Check for identical solutions
    if e21 == 0 and e32 == 0:
        return {
            "inputs": {
                "spacings": [h1, h2, h3],
                "values": [f1, f2, f3],
                "safety_factor": safety_factor,
            },
            "results": {
                "observed_order": None,
                "gci_fine": 0.0,
                "gci_coarse": 0.0,
                "asymptotic_ratio": None,
                "in_asymptotic_range": None,
                "refinement_ratio_21": r21,
                "refinement_ratio_32": r32,
                "extrapolated_value": f1,
            },
        }

    # Observed order
    if e21 == 0:
        raise ValueError("Fine and medium solutions are identical but coarse differs; cannot compute order")
    if e32 == 0:
        raise ValueError("Medium and coarse solutions are identical but fine differs; cannot compute order")

    p = abs(math.log(abs(e32 / e21))) / math.log(r21)

    # GCI fine and coarse
    gci_fine = safety_factor * abs(e21 / f1) / (r21 ** p - 1) if f1 != 0 else float("inf")
    gci_coarse = safety_factor * abs(e32 / f2) / (r32 ** p - 1) if f2 != 0 else float("inf")

    # Asymptotic ratio
    if gci_fine != 0 and math.isfinite(gci_fine) and math.isfinite(gci_coarse):
        asymptotic_ratio = gci_coarse / (r21 ** p * gci_fine)
    else:
        asymptotic_ratio = None

    in_asymptotic_range = None
    if asymptotic_ratio is not None:
        in_asymptotic_range = abs(asymptotic_ratio - 1.0) < 0.1

    # Richardson extrapolated value
    extrapolated_value = f1 + (f1 - f2) / (r21 ** p - 1)

    return {
        "inputs": {
            "spacings": [h1, h2, h3],
            "values": [f1, f2, f3],
            "safety_factor": safety_factor,
        },
        "results": {
            "observed_order": p,
            "gci_fine": gci_fine,
            "gci_coarse": gci_coarse,
            "asymptotic_ratio": asymptotic_ratio,
            "in_asymptotic_range": in_asymptotic_range,
            "refinement_ratio_21": r21,
            "refinement_ratio_32": r32,
            "extrapolated_value": extrapolated_value,
        },
    }


def main():
    args = parse_args()
    try:
        spacings = [float(x) for x in args.spacings.split(",")]
        values = [float(x) for x in args.values.split(",")]
    except ValueError:
        print("Error: spacings and values must be comma-separated numbers", file=sys.stderr)
        sys.exit(2)

    try:
        result = compute_gci(spacings, values, args.safety_factor)
    except ValueError as exc:
        print("Error: %s" % exc, file=sys.stderr)
        sys.exit(2)

    if args.json:
        print(json.dumps(result, indent=2))
    else:
        r = result["results"]
        print("Grid Convergence Index (GCI)")
        print("  refinement ratio r21: %.4f" % r["refinement_ratio_21"])
        print("  refinement ratio r32: %.4f" % r["refinement_ratio_32"])
        if r["observed_order"] is not None:
            print("  observed order: %.4f" % r["observed_order"])
        print("  GCI fine: %.6g" % r["gci_fine"])
        print("  GCI coarse: %.6g" % r["gci_coarse"])
        if r["asymptotic_ratio"] is not None:
            print("  asymptotic ratio: %.4f" % r["asymptotic_ratio"])
        if r["in_asymptotic_range"] is not None:
            print("  in asymptotic range: %s" % r["in_asymptotic_range"])
        print("  extrapolated value: %.6g" % r["extrapolated_value"])


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Spatial convergence study via h-refinement analysis."""
import argparse
import json
import math
import sys


def parse_args():
    parser = argparse.ArgumentParser(
        description="Analyze spatial convergence by computing observed order from grid refinement data.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--spacings", type=str, required=True,
        help="Comma-separated grid spacings (e.g. 0.4,0.2,0.1)",
    )
    parser.add_argument(
        "--values", type=str, required=True,
        help="Comma-separated solution values at each spacing level",
    )
    parser.add_argument(
        "--expected-order", type=float, default=None,
        help="Expected convergence order for assessment",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def compute_h_refinement(spacings, values, expected_order=None):
    """Compute spatial convergence order from grid refinement data.

    Parameters
    ----------
    spacings : list of float
        Grid spacings (must be positive).
    values : list of float
        Solution values corresponding to each spacing.
    expected_order : float or None
        Expected convergence order for assessment.

    Returns
    -------
    dict
        Results with observed_orders, mean_order, richardson_extrapolated_value,
        in_asymptotic_range, convergence_assessment, and notes.
    """
    if len(spacings) != len(values):
        raise ValueError("spacings and values must have the same length")
    if len(spacings) < 2:
        raise ValueError("At least 2 refinement levels required")
    for s in spacings:
        if not math.isfinite(s) or s <= 0:
            raise ValueError("All spacings must be positive finite numbers")
    for v in values:
        if not math.isfinite(v):
            raise ValueError("All values must be finite numbers")

    # Sort by spacing descending (coarsest first)
    paired = sorted(zip(spacings, values), key=lambda x: -x[0])
    spacings_sorted = [p[0] for p in paired]
    values_sorted = [p[1] for p in paired]

    notes = []
    observed_orders = []

    # Compute observed orders from consecutive triplets using log-ratio
    if len(spacings_sorted) >= 3:
        for i in range(len(spacings_sorted) - 2):
            h_coarse = spacings_sorted[i]
            h_mid = spacings_sorted[i + 1]
            h_fine = spacings_sorted[i + 2]
            f_coarse = values_sorted[i]
            f_mid = values_sorted[i + 1]
            f_fine = values_sorted[i + 2]

            e_coarse = abs(f_coarse - f_mid)
            e_fine = abs(f_mid - f_fine)

            if e_fine == 0 or e_coarse == 0:
                notes.append(
                    "Zero error difference at levels %d-%d; "
                    "cannot compute order for this pair" % (i, i + 2)
                )
                continue

            r_coarse = h_coarse / h_mid
            r_fine = h_mid / h_fine

            if r_coarse <= 0 or r_fine <= 0:
                continue

            p = math.log(e_coarse / e_fine) / math.log(r_coarse)
            observed_orders.append(p)

            if p < 0:
                notes.append(
                    "Negative observed order (%.2f) at levels %d-%d: "
                    "possible divergence" % (p, i, i + 2)
                )

    # Mean order
    mean_order = None
    if observed_orders:
        mean_order = sum(observed_orders) / len(observed_orders)

    # Check for pre-asymptotic behavior
    in_asymptotic_range = True
    if len(observed_orders) >= 2:
        for i in range(len(observed_orders) - 1):
            p1 = observed_orders[i]
            p2 = observed_orders[i + 1]
            avg = (abs(p1) + abs(p2)) / 2.0
            if avg > 0 and abs(p1 - p2) / avg > 0.5:
                in_asymptotic_range = False
                notes.append("Pre-asymptotic behavior detected: observed order varies >50%% between consecutive pairs")
                break
    elif len(observed_orders) == 1:
        in_asymptotic_range = True
    else:
        in_asymptotic_range = None

    # Richardson extrapolation from finest two levels
    richardson_extrapolated_value = None
    h_fine = spacings_sorted[-1]
    h_next = spacings_sorted[-2]
    f_fine = values_sorted[-1]
    f_next = values_sorted[-2]
    r = h_next / h_fine

    if mean_order is not None and mean_order > 0:
        richardson_extrapolated_value = f_fine + (f_fine - f_next) / (r ** mean_order - 1)
    elif expected_order is not None and expected_order > 0:
        richardson_extrapolated_value = f_fine + (f_fine - f_next) / (r ** expected_order - 1)
        notes.append("Richardson extrapolation used expected order (no observed order available)")

    # Convergence assessment
    convergence_assessment = "unknown"
    if expected_order is not None and mean_order is not None:
        if abs(mean_order - expected_order) / expected_order <= 0.1:
            convergence_assessment = "PASS: observed order (%.2f) within 10%% of expected (%.2f)" % (
                mean_order, expected_order,
            )
        else:
            convergence_assessment = "FAIL: observed order (%.2f) differs from expected (%.2f) by >10%%" % (
                mean_order, expected_order,
            )
    elif mean_order is not None:
        convergence_assessment = "Observed order: %.2f (no expected order given for comparison)" % mean_order
    elif len(spacings_sorted) < 3:
        convergence_assessment = "Insufficient levels to compute observed order (need >= 3)"

    return {
        "inputs": {
            "spacings": spacings_sorted,
            "values": values_sorted,
            "expected_order": expected_order,
        },
        "results": {
            "observed_orders": observed_orders,
            "mean_order": mean_order,
            "richardson_extrapolated_value": richardson_extrapolated_value,
            "in_asymptotic_range": in_asymptotic_range,
            "convergence_assessment": convergence_assessment,
            "notes": notes,
        },
    }


def main():
    args = parse_args()
    try:
        spacings = [float(x) for x in args.spacings.split(",")]
        values = [float(x) for x in args.values.split(",")]
    except ValueError:
        print("Error: spacings and values must be comma-separated numbers", file=sys.stderr)
        sys.exit(2)

    try:
        result = compute_h_refinement(spacings, values, args.expected_order)
    except ValueError as exc:
        print("Error: %s" % exc, file=sys.stderr)
        sys.exit(2)

    if args.json:
        print(json.dumps(result, indent=2))
    else:
        r = result["results"]
        print("Spatial Convergence Study (h-refinement)")
        print("  levels: %d" % len(spacings))
        if r["observed_orders"]:
            print("  observed orders: %s" % ", ".join("%.4f" % o for o in r["observed_orders"]))
        if r["mean_order"] is not None:
            print("  mean order: %.4f" % r["mean_order"])
        if r["richardson_extrapolated_value"] is not None:
            print("  Richardson extrapolated value: %.6g" % r["richardson_extrapolated_value"])
        if r["in_asymptotic_range"] is not None:
            print("  in asymptotic range: %s" % r["in_asymptotic_range"])
        print("  assessment: %s" % r["convergence_assessment"])
        for note in r["notes"]:
            print("  note: %s" % note)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Richardson extrapolation for solution verification."""
import argparse
import json
import math
import sys


def parse_args():
    parser = argparse.ArgumentParser(
        description="Perform Richardson extrapolation to estimate grid/timestep-independent solution.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--spacings", type=str, required=True,
        help="Comma-separated spacings (e.g. 0.02,0.01)",
    )
    parser.add_argument(
        "--values", type=str, required=True,
        help="Comma-separated solution values at each spacing level",
    )
    parser.add_argument(
        "--order", type=float, required=True,
        help="Assumed convergence order",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def compute_richardson_extrapolation(spacings, values, order):
    """Perform Richardson extrapolation.

    Parameters
    ----------
    spacings : list of float
        Grid spacings or timestep sizes (positive).
    values : list of float
        Solution values at each level.
    order : float
        Assumed convergence order (positive).

    Returns
    -------
    dict
        Results with extrapolated_value, error_estimate, relative_error_estimate,
        and optionally observed_order and order_consistent.
    """
    if len(spacings) != len(values):
        raise ValueError("spacings and values must have the same length")
    if len(spacings) < 2:
        raise ValueError("At least 2 refinement levels required")
    if order <= 0:
        raise ValueError("order must be positive")
    for s in spacings:
        if not math.isfinite(s) or s <= 0:
            raise ValueError("All spacings must be positive finite numbers")
    for v in values:
        if not math.isfinite(v):
            raise ValueError("All values must be finite numbers")

    # Sort by spacing descending (coarsest first)
    paired = sorted(zip(spacings, values), key=lambda x: -x[0])
    spacings_sorted = [p[0] for p in paired]
    values_sorted = [p[1] for p in paired]

    # Richardson extrapolation from finest two levels
    h_fine = spacings_sorted[-1]
    h_coarse = spacings_sorted[-2]
    f_fine = values_sorted[-1]
    f_coarse = values_sorted[-2]
    r = h_coarse / h_fine

    extrapolated_value = f_fine + (f_fine - f_coarse) / (r ** order - 1)
    error_estimate = abs(f_fine - f_coarse) / (r ** order - 1)
    relative_error_estimate = error_estimate / abs(extrapolated_value) if extrapolated_value != 0 else float("inf")

    result = {
        "extrapolated_value": extrapolated_value,
        "error_estimate": error_estimate,
        "relative_error_estimate": relative_error_estimate,
    }

    # With 3+ levels, compute observed order and check consistency
    if len(spacings_sorted) >= 3:
        h1 = spacings_sorted[-1]  # finest
        h2 = spacings_sorted[-2]
        h3 = spacings_sorted[-3]  # coarsest of three
        f1 = values_sorted[-1]
        f2 = values_sorted[-2]
        f3 = values_sorted[-3]

        e_fine = abs(f2 - f1)
        e_coarse = abs(f3 - f2)

        if e_fine > 0 and e_coarse > 0:
            r_ratio = h2 / h1
            observed_order = math.log(e_coarse / e_fine) / math.log(r_ratio)
            result["observed_order"] = observed_order
            result["order_consistent"] = abs(observed_order - order) / order <= 0.1
        else:
            result["observed_order"] = None
            result["order_consistent"] = None

    return {
        "inputs": {
            "spacings": spacings_sorted,
            "values": values_sorted,
            "order": order,
        },
        "results": result,
    }


def main():
    args = parse_args()
    try:
        spacings = [float(x) for x in args.spacings.split(",")]
        values = [float(x) for x in args.values.split(",")]
    except ValueError:
        print("Error: spacings and values must be comma-separated numbers", file=sys.stderr)
        sys.exit(2)

    try:
        result = compute_richardson_extrapolation(spacings, values, args.order)
    except ValueError as exc:
        print("Error: %s" % exc, file=sys.stderr)
        sys.exit(2)

    if args.json:
        print(json.dumps(result, indent=2))
    else:
        r = result["results"]
        print("Richardson Extrapolation")
        print("  assumed order: %.2f" % args.order)
        print("  extrapolated value: %.6g" % r["extrapolated_value"])
        print("  error estimate: %.6g" % r["error_estimate"])
        print("  relative error estimate: %.6g" % r["relative_error_estimate"])
        if "observed_order" in r and r["observed_order"] is not None:
            print("  observed order: %.4f" % r["observed_order"])
            print("  order consistent: %s" % r["order_consistent"])


if __name__ == "__main__":
    main()