USP

Gives AI coding agents deep domain expertise in numerical methods and scientific computing, eliminating repetitive guidance. Features 17 skills, 67 scripts, and robust security, all adhering to the open Agent Skills standard for portabilit…

Use cases

• 01Checking numerical stability (CFL, von Neumann) for simulations.
• 02Performing grid and time convergence studies with Richardson extrapolation.
• 03Selecting appropriate linear and nonlinear solvers for Ax=b systems.
• 04Generating and evaluating meshes for PDE simulations based on physics scales.
• 05Creating SLURM job scripts for HPC deployments.

Detected files (8)

• skills/core-numerical/time-stepping/SKILL.mdskill
  Show content (6511 bytes)

  ---
  name: time-stepping
  description: >
    Plan and control time-step policies for transient simulations — couple
    CFL and physics-based stability limits with adaptive stepping, ramp initial
    transients through sharp gradients or phase changes, schedule output intervals
    and checkpoint cadence, and plan restart strategies for long-running jobs.
    Use when choosing dt for a new simulation, diagnosing adaptive time-step
    oscillations, deciding checkpoint frequency to minimize lost work, or
    setting up output schedules aligned with physical time scales, even if
    the user only says "my run is too slow" or "how often should I save."
  allowed-tools: Read, Bash, Write, Grep, Glob
  metadata:
    author: HeshamFS
    version: "1.1.0"
    security_tier: high
    security_reviewed: true
    tested_with:
      - claude-code
      - gemini-cli
      - vs-code-copilot
    eval_cases: 2
    last_reviewed: "2026-03-26"
  ---
  
  # Time Stepping
  
  ## Goal
  
  Provide a reliable workflow for choosing, ramping, and monitoring time steps plus output/checkpoint cadence.
  
  ## Requirements
  
  - Python 3.8+
  - No external dependencies (uses stdlib)
  
  ## Inputs to Gather
  
  | Input | Description | Example |
  |-------|-------------|---------|
  | Stability limits | CFL/Fourier/reaction limits | `dt_max = 1e-4` |
  | Target dt | Desired time step | `1e-5` |
  | Total run time | Simulation duration | `10 s` |
  | Output interval | Time between outputs | `0.1 s` |
  | Checkpoint cost | Time to write checkpoint | `120 s` |
  
  ## Decision Guidance
  
  ### Time Step Selection
  
  ```
  Is stability limit known?
  ├── YES → Use min(dt_target, dt_limit × safety)
  └── NO → Start conservative, increase adaptively
  
  Need ramping for startup?
  ├── YES → Start at dt_init, ramp to dt_target over N steps
  └── NO → Use dt_target from start
  ```
  
  ### Ramping Strategy
  
  | Problem Type | Ramp Steps | Initial dt |
  |--------------|------------|------------|
  | Smooth IC | None needed | Full dt |
  | Sharp gradients | 5-10 | 0.1 × dt |
  | Phase change | 10-20 | 0.01 × dt |
  | Cold start | 10-50 | 0.001 × dt |
  
  ## Script Outputs (JSON Fields)
  
  | Script | Key Outputs |
  |--------|-------------|
  | `scripts/timestep_planner.py` | `dt_limit`, `dt_recommended`, `ramp_schedule` |
  | `scripts/output_schedule.py` | `output_times`, `interval`, `count` |
  | `scripts/checkpoint_planner.py` | `checkpoint_interval`, `checkpoints`, `overhead_fraction` |
  
  ## Workflow
  
  1. **Get stability limits** - Use numerical-stability skill
  2. **Plan time stepping** - Run `scripts/timestep_planner.py`
  3. **Schedule outputs** - Run `scripts/output_schedule.py`
  4. **Plan checkpoints** - Run `scripts/checkpoint_planner.py`
  5. **Monitor during run** - Adjust dt if limits change
  
  ## Conversational Workflow Example
  
  **User**: I'm running a 10-hour phase-field simulation. How often should I checkpoint?
  
  **Agent workflow**:
  1. Plan checkpoints based on acceptable lost work:
     ```bash
     python3 scripts/checkpoint_planner.py --run-time 36000 --checkpoint-cost 120 --max-lost-time 1800 --json
     ```
  2. Interpret: Checkpoint every 30 minutes, overhead ~0.7%, max 30 min lost work on crash.
  
  ## Pre-Run Checklist
  
  - [ ] Confirm dt limits from stability analysis
  - [ ] Define ramping strategy for transient startup
  - [ ] Choose output interval consistent with physics time scales
  - [ ] Plan checkpoints based on restart risk
  - [ ] Re-evaluate dt after parameter changes
  
  ## CLI Examples
  
  ```bash
  # Plan time stepping with ramping
  python3 scripts/timestep_planner.py --dt-target 1e-4 --dt-limit 2e-4 --safety 0.8 --ramp-steps 10 --json
  
  # Schedule output times
  python3 scripts/output_schedule.py --t-start 0 --t-end 10 --interval 0.1 --json
  
  # Plan checkpoints for long run
  python3 scripts/checkpoint_planner.py --run-time 36000 --checkpoint-cost 120 --max-lost-time 1800 --json
  ```
  
  ## Error Handling
  
  | Error | Cause | Resolution |
  |-------|-------|------------|
  | `dt-target must be positive` | Invalid time step | Use positive value |
  | `t-end must be > t-start` | Invalid time range | Check time bounds |
  | `checkpoint-cost must be < run-time` | Checkpoint too expensive | Reduce checkpoint size |
  
  ## Interpretation Guidance
  
  ### dt Behavior
  
  | Observation | Meaning | Action |
  |-------------|---------|--------|
  | dt stable at target | Good | Continue |
  | dt shrinking | Stability issue | Check CFL, reduce target |
  | dt oscillating | Borderline stability | Add safety factor |
  
  ### Checkpoint Overhead
  
  | Overhead | Acceptability |
  |----------|---------------|
  | < 1% | Excellent |
  | 1-5% | Good |
  | 5-10% | Acceptable |
  | > 10% | Too frequent, increase interval |
  
  ## Security
  
  ### Input Validation
  - All numeric parameters (`dt-target`, `dt-limit`, `safety`, `t-start`, `t-end`, `interval`, `run-time`, `checkpoint-cost`, `max-lost-time`) are validated as finite positive numbers
  - `ramp-steps` is validated as a non-negative integer with an upper bound
  - Time range consistency is enforced (`t-end` must exceed `t-start`; `checkpoint-cost` must be less than `run-time`)
  
  ### File Access
  - Scripts read no external files; all inputs are provided via CLI arguments
  - Scripts write only to stdout (JSON output); no files are created unless the agent explicitly uses the Write tool
  
  ### Tool Restrictions
  - **Read**: Used to inspect script source, references, and user configuration files
  - **Bash**: Used to execute the three Python planning scripts (`timestep_planner.py`, `output_schedule.py`, `checkpoint_planner.py`) with explicit argument lists
  - **Write**: Used to save generated time-step plans or checkpoint schedules; writes are scoped to the user's working directory
  - **Grep/Glob**: Used to locate relevant files and search references
  
  ### Safety Measures
  - No `eval()`, `exec()`, or dynamic code generation
  - All subprocess calls use explicit argument lists (no `shell=True`)
  - Scripts use only Python standard library; no pickle loading or deserialization of untrusted data
  - All output is deterministic JSON with no shell-interpretable content
  
  ## Limitations
  
  - **Not adaptive control**: Plans static schedules, not runtime adaptation
  - **Assumes constant physics**: If parameters change, re-plan
  
  ## References
  
  - `references/cfl_coupling.md` - Combining multiple stability limits
  - `references/ramping_strategies.md` - Startup policies
  - `references/output_checkpoint_guidelines.md` - Cadence rules
  
  ## Version History
  
  - **v1.1.0** (2024-12-24): Enhanced documentation, decision guidance, examples
  - **v1.0.0**: Initial release with 3 planning scripts
  

