wigging
commented
3 years ago The comments below are reposted from an Issue on the SciPy repository about the TR-BDF2 algorithm. The comments are from @laurent90git and are related to code in this repo.
I've quickly tested it. adding printouts of the time and time step in solve_ivp shows that the solution is actually advancing. I ran out of patience after a few minutes, having reached t~1ms. As your system is likely very stiff, at least during the initial transient, explicit methods perform very poorly and will explode/overflow if you don't use tight integration tolerances. They will anyway be very inefficient.
By the way, I see that you just call solve_ivp without any integration-related parameters (rtol, atol). I think you should take a closer look at these parameters and.their impact on the solution quality.
With implicit methods, every integration step is very long, because, as your system is quite nonlinear, the Jacobian of your ODE system must be computed anew very often. As your system is large (1500 ODEs) and scipy does not have any information regarding its sparsity pattern, the solver calls your ODE function dydt at least 1500 times for each Jacobian update... This results in very slow performance.
One thing I noticed is that, in your vector of unknowns, you have sorted the variables by type (all discrete rho_xx are contiguous for example). As the different fields interact with one another, this results in a Jacobian sparsity pattern that has components way off the diagonal. For one-dimensional models, it is usually best to sort all the variables according to the index of the mesh cell (or spatial point) that are related to. Then within each cell, the variables may be ordered as you wish, as long as the ordering is constant from one cell to the next. This way, the Jacobian will only have nonzero components very close to the diagonal. Consequently, much more efficient approaches can be used to compute this Jacobian without performing s as many function calls. You have something like 15 fields if I remember correctly, all discretise on the same mesh. In the ideal case, with the previous reordering of the unknowns, you could compute the Jacobian in 3*15=45 calls only, which is way better (I assume your spatial schemes only use a 3-points stencil, for exemple as in the second order centered finite difference for diffusive terms). This assumes that the Jacobian is not formed analytically but by finite-differences. A banded Jacobian may also enable the use of optimized linear solvers for the Newton loop, even though solve_ivp does not offer that level of customisation yet, as far as I am aware. I don't have time this weekend to edit your code to show you what I mean, but I'll try to do that next week!
Overall, as far as I've investigated, the issue is not linked to the lack of a performant integration method for stiff equations. I think it's just a combination of stiffness (requiring implicit integration), non-linearity (requiring repeated Jacobian updates) and large system size (inefficient Jacobian evaluation), all of which, added to Python's overall lack of performance (compared to other codes like Fortran or others) result in long computation times.
In the article you cited, I don't see any comment on the bad performance of other solvers (compared to ode23tb). Have you read that somewhere else?
laurent90git
The current implementation of the model using the Radau method with default tolerances does not converge to a solution. The model will run for several minutes and eventually output the warnings shown below. I terminated the program after 26 minutes because the solver had not obtained a solution. I'm not sure what's going on but suggestions on how to use SciPy's
solve_ivpfunction for this gasifier model would be helpful.