VEP investigation: bug fixes, benchmark suite, TI-VEP fault stability

[lmoresi]

Summary

Cumulative work from the VEP loading/unloading investigation. Three groups of changes, all preserving default behaviour for existing callers.

Bug fixes

• MathematicalMixin.__getattr__ infinite recursion on sym — when a subclass's sym @Property raises AttributeError during partial __init__ (e.g. SwarmVariable.sym before self._meshVar is set), Python falls back to __getattr__, which calls _validate_sym, which accesses self.sym, which re-enters __getattr__ for name="sym" — unbounded loop, 100% CPU, no output. Reproduces in any test that creates a VE solver after a previous one in the same process. Added "sym" to the existing guard.

• VEP yield-surface drift under variable dt — projection-snapshot fix in SemiLagrangian DDt that prevents σ from breaching ±τ_y by ~30% under repeated halve/double-dt changes. Two new regression tests (test_1052_VEP_stability_regression).

• TI-VEP BDF-2 instability with spatial yield_stress — newly diagnosed in this session. At order=2 with an influence_function-localised yield_stress field (i.e. realistic fault simulations), the simulation drifts unstably: |σ_xy| → 10⁸ over ~10 t_r at θ=±15°. Iso VEP with the same spatial τ_y was masked stable by yield clipping; pure VE on the same mesh+BC+IC also drifts. Recommended workaround: use order=1 for TI-VEP fault problems (the class default). The bdf_blend knob is exposed for users who must use order=2 on a uniform-yield_stress problem and want explicit damping. Default bdf_blend=1.0 (no silent default-degradation). Three planning items filed for the deeper fix (flux_with_history_order(k) API, UWexpression mutation audit, exponential-integration alternative).

• Restore timestep=dt in run_ve_vep_oscillatory_plot.py — one-line test-script fix.

• Smooth yield-mode retired — under-clipped the yield surface by ~50% with no compensating benefit. Setter now redirects users to softmin.

New solver APIs

• DDt.set_initial_history(values, dt=None) on SemiLagrangian and Eulerian — public API for planting BDF history at the start of a run. Two use cases: (a) analytical IC for benchmarks (no startup transient), (b) checkpoint restart (resume multistep history without re-ramping effective_order from BDF-1). Replaces a four-private-attribute pattern in bench_ve_harmonic.py. 6 unit tests in test_1052_ddt_set_initial_history.

• SaddlePt.set_jacobian_F1_source(F1_alt, linesearch=\"cp\") — optional override that autodiffs an alternate F1 expression for the velocity-block Jacobian while the residual F1 stays unchanged. Enables inexact-Newton tricks (smooth-Jacobian / sharp-residual). Default-installs cp linesearch because bt (PETSc default) rejects useful steps when the Jacobian is inexact.

• divergence_retries kwarg on all SNES solvers (default 0). Warm-start re-solve helps Newton step off transient kinks.

Benchmark suite

A self-contained run / log / replot architecture for VE/VEP analytical benchmarks (in docs/advanced/benchmarks/):

• VE harmonic with peak-start IC — eliminates startup transient. BDF-2 max|err| = 1.34e-3.
• VE / VEP square-wave — BDF history at BC discontinuities, projection-snapshot fix verification.
• Variable-dt VEP at fixed and variable Δt with both Min and softmin yield modes. Demonstrated: switching snes_linesearch_type from bt to cp for inexact-Jacobian problems eliminates 196/413 spurious line-search rejections without losing accuracy.
• TI-VEP harmonic angled-fault at θ ∈ {0°, ±15°} × τ_y ∈ {0.15, 0.30}. Dedicated plot script with V_top driving overlay.

The benchmark npz files are self-contained (parameters + traces + analytical reference); replotting is decoupled from the slow simulation.

Empirical impact (from the benchmark traces, all on this branch)

Problem	BDF-1 max|err|	BDF-2 max|err|	Notes
VE harmonic	1.82e-2	1.34e-3	BDF-2 wins — pure VE benefits from second-order
VEP square-wave (Min, BDF-2 default bt)	—	9.4e-2	95% of steps are 1-iter — exits early on rel-tol
VEP square-wave (softmin Jac + cp linesearch)	—	1.14e-2	0 line-search rejections, mean 0.72 iter/step
TI-VEP fault at order=1	clean	n/a	Recommended for fault simulations
TI-VEP fault at order=2 (spatial τ_y)	n/a	blow-up to 10⁸	Use order=1 instead

Test plan

• All 20 VE/VEP regression tests pass at HEAD
• ./uw build succeeds
• TI-VEP harmonic benchmark produces stable traces at order=1 for all (θ, τ_y) combinations
• CI smoke (level_1 + tier_a)

Out of scope / follow-ups

Filed in ~/Library/CloudStorage/Box-Box/Planning/underworld.md (Bugs and Nice to Have):

1. TI-VEP BDF-2 + spatial yield_stress — workaround shipped; deeper fix needs cm.flux_with_history_order(k) API and the inexact-Newton path investigated
2. UWexpression mutation pattern in constitutive_models.py (5 sites) — silent state aliasing, latent disaster
3. sympy.simplify=True defaults in evaluate* and visualisation — flip them
4. Exponential integration for the VE/VEP constitutive update — structural rethink that would sidestep BDF stability altogether (Hochbruck-Ostermann style); lead from offline review

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

This PR bundles three fixes discovered during the VEP loading/unloading investigation, aiming to preserve legacy defaults while improving robustness for VE/VEP workflows and SNES divergence handling.

Changes:

• Restores the required timestep= kwarg in tests/run_ve_vep_oscillatory_plot.py so the VE/VEP oscillatory script runs end-to-end again.
• Prevents infinite recursion in MathematicalMixin.__getattr__ when sym property access fails during partial initialization.
• Adds an opt-in divergence_retries kwarg to SNES-based solve() methods and implements warm-start retries in the Cython SNES base layer.

Reviewed changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated 2 comments.

File	Description
tests/run_ve_vep_oscillatory_plot.py	Fixes test-script invocation by passing timestep=dt to VE_Stokes.solve().
src/underworld3/utilities/mathematical_mixin.py	Guards __getattr__('sym') to avoid recursion when sym raises AttributeError during early init.
src/underworld3/systems/solvers.py	Threads divergence_retries through high-level solver wrappers and documents the new kwarg.
src/underworld3/cython/petsc_generic_snes_solvers.pyx	Implements _snes_solve_with_retries and wires it into SNES solve paths (Scalar/Vector/MultiComponent/Stokes).

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.