• skills/core-numerical/numerical-integration/SKILL.mdskill
  Show content (8556 bytes)

  ---
  name: numerical-integration
  description: >
    Select and configure time integration methods for ODE and PDE simulations —
    choose among explicit Runge-Kutta, BDF, Rosenbrock, and Adams families,
    set relative and absolute error tolerances, implement adaptive step-size
    control with I/PI/PID controllers, plan IMEX operator splitting for mixed
    stiff and non-stiff terms, and estimate splitting errors. Use when picking
    an integrator for a new simulation, diagnosing step rejections or tolerance
    failures, setting up operator splitting for phase-field or reaction-diffusion
    problems, or deciding between explicit and implicit time marching, even if
    the user only says "my solver keeps rejecting steps" or "which ODE method
    should I use."
  allowed-tools: Read, Write, Grep, Glob
  metadata:
    author: HeshamFS
    version: "1.1.0"
    security_tier: medium
    security_reviewed: true
    tested_with:
      - claude-code
      - gemini-cli
      - vs-code-copilot
    eval_cases: 2
    last_reviewed: "2026-03-26"
  ---
  
  # Numerical Integration
  
  ## Goal
  
  Provide a reliable workflow to select integrators, set tolerances, and manage adaptive time stepping for time-dependent simulations.
  
  ## Requirements
  
  - Python 3.8+
  - NumPy (for some scripts)
  - No heavy dependencies for core functionality
  
  ## Inputs to Gather
  
  | Input | Description | Example |
  |-------|-------------|---------|
  | Problem type | ODE/PDE, stiff/non-stiff | `stiff PDE` |
  | Jacobian available | Can compute ∂f/∂u? | `yes` |
  | Target accuracy | Desired error level | `1e-6` |
  | Constraints | Memory, implicit allowed? | `implicit OK` |
  | Time scale | Characteristic time | `1e-3 s` |
  
  ## Decision Guidance
  
  ### Choosing an Integrator
  
  ```
  Is the problem stiff?
  ├── YES → Is Jacobian available?
  │   ├── YES → Use Rosenbrock or BDF
  │   └── NO → Use BDF with numerical Jacobian
  └── NO → Is high accuracy needed?
      ├── YES → Use RK45 or DOP853
      └── NO → Use RK4 or Adams-Bashforth
  ```
  
  ### Stiff vs Non-Stiff Detection
  
  | Symptom | Likely Stiff | Action |
  |---------|--------------|--------|
  | dt shrinks to tiny values | Yes | Switch to implicit |
  | Eigenvalues span many decades | Yes | Use BDF/Radau |
  | Smooth solution, reasonable dt | No | Stay explicit |
  
  ## Script Outputs (JSON Fields)
  
  | Script | Key Outputs |
  |--------|-------------|
  | `scripts/error_norm.py` | `error_norm`, `scale_min`, `scale_max` |
  | `scripts/adaptive_step_controller.py` | `accept`, `dt_next`, `factor` |
  | `scripts/integrator_selector.py` | `recommended`, `alternatives`, `notes` |
  | `scripts/imex_split_planner.py` | `implicit_terms`, `explicit_terms`, `splitting_strategy` |
  | `scripts/splitting_error_estimator.py` | `error_estimate`, `substeps` |
  
  ## Workflow
  
  1. **Classify stiffness** - Check eigenvalue spread or use stiffness_detector
  2. **Choose tolerances** - See `references/tolerance_guidelines.md`
  3. **Select integrator** - Run `scripts/integrator_selector.py`
  4. **Compute error norms** - Use `scripts/error_norm.py` for step acceptance
  5. **Adapt step size** - Use `scripts/adaptive_step_controller.py`
  6. **Plan IMEX/splitting** - If mixed stiff/nonstiff, use `scripts/imex_split_planner.py`
  7. **Validate convergence** - Repeat with tighter tolerances
  
  ## Conversational Workflow Example
  
  **User**: I'm solving the Allen-Cahn equation with a stiff double-well potential. What integrator should I use?
  
  **Agent workflow**:
  1. Check integrator options:
     ```bash
     python3 scripts/integrator_selector.py --stiff --jacobian-available --accuracy high --json
     ```
  2. Plan the IMEX splitting (diffusion implicit, reaction explicit):
     ```bash
     python3 scripts/imex_split_planner.py --stiff-terms diffusion --nonstiff-terms reaction --coupling weak --json
     ```
  3. Recommend: Use IMEX-BDF2 with diffusion term implicit, double-well reaction explicit.
  
  ## Pre-Integration Checklist
  
  - [ ] Identify stiffness and dominant time scales
  - [ ] Set `rtol`/`atol` consistent with physics and units
  - [ ] Confirm integrator compatibility with stiffness
  - [ ] Use error norm to accept/reject steps
  - [ ] Verify convergence with tighter tolerance run
  
  ## CLI Examples
  
  ```bash
  # Select integrator for stiff problem with Jacobian
  python3 scripts/integrator_selector.py --stiff --jacobian-available --accuracy high --json
  
  # Compute scaled error norm
  python3 scripts/error_norm.py --error 0.01,0.02 --solution 1.0,2.0 --rtol 1e-3 --atol 1e-6 --json
  
  # Adaptive step control with PI controller
  python3 scripts/adaptive_step_controller.py --dt 1e-2 --error-norm 0.8 --order 4 --controller pi --json
  
  # Plan IMEX splitting
  python3 scripts/imex_split_planner.py --stiff-terms diffusion,elastic --nonstiff-terms reaction --coupling strong --json
  
  # Estimate splitting error
  python3 scripts/splitting_error_estimator.py --dt 1e-4 --scheme strang --commutator-norm 50 --target-error 1e-6 --json
  ```
  
  ## Error Handling
  
  | Error | Cause | Resolution |
  |-------|-------|------------|
  | `rtol and atol must be positive` | Invalid tolerances | Use positive values |
  | `error-norm must be positive` | Negative error norm | Check error computation |
  | `Unknown controller` | Invalid controller type | Use `i`, `pi`, or `pid` |
  | `Splitting requires at least one term` | Empty term list | Specify stiff or nonstiff terms |
  
  ## Interpretation Guidance
  
  ### Error Norm Values
  
  | Error Norm | Meaning | Action |
  |------------|---------|--------|
  | < 1.0 | Step acceptable | Accept, maybe increase dt |
  | ≈ 1.0 | At tolerance boundary | Accept with current dt |
  | > 1.0 | Step rejected | Reject, reduce dt |
  
  ### Controller Selection
  
  | Controller | Properties | Best For |
  |------------|------------|----------|
  | I (integral) | Simple, some overshoot | Non-stiff, moderate accuracy |
  | PI (proportional-integral) | Smooth, robust | General use |
  | PID | Aggressive adaptation | Rapidly varying dynamics |
  
  ### IMEX Strategy
  
  | Coupling | Strategy |
  |----------|----------|
  | Weak | Simple operator splitting |
  | Moderate | Strang splitting |
  | Strong | Fully coupled IMEX-RK |
  
  ## Security
  
  ### Input Validation
  - All numeric inputs (`dt`, `rtol`, `atol`, `error_norm`, `stiffness_ratio`, `commutator_norm`, etc.) are validated as finite numbers at the function boundary
  - `imex_split_planner.py` validates term names against `[a-zA-Z_][a-zA-Z0-9_ -]*` with length and count limits, preventing injection payloads in user-supplied term lists
  - Comma-separated value lists are capped at 100,000 entries to prevent resource exhaustion
  - Numeric bounds enforced: `dimension` capped at 10 billion, `order` at 20, `stiffness_ratio` at 1e30
  - `--controller` is validated against a fixed allowlist (`i`, `pi`, `pid`)
  - `--scheme` is validated against known splitting schemes (`lie`, `strang`)
  
  ### File Access
  - Scripts read no external files; all inputs are provided via CLI arguments
  - Scripts write only to stdout (JSON output); no files are created unless the agent explicitly uses the Write tool
  
  ### Tool Restrictions
  - **Read**: Used to inspect script source, references, and user configuration files
  - **Write**: Used to save integrator recommendations or splitting plans; writes are scoped to the user's working directory
  - **Grep/Glob**: Used to locate relevant files and search references
  - The skill's `allowed-tools` excludes `Bash` to prevent the agent from executing arbitrary commands when processing user-provided inputs
  
  ### Safety Measures
  - No `eval()`, `exec()`, or dynamic code generation
  - All subprocess calls use explicit argument lists (no `shell=True`)
  - Reduced tool surface (no Bash) limits the agent to read/write operations only
  - Term names are sanitized before use, preventing shell metacharacter injection
  
  ## Limitations
  
  - **No automatic stiffness detection**: Use stiffness_detector from numerical-stability
  - **Splitting assumes separability**: Terms must be cleanly separable
  - **Jacobian requirement**: Some methods need analytical or numerical Jacobian
  
  ## References
  
  - `references/method_catalog.md` - Integrator options and properties
  - `references/tolerance_guidelines.md` - Choosing rtol/atol
  - `references/error_control.md` - Error norm and adaptation formulas
  - `references/imex_guidelines.md` - Stiff/non-stiff splitting
  - `references/splitting_catalog.md` - Operator splitting patterns
  - `references/multiphase_field_patterns.md` - Phase-field specific splits
  
  ## Version History
  
  - **v1.1.0** (2024-12-24): Enhanced documentation, decision guidance, examples
  - **v1.0.0**: Initial release with 5 integration scripts
  

