Benchmark script and performance-contract tests

[janetyq]

Two pieces from the performance audit's "how do we measure progress" question: a benchmark script for the timing side, run by hand, and count-based contract tests for the regression side, run in CI.

benchmarks/bench.py

Phase timings of the solver's main workloads: 3D linear elastic (build, direct, AMG-CG, both stress recoveries), 2D Poisson P1 and P2 (build, both backends, both gradient recoveries, both estimators), Newton St-VK P1 and P2, ThetaMethod heat with constant and with TimeDependent data, AdaptiveRefinement with each estimator, SIMP, Ruppert at two area caps, RedGreenRefiner, and ModalAnalysis on a P2 beam. The default sizes are those of attic/performance-audit-2026-09-03.md section 1 (a few minutes); --quick runs every workload in about 25 s. --only picks workloads and --json PATH appends one record per phase (commit, timestamp, numpy/scipy/pyamg versions, workload, phase, seconds) as JSON lines.

uv run python benchmarks/bench.py --quick
uv run python benchmarks/bench.py --only poisson2d newton

It exists to produce a before/after table to paste into a PR that touches a hot path. Ratios within one run on one machine are meaningful; absolutes are not (one laptop drifted 3x in a day; CI runners differ), so it is not in CI and not a threshold test. benchmarks is added to pyright's include list so the script stays type-checked; ruff already covers it.

tests/test_perf_contracts.py

What makes a path fast is countable and deterministic, so these assert counts through a CountingBackend (new in tests/helpers.py, wrapping any backend and counting prepare calls and factorization solves):

• LinearSolve factors exactly once and back-substitutes exactly once.
• ThetaMethod over N steps factors once and solves N times.
• NewmarkMethod over N steps factors exactly twice (the mass for the initial acceleration, then the effective operator) and solves N + 1 times.
• SensitivityAnalysis shares one factorization between forward and adjoint; a self-adjoint QoI (Compliance under homogeneous supports) makes no adjoint solve, PointValue exactly one.
• DirectBackend's ordering keeps L + U fill under 8x nnz(A) on a 100x100 P1 Poisson block (5.3x today), and the test checks that scipy's default COLAMD ordering exceeds the same bound (12.2x), so the bound really separates the two.
• AMG-CG matvecs on Poisson grow by at most 1.5x on a mesh doubling (measured 14 at 60x60, 17 at 120x120) and stay under 40; the bound is loose because the exact count is pyamg's aggregation, not this library's contract. The count is read by swapping a counting LinearOperator into the CG solver's operator, a private attribute, noted in the test.
• The constant-Source contract from Source: integrate a constant load element-wise, optimize the sampled contraction #172 (tests/test_loads.py) is referenced as part of the set rather than duplicated.

A last test checks the counting wrapper changes nothing but the counts.

Verification

Full suite 1122 passed (including the 7 new tests); ruff check and pyright clean over the whole repo with the new include. Quick preset on this branch (dcd55fa base, one laptop):

== elastic3d ==
    box_mesh (5000 tets)             0.006 s
    equation.problem (3993 dofs)     0.045 s
    solve direct                     0.119 s
    solve AMG-CG                     0.068 s
    nodal_stress('average')          0.040 s
    nodal_stress('l2')               0.082 s
== poisson2d ==
    P1 equation.problem (44402 tris, 22500 dofs)     0.240 s
    P1 solve direct                                  0.191 s
    P1 solve AMG-CG                                  0.272 s
    P1 nodal_gradient('average')                     0.179 s
    P1 nodal_gradient('l2')                          0.425 s
    P1 RecoveryEstimator                             0.906 s
    P1 ResidualEstimator                             0.697 s
    P2 equation.problem (10952 tris, 22201 dofs)     0.210 s
    P2 solve direct                                  0.511 s
    P2 solve AMG-CG                                  0.700 s
    P2 nodal_gradient('average')                     0.271 s
    P2 nodal_gradient('l2')                          1.306 s
    P2 RecoveryEstimator                             1.432 s
    P2 ResidualEstimator                             0.243 s
== newton ==
    P1 Newton solve (3600 dofs)     1.025 s
    P2 Newton solve (3422 dofs)     1.481 s
== heat ==
    100 steps, constant data (10000 dofs)                          0.683 s
    100 steps, TimeDependent source and Dirichlet (10000 dofs)     2.124 s
== adapt ==
    ResidualEstimator to 6374 tris     1.163 s
    RecoveryEstimator to 4143 tris     4.589 s
== simp ==
    SIMPModel + filter (2400 tris)     0.017 s
    5 iterations (2562 dofs)           0.322 s
== ruppert ==
    max_area=0.05: 627 tris      0.647 s
    max_area=0.02: 1413 tris     1.306 s
== redgreen ==
    RedGreenRefiner.__init__ (6962 tris)     0.089 s
    refine all -> 27848                      1.106 s
    refine 1/7 -> 56034                      2.292 s
== modal ==
    6 modes (3402 dofs)     0.259 s
total 25.8 s