[mleot]

Description

Fixes #5798
Fixes #5799

A differentiated 1D Interpolant now works in the solvers for every interpolator. The
heat of mixing differentiates the OCP, so "heat of mixing": "true" failed with any
tabulated OCP or entropic change: a "cubic" table crashed at solve time and a
"linear" one (the default interpolator) at build time. Ai2020 as shipped hit the
cubic case.

Cubic (#5798). scipy's BSpline.derivative() keeps a coefficient array longer than
the spline's basis (52 coefficients for 49 basis functions, from 50 points) and ignores
the tail. _cubic_to_casadi passed the whole array. CasADi evaluates that, but its own
derivative of the spline, which the solver's Jacobian needs, rejects it:

reshape: number of elements must remain the same. Old shape is 52x1. New shape is 49x1=49.

The #5582 tests compare values only, so they passed. Only the basis coefficients are
now passed; values are unchanged.

Linear (#5799). interp1d has no .derivative(), so .diff raised
AttributeError. The derivative is now the equivalent degree-1 spline differentiated
once (the slope of the interval containing the point, right-continuous at knots, the end
intervals' slope when extrapolating), converted to CasADi as a slope lookup in the same
way as the pchip conversion looks up its coefficients. Higher derivatives are zero.

Unchanged: every interpolant value and derivative that worked before is
bit-identical in CasADi inside the data range; undifferentiated linear interpolants
still use interp1d and interpn_linear.

With a 201-point table of Chen2020's graphite OCP, the heat of mixing now matches the
analytic OCP's to 1e-3 (linear), 1e-4 (pchip) and 3e-6 (cubic), for SPM and DFN.

Not addressed here: a cubic Interpolant's CasADi conversion is zero outside its
data range (scipy extrapolates), for values as well as derivatives. That predates this
PR and will get its own issue.

Type of change

Bug fix: an entry under ## Bug fixes in CHANGELOG.md.

Important checks:

• No style issues: nox -s pre-commit (every hook passes on the branch diff)

• All tests pass: nox -s tests (locally 4558 unit, 515 integration and 63
  pybammsolvers tests pass on this branch)

• The documentation builds: nox -s doctests (not run; no docs changed)

• Code is commented for hard-to-understand areas

• Tests added that prove fix is effective or that feature works:

  • CasADi's Jacobian of a differentiated interpolant must equal its second derivative,
    for linear, cubic and pchip on uniform and non-uniform grids
  • a linear interpolant's derivative is the interval slope at knots, inside the data and
    when extrapolating, for vector data, and its second derivative is zero
  • the heat of mixing with a tabulated OCP (linear, pchip, cubic) matches the analytic
    OCP for SPM and DFN

  Nine unit and four integration cases fail without the fix.

🤖 Generated with Claude Code

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.00%. Comparing base (aa8a48b) to head (8bcae57).

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #5804   +/-   ##
=======================================
  Coverage   98.00%   98.00%           
=======================================
  Files         346      346           
  Lines       35259    35280   +21     
=======================================
+ Hits        34554    34575   +21     
  Misses        705      705           

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[martinjrobins]

thanks @mleot, can you also check if this still works with serialising the interpolant (packages/pybamm/src/pybamm/expression_tree/operations/serialise.py:105), and diffsl conversion of the interpolant (packages/pybamm/src/pybamm/expression_tree/operations/diffsl.py:1280)?

[mleot]

Thanks @martinjrobins, good call on both of these.

For serialisation, the JSON round trip was already fine, but to_source was dropping _num_derivatives, so a differentiated interpolant would come back as the interpolant itself. I added it to the generated constructor, with tests for linear and cubic, first and second derivatives.

DiffSL was silently wrong too. interp1d just interpolates the table values, so a derivative was being exported as U instead of dU/dx (0.26 instead of 1.0 for U = x² at 0.5). I don't think there's a clean way to do a linear derivative with interp1d since it's piecewise constant, so for now it raises a NotImplementedError. I can open an issue for proper support if that would be useful.

I also rebased onto main. Brady's #5955 changed the CasADi interpolants to find which interval a point is in by binary search, since the old shortcut could pick the wrong interval right at a knot, so I updated the new linear derivative to do the same. I think it's ready for another look when you get the chance.