• skills/core-numerical/convergence-study/SKILL.mdskill
  Show content (6265 bytes)

  ---
  name: convergence-study
  description: >
    Perform spatial and temporal convergence analysis for solution verification —
    compute observed convergence orders from grid or timestep refinement studies,
    apply Richardson extrapolation to estimate discretization error, and calculate
    the Grid Convergence Index (GCI) per ASME V&V 20 standards. Use when verifying
    that a numerical solution converges at the expected rate, estimating the
    error on the finest mesh, checking whether grids are in the asymptotic range,
    or preparing formal verification reports, even if the user only asks "is my
    mesh fine enough" or "how accurate is my solution."
  allowed-tools:
    - Bash
    - Read
  metadata:
    author: HeshamFS
    version: "1.1.0"
    security_tier: high
    security_reviewed: true
    tested_with:
      - claude-code
      - gemini-cli
      - vs-code-copilot
    eval_cases: 2
    last_reviewed: "2026-03-26"
  ---
  
  # Convergence Study
  
  ## Goal
  
  Provide script-driven convergence analysis for verifying that numerical solutions converge at the expected rate as the mesh or timestep is refined.
  
  ## Requirements
  
  - Python 3.8+
  - NumPy (not required; scripts use only math stdlib)
  
  ## Inputs to Gather
  
  | Input | Description | Example |
  |-------|-------------|---------|
  | Grid spacings | Sequence of mesh sizes (coarse to fine) | `0.4,0.2,0.1,0.05` |
  | Timestep sizes | Sequence of dt values | `0.04,0.02,0.01` |
  | Solution values | QoI at each refinement level | `1.16,1.04,1.01,1.0025` |
  | Expected order | Formal order of the numerical scheme | `2.0` |
  | Safety factor | GCI safety factor (1.25 default) | `1.25` |
  
  ## Script Outputs (JSON Fields)
  
  | Script | Key Outputs |
  |--------|-------------|
  | `scripts/h_refinement.py` | `results.observed_orders`, `results.mean_order`, `results.richardson_extrapolated_value`, `results.convergence_assessment` |
  | `scripts/dt_refinement.py` | Same as h_refinement but for temporal convergence |
  | `scripts/richardson_extrapolation.py` | `results.extrapolated_value`, `results.error_estimate`, `results.observed_order` |
  | `scripts/gci_calculator.py` | `results.observed_order`, `results.gci_fine`, `results.gci_coarse`, `results.asymptotic_ratio`, `results.in_asymptotic_range` |
  
  ## Workflow
  
  1. **Run grid/timestep refinement study** with at least 3 levels
  2. **Compute observed convergence order** with `h_refinement.py` or `dt_refinement.py`
  3. **Compare** observed order to expected order of the scheme
  4. **Estimate discretization error** via Richardson extrapolation
  5. **Report GCI** for formal solution verification using `gci_calculator.py`
  6. **Document** convergence results and any anomalies
  
  ## Decision Guidance
  
  ```
  Do you have 3+ refinement levels?
  +-- YES --> Run h_refinement.py or dt_refinement.py
  |           +-- Observed order matches expected? --> Solution verified
  |           +-- Order too low? --> Check: pre-asymptotic, coding error, insufficient resolution
  |           +-- Order too high? --> Check: superconvergence or cancellation effects
  +-- NO (only 2 levels) --> Use richardson_extrapolation.py with assumed order
                             (less reliable without order verification)
  ```
  
  ## CLI Examples
  
  ```bash
  # Spatial convergence with 4 grid levels
  python3 scripts/h_refinement.py --spacings 0.4,0.2,0.1,0.05 --values 1.16,1.04,1.01,1.0025 --expected-order 2.0 --json
  
  # Temporal convergence with 3 timestep levels
  python3 scripts/dt_refinement.py --timesteps 0.04,0.02,0.01 --values 2.12,2.03,2.0075 --expected-order 2.0 --json
  
  # Richardson extrapolation with assumed 2nd-order
  python3 scripts/richardson_extrapolation.py --spacings 0.02,0.01 --values 1.0032,1.0008 --order 2.0 --json
  
  # GCI for 3-mesh verification
  python3 scripts/gci_calculator.py --spacings 0.04,0.02,0.01 --values 1.0128,1.0032,1.0008 --json
  ```
  
  ## Error Handling
  
  | Error | Cause | Resolution |
  |-------|-------|------------|
  | `spacings and values must have the same length` | Mismatched input arrays | Provide equal-length lists |
  | `At least 2 refinement levels required` | Too few data points | Add more refinement levels |
  | `Exactly 3 refinement levels required` | GCI needs 3 levels | Provide fine/medium/coarse |
  | `Oscillatory convergence detected` | Non-monotone convergence | Check mesh quality or scheme |
  
  ## Interpretation Guidance
  
  | Scenario | Meaning | Action |
  |----------|---------|--------|
  | Observed order matches expected | Solution in asymptotic range | Report GCI, extrapolate |
  | Observed order < expected | Pre-asymptotic or coding bug | Refine further or debug |
  | Negative observed order | Solution diverging | Check implementation |
  | GCI asymptotic ratio near 1.0 | Grids in asymptotic range | Results are reliable |
  | GCI asymptotic ratio far from 1.0 | Not in asymptotic range | Refine further |
  
  ## Security
  
  ### Input Validation
  - All numeric parameters (`spacings`, `timesteps`, `values`, `expected-order`, `order`) are validated as finite positive numbers
  - Comma-separated value lists are length-matched (spacings and values must have equal length) and capped at 10,000 entries
  - GCI calculator enforces exactly 3 refinement levels; Richardson extrapolation requires at least 2
  - Safety factor is validated as a finite number greater than 1.0
  
  ### File Access
  - Scripts read no external files; all inputs are provided via CLI arguments
  - Scripts write only to stdout (JSON output); no files are created unless the agent explicitly uses the Write tool
  
  ### Tool Restrictions
  - **Bash**: Used to execute the four Python analysis scripts (`h_refinement.py`, `dt_refinement.py`, `richardson_extrapolation.py`, `gci_calculator.py`) with explicit argument lists
  - **Read**: Used to inspect script source and reference documentation
  
  ### Safety Measures
  - No `eval()`, `exec()`, or dynamic code generation
  - All subprocess calls use explicit argument lists (no `shell=True`)
  - Scripts use only Python standard library (`math` module); no pickle loading or deserialization of untrusted data
  - Minimal tool surface (Bash and Read only) limits the agent's ability to modify the filesystem
  
  ## References
  
  - `references/convergence_theory.md` - Formal convergence order, log-log analysis, asymptotic range
  - `references/gci_guidelines.md` - Roache's GCI method, ASME V&V 20, safety factors
  

