Solver differences with sinusoidal currents

[bardsleypt]

Hello,

I am using PyBaMM to run EIS simulations, by applying a sinusoidal current, and measuring the responding voltage, from which I can extract impedance measurements. The issue I'm having, is I can get fairly different voltage results by using different solvers (Casadi with various options, and Scikit).

I've attached a template script to this discussion, as well as the output figure that shows the resulting voltage curves. I tend to believe the 'scikit' and 'safe w/o grid' solutions (as indicated in my attached plot, they lie on top of one another), but would like to understand a few things more clearly. Namely,

• Why does the 'fast' solver introduce voltage discontinuities/higher-modes in the voltage curve? Note, this behavior is identical to using the casadi mode='safe', so the mode='safe without grid' seems to be a better option here.
• Why does the 'fast' solver create an entirely different voltage curve when 'stepping' through the solution per each time step (as indicated by the 'step' solution in my attached plot)?
• I've also tested this with casadi mode='safe' while also explicitly controlling dt_max, this gives similar behavior as the 'safe without grid'. This leads me to believe pybamm+casadi determines its own grid, and simply interpolates to the desired 't_eval' points when returning the solution. Is this correct?

Hopefully these questions are clear enough, thanks to anyone who can provide insight. Happy to clarify if needed.

Best,
Patrick

import pybamm
import numpy as np
from functools import partial
import matplotlib.pyplot as plt

model = pybamm.lithium_ion.DFN(options={'surface form': 'differential'})
params = pybamm.ParameterValues('Ramadass2004')


def current(t, amp=100e-3, freq=1):
    return amp * pybamm.sin(2 * np.pi * freq * t)


# set up frequency params
freq = 2e3
samples_per_period = 128
num_periods = 3

# set up t_eval
dt = 1 / freq / samples_per_period
t_eval = np.array(range(0, samples_per_period*num_periods)) * dt

# Update current function with 2 kHz sinusoid
params['Current function [A]'] = partial(current, amp=100e-3, freq=freq)

# Set up different solvers
solver_fast = pybamm.CasadiSolver(mode='fast')
solver_safe = pybamm.CasadiSolver(mode='safe without grid')
solver_scik = pybamm.ScikitsDaeSolver()

# Get different solutions
sol_fast = pybamm.Simulation(model=model, parameter_values=params).\
    solve(t_eval=t_eval, solver=solver_fast)
sol_safe = pybamm.Simulation(model=model, parameter_values=params).\
    solve(t_eval=t_eval, solver=solver_safe)
sol_scik = pybamm.Simulation(model=model, parameter_values=params).\
    solve(t_eval=t_eval, solver=solver_scik)
sim_step = pybamm.Simulation(model=model, parameter_values=params, solver=solver_fast)
t, t_max, dt = 0, t_eval[-1], t_eval[1]
while t < t_max:
    sol_step = sim_step.step(dt)
    t += dt

fig, axs = plt.subplots(2, 1, figsize=(10, 10/1.6), sharex=True,
                        constrained_layout=True)
for y_key, ax in zip(['Current [A]', 'Terminal voltage [V]'], axs):
    for sol, label in zip([sol_fast, sol_safe, sol_scik, sol_step],
                          ['fast', 'safe', 'scik', 'step']):
        ax.plot(sol['Time [s]'].data, sol[y_key].data, lw=2, label=label)
    ax.legend(bbox_to_anchor=(1, 0.5), loc='center left', fontsize=16)
    ax.set_ylabel(y_key, fontsize=18)
    ax.tick_params(labelsize=16)
    ax.grid('on')
axs[1].set_xlabel('Time [s]', fontsize=18)
fig.suptitle('Solver diffs', fontsize=20)
fig.show()

[valentinsulzer]

Yeah this was reported before and I couldn't get to the bottom of it #1482

I don't know what's happening. Does it happen with different models?

[rtimms]

It works ok for SPM and SPMe, but not Newman-Tobias. It also works fine at lower frequencies. Some strange combination of model + solver choice...

[valentinsulzer]

Is it a ODE vs DAE thing?

[rtimms]

I thought that but formulated with capacitance it's an ODE. I was hoping it would be that.

It's weird that the solver handles it OK with SPM, which seems to suggest the way it's constructing the grid isn't the issue.

I've also tried with finer meshes which doesn't seem to help either.

[bardsleypt]

I did just now check it for 'surface form' both None and 'differential', which I believe changes the system from being DAE-based to purely DE based, respectively (I might be wrong, you both probably know better). In any case, the discretization/jumping effect is present in both cases.

[bardsleypt]

Yes, to confirm with Robert, here is the same plot I generated before but for SPMe; the discretization effects are not present here. I may have some time today to step into the code at a deeper level - I'll report anything else I may find out. It is good to know this is a known issue and it wasn't me doing something entirely wrong.

[valentinsulzer]

The "fast" solution is correct after each jump, while between the jumps the gradient of the fast and step solutions are the same. That might help elucidate what is going on.

[bardsleypt]

Is there any chance it is some type of indexing problem? When the Casadi solution points are put into the pybamm.Solution() object, perhaps they are out of order in some way so that upon retrieving the variables from the solution object they give the jumping effect?

I say that, but looking at the plots, the flat parts of the jumps don't seem to correspond with any of the correct solution points, which I would expect if it were simply an indexing problem.

[rtimms]

Also, the slope of the "fast" solution agrees with the slope of the stepped one, but it somehow gets back on track with the jumps where the stepped one doesn't. It's odd, particularly that it's not for all models.

[valentinsulzer]

We should check if we can reproduce this with pure casadi code. If so, we can ask the casadi people for help. If not, it's a bug in our casadi solver

[valentinsulzer]

Would also be interesting to see whether dt_max has any effect on the results