fix: improve Poisson solver IVP robustness (#162)

[alok-108]

Description

This PR addresses the numerical instabilities in the Poisson Solver Initial Value Problem (IVP) implementation as described in Issue #162.

When solving Poisson's equation as an IVP, inward integration for higher angular momentum modes ($l > 0$) is notoriously susceptible to catastrophic amplification of numerical noise from the starting boundary. To make solve_poisson_ivp robust across both compact and diffuse charge densities, this PR introduces several mathematical and programmatic safeguards:

Key Changes:

1. Adaptive Integration Bounds (r_max):
   Instead of using a fixed or excessively large starting radius, the integration bound is now adaptively scaled based on the density's decay. By calculating the cumulative integral of the radial density, we estimate r_99 (the radius enclosing 99% of the charge) and cap r_max to prevent integrating over unnecessarily massive distances (e.g., from r=1000) for highly compact densities.

2. Fallback ODE Solver Chain:
   The highly-accurate DOP853 is maintained as the default solver. However, if it fails to converge due to stiffness or step-size limits, solve_ode_ivp now implements a graceful fallback chain (DOP853 → RK45 → BDF → Radau → LSODA) and raises a warning instead of immediately crashing with a ValueError.

3. Filtering Numerical Noise ($l > 0$):
   For spherically symmetric charge densities, numerical integration on the Lebedev grid leaves tiny precision errors (e.g., ~1e-16) in the $l > 0$ radial components. The solver previously amplified this tiny noise to $10^{24}$! We now threshold components where the maximum absolute value is < 1e-12 and skip the ODE integration entirely, returning a zero spline to avoid artificial blowups.

4. Safe Boundary Extrapolation (make_safe_spline):
   Added a safeguard to handle evaluations outside the transform domain, returning the correct analytical behavior (Q/r for the monopole $l=0$, and 0 for $l>0$) beyond the boundary, rather than letting the CubicSpline extrapolate blindly.

5. Robustness Test Suite:
   Added tests/test_poisson_ivp_robustness.py to explicitly verify that the method is now stable for both compact and diffuse densities, and to ensure the fallback warning mechanism works correctly.

Closes #162

[Copilot]

Pull request overview

Addresses numerical instability in the Poisson IVP solver (#162) by adding adaptive integration bounds, a multi-method ODE fallback chain, noise filtering for higher-l components, and safe analytic extrapolation outside the radial domain.

Changes:

• poisson.py: Adaptive r_max from cumulative |density|, skip ODE integration for negligible l>0 components, wrap returned splines with make_safe_spline for r > r_max analytic continuation.
• ode.py: Replace single-method solve_ivp call with a fallback chain (DOP853 → RK45 → BDF → Radau → LSODA), warning on each failure and raising only if all fail.
• tests/test_poisson_ivp_robustness.py: New tests for IVP-vs-BVP agreement, fallback warning behavior, and adaptive r_max heuristic.

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated 7 comments.

File	Description
src/grid/poisson.py	Adaptive r_max, l>0 zero-source skip, safe-spline boundary wrapper
src/grid/ode.py	Multi-method ODE solver fallback chain with warnings
src/grid/tests/test_poisson_ivp_robustness.py	New robustness tests for IVP solver

---

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

[alok-108]

Hi @marco-2023, here's a visual summary of the changes in this PR. All Copilot review comments have been addressed. Please take a look when you can. Thanks!

🔍 What this PR adds (flow)

graph TD
    A[Poisson IVP call] --> B[Adaptive r_max<br/>from density decay]
    B --> C[For each l,m component]
    C --> D{Source term<br/>negligible?}
    D -- yes --> E[Skip ODE,<br/>return zero spline]
    D -- no --> F[Try DOP853]
    F -- success --> G[Wrap in safe spline]
    F -- fail --> H[Fallback chain:<br/>RK45 → BDF → Radau → LSODA]
    H -- success --> G
    H -- all fail --> I[Raise RuntimeError]
    G --> J[Interpolate safely<br/>even outside domain]
    J --> K[Return potential]

Loading