• skills/core-numerical/differentiation-schemes/SKILL.mdskill
  Show content (7597 bytes)

  ---
  name: differentiation-schemes
  description: >
    Select and apply numerical differentiation schemes for PDE and ODE
    discretization — generate finite-difference stencils at arbitrary order and
    accuracy, choose between central, upwind, compact (Pade), and spectral
    methods, handle boundary stencils, and estimate truncation error scaling.
    Use when discretizing spatial derivatives, picking a scheme for advection-
    or diffusion-dominated problems, building custom stencils for nonstandard
    operators, or comparing dispersion and dissipation properties of candidate
    schemes, even if the user just says "how do I approximate this derivative"
    or "my solution is too diffusive."
  allowed-tools: Read, Bash, Write, Grep, Glob
  metadata:
    author: HeshamFS
    version: "1.1.0"
    security_tier: high
    security_reviewed: true
    tested_with:
      - claude-code
      - gemini-cli
      - vs-code-copilot
    eval_cases: 2
    last_reviewed: "2026-03-26"
  ---
  
  # Differentiation Schemes
  
  ## Goal
  
  Provide a reliable workflow to select a differentiation scheme, generate stencils, and assess accuracy for simulation discretization.
  
  ## Requirements
  
  - Python 3.8+
  - NumPy (for stencil computations)
  - No heavy dependencies
  
  ## Inputs to Gather
  
  | Input | Description | Example |
  |-------|-------------|---------|
  | Derivative order | First, second, etc. | `1` or `2` |
  | Target accuracy | Order of truncation error | `2` or `4` |
  | Grid type | Uniform, nonuniform | `uniform` |
  | Boundary type | Periodic, Dirichlet, Neumann | `periodic` |
  | Smoothness | Smooth or discontinuous | `smooth` |
  
  ## Decision Guidance
  
  ### Scheme Selection Flowchart
  
  ```
  Is the field smooth?
  ├── YES → Is domain periodic?
  │   ├── YES → Use central differences or spectral
  │   └── NO → Use central interior + one-sided at boundaries
  └── NO → Are there shocks/discontinuities?
      ├── YES → Use upwind, TVD, or WENO
      └── NO → Use central with limiters
  ```
  
  ### Quick Reference
  
  | Situation | Recommended Scheme |
  |-----------|-------------------|
  | Smooth, periodic | Central, spectral |
  | Smooth, bounded | Central + one-sided BCs |
  | Advection-dominated | Upwind |
  | Shocks/fronts | TVD, WENO |
  | High accuracy needed | Compact (Padé), spectral |
  
  ## Script Outputs (JSON Fields)
  
  | Script | Key Outputs |
  |--------|-------------|
  | `scripts/stencil_generator.py` | `offsets`, `coefficients`, `order`, `accuracy` |
  | `scripts/scheme_selector.py` | `recommended`, `alternatives`, `notes` |
  | `scripts/truncation_error.py` | `error_scale`, `order`, `notes` |
  
  ## Workflow
  
  1. **Identify requirements** - derivative order, accuracy, smoothness
  2. **Select scheme** - Run `scripts/scheme_selector.py`
  3. **Generate stencils** - Run `scripts/stencil_generator.py`
  4. **Estimate error** - Run `scripts/truncation_error.py`
  5. **Validate** - Test with manufactured solutions or grid refinement
  
  ## Conversational Workflow Example
  
  **User**: I need to discretize a second derivative for a diffusion equation on a uniform grid. I want 4th-order accuracy.
  
  **Agent workflow**:
  1. Select appropriate scheme:
     ```bash
     python3 scripts/scheme_selector.py --smooth --periodic --order 2 --accuracy 4 --json
     ```
  2. Generate the stencil:
     ```bash
     python3 scripts/stencil_generator.py --order 2 --accuracy 4 --scheme central --json
     ```
  3. Result: 5-point stencil with coefficients `[-1/12, 4/3, -5/2, 4/3, -1/12]` / dx².
  
  ## Pre-Discretization Checklist
  
  - [ ] Confirm derivative order and target accuracy
  - [ ] Choose scheme appropriate to smoothness and boundaries
  - [ ] Generate and inspect stencils at boundaries
  - [ ] Estimate truncation error vs physics scales
  - [ ] Verify with grid refinement study
  
  ## CLI Examples
  
  ```bash
  # Select scheme for smooth periodic problem
  python3 scripts/scheme_selector.py --smooth --periodic --order 1 --accuracy 4 --json
  
  # Generate central difference stencil for first derivative
  python3 scripts/stencil_generator.py --order 1 --accuracy 2 --scheme central --json
  
  # Generate 4th-order second derivative stencil
  python3 scripts/stencil_generator.py --order 2 --accuracy 4 --scheme central --json
  
  # Estimate truncation error
  python3 scripts/truncation_error.py --dx 0.01 --order 2 --accuracy 2 --scale 1.0 --json
  ```
  
  ## Error Handling
  
  | Error | Cause | Resolution |
  |-------|-------|------------|
  | `order must be positive` | Invalid derivative order | Use 1, 2, 3, ... |
  | `accuracy must be even for central` | Odd accuracy requested | Use 2, 4, 6, ... |
  | `Unknown scheme` | Invalid scheme type | Use central, upwind, compact |
  
  ## Interpretation Guidance
  
  ### Stencil Properties
  
  | Property | Meaning |
  |----------|---------|
  | Symmetric offsets | Central scheme (no directional bias) |
  | Asymmetric offsets | One-sided or upwind scheme |
  | More points | Higher accuracy but wider stencil |
  
  ### Truncation Error Scaling
  
  | Accuracy Order | Error Scales As | Refinement Factor |
  |----------------|-----------------|-------------------|
  | 2nd order | O(dx²) | 2× refinement → 4× error reduction |
  | 4th order | O(dx⁴) | 2× refinement → 16× error reduction |
  | 6th order | O(dx⁶) | 2× refinement → 64× error reduction |
  
  ### Common Stencils
  
  | Derivative | Accuracy | Points | Coefficients (× 1/dx or 1/dx²) |
  |------------|----------|--------|-------------------------------|
  | 1st | 2 | 3 | [-1/2, 0, 1/2] |
  | 1st | 4 | 5 | [1/12, -2/3, 0, 2/3, -1/12] |
  | 2nd | 2 | 3 | [1, -2, 1] |
  | 2nd | 4 | 5 | [-1/12, 4/3, -5/2, 4/3, -1/12] |
  
  ## Security
  
  ### Input Validation
  - `--order` (derivative order) is validated as a positive integer with an upper bound
  - `--accuracy` is validated as a positive even integer for central schemes
  - `--scheme` is validated against a fixed allowlist (`central`, `upwind`, `compact`)
  - `--dx` and `--scale` are validated as finite positive numbers
  - No user-supplied strings are interpolated into code paths or shell commands
  
  ### File Access
  - Scripts read no external files; all inputs are provided via CLI arguments
  - Scripts write only to stdout (JSON output); no files are created unless the agent explicitly uses the Write tool
  
  ### Tool Restrictions
  - **Read**: Used to inspect script source, references, and user configuration files
  - **Bash**: Used to execute the three Python scripts (`stencil_generator.py`, `scheme_selector.py`, `truncation_error.py`) with explicit argument lists
  - **Write**: Used to save generated stencil coefficients or scheme recommendations; writes are scoped to the user's working directory
  - **Grep/Glob**: Used to locate relevant files and search references
  
  ### Safety Measures
  - No `eval()`, `exec()`, or dynamic code generation
  - All subprocess calls use explicit argument lists (no `shell=True`)
  - Stencil computation uses only NumPy linear algebra on small, bounded matrices (stencil width limited by accuracy order)
  - All output is deterministic JSON with no shell-interpretable content
  
  ## Limitations
  
  - **Boundary handling**: Stencil generator provides interior stencils; boundaries need special treatment
  - **Nonuniform grids**: Standard stencils assume uniform spacing
  - **Spectral**: Not covered by stencil generator
  
  ## References
  
  - `references/stencil_catalog.md` - Common stencils
  - `references/boundary_handling.md` - One-sided schemes
  - `references/scheme_selection.md` - FD/FV/spectral comparison
  - `references/error_guidance.md` - Truncation error scaling
  
  ## Version History
  
  - **v1.1.0** (2024-12-24): Enhanced documentation, decision guidance, examples
  - **v1.0.0**: Initial release with 3 differentiation scripts
  

