#!/usr/bin/env python3
"""Tier-2 item 11 -- turn the tangent stability check back on.

This one is nearly free and closes a real gap, so it is worth running alongside
the Tier-1 batch.

`NonlinearSolveOptions.stability_modes` defaults to 6 in the solver and
`_solve_tangent_stability_modes` is fully implemented, but the two-group runner
sets `stability_modes=0`, so no accepted equilibrium has ever been checked for
stability. Newton convergence plus a positive Jacobian does not imply a stable
equilibrium, and this model has two reasons to care:

  * the shell is 0.12 um thick and the hyperosmotic compression area reaches
    62.6 % (Fig. S4d), so much of it is in minimum-principal compression, where
    a real lamina would wrinkle;
  * the hyperosmotic residual climbs about two orders of magnitude around load
    factor 0.6 before recovering (Fig. S1), which is what passing close to a
    critical point looks like.

A smooth Newton solve on a symmetric mesh will happily walk past a bifurcation
without reporting it. Report the smallest eigenvalues of the projected tangent
for every accepted equilibrium, and rule any case whose global mean resultant is
non-positive out of the model's domain of validity.
"""

from __future__ import annotations

import argparse
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent))
from _common import add_common_arguments, import_runner, patch_options  # noqa: E402


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__)
    add_common_arguments(parser)
    parser.add_argument(
        "--stability-modes",
        type=int,
        default=6,
        help="Number of smallest tangent eigenvalues to extract per equilibrium.",
    )
    parser.add_argument(
        "--element-order",
        choices=("p1", "p2"),
        default="p2",
    )
    args = parser.parse_args()

    print("Tangent stability check -- planned run")
    print(f"  project root    : {args.project_root}")
    print(f"  output root     : {args.out_root}")
    print(f"  primary mesh    : {args.primary_mesh}")
    print(f"  stability modes : {args.stability_modes}")
    print(
        "  Acceptance: every reported equilibrium must have a strictly positive "
        "smallest eigenvalue of the projected tangent."
    )
    if args.dry_run:
        print("\n--dry-run: nothing solved.")
        return 0

    runner = import_runner(args.project_root)
    patch_options(runner, stability_modes=args.stability_modes)
    args.out_root.mkdir(parents=True, exist_ok=True)

    out_dir = args.out_root / "stability_check"
    print(f"\n=== stability check -> {out_dir} ===", flush=True)
    runner.run_analysis(
        out_dir,
        primary_mesh=args.primary_mesh,
        run_intermediate_p2=(args.element_order == "p2"),
        resume=args.resume,
        sparse_backend=args.sparse_backend,
        rigid_constraint_linear_solver=args.rigid_constraint_linear_solver,
    )
    print(
        "\nDone. The per-case JSON now carries a stability block "
        "(stability_analysis_requested = true). Pull the smallest eigenvalue "
        "for each accepted equilibrium into Figure S1."
    )
    return 0


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