heshamfs/materials-simulation-skills/nonlinear-solvers

Nonlinear Solvers

Goal

Provide a universal workflow to select a nonlinear solver, configure globalization strategies, and diagnose convergence for root-finding, optimization, and least-squares problems.

Requirements

• Python 3.10+
• NumPy (for Jacobian diagnostics)
• SciPy (optional, for advanced analysis)

Inputs to Gather

Decision Guidance

Solver Selection Flowchart

Is Jacobian available and cheap?
├── YES → Problem size?
│   ├── Small (n < 1000) → Newton (full)
│   └── Large (n ≥ 1000) → Newton-Krylov
└── NO → Is objective smooth?
    ├── YES → Memory limited?
    │   ├── YES → L-BFGS or Broyden
    │   └── NO → BFGS
    └── NO → Anderson acceleration or Picard

Quick Reference

Script Outputs (JSON Fields)

Workflow

1. Characterize problem - Identify type, size, Jacobian availability
2. Select solver - Run scripts/solver_selector.py
3. Choose globalization - Run scripts/globalization_advisor.py
4. Analyze Jacobian - If available, run scripts/jacobian_diagnostics.py
5. Monitor residuals - During solve, use scripts/residual_monitor.py
6. Analyze convergence - Run scripts/convergence_analyzer.py
7. Evaluate steps - For trust region, use scripts/step_quality.py

Conversational Workflow Example

User: My Newton solver for a phase-field simulation is converging very slowly. After 50 iterations, the residual only dropped from 1 to 0.1.

Agent workflow:

1. Analyze convergence:

``bash python3 scripts/convergence_analyzer.py --residuals 1,0.8,0.6,0.5,0.4,0.3,0.2,0.15,0.12,0.1 --json ``

1. Check globalization strategy:

``bash python3 scripts/globalization_advisor.py --problem-type root-finding --jacobian-quality ill-conditioned --previous-failures 0 --json ``

1. Recommend: Switch to trust region with Levenberg-Marquardt regularization, or use Newton-Krylov with better preconditioning.

Pre-Solve Checklist

• [ ] Confirm problem type (root-finding, optimization, least-squares)
• [ ] Assess Jacobian availability and cost
• [ ] Check initial guess quality
• [ ] Set appropriate tolerances
• [ ] Choose globalization strategy
• [ ] Prepare to monitor convergence

CLI Examples

# Select solver for large unconstrained optimization
python3 scripts/solver_selector.py --size 50000 --smooth --memory-limited --json
 
# Select solver for a small nonlinear least-squares (data-fitting) problem
python3 scripts/solver_selector.py --problem-type least-squares --size 6 --jacobian-available --smooth --json
 
# Analyze convergence from residual history
python3 scripts/convergence_analyzer.py --residuals 1,0.1,0.01,0.001,0.0001 --tolerance 1e-6 --json
 
# Diagnose Jacobian quality
python3 scripts/jacobian_diagnostics.py --matrix jacobian.txt --json
 
# Get globalization recommendation
python3 scripts/globalization_advisor.py --problem-type optimization --jacobian-quality good --json
 
# Globalization for a distant initial guess (favors trust region)
python3 scripts/globalization_advisor.py --problem-type root-finding --jacobian-quality good --far-from-solution --json
 
# Monitor residual patterns
python3 scripts/residual_monitor.py --residuals 1,0.8,0.9,0.7,0.75,0.6 --target-tolerance 1e-8 --json
 
# Evaluate step quality for trust region
python3 scripts/step_quality.py --predicted-reduction 0.5 --actual-reduction 0.4 --step-norm 0.8 --gradient-norm 1.0 --trust-radius 1.0 --json

Error Handling

Interpretation Guidance

Convergence Type

Jacobian Quality

Step Quality (Trust Region)

Verification checklist

Do not trust a "solved" claim until these concrete artifacts are recorded:

• [ ] Logged the full residual norm history and ran convergence_analyzer.py --residuals <history>; recorded convergence_type and estimated_rate, and confirmed converged: true against the actual solver tolerance (not the default 1e-10).
• [ ] Confirmed the residual sequence is monotone-decreasing or fed it to residual_monitor.py; recorded patterns_detected and verified it does NOT include diverging, oscillating, plateau, or slow_convergence while still above tolerance.
• [ ] If a Jacobian is available, ran jacobian_diagnostics.py --matrix J.txt and recorded condition_number and jacobian_quality; for an analytic Jacobian, passed --finite-diff-matrix and confirmed finite_diff_error is below ~1e-2 (no "Large discrepancy" note).
• [ ] Checked rank_deficient from jacobian_diagnostics.py is false (or documented why a rank-deficient/near-singular Jacobian is expected and that Levenberg-Marquardt regularization is in use).
• [ ] For a trust-region solve, evaluated accepted steps with step_quality.py and recorded the reduction ratio; confirmed accepted steps have ratio >= 0.25 (not very_poor/poor) and that the trust_radius_action matches the recorded ρ.
• [ ] Recorded the solver and globalization actually used and confirmed they match solver_selector.py and globalization_advisor.py recommendations for the stated problem type, size, and Jacobian quality (e.g., large/expensive-Jacobian → Newton-Krylov; least-squares → Levenberg-Marquardt trust region).
• [ ] Re-confirmed convergence after any change to tolerance, initial guess, or preconditioner — the convergence type can flip (e.g., quadratic → linear/stagnated) and must be re-classified, not assumed.

Common pitfalls & rationalizations

Security