Adaptive (error-controlled) timestep for the MetSim erosion model (on by default)

[dvida]

Adaptive (error-controlled) timestep for the MetSim erosion model — on by default

Why

The fixed RK4 step const.dt = 0.005 s has no convergence basis and is badly under-resolved:

• On a fast Orionid single body, the dt=0.005 velocity error vs a converged run is ~6 km/s, converging only ~1st order (dt=1e-4 still ~0.56 km/s).
• The light-curve shape depends on the step size — historically you had to crank dt way down (and eat the runtime) to resolve the erosion tail.
• Small grains near erosion_mass_min overshoot their true velocity by tens of percent in a single step.

This PR makes accuracy a tolerance rather than a magic constant, turns it on by default, and keeps the fixed path bit-for-bit for anyone who wants the old behaviour.

What

• Adaptive steppers (MetSimErosionCyTools.pyx), per fragment, within a fixed output cadence:
  • adaptiveDP45Step() — the default: an embedded Dormand-Prince RK45 pair on the coupled [mass, velocity] system. One step yields the 5th-order solution and a 4th-order error estimate in ~7 RHS evals; the higher order takes larger steps (~3–4× faster than step-doubling).
  • adaptiveSingleBodyStep() — step-doubling on the original operator-split RK4, kept as a fallback (adaptive_high_order=False).
  • Both: PI step-size controller, atmosphere refreshed each sub-step, drag mass floored at m_kill so near-exhaustion grains stay stable, and erosion grains shed per sub-step so the light-curve shape no longer depends on dt.
• On by default. Constants.adaptive_dt=True, adaptive_rtol=1e-5 (targets numerical error below measurement noise ~0.1 mag / ~0.1 km/s at ~25 FPS, not machine precision). ablateAll() dispatches on const.adaptive_dt; False → the original fixed path, untouched. The macro step stays the fixed output cadence (one results row per dt), so no downstream consumer (GUI residuals, Dynesty likelihood, wake) changes.
• Grain-release cadence is a tunable physical resolution, decoupled from dt:
  • erosion_release_length (default 50 m) — an eroding fragment's sub-step is capped so grains are shed ~every this many metres of flight → grain births independent of dt/rtol. Now editable in the GUI.
  • erosion_release_vref (default 30000 m/s) — above this speed the cadence is frozen in time instead of distance, so the sub-step count stays bounded for fast meteors (a distance cadence makes cost grow with v). Slower fragments keep the exact distance cadence; vref ≤ 0 disables the cap. JSON-level knob.
• GUI: an "adaptive dt" checkbox (checked by default) that disables the fixed-dt box when on; rtol and grain-release-length boxes alongside, enabled with the checkbox (tooltips explain tuning). Fixed-dt box is used only when adaptive is off.
• Constants gain adaptive_dt / adaptive_high_order / adaptive_rtol / adaptive_atol_m / adaptive_atol_v / adaptive_dt_min / adaptive_dt_max / adaptive_max_substeps / erosion_release_length / erosion_release_vref; runSimulation/loadConstants back-fill them for older JSONs (which now pick up adaptive-on).
• Tests: wmpl/MetSim/Tests is an importable package; added adaptive-mode regression tests (no-NaN / mass-non-negative / grain stability, radiated-energy conservation vs fixed). Suite passes 8/8.

Performance (measured on the 13 fitted 2019 Orionids)

• Adaptive on by default is the win. Against a fixed dt fine enough to reach the same accuracy, adaptive rtol=1e-4 matches fixed dt≈0.001 at ~3.4× less time, and ~11× less than the fixed dt that beats it. So you get resolved erosion without hand-tuning dt down — the old workaround.
• Velocity-scaled cadence (erosion_release_vref=30000) gives a further mean 1.81× over a pure-distance cadence on these fast (~67 km/s) events, with a max peak-magnitude shift of 0.020 mag (≪ measurement noise; the change is confined to the faint erosion tail).
• Against the coarse legacy dt=0.005, adaptive costs more wall time (it actually resolves the physics), but that step was never accurate — see Why.
• rtol=1e-4 vs 1e-5 is a wash on eroding meteors (the 50 m grain cap, not rtol, sets the sub-step count), so the default stays at the tighter, safe 1e-5.

Validation