• skills/core-numerical/linear-solvers/SKILL.mdskill
  Show content (8240 bytes)

  ---
  name: linear-solvers
  description: >
    Select and configure linear solvers for Ax=b systems arising in numerical
    simulations — choose between direct (LU, Cholesky) and iterative (CG, GMRES,
    BiCGSTAB, MINRES) methods, analyze sparsity patterns and matrix conditioning,
    recommend preconditioners (AMG, ILU, IC), apply row/column scaling, and
    diagnose convergence stagnation from residual histories. Use when setting up
    a linear solve for FEM/FVM assembly, debugging slow or stalled Krylov
    iterations, choosing a preconditioner for SPD or nonsymmetric systems, or
    investigating ill-conditioning, even if the user only says "my solver is
    slow" or "GMRES won't converge."
  allowed-tools: Read, Write, Grep, Glob
  metadata:
    author: HeshamFS
    version: "1.1.0"
    security_tier: medium
    security_reviewed: true
    tested_with:
      - claude-code
      - gemini-cli
      - vs-code-copilot
    eval_cases: 2
    last_reviewed: "2026-03-26"
  ---
  
  # Linear Solvers
  
  ## Goal
  
  Provide a universal workflow to select a solver, assess conditioning, and diagnose convergence for linear systems arising in numerical simulations.
  
  ## Requirements
  
  - Python 3.8+
  - NumPy, SciPy (for matrix operations)
  - See individual scripts for dependencies
  
  ## Inputs to Gather
  
  | Input | Description | Example |
  |-------|-------------|---------|
  | Matrix size | Dimension of system | `n = 1000000` |
  | Sparsity | Fraction of nonzeros | `0.01%` |
  | Symmetry | Is A = Aᵀ? | `yes` |
  | Definiteness | Is A positive definite? | `yes (SPD)` |
  | Conditioning | Estimated condition number | `10⁶` |
  
  ## Decision Guidance
  
  ### Solver Selection Flowchart
  
  ```
  Is matrix small (n < 5000) and dense?
  ├── YES → Use direct solver (LU, Cholesky)
  └── NO → Is matrix symmetric?
      ├── YES → Is it positive definite?
      │   ├── YES → Use CG with AMG/IC preconditioner
      │   └── NO → Use MINRES
      └── NO → Is it nearly symmetric?
          ├── YES → Use BiCGSTAB
          └── NO → Use GMRES with ILU/AMG
  ```
  
  ### Quick Reference
  
  | Matrix Type | Solver | Preconditioner |
  |-------------|--------|----------------|
  | SPD, sparse | CG | AMG, IC |
  | Symmetric indefinite | MINRES | ILU |
  | Nonsymmetric | GMRES, BiCGSTAB | ILU, AMG |
  | Dense | LU, Cholesky | None |
  | Saddle point | Schur complement, Uzawa | Block preconditioner |
  
  ## Script Outputs (JSON Fields)
  
  | Script | Key Outputs |
  |--------|-------------|
  | `scripts/solver_selector.py` | `recommended`, `alternatives`, `notes` |
  | `scripts/convergence_diagnostics.py` | `rate`, `stagnation`, `recommended_action` |
  | `scripts/sparsity_stats.py` | `nnz`, `density`, `bandwidth`, `symmetry` |
  | `scripts/preconditioner_advisor.py` | `suggested`, `notes` |
  | `scripts/scaling_equilibration.py` | `row_scale`, `col_scale`, `notes` |
  | `scripts/residual_norms.py` | `residual_norms`, `relative_norms`, `converged` |
  
  ## Workflow
  
  1. **Characterize matrix** - symmetry, definiteness, sparsity
  2. **Analyze sparsity** - Run `scripts/sparsity_stats.py`
  3. **Select solver** - Run `scripts/solver_selector.py`
  4. **Choose preconditioner** - Run `scripts/preconditioner_advisor.py`
  5. **Apply scaling** - If ill-conditioned, use `scripts/scaling_equilibration.py`
  6. **Monitor convergence** - Use `scripts/convergence_diagnostics.py`
  7. **Diagnose issues** - Check residual history with `scripts/residual_norms.py`
  
  ## Conversational Workflow Example
  
  **User**: My GMRES solver is stagnating after 50 iterations. The residual drops to 1e-3 then stops improving.
  
  **Agent workflow**:
  1. Diagnose convergence:
     ```bash
     python3 scripts/convergence_diagnostics.py --residuals 1,0.1,0.01,0.005,0.003,0.002,0.002,0.002 --json
     ```
  2. Check for preconditioning advice:
     ```bash
     python3 scripts/preconditioner_advisor.py --matrix-type nonsymmetric --sparse --stagnation --json
     ```
  3. Recommend: Increase restart parameter, try ILU(k) with higher k, or switch to AMG.
  
  ## Pre-Solve Checklist
  
  - [ ] Confirm matrix symmetry/definiteness
  - [ ] Decide direct vs iterative based on size and sparsity
  - [ ] Set residual tolerance relative to physics scale
  - [ ] Choose preconditioner appropriate to matrix structure
  - [ ] Apply scaling/equilibration if needed
  - [ ] Track convergence and adjust if stagnation occurs
  
  ## CLI Examples
  
  ```bash
  # Analyze sparsity pattern
  python3 scripts/sparsity_stats.py --matrix A.npy --json
  
  # Select solver for SPD sparse system
  python3 scripts/solver_selector.py --symmetric --positive-definite --sparse --size 1000000 --json
  
  # Get preconditioner recommendation
  python3 scripts/preconditioner_advisor.py --matrix-type spd --sparse --json
  
  # Diagnose convergence from residual history
  python3 scripts/convergence_diagnostics.py --residuals 1,0.2,0.05,0.01 --json
  
  # Apply scaling
  python3 scripts/scaling_equilibration.py --matrix A.npy --symmetric --json
  
  # Compute residual norms
  python3 scripts/residual_norms.py --residual 1,0.1,0.01 --rhs 1,0,0 --json
  ```
  
  ## Error Handling
  
  | Error | Cause | Resolution |
  |-------|-------|------------|
  | `Matrix file not found` | Invalid path | Check file exists |
  | `Matrix must be square` | Non-square input | Verify matrix dimensions |
  | `Residuals must be positive` | Invalid residual data | Check input format |
  
  ## Interpretation Guidance
  
  ### Convergence Rate
  
  | Rate | Meaning | Action |
  |------|---------|--------|
  | < 0.1 | Excellent | Current setup optimal |
  | 0.1 - 0.5 | Good | Acceptable for most problems |
  | 0.5 - 0.9 | Slow | Consider better preconditioner |
  | > 0.9 | Stagnation | Change solver or preconditioner |
  
  ### Stagnation Diagnosis
  
  | Pattern | Likely Cause | Fix |
  |---------|--------------|-----|
  | Flat residual | Poor preconditioner | Improve preconditioner |
  | Oscillating | Near-singular or indefinite | Check matrix, try different solver |
  | Very slow decay | Ill-conditioned | Apply scaling, use AMG |
  
  ## Security
  
  ### Input Validation
  - All numeric inputs (residuals, tolerances, matrix entries) are validated as finite numbers
  - Comma-separated residual/vector inputs are capped at 100,000 entries
  - The `solver_selector.py` `--size` parameter is bounded at 10 billion
  - `--matrix-type` is validated against a fixed allowlist (`spd`, `symmetric`, `nonsymmetric`)
  - Boolean flags (`--symmetric`, `--positive-definite`, `--sparse`, `--stagnation`) are type-safe argparse flags
  
  ### File Access
  - `sparsity_stats.py` and `scaling_equilibration.py` read a single matrix file (`.npy` format) specified by `--matrix`
  - `np.load()` is called with `allow_pickle=False` to prevent arbitrary code execution via crafted `.npy` files
  - Matrix files are rejected if they exceed 500 MB before any parsing occurs
  - Matrix dimension limits (100,000 per dimension) prevent memory exhaustion
  - All other scripts read no external files; inputs are provided via CLI arguments
  
  ### Tool Restrictions
  - **Read**: Used to inspect script source, references, and matrix files
  - **Write**: Used to save analysis results or solver recommendations; writes are scoped to the user's working directory
  - **Grep/Glob**: Used to locate relevant files and search references
  - The skill's `allowed-tools` excludes `Bash` to prevent the agent from executing arbitrary commands when processing untrusted matrix files or numeric inputs
  
  ### Safety Measures
  - No `eval()`, `exec()`, or dynamic code generation
  - All subprocess calls use explicit argument lists (no `shell=True`)
  - Reduced tool surface (no Bash) limits the agent to read/write operations only
  - JSON output mode produces structured, parseable results without shell-interpretable content
  
  ## Limitations
  
  - **Large dense matrices**: Direct solvers may run out of memory
  - **Highly indefinite**: Standard preconditioners may fail
  - **Saddle-point**: Requires specialized block preconditioners
  
  ## References
  
  - `references/solver_decision_tree.md` - Selection logic
  - `references/preconditioner_catalog.md` - Preconditioner options
  - `references/convergence_patterns.md` - Diagnosing failures
  - `references/scaling_guidelines.md` - Equilibration guidance
  
  ## Version History
  
  - **v1.1.0** (2024-12-24): Enhanced documentation, decision guidance, examples
  - **v1.0.0**: Initial release with 6 solver analysis scripts
  

