pybamm-team
Sundial Solver Exposing Frequency Of Jacobian Update
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
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?
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.
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)
@martinjrobins I have sent an invite, the model that is problematic is the one called Cornish2021.py in the models folder.
@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.
@martinjrobins this is what I get when running with my bespoke solver btw. The solver uses numpy and always recalculates the jacobian at every time step.
Thanks @Dharshannan, what method are you using, and does your solver have adaptive timestepping? as you don't seem to be resolving the discontinuity near 2 Ah very well. I am able to get sundials to get over this discontinuity by reducing the max order of the solver to 2, but the simulation is still quite unstable, so I'm not sure how useful this is. Are you able to re-formulate your equations to remove the singularity?
@martinjrobins I am using a modified newton raphson approach with adaptive time stepping using linesearch (and another ADHoc time stepping method that keeps track of the determinant magnitude of the jacobian matrix). For Li-S the discontinuity is meant to happen, there is meant to be a dip and rise in the voltage which is characterized by the precipitation reaction of sulphide species. By adjusting the equilibrium potential the dip and rise can happen at either a lower or higher potential.
yea, its just that the BDF solver in sundials is a multi-step method and does not typically deal well with discontinuities. This is why we need to stop and restart the solver whenever there are any discontinuities in the input or model.
@martinjrobins what would you advice in terms of tuning the IDA solver parameters, or just to deal with this instability in general?