IDA adaptive time stepping

[MarcBerliner]

Description

This PR makes several changes to adaptive time-stepping in IDA. Previously, the t_eval decided both the points at which we store the solution and the points to explicitly stop the IDA solver. Now, the time-stepping has been broken down into t_eval (which retains the same behavior) and t_interp, which uses IDA's internal Hermite interpolator to efficiently store the solution. Additionally, we default to saving the adaptive time steps, which gives us full-order solution accuracy without negatively impacting the solve speed with periods or dense t_eval vectors.

This PR cleans up the main solve loop of the IDA solver and clearly delineates saving the full solution y vs. specific output variables.

Example: adaptive time outputs (new) vs. fixed outputs

The above picture shows the difference in resolution between the new adaptive time-stepping scheme and the previous approach without affecting the performance. The markers show every data point available in sol.t

Period/drive cycles improvements

parameter_values = pybamm.ParameterValues("Chen2020")
parameter_values.set_initial_stoichiometries(1)

experiment = pybamm.step.CRate(0.1, period = period, duration = 36000)
sim = pybamm.Simulation(
    model,
    solver=solver,
    parameter_values=parameter_values,
    experiment=experiment,
)

sim.solve()

The above figure shows the solve time for a full C/10 discharge. The new IDA solver is strictly faster than the old IDA and the Casadi solver. For the SPM, the IDA vs. Casadi timings are closer because IDA is a DAE solver while Casadi is using an ODE solver. Adding CVODE to pybamm could be future work to further improve our ODE speed.

For small periods, the new IDA solver is about an order of magnitude faster than the old IDA. For large periods, it's about half an order of magnitude faster than Casadi.

In a parameter estimation context, you can also post-interpolate the adaptive time-stepping solution without specifying a t_interp or period for just the variable of interest (such as voltage) which will dramatically improve the performance.

Adaptive time-stepping (no period)

parameter_values = pybamm.ParameterValues("Chen2020")
parameter_values.set_initial_stoichiometries(0)

experiment = pybamm.Experiment([f"Charge at {C_rate}C until 4.2 V"])
sim = pybamm.Simulation(
    model,
    solver=solver,
    parameter_values=parameter_values,
    experiment=experiment,
)

sim.solve()

When t_eval or a period is not set for an experiment, the new IDA is strictly faster than the previous IDA. Additionally, adaptive time-stepping has much better performance for slower simulations that previously used a dense t_eval vector by default.

Fixes #4312

Type of change

Please add a line in the relevant section of CHANGELOG.md to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

• Optimization (back-end change that speeds up the code)

Key checklist:

• No style issues: $ pre-commit run (or $ nox -s pre-commit) (see CONTRIBUTING.md for how to set this up to run automatically when committing locally, in just two lines of code)
• All tests pass: $ python run-tests.py --all (or $ nox -s tests)
• The documentation builds: $ python run-tests.py --doctest (or $ nox -s doctests)

You can run integration tests, unit tests, and doctests together at once, using $ python run-tests.py --quick (or $ nox -s quick).

Further checks:

• Code is commented, particularly in hard-to-understand areas
• Tests added that prove fix is effective or that feature works

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 99.46%. Comparing base (fcab586) to head (0a9e788).
Report is 195 commits behind head on develop.

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #4351      +/-   ##
===========================================
- Coverage    99.50%   99.46%   -0.05%     
===========================================
  Files          289      289              
  Lines        22146    22200      +54     
===========================================
+ Hits         22037    22081      +44     
- Misses         109      119      +10     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

---

🚨 Try these New Features:

• Flaky Tests Detection - Detect and resolve failed and flaky tests

[valentinsulzer]

Looks good to me other than some very minor comments. I have a feeling there are several bits of code we'll be able to refactor once this is in as the default everywhere but that can be a job for later. As usual will leave the C++ code to others to review

[martinjrobins]

Really great work @MarcBerliner, this looks excellent and some really great speed-ups! :) I've suggested a few changes below that need to be fixed but mostly minor. One general question I had: many of the idaklu tests use a dense t_eval and an identical t_interp, which is now the "slow" way of running the solver. When not specificly testing the interpolation code, I would use the most efficient / standard way of running solver.solve, which I believe would be t_eval = [start, stop] and t_interp = None or t_interp = np.linspace(start, stop, N)?

[MarcBerliner]

One general question I had: many of the idaklu tests use a dense t_eval and an identical t_interp, which is now the "slow" way of running the solver. When not specificly testing the interpolation code, I would use the most efficient / standard way of running solver.solve, which I believe would be t_eval = [start, stop] and t_interp = None or t_interp = np.linspace(start, stop, N)?

@martinjrobins good point. The most efficient way to run it is using the t_eval = [start, stop] and t_interp = None. A current limitation with this PR is that we do not have the Hermite interpolator available in python (adding it is difficult), so all the post-interpolation tests like sol[v](t) use linear interpolation, which is less accurate. I changed most tests to use t_interp = np.linspace(start, stop, N), which uses the Hermite interpolator.

Once this PR merges, we can fully migrate all the tests to using t_eval = [start, stop] and t_interp = np.linspace(start, stop, N). When we have more accurate post-interpolation, we can then change it to t_eval = [start, stop] and t_interp = None

[martinjrobins]

thanks @MarcBerliner, just a few more minor changes below

[martinjrobins]

Looks great, thanks @MarcBerliner