remove explicit sensitivity calculation (casadi solver)

[martinjrobins]

In order to calculate sensitivities, the pybamm casadi solver solves the full sensitivity equations due to the fact that casadi does not currently calculate sensitivities natively. This is inefficient compared with the staggard approach taken by libraries like Sundials (which we use in the IDAKLU solver), and will soon be unneccesary once Casadi completes current work on adding sensitivity calculations (casadi/casadi#3682). More importantly, the need to support explicit sensitivity calculations complicates a lot of BaseSolver, ProcessedVariable and Solution code, and forces extra complexity/inefficiency on the other solvers. Eg. the IDAKLU solver currently reshapes and post-processes its sensitivity calculations to fit in with the casadi solver's sensitivity output.

This issue proposes to remove the explicit sensitivity calculation support in PyBaMM, leaving the casadi solver unable to calculate sensitivities, at least until Casadi completes casadi/casadi#3682. The goal is to simplify on-going work on the solvers in #3910. If the user attempts to calculate sensivities with the casadi solver an error will be raised recommending the use of the IDAKLU solver instead.

[martinjrobins]

this was discussed at the developer meeting 28/10/2024 and everyone present was happy with removing the explicit sensitivity calculation for the casadi solver

[BradyPlanden]

Hi @martinjrobins, do you have a target release for this issue to land in? It would be great if this was packaged alongside a release with the IDAKLU becoming the default solver, as I think we want the default solver to support sensitivities, right? A user error is helpful, but having the default solver supporting sensitivities seems like something we would want to maintain.

[martinjrobins]

I agree that changing IDAKLU to be the default solver at the same time is a good idea. We can do this once we've got the solver separated out in #4487

[Azure-Dreamer]

it seems now casadi can calculate sensitivities #4678, though failed in solving my P2D sensitivity. I tried a simple ode in dae form and it gave correct result.

import casadi
x = casadi.SX.sym('x')
z = casadi.SX.sym('z')
p = casadi.SX.sym('p')

ode = x + p
alg = z
quad = x
dae = {'x':x, 'z':z, 'p':p, 'ode':ode, 'alg':alg, 'quad':quad}
x0 = 1.
p0 = 1.
# casadi.exp(1) = 2.718281828459045

opts = {}
I_01 = casadi.integrator('I', 'idas', dae, 0, 1, opts)
res_01 = I_01(x0=x0, p=p0)
xf_01 = res_01['xf']
zf_01 = res_01['zf']
qf_01 = res_01['qf']
print(f"xf_01 = {xf_01}, zf_01 = {zf_01}, qf_01 = {qf_01}")
# xf_01 = 4.43659 = 2 * exp(1) - 1
# zf_01 = 0, qf_01 = 2.43659

I_fwd_p_01 = I_01.factory('I_fwd', ['x0', 'z0', 'p', 'fwd:p'], ['fwd:xf', 'fwd:zf', 'fwd:qf'])
res_fwd_p_01 = I_fwd_p_01(x0=x0, p=p0, fwd_p=1)
fwd_xf_p_01 = res_fwd_p_01['fwd_xf']
fwd_zf_p_01 = res_fwd_p_01['fwd_zf']
fwd_qf_p_01 = res_fwd_p_01['fwd_qf']
print(f"fwd_xf_p_01 = {fwd_xf_p_01}, fwd_zf_p_01 = {fwd_zf_p_01}, fwd_qf_p_01 = {fwd_qf_p_01}")
# fwd_xf_p_01 = 1.71828 = exp(1) - 1
# fwd_zf_p_01 = 0, fwd_qf_p_01 = 0.718284

[martinjrobins]

it seems now casadi can calculate sensitivities #4678, though failed in solving my P2D sensitivity. I tried a simple ode in dae form and it gave correct result.

The explicit sensitivity calculation occur in pybamm rather than casadi. I was incorrect to say above that casadi does not have native forward sens calculation, it does, the main issue was that casadi did not support events (until casadi/casadi#3682). Either way we should still remove the explicit sens calculations (in the pybamm casadi solver wrapper)

[Azure-Dreamer]

@martinjrobins I’m sorry, but I’m not quite sure what is specifically meant by 'explicit sens calculations.' However, CasADi can calculate the sensitivity of xf and zf to p step by step, which is not exactly the same with Park's method, although already failed with my experiment in #4678.

d(xf)/d(p) = partial(xf)/partial(x0) + partial(xf)/partial(z0) + partial(xf)/partial(p)
d(zf)/d(p) = partial(zf)/partial(x0) + partial(zf)/partial(z0) + partial(zf)/partial(p)

where x0 and z0 are initial value, xf and zf are output. d(xf)/d(p) means the total derivative of xf with respect to p, and partial(xf)/partial(p) means the partial derivative of xf with respect to p at current step.
Based on my understanding, I think this can support events as well. I'm not sure if my idea is correct.

[kratman]

I marked myself as assigned for this one since I am planning to start switching the default solver to IDAKLU