• Backward compat: adaptive_dt=False reproduces master bit-for-bit (np.array_equal, equal_nan=True, all output columns) on the synthetic scenario + 13 Orionids. Fixed mode is also bit-for-bit invariant to the new erosion_release_vref.
• Convergence (single body): adaptive cuts the velocity error from ~6 km/s (fixed dt=0.005) to well below measurement precision; DP45 agrees with step-doubling to ~0.06 mag.
• Erosion: light curves converge to the fine-fixed reference; radiated energy conserved (luminosity ratio ~0.997–1.005).

Tuning map

rtol = dynamics accuracy · erosion_release_length = grain-birth resolution · erosion_release_vref = fast-meteor cost cap · dt = output sampling cadence.

🤖 Generated with Claude Code

[dvida]

Adaptive timestep now on by default (adaptive_dt=True, rtol=1e-4)

Latest commit flips the defaults so the error-controlled adaptive integrator (DP45) is used out of the box, with adaptive_rtol=1e-4. The fixed-step path is untouched and remains bit-for-bit identical to master whenever adaptive_dt=False is set in the JSON/GUI, so legacy runs are fully reproducible.

Benchmarked on the 13 fitted 2019 Orionids (all strongly eroding, JSON dt=0.005). Two separate questions, two separate answers — worth stating plainly because they're easy to conflate.

1. Turning adaptive on is the real win (vs an equally-accurate fixed dt)

Getting well-resolved erosion out of the fixed scheme means cranking dt down by hand — the old workaround. Adaptive gives that resolution automatically and far cheaper. Cost vs light-curve accuracy for one representative event (20191023_091310), LC RMS measured against a fine dt=5e-5 fixed reference:

config	wall time	LC RMS vs fine ref
fixed dt=0.005 (old default)	0.05 s	0.554 mag
fixed dt=0.001	0.90 s	0.516 mag
fixed dt=0.0005	3.04 s	0.496 mag
fixed dt=0.0002	15.5 s	0.433 mag
adaptive rtol=1e-4	0.26 s	0.539 mag

Adaptive lands at roughly the accuracy of fixed dt=0.001 but at ~3.4× less time, and ~11× less than the fixed dt (0.0005) that actually beats it. That's the case for making it the default: good resolution without per-event dt tuning.

2. rtol=1e-4 vs 1e-5 is a wash on eroding meteors — but safe

Across all 13 Orionids, loosening rtol from 1e-5 to 1e-4 gave mean 1.0× speedup (1e-4 uses ~97% of the sub-steps) and a light curve identical to <0.0005 mag (max 0.0002 mag, max 0.03 km/s). The reason: for an eroding meteor the sub-step count is set by the erosion_release_length=50 m grain-release cap, not by the error tolerance — so rtol barely moves the cost. Same holds for single-body cases (313 → 302 sub-steps over 1e-5 → 1e-3), where the dt_max=0.005 ceiling binds instead.

So rtol=1e-4 is a slightly-looser, zero-downside default (identical output), but the speedup does not come from loosening rtol on these events. If a faster-with-1e-4 result shows up elsewhere, it'd be an event class where the tolerance actually binds (few grains, short flight) rather than the grain cap — none of the fitted Orionids are in that regime.

Backward compatibility & tests

• adaptive_dt=False → fixed path, bit-for-bit vs master (checked with np.array_equal, equal_nan=True, on the 13 Orionids + synthetic single-body).
• Older JSONs without the adaptive fields now back-fill to adaptive-on; set adaptive_dt=False for legacy fixed-step behaviour.
• Regression suite wmpl/MetSim/Tests/test_MetSimErosion.py: 8/8 pass.

Figures saved under 2019orionids/pr86_adaptive_benchmark/ (pr86_defaults_benchmark.png, default_rtol_1e4_vs_1e5.png).

[dvida]

Follow-up — reverted the rtol default to 1e-5, kept adaptive_dt=True.

The benchmark above showed rtol=1e-4 buys no speedup over 1e-5 on these events (sub-step count is set by the 50 m erosion_release_length grain cap, not the tolerance; single-body cases are bounded by dt_max). Since 1e-5 is tighter — more accuracy headroom in the regime where rtol actually binds — at negligible extra cost, it stays the default. So the net change is just adaptive on by default (adaptive_dt=True), which is the real win (§1: ~3.4–11× faster than an equally-accurate fixed dt). Commit 267cf47.