Sundial Solver Exposing Frequency Of Jacobian Update

[Dharshannan]

Description

Current IDAKLU solver uses the Sundials solver however, the Jacobian update frequency is not exposed as a user input. This causes issues with some models (like the Li-S multi-reaction phases) that rely on accurate Jacobian information since the solver uses an outdated Jacobian most of the time.

Motivation

This exposing of the functionality for the frequency of Jacobian update would help in the development of more complex models, specifically those that rely of multi-reaction stages like the Li-S models, as well as serve as a good tuning baseline for the solver to validate accuracy and performance.

Possible Implementation

No response

Additional context

No response

[martinjrobins]

Hi @Dharshannan, do you know which particular sundials options you would like to expose? My first suggestion would be to be reduce the max number of iterations of the non-linear solver, so that the solver would fail more often and so the jacobian would be updated more frequently. You can see all the sundials option functions here: https://sundials.readthedocs.io/en/latest/idas/Usage/SIM.html#optional-input-functions

However, I would also note that this jacobian is only used internally by sundials for the non-linear solver newton step, and this jacobian would not be used (directly) in the evaluation of the model rhs and algebraic equations. Can you please define what you mean by "some models (like the Li-S multi-reaction phases) that rely on accurate Jacobian information". Exactly how do these models use the jacobian?

[Dharshannan]

Hi @martinjrobins I am working on a separate repository to implement Li-S models into PyBaMM (with @rtimms and @brosaplanella), these models are very unstable, the more recent model, with 3 reaction stages is not able to solve a full discharge/charge. I can invite you to the repo for this project if that would help?
And we can work from there, perhaps organise a meeting to explain further would be helpful. Thanks.

[martinjrobins]

sure, please invite me to the repo it would be useful to see the equations you are trying to solve (if you have them written down that would be good too)

[Dharshannan]

@martinjrobins I have sent an invite, the model that is problematic is the one called Cornish2021.py in the models folder.

[martinjrobins]

@Dharshannan: I had a quick look at this. The numerical issues with your model seem to be occuring when S8 and S4 both go towards zero. When they are both under the level of numerical precision the ratio in your model (eq 4a in the paper) becomes quite unstable, see attached plot of the High plateau potential, and I think this might be causing the model to fail. I'm unsure how increasing the frequency of jacobian updates will help with this though, as once the newton solver starts failing it will update the jacobian anyway before trying again.