Tip

An interactive online version of this notebook is available, which can be accessed via Open this notebook in Google Colab

Alternatively, you may download this notebook and run it offline.

Attention

You are viewing this notebook on the latest version of the documentation, where these notebooks may not be compatible with the stable release of PyBaMM since they can contain features that are not yet released. We recommend viewing these notebooks from the stable version of the documentation. To install the latest version of PyBaMM that is compatible with the latest notebooks, build PyBaMM from source.

Interpolation and evaluation points#

In this notebook we will discuss the use of interpolation and evaluation points in PyBaMM. These points are used to control the solver and the output of the simulation, and can be used to improve the performance of the solver, or if they are used incorrectly, can lead to decreased performance.

[1]:

%pip install "pybamm[plot,cite]" -q    # install PyBaMM if it is not installed
import pybamm
import numpy as np
import time

Note: you may need to restart the kernel to use updated packages.

Solver timestepping#

Every PyBaMM solver solves the equations by starting from a set of initial conditions and then “stepping” forward in time. At each step, the solver solves the equations, including any algebraic conditions, and evaluates the local error to determine how long each time-step should be in order to satisfy the tolerances supplied by the user. Therefore, the progress of the solver is defined by a sequence of internal time-points at which the solver solved the equations. For example if we wanted to solve the equations from t=0 to t=10s, the solver might choose to solve the equations at the following time-points:

|-------|-------------------------------|---------------------------|----------|------|-----|
0      0.9                             4.9                         7.3        8.5     9.1   10

At each time point \(i\), the solver will store the solution to the equations \(\mathbf{y}_i\) and perhaps some other relvent information about the timestep. At the end of the solve, the solver will return the solution comprised of the \(n\) solutions \(\mathbf{y}_i\) at each time-point \(t_i\). In our example above the solution would consist of \(n=7\) solutions. This data can be later used to either plot the solution, or to evaluate any output variables that the user has requested.

[2]:

sim = pybamm.Simulation(pybamm.lithium_ion.SPM(), solver=pybamm.IDAKLUSolver())
sol = sim.solve([0, 10])
print("solution was generated at times", sol.t)

solution was generated at times [ 0.    0.01  0.02  0.04  0.08  0.16  0.32  0.64  1.28  2.56  3.84  5.12
  7.68 10.  ]

You can see that the initial timestep used by the solver starts out quite small, but then exponentially increases as the solver becomes more confident in the solution. The method in this case is a multi-step method with variable order, so the solver is able to take longer steps with an increased order of accuracy as it builds up information over multiple steps.

Evaluation points#

Thus far, we have described how the solver chooses its own internal timepoints in order to solve the equations at. In some cases, however, these timepoints need to be supplied by the user to ensure an accurate solution. A trivial example of this is the start and end timepoints that you wish to perform the solve. Other points include the points at which the solution is discontinuous, which the solver cannot know in advance or detect. For this reason PyBaMM also allows the user to supply a set of evaluation points to the solver using the `t_eval <https://docs.pybamm.org/en/stable/source/api/solvers/base_solver.html#pybamm.BaseSolver.solve>`__ argument to the solve function. The solver will then make sure to stop at each evaluation point and return the solution at these points. For example, say if t_eval = [0, 5, 10], the solution returned by the solver might be:

|-------|-----------------------------|--|--------|-----------------|-----------------|-----|
0      0.9                           4.9 5        5.9               7.3               9     10

In this case the solver has returned not only the evaluation points requested, but also all the internal time points that it used.

Normally, a PyBaMM user would not be required to supply any evaluation points other than the start and end points of the simulation. PyBaMM itself detects any discontinuous events in the model and adds these to the evaluation points before passing them to the solver. However, in some cases the user may know in advance that the solution is discontinuous at a certain point, and it is more accurate and numerically stable to supply this point to the solver. Note, however, that every additional evaluation point will increase the time taken to solve the equations as the solver will be required to take more steps to solve the equations, so only add additional evaluation points if required.

Below we add an additional evaluation point at \(t=5s\) to the example above.

[3]:

sim = pybamm.Simulation(pybamm.lithium_ion.SPM(), solver=pybamm.IDAKLUSolver())
sol = sim.solve(t_eval=np.array([0, 5, 10]))
print("solution was generated at times", sol.t)

solution was generated at times [ 0.     0.01   0.02   0.04   0.08   0.16   0.32   0.64   1.28   2.56
  3.84   5.     5.005  5.01   5.02   5.04   5.08   5.16   5.32   5.64
  6.28   7.56  10.   ]

Adding the evaluation point at \(t=5s\) will force the solver to stop and restart integration at this point. You can see after this point the solver takes smaller steps as it restarts, and then gradually increases the step size again as it becomes more confident in the solution.

Interpolation points#

When evaluating output variables using the solution, if no interpolation points are provided PyBaMM will interpolate between the internal time-points to get the solution at the time-points requested by the user. For example, if the user requested the solution at \(t=0.1s\), PyBaMM would interpolate the solution between \(t=0\) and \(t=0.9\) to get the solution at \(t=0.1s\). However, often a user will know in advance the time-points at which they want to evaluate the solution, and it can be more efficient to simply provide these time-points to the solver in advance so that it can do this interpolation during the solve.

The IDAKLU solver allows the user to provide a set of interpolation points to the solver using the `t_interp <https://docs.pybamm.org/en/stable/source/api/solvers/base_solver.html#pybamm.BaseSolver.solve>`__ argument to the solve function. The solver will then interpolate the solution on-the-fly during the solve, and return the solution at the requested time-points. For example, say if t_interp = [2, 4, 6, 8], and the solver takes the same internal time-points as in the example above, the solution returned by the solver would be:

|-------|---------*----------------*----|--------------------*------|------*---|------|-----|
0                 2                4                         6             8                10

where the * represent the solution at the requested time-points. The solver is still stepping to the same internal time-points in order to solve the equations, but it is also interpolating the solution at the requested time-points * and storing only these interpolated solutions. Therefore, use interpolation points if you know in advance the time-points at which you want to evaluate the solution, as interpolating to any other time-point post-solve will be much less accurate.

[4]:

import pybamm
import numpy as np

sim = pybamm.Simulation(pybamm.lithium_ion.SPM(), solver=pybamm.IDAKLUSolver())
sol = sim.solve(t_eval=[0, 10], t_interp=[2, 4, 6])
print("solution was generated at times", sol.t)

solution was generated at times [ 0.  2.  4.  6. 10.]

You can see that once we provide interpolation points to the solver, the solution will no longer store the solution at every internal time-point, but only at the interpolation and evaluation points provided.

Performance considerations#

We have already mentioned that adding additional evaluation points will increase the time taken to solve the equations, as the solver will be required to take more steps and will restart itself at each evaluation point, which can be computationally expensive.

Using interpolation points can help to decrease the time taken by the solver, but only if the number of internal timesteps is much greater than the number of interpolation points. For example, the simulation below uses >3000 internal timesteps to compute the whole solution, but we know in advance that we only want to evaluate the solution at 1000 known time-points. In this case, providing interpolation points to the solver will be more efficient than storing all the internal time-points and then interpolating the solution post-solve.

[5]:

parameter_values = pybamm.ParameterValues("Chen2020")
parameter_values.set_initial_stoichiometries(1)
experiment = pybamm.step.CRate(0.1, period=10, duration=36000)
sim = pybamm.Simulation(
    pybamm.lithium_ion.DFN(),
    solver=pybamm.IDAKLUSolver(),
    parameter_values=parameter_values,
    experiment=experiment,
)
sol = sim.solve()
print(f"Number of internal time steps: {len(sol.t)}")
t_final = sol["Time [h]"].entries[-1]
t_data = np.linspace(0, t_final, 1000)

start_time = time.perf_counter()
sol = sim.solve()
voltage = sol["Terminal voltage [V]"](t_data)
end_time = time.perf_counter()
print(f"Time to solve (no t_interp): {end_time - start_time}s")

start_time = time.perf_counter()
sol = sim.solve(t_interp=t_data)
voltage = sol["Terminal voltage [V]"].data
end_time = time.perf_counter()
print(f"Time to solve (with t_interp): {end_time - start_time}s")

Number of internal time steps: 3601
Time to solve (no t_interp): 0.07621171302162111s
Time to solve (with t_interp): 0.030624893959611654s

In the opposite case, lets increase the number of interpolation points so that they are greater than the number of internal time-points. In this case, more work will be required in both cases due to the cost of interpolating at more points. However, storing these points instead of only the internal time-points will be costly, so using interpolation points is slower in this case.

[6]:

t_data = np.linspace(0, t_final, 10000)
start_time = time.perf_counter()
sol = sim.solve()
voltage = sol["Terminal voltage [V]"](t_data)
end_time = time.perf_counter()
print(f"Time to solve (no t_interp): {end_time - start_time}s")

start_time = time.perf_counter()
sol = sim.solve(t_interp=t_data)
voltage = sol["Terminal voltage [V]"].data
end_time = time.perf_counter()
print(f"Time to solve (with t_interp): {end_time - start_time}s")

Time to solve (no t_interp): 0.06396568496711552s
Time to solve (with t_interp): 0.12752164108678699s