Have we leveraged all capabilities of the cvode solver?

[StephanSchaller]

Now cvode seems to be a reliable solver for stiff problems, but we still run into numerical issues from time to time (well we do go complex, but still).

I am just wondering if we are sqeezing everything there is from the solver, or maybe if there are more sophisticated solvers out there by now?!

Most of the time our issue is: Simulation run failed: some variables became negative

Now, sometimes things can be improved by "thoughtful" implementations like in #919, or by using the "Scale Factor" property of setSpeciesInitialValue, but that's still not enough.

Any thoughts appreciated. Cheers, Stephan

[Yuri05]

@StephanSchaller First of all, please check CVODE usage notes (especially the section On controlling unphysical solutions and integrating over discontinuities)

I am just wondering if we are sqeezing everything there is from the solver

We don't. Some things could be improved.

• E.g. we use just one "global" absolute tolerance per simulation instead of per variable definition (which actually can be compensated by using scale factors).
• E.g. we could also use negative values check for automated reducing of the internal step size used by the solver. In some situations this would help to avoid negative values. However if the problem is not a numerical issue but the underlying model is not properly defined, it would be harder to estimate which exactly model components cause the error.

maybe if there are more sophisticated solvers out there by now

Maybe. CVODE is one of the best solvers for stiff ODEs. Which does not mean there were no alternatives.

@HDiedam FYI

[PavelBal]

A Comparison Between Differential Equation Solver Suites In MATLAB, R, Julia, Python, C, Mathematica, Maple, and Fortran.

According to the article, SVODE is indeed one of the best. DifferentialEquations.jl (Julia) seems to offer more feautres, better event handling, and can solve equations on GPU, which may be an advantage for feature tasks.

Also, just a quote from the comments:

Yeah, SUNDIALS is amazing and we all need to continuously learn from them. But, I wouldn't say it's far superior. Multistep methods can fail quite often (showed up just the other day in a StackOverflow post) especially at higher tolerances, and they don't work well when you have lots of events.

At some day, I would like to see the OSPS simulation core running on GPU. I must admit I have no idea how well ODE can be solved on GPU, but upgrading a GPU is always way simpler than upgrading the CPU(s).

[hdiedam]

There is definitely room for improvement in the solver interface, e.g. the proper treatment of implicit switches/events, but there is no simple solution for the stability of states that vary over multiple orders of magnitudes. Regarding GPUs: Solving stiff ODEs is a quite complex task, involving many evaluations of the custom right-hand-side function and for most algorithms its derivative, and therefore hard to parallelize on a GPU. There are very few papers that tried to bring implicit algorithms to GPUs at all and none I have seen so far suggested a performance increase for our problem sizes. Actually, even my limited tests with CPU-parallel implementations (which is easier to achieve than GPU-parallelism) have not let to any speed-ups for our typical PBPK models. In contrast, a proper sparse matrix handling yields a significant speed-up and is, outside of CVODES, equally hard to find in my experience. Regarding DifferentialEquations.jl: They offer quite a collection of algorithms, but we don't need most of them. They more or less suggest to use CVODES for our class of problems themselves (https://docs.juliadiffeq.org/latest/solvers/ode_solve.html -> stiff problems), but as they offer an interface, it might be quite easy to benchmark alternative algorithms for any problems you regard as too unstable. But also please keep in mind that apart from the algorithm, taking another programming language into the fold makes maintaining OSP harder.

[StephanSchaller]

@hdiedam

but there is no simple solution for the stability of states that vary over multiple orders of magnitudes.

We have the scale factors, but they have their limitations, as, especially for QSP models, if your simulation does not "engage" all system structures (when you have a model with multiple inputs but do not use all in one single experiment), only parts of the factors get set correctly, so it's more or less required to manually set all factors.