• skills/core-numerical/mesh-generation/SKILL.mdskill
  Show content (6612 bytes)

  ---
  name: mesh-generation
  description: >
    Plan and evaluate mesh generation for numerical simulations — estimate grid
    resolution from physics scales (interface width, boundary layers, wavelengths),
    check aspect ratios and skewness against quality thresholds, choose between
    structured, unstructured, and adaptive mesh refinement strategies, and compute
    grid sizing for 1D/2D/3D domains. Use when setting up a new mesh, diagnosing
    poor solver convergence caused by mesh quality, deciding how many points to
    place across a phase-field interface or boundary layer, or preparing a mesh
    convergence study, even if the user only asks "what resolution do I need"
    or "why is my solver failing."
  allowed-tools: Read, Write, Grep, Glob
  metadata:
    author: HeshamFS
    version: "1.1.0"
    security_tier: medium
    security_reviewed: true
    tested_with:
      - claude-code
      - gemini-cli
      - vs-code-copilot
    eval_cases: 2
    last_reviewed: "2026-03-26"
  ---
  
  # Mesh Generation
  
  ## Goal
  
  Provide a consistent workflow for selecting mesh resolution and checking mesh quality for PDE simulations.
  
  ## Requirements
  
  - Python 3.8+
  - No external dependencies (uses stdlib)
  
  ## Inputs to Gather
  
  | Input | Description | Example |
  |-------|-------------|---------|
  | Domain size | Physical dimensions | `1.0 × 1.0 m` |
  | Feature size | Smallest feature to resolve | `0.01 m` |
  | Points per feature | Resolution requirement | `10 points` |
  | Aspect ratio limit | Maximum dx/dy ratio | `5:1` |
  | Quality threshold | Skewness limit | `< 0.8` |
  
  ## Decision Guidance
  
  ### Resolution Selection
  
  ```
  What is the smallest feature size?
  ├── Interface width → dx ≤ width / 5
  ├── Boundary layer → dx ≤ layer_thickness / 10
  ├── Wave length → dx ≤ lambda / 20
  └── Diffusion length → dx ≤ sqrt(D × dt) / 2
  ```
  
  ### Mesh Type Selection
  
  | Problem | Recommended Mesh |
  |---------|------------------|
  | Simple geometry, uniform | Structured Cartesian |
  | Complex geometry | Unstructured triangular/tetrahedral |
  | Boundary layers | Hybrid (structured near walls) |
  | Adaptive refinement | Quadtree/Octree or AMR |
  
  ## Script Outputs (JSON Fields)
  
  | Script | Key Outputs |
  |--------|-------------|
  | `scripts/grid_sizing.py` | `dx`, `nx`, `ny`, `nz`, `notes` |
  | `scripts/mesh_quality.py` | `aspect_ratio`, `skewness`, `quality_flags` |
  
  ## Workflow
  
  1. **Estimate resolution** - From physics scales
  2. **Compute grid sizing** - Run `scripts/grid_sizing.py`
  3. **Check quality metrics** - Run `scripts/mesh_quality.py`
  4. **Adjust if needed** - Fix aspect ratios, reduce skewness
  5. **Validate** - Mesh convergence study
  
  ## Conversational Workflow Example
  
  **User**: I need to mesh a 1mm × 1mm domain for a phase-field simulation with interface width of 10 μm.
  
  **Agent workflow**:
  1. Compute grid sizing:
     ```bash
     python3 scripts/grid_sizing.py --length 0.001 --resolution 200 --json
     ```
  2. Verify interface is resolved: dx = 5 μm, interface width = 10 μm → 2 points per interface width.
  3. Recommend: Increase to 500 points (dx = 2 μm) for 5 points across interface.
  
  ## Pre-Mesh Checklist
  
  - [ ] Define target resolution per feature/interface
  - [ ] Ensure dx meets stability constraints (see numerical-stability)
  - [ ] Check aspect ratio < limit (typically 5:1)
  - [ ] Check skewness < threshold (typically 0.8)
  - [ ] Validate mesh convergence with refinement study
  
  ## CLI Examples
  
  ```bash
  # Compute grid sizing for 1D domain
  python3 scripts/grid_sizing.py --length 1.0 --resolution 200 --json
  
  # Check mesh quality
  python3 scripts/mesh_quality.py --dx 1.0 --dy 0.5 --dz 0.5 --json
  
  # High aspect ratio check
  python3 scripts/mesh_quality.py --dx 1.0 --dy 0.1 --json
  ```
  
  ## Error Handling
  
  | Error | Cause | Resolution |
  |-------|-------|------------|
  | `length must be positive` | Invalid domain size | Use positive value |
  | `resolution must be > 1` | Insufficient points | Use at least 2 |
  | `dx, dy must be positive` | Invalid spacing | Use positive values |
  
  ## Interpretation Guidance
  
  ### Aspect Ratio
  
  | Aspect Ratio | Quality | Impact |
  |--------------|---------|--------|
  | 1:1 | Excellent | Optimal accuracy |
  | 1:1 - 3:1 | Good | Acceptable |
  | 3:1 - 5:1 | Fair | May affect accuracy |
  | > 5:1 | Poor | Solver issues likely |
  
  ### Skewness
  
  | Skewness | Quality | Impact |
  |----------|---------|--------|
  | 0 - 0.25 | Excellent | Optimal |
  | 0.25 - 0.50 | Good | Acceptable |
  | 0.50 - 0.80 | Fair | May affect accuracy |
  | > 0.80 | Poor | Likely problems |
  
  ### Resolution Guidelines
  
  | Application | Points per Feature |
  |-------------|-------------------|
  | Phase-field interface | 5-10 |
  | Boundary layer | 10-20 |
  | Shock | 3-5 (with capturing) |
  | Wave propagation | 10-20 per wavelength |
  | Smooth gradients | 5-10 |
  
  ## Security
  
  ### Input Validation
  - All inputs (`length`, `resolution`, `dx`, `dy`, `dz`) are validated as finite positive numbers with upper bounds to prevent resource exhaustion
  - `dims` is restricted to `{1, 2, 3}`
  - `argparse` type parameters reject non-numeric input at the CLI boundary before any processing occurs
  
  ### File Access
  - Scripts read no external files; all inputs are provided via CLI arguments
  - Scripts write only to stdout (JSON output); no files are created unless the agent explicitly uses the Write tool
  
  ### Tool Restrictions
  - **Read**: Used to inspect script source, references, and user configuration files
  - **Write**: Used to save grid sizing results or mesh quality reports; writes are scoped to the user's working directory
  - **Grep/Glob**: Used to locate relevant files and search references
  - The skill's `allowed-tools` excludes `Bash` to prevent the agent from executing arbitrary commands when processing user-provided inputs
  
  ### Safety Measures
  - No `eval()`, `exec()`, or dynamic code generation
  - All subprocess calls use explicit argument lists (no `shell=True`)
  - Reduced tool surface (no Bash) means the agent should use `Read` and `Write` to prepare inputs and capture outputs rather than constructing shell commands from user text
  - All output is deterministic JSON with no shell-interpretable content
  
  ## Limitations
  
  - **2D/3D only**: No unstructured mesh generation
  - **Quality metrics**: Basic aspect ratio and skewness only
  - **No mesh generation**: Sizing recommendations only
  
  ## References
  
  - `references/mesh_types.md` - Structured vs unstructured
  - `references/quality_metrics.md` - Aspect ratio/skewness thresholds
  
  ## Version History
  
  - **v1.1.0** (2024-12-24): Enhanced documentation, decision guidance, examples
  - **v1.0.0**: Initial release with 2 mesh quality scripts
  

• skills/core-numerical/nonlinear-solvers/SKILL.mdskill
  Show content (9784 bytes)

  ---
  name: nonlinear-solvers
  description: >
    Select and configure nonlinear solvers for root-finding f(x)=0, optimization
    min F(x), and least-squares problems — choose among Newton, Newton-Krylov,
    quasi-Newton (BFGS, L-BFGS), Broyden, Anderson acceleration, and
    Levenberg-Marquardt methods, configure line search or trust-region
    globalization, diagnose convergence rate (quadratic, linear, stagnated),
    and assess Jacobian quality and conditioning. Use when a Newton solver
    converges slowly or diverges, choosing between line search and trust region,
    debugging nonlinear iteration failures in FEM or phase-field codes, or
    selecting a solver for large-scale unconstrained optimization, even if
    the user only says "my Newton iterations aren't converging."
  allowed-tools: Read, Bash, Write, Grep, Glob
  metadata:
    author: HeshamFS
    version: "1.1.0"
    security_tier: high
    security_reviewed: true
    tested_with:
      - claude-code
      - gemini-cli
      - vs-code-copilot
    eval_cases: 2
    last_reviewed: "2026-03-26"
  ---
  
  # 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.8+
  - NumPy (for Jacobian diagnostics)
  - SciPy (optional, for advanced analysis)
  
  ## Inputs to Gather
  
  | Input | Description | Example |
  |-------|-------------|---------|
  | Problem type | Root-finding, optimization, least-squares | `root-finding` |
  | Problem size | Number of unknowns | `n = 10000` |
  | Jacobian availability | Analytic, finite-diff, unavailable | `analytic` |
  | Jacobian cost | Cheap or expensive to compute | `expensive` |
  | Constraints | None, bounds, equality, inequality | `none` |
  | Smoothness | Is objective/residual smooth? | `yes` |
  | Residual history | Sequence of residual norms | `1,0.1,0.01,...` |
  
  ## 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
  
  | Problem Type | First Choice | Alternative | Globalization |
  |--------------|--------------|-------------|---------------|
  | Small root-finding | Newton | Broyden | Line search |
  | Large root-finding | Newton-Krylov | Anderson | Trust region |
  | Optimization | L-BFGS | BFGS | Wolfe line search |
  | Least-squares | Levenberg-Marquardt | Gauss-Newton | Trust region |
  | Bound constrained | L-BFGS-B | Trust-region reflective | Projected |
  
  ## Script Outputs (JSON Fields)
  
  | Script | Key Outputs |
  |--------|-------------|
  | `scripts/solver_selector.py` | `recommended`, `alternatives`, `notes` |
  | `scripts/convergence_analyzer.py` | `converged`, `convergence_type`, `estimated_rate`, `diagnosis` |
  | `scripts/jacobian_diagnostics.py` | `condition_number`, `jacobian_quality`, `rank_deficient` |
  | `scripts/globalization_advisor.py` | `strategy`, `line_search_type`, `trust_region_type`, `parameters` |
  | `scripts/residual_monitor.py` | `patterns_detected`, `alerts`, `recommendations` |
  | `scripts/step_quality.py` | `ratio`, `step_quality`, `accept_step`, `trust_radius_action` |
  
  ## 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
     ```
  2. Check globalization strategy:
     ```bash
     python3 scripts/globalization_advisor.py --problem-type root-finding --jacobian-quality ill-conditioned --previous-failures 0 --json
     ```
  3. 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
  
  ```bash
  # Select solver for large unconstrained optimization
  python3 scripts/solver_selector.py --size 50000 --smooth --memory-limited --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
  
  # 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
  
  | Error | Cause | Resolution |
  |-------|-------|------------|
  | `problem_size must be positive` | Invalid size | Check problem dimension |
  | `constraint_type must be one of...` | Unknown constraint | Use: none, bound, equality, inequality |
  | `residuals must be non-negative` | Invalid residual data | Check residual computation |
  | `Matrix file not found` | Invalid path | Verify Jacobian file exists |
  
  ## Interpretation Guidance
  
  ### Convergence Type
  
  | Type | Meaning | Action |
  |------|---------|--------|
  | quadratic | Optimal Newton | Continue, near solution |
  | superlinear | Quasi-Newton working | Monitor for stagnation |
  | linear | Acceptable | May improve with preconditioner |
  | sublinear | Too slow | Change method or formulation |
  | stagnated | No progress | Check Jacobian, preconditioner |
  | diverged | Increasing residual | Add globalization, check Jacobian |
  
  ### Jacobian Quality
  
  | Quality | Condition Number | Action |
  |---------|------------------|--------|
  | good | < 10⁶ | Standard Newton works |
  | moderately-conditioned | 10⁶ - 10¹⁰ | Consider scaling |
  | ill-conditioned | > 10¹⁰ | Use regularization |
  | near-singular | ∞ | Reformulate or use LM |
  
  ### Step Quality (Trust Region)
  
  | Ratio ρ | Quality | Trust Radius |
  |---------|---------|--------------|
  | ρ < 0 | very_poor | Shrink aggressively |
  | ρ < 0.25 | marginal | Shrink |
  | 0.25 ≤ ρ < 0.75 | good | Maintain |
  | ρ ≥ 0.75 | excellent | Expand if at boundary |
  
  ## Security
  
  ### Input Validation
  - `--size` (problem size) is validated as a positive integer, bounded at 10 billion
  - `--residuals` are validated as finite non-negative numbers, capped at 100,000 entries
  - `--tolerance` and `--target-tolerance` are validated as finite positive numbers
  - `--problem-type` and `--constraint-type` are validated against fixed allowlists
  - `--jacobian-quality` is validated against a fixed allowlist (`good`, `ill-conditioned`, etc.)
  - Step quality parameters (`predicted-reduction`, `actual-reduction`, `step-norm`, `gradient-norm`, `trust-radius`) are validated as finite numbers
  
  ### File Access
  - `jacobian_diagnostics.py` reads a single matrix file specified by `--matrix`; no directory traversal beyond the given path
  - Matrix files are size-limited and loaded with `allow_pickle=False` to prevent code execution
  - All other scripts read no external files; inputs are provided via CLI arguments
  - Scripts write only to stdout (JSON output)
  
  ### Tool Restrictions
  - **Read**: Used to inspect script source, references, and user configuration files
  - **Bash**: Used to execute the six Python analysis scripts (`solver_selector.py`, `convergence_analyzer.py`, `jacobian_diagnostics.py`, `globalization_advisor.py`, `residual_monitor.py`, `step_quality.py`) with explicit argument lists
  - **Write**: Used to save analysis results or solver recommendations; writes are scoped to the user's working directory
  - **Grep/Glob**: Used to locate relevant files and search references
  
  ### Safety Measures
  - No `eval()`, `exec()`, or dynamic code generation
  - All subprocess calls use explicit argument lists (no `shell=True`)
  - Matrix dimension limits prevent memory exhaustion when loading Jacobian files
  - Residual history analysis operates on bounded-length numeric arrays only
  
  ## Limitations
  
  - **No global convergence guarantee**: All methods may fail for pathological problems
  - **Jacobian accuracy**: Finite-difference Jacobian may be inaccurate near discontinuities
  - **Large dense problems**: May require specialized solvers not covered here
  - **Constrained optimization**: Complex constraints need SQP or interior point methods
  
  ## References
  
  - `references/solver_decision_tree.md` - Problem-based solver selection
  - `references/method_catalog.md` - Method details and parameters
  - `references/convergence_diagnostics.md` - Diagnosing convergence issues
  - `references/globalization_strategies.md` - Line search and trust region
  
  ## Version History
  
  - **v1.0.0** : Initial release with 6 analysis scripts
  

• skills/core-numerical/numerical-stability/SKILL.mdskill
  Show content (8094 bytes)

  ---
  name: numerical-stability
  description: >
    Analyze numerical stability for time-dependent PDE simulations — check CFL
    and Fourier criteria, perform von Neumann stability analysis, detect stiffness,
    evaluate matrix conditioning, and recommend explicit vs implicit time-stepping
    schemes. Use when selecting time steps, diagnosing numerical blow-up or solver
    divergence, checking convergence criteria, or evaluating scheme stability for
    advection, diffusion, or reaction problems, even if the user doesn't explicitly
    mention "stability" or "CFL."
  allowed-tools: Read, Bash, Write, Grep, Glob
  metadata:
    author: HeshamFS
    version: "1.1.0"
    security_tier: high
    security_reviewed: true
    tested_with:
      - claude-code
      - gemini-cli
      - vs-code-copilot
    eval_cases: 2
    last_reviewed: "2026-03-26"
  ---
  
  # Numerical Stability
  
  ## Goal
  
  Provide a repeatable checklist and script-driven checks to keep time-dependent simulations stable and defensible.
  
  ## Requirements
  
  - Python 3.8+
  - NumPy (for matrix_condition.py and von_neumann_analyzer.py)
  - See `scripts/requirements.txt` for dependencies
  
  ## Inputs to Gather
  
  | Input | Description | Example |
  |-------|-------------|---------|
  | Grid spacing `dx` | Spatial discretization | `0.01 m` |
  | Time step `dt` | Temporal discretization | `1e-4 s` |
  | Velocity `v` | Advection speed | `1.0 m/s` |
  | Diffusivity `D` | Thermal/mass diffusivity | `1e-5 m²/s` |
  | Reaction rate `k` | First-order rate constant | `100 s⁻¹` |
  | Dimensions | 1D, 2D, or 3D | `2` |
  | Scheme type | Explicit or implicit | `explicit` |
  
  ## Decision Guidance
  
  ### Choosing Explicit vs Implicit
  
  ```
  Is the problem stiff (fast + slow dynamics)?
  ├── YES → Use implicit or IMEX scheme
  │         └── Check conditioning with matrix_condition.py
  └── NO → Is CFL/Fourier satisfied with reasonable dt?
      ├── YES → Use explicit scheme (cheaper per step)
      └── NO → Consider implicit or reduce dx
  ```
  
  ### Stability Limit Quick Reference
  
  | Physics | Number | Explicit Limit (1D) | Formula |
  |---------|--------|---------------------|---------|
  | Advection | CFL | C ≤ 1 | `C = v·dt/dx` |
  | Diffusion | Fourier | Fo ≤ 0.5 | `Fo = D·dt/dx²` |
  | Reaction | Reaction | R ≤ 1 | `R = k·dt` |
  
  **Multi-dimensional correction**: For d dimensions, diffusion limit is `Fo ≤ 1/(2d)`.
  
  ## Script Outputs (JSON Fields)
  
  | Script | Key Outputs |
  |--------|-------------|
  | `scripts/cfl_checker.py` | `metrics.cfl`, `metrics.fourier`, `recommended_dt`, `stable` |
  | `scripts/von_neumann_analyzer.py` | `results.max_amplification`, `results.stable` |
  | `scripts/matrix_condition.py` | `results.condition_number`, `results.is_symmetric` |
  | `scripts/stiffness_detector.py` | `results.stiffness_ratio`, `results.stiff`, `results.recommendation` |
  
  ## Workflow
  
  1. **Identify dominant physics** (advection vs diffusion vs reaction)
  2. **Run CFL checker** with `scripts/cfl_checker.py`
  3. **Compare to limits** and adjust `dt` if needed
  4. **Check stiffness** with `scripts/stiffness_detector.py` if multiple scales
  5. **Analyze custom schemes** with `scripts/von_neumann_analyzer.py`
  6. **Check conditioning** with `scripts/matrix_condition.py` for implicit solves
  7. **Document** the stability verdict and recommended time step
  
  ## Conversational Workflow Example
  
  **User**: My phase-field simulation is blowing up after 100 steps. I'm using explicit Euler with dx=0.01, dt=1e-4, and diffusivity D=1e-3.
  
  **Agent workflow**:
  1. Check stability criteria:
     ```bash
     python3 scripts/cfl_checker.py --dx 0.01 --dt 1e-4 --diffusivity 1e-3 --dimensions 2 --json
     ```
  2. Interpret results:
     - Fourier number: `Fo = 1e-3 × 1e-4 / (0.01)² = 1.0`
     - 2D limit: `Fo ≤ 0.25`
     - **Violation**: Fo = 1.0 > 0.25, unstable!
  3. Recommend fix:
     - Reduce dt to `2.5e-5` (to get Fo = 0.25)
     - Or increase dx, or switch to implicit
  
  ## Pre-Simulation Stability Checklist
  
  - [ ] Identify dominant physics and nondimensional groups
  - [ ] Compute CFL/Fourier/Reaction numbers with `cfl_checker.py`
  - [ ] If explicit and limit violated, reduce `dt` or change scheme
  - [ ] If stiffness ratio > 1000, select implicit/stiff integrator
  - [ ] For custom schemes, verify amplification factor ≤ 1
  - [ ] Document stability reasoning with inputs and outputs
  
  ## CLI Examples
  
  ```bash
  # Check CFL/Fourier for 2D diffusion-advection
  python3 scripts/cfl_checker.py --dx 0.1 --dt 0.01 --velocity 1.0 --diffusivity 0.1 --dimensions 2 --json
  
  # Von Neumann analysis for custom 3-point stencil
  python3 scripts/von_neumann_analyzer.py --coeffs 0.2,0.6,0.2 --dx 1.0 --nk 128 --json
  
  # Detect stiffness from eigenvalue estimates
  python3 scripts/stiffness_detector.py --eigs=-1,-1000 --json
  
  # Check matrix conditioning for implicit system
  python3 scripts/matrix_condition.py --matrix A.npy --norm 2 --json
  ```
  
  ## Error Handling
  
  | Error | Cause | Resolution |
  |-------|-------|------------|
  | `dx and dt must be positive` | Zero or negative values | Provide valid positive numbers |
  | `No stability criteria applied` | Missing velocity/diffusivity | Provide at least one physics parameter |
  | `Matrix file not found` | Invalid path | Check matrix file exists |
  | `Could not compute eigenvalues` | Singular or ill-formed matrix | Check matrix validity |
  
  ## Interpretation Guidance
  
  | Scenario | Meaning | Action |
  |----------|---------|--------|
  | `stable: true` | All checked criteria satisfied | Proceed with simulation |
  | `stable: false` | At least one limit violated | Reduce dt or change scheme |
  | `stable: null` | No criteria could be applied | Provide more physics inputs |
  | Stiffness ratio > 1000 | Problem is stiff | Use implicit integrator |
  | Condition number > 10⁶ | Ill-conditioned | Use scaling/preconditioning |
  
  ## Security
  
  ### Input Validation
  - All numeric parameters (`dx`, `dt`, `velocity`, `diffusivity`, `dimensions`) are validated as finite positive numbers before any computation
  - `--dimensions` is restricted to `{1, 2, 3}`
  - Comma-separated eigenvalue lists (`--eigs`) are capped at 10,000 entries and validated as finite numbers
  - Stencil coefficient lists (`--coeffs`) are length-limited and validated as finite floats
  
  ### File Access
  - `matrix_condition.py` reads a single matrix file (`.npy` format) specified by `--matrix`; no directory traversal beyond the given path
  - Matrix files are rejected if they exceed 500 MB before parsing
  - `np.load()` is called with `allow_pickle=False` to prevent arbitrary code execution via crafted `.npy` files
  - Scripts write only to stdout (JSON output); no files are created unless the agent explicitly uses the Write tool
  
  ### Tool Restrictions
  - **Read**: Used to inspect script source, references, and user configuration files
  - **Bash**: Used to execute the four Python analysis scripts (`cfl_checker.py`, `von_neumann_analyzer.py`, `matrix_condition.py`, `stiffness_detector.py`) with explicit argument lists
  - **Write**: Used to save analysis results or generated reports; writes are scoped to the user's working directory
  - **Grep/Glob**: Used to locate relevant files and search references
  
  ### Safety Measures
  - No `eval()`, `exec()`, or dynamic code generation
  - All subprocess calls use explicit argument lists (no `shell=True`)
  - Matrix dimension limits (100,000 per dimension) prevent memory exhaustion
  - JSON output mode (`--json`) produces structured, parseable results without shell-interpretable content
  
  ## Limitations
  
  - **Explicit schemes only** for CFL/Fourier checks (implicit is unconditionally stable)
  - **Von Neumann analysis** assumes linear, constant-coefficient, periodic BCs
  - **Stiffness detection** requires eigenvalue estimates from user
  
  ## References
  
  - `references/stability_criteria.md` - Decision thresholds and formulas
  - `references/common_pitfalls.md` - Frequent failure modes and fixes
  - `references/scheme_catalog.md` - Stability properties of common schemes
  
  ## Version History
  
  - **v1.1.0** (2024-12-24): Enhanced documentation, decision guidance, examples
  - **v1.0.0**: Initial release with 4 stability analysis scripts
  

README

You: Check if dt=0.001 is stable for my advection problem with v=2.0 m/s and dx=0.01.

Agent: I'll use the numerical-stability skill to check this.
       [runs cfl_checker.py --dx 0.01 --dt 0.001 --velocity 2.0 --json]

       CFL number = 0.2 (limit: 1.0) -- stable.
       Recommended max dt = 0.005 with safety factor 1.0.

skills/core-numerical/numerical-stability/
    SKILL.md              # Instructions + YAML metadata (loaded when skill triggers)
    scripts/              # Python CLI tools (executed for reproducible results)
        cfl_checker.py
        von_neumann_analyzer.py
        matrix_condition.py
        stiffness_detector.py
    references/           # Deep domain knowledge (loaded only when needed)
        stability_criteria.md
        common_pitfalls.md
        scheme_catalog.md
    evals/                # Evaluation suite per agentskills.io spec
        evals.json        # Test cases with prompts, assertions, expected outputs
    CHANGELOG.md          # Version history

git clone https://github.com/heshamfs/materials-simulation-skills.git
cd materials-simulation-skills
pip install -r requirements-dev.txt

python -m pytest tests/ -v --tb=short          # All 932 tests
python -m pytest tests/unit -v --tb=short       # Unit tests only
python -m pytest tests/integration -v           # Integration tests only

# Personal (available across all projects)
cp -r skills/core-numerical/numerical-stability ~/.claude/skills/numerical-stability

# Project-level (committed to version control)
cp -r skills/core-numerical/numerical-stability .claude/skills/numerical-stability

claude --add-dir /path/to/materials-simulation-skills/skills

# User-scoped (available across all workspaces)
cp -r skills/core-numerical/numerical-stability ~/.gemini/skills/numerical-stability

# Workspace-scoped (project-specific)
cp -r skills/core-numerical/numerical-stability .gemini/skills/numerical-stability

# User-scoped
cp -r skills/core-numerical/numerical-stability ~/.agents/skills/numerical-stability

# Repository-scoped
cp -r skills/core-numerical/numerical-stability .agents/skills/numerical-stability

# Workspace (committed to version control)
cp -r skills/core-numerical/numerical-stability .github/skills/numerical-stability

# Personal (across all workspaces)
cp -r skills/core-numerical/numerical-stability ~/.copilot/skills/numerical-stability

cp -r skills/core-numerical/numerical-stability skills/numerical-stability

Use numerical-stability to check a proposed dt for my phase-field run.

skills/
    core-numerical/          # 8 skills: stability, solvers, meshing, convergence, ...
    simulation-workflow/     # 5 skills: validation, optimization, orchestration, ...
    hpc-deployment/          # 1 skill: SLURM job script generation
    ontology/                # 3 skills: ontology exploration, mapping, validation
    <each-skill>/
        SKILL.md             # Instructions + YAML frontmatter (with metadata block)
        scripts/             # Python CLI tools with --json output
        references/          # Domain knowledge documents
        evals/evals.json     # Evaluation suite (prompts, assertions)
        CHANGELOG.md         # Version history
tests/
    unit/                    # Pure-function tests via load_module()
    integration/             # Subprocess + JSON schema validation
    fixtures/                # Sample data files for CI smoke tests
.github/
    workflows/ci.yml         # Cross-platform CI + quality validation
    ISSUE_TEMPLATE/          # Bug reports, skill proposals
    PULL_REQUEST_TEMPLATE.md # PR checklist