convergence-study

Files (10)

• CHANGELOG.md
• SKILL.md
• evals.json · evals/
• convergence_theory.md · references/
• gci_guidelines.md · references/
• refinement_order.py · scripts/_support/
• dt_refinement.py · scripts/
• gci_calculator.py · scripts/
• h_refinement.py · scripts/
• richardson_extrapolation.py · scripts/

SKILL.md (first lines)

# 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.10+
- No third-party packages — scripts use only the Python standard library (`math`).

## 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`, `results.extrapolated_value`, `results.notes` |

## 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? --> Order agreement; check persistence and error assumptions
|           +-- 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 | Evidence consistent with the expected error law | Check persistence across levels; label estimates conditional |
| 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 | See caveat below | Confirm with order comparison |
| GCI asymptotic ratio far from 1.0 | Not in asymptotic range | Refine further |

> **Asymptotic-ratio caveat (constant refinement ratios).** When the refinement
> ratios are equal (`r21 == r32`, the common case), the asymptotic ratio
> `AR = GCI_coarse / (r^p * GCI_fine)` reduces algebraically to `f1/f2`. It then
> only measures the relative gap between the two finest QoI values — **not** whether
> the data follow the assumed power law — so an `AR` near 1.0 can give false
> reassurance even when the observed order is far from expected. The
> `gci_calculator.py` JSON emits a `notes` entry flagging this. For a real
> asymptotic-range determination:
> 1. Compare the observed order `p` to the scheme's theoretical/expected order — a
>    match supports the expected error law but does not establish an asymptotic range.
> 2. For stronger verification, use 4+ systematically refined grids and check that
>    the observed order is consistent across successive grid triplets
>    (`h_refinement.py` reports one order per triplet plus `mean_order`).

## Verification checklist

When explaining a refinement result, compare only the dimensionless
`results.mean_order` with the expected order. The
`results.richardson_extrapolated_value` estimates the observable at zero spacing
and has the same units as the input observable; it is never the convergence
order. Use `order_comparison`, `extrapolation_evidence`, and
`scientific_interpretation` from `h_refinement.py` to keep these roles explicit.
For h=(0.2,0.1,0.05) and values=(1.04,1.01,1.0025), p=2 and the extrapolated
observable=1. Compare 2 with the expected order, not 1. Keep this guidance in
the final synthesis context. A single grid triplet cannot establish stable
order across triplets: `in_asymptotic_range` is null for a positive single
estimate. More grids, residual convergence and model validation remain needed.

Use the [deterministic explanation renderer](../../GROUNDED_EXPLANATIONS.md)
for final scientific text from trusted, input-verified `h_refinement` JSON.
An order agreement decision is true when the measured order is within the
stated tolerance, even though physical accuracy and asymptotic behavior remain
unproven. The relative tolerance 0.1 means 10%, not 0.1%.

From HeshamFS/materials-simulation-skills, Apache-2.0.

More from HeshamFS/materials-simulation-skills