PyBaMM
PyBaMM copied to clipboard
pybamm-team
[Bug]: Unphysical solutions from IDAKLU solver
PyBaMM Version
25.6.0
Python Version
3.10.12
Describe the bug
I made a Q&A post on this on the discourse channel, but I think it is actually more appropriate here:
I’ve encountered some concerning problems using the IDAKLU solver recently. I am coding a new model, as described in this paper. When I solve it using the Casadi solver, results look great. However, when I change to the IDAKLU solver, a strange discontinuity is observed in the particle concentration and the cell voltage.
You can reproduce this using the MWE below, by changing the solver choice at the bottom of the code. In Yang's model, there is a discontinuity in the time derivatives of surface concentration - I've highlighted this in the plot in my MWE. I don't see why this should be problematic, and I have solved other models which have discontinuities like this using IDAKLU without a problem. Nonetheless, the error occurs at the point of this discontinuity.
Some details on the model: 3 out of the 4 conservation equations in Yang’s model are identical to the DFN, whilst the particle mass balance is solved on the macroscale rather than the microscale (as in the DFN). They use a boundary layer theory to map between the average particle concentration and the surface concentration, for use in the Butler Volmer. I have made this version of Yang’s model from the Basic DFN framework.
Steps to Reproduce
#
# Yang's model, built from Basic Doyle-Fuller-Newman (DFN) Half Cell Model
#
import pybamm
import matplotlib.pyplot as plt
from pybamm.models.full_battery_models.lithium_ion.base_lithium_ion_model import BaseModel
class BasicYangModel(BaseModel):
"""
Rebuilt from basic DFN
Parameters
----------
options : dict
A dictionary of options to be passed to the model. For the half cell it should
include which is the working electrode.
name : str, optional
The name of the model.
"""
def __init__(self, options=None, name="Yang Model", my_params=None, model_type="averaged"):
options = {"working electrode": "positive"}
super().__init__(options, name)
pybamm.citations.register("Marquis2019")
# `param` is a class containing all the relevant parameters and functions for
# this model. These are purely symbolic at this stage, and will be set by the
# `ParameterValues` class when the model is processed.
self.my_params = my_params
######################
# Variables
######################
# Variables that depend on time only are created without a domain
Q = pybamm.Variable("Discharge capacity [A.h]")
# Variables that vary spatially are created with a domain.
c_e_s = pybamm.Variable(
"Separator electrolyte concentration [mol.m-3]", domain="separator"
)
c_e_w = pybamm.Variable(
"Positive electrolyte concentration [mol.m-3]", domain="positive electrode"
)
c_s_surf_w = pybamm.Variable("Positive particle surface concentration [mol.m-3]", domain="positive electrode")
c_e = pybamm.concatenation(c_e_s, c_e_w)
c_s_w = pybamm.Variable(
"Positive particle concentration [mol.m-3]", domain="positive electrode",
)
phi_s_w = pybamm.Variable(
"Positive electrode potential [V]", domain="positive electrode"
)
phi_e_s = pybamm.Variable(
"Separator electrolyte potential [V]", domain="separator"
)
phi_e_w = pybamm.Variable(
"Positive electrolyte potential [V]", domain="positive electrode"
)
phi_e = pybamm.concatenation(phi_e_s, phi_e_w)
# Constant temperature
T = self.param.T_init
######################
# Other set-up
######################
# Current density
i_cell = self.param.current_density_with_time
# Porosity and Transport_efficiency
eps_s = pybamm.PrimaryBroadcast(
pybamm.Parameter("Separator porosity"), "separator"
)
eps_w = pybamm.PrimaryBroadcast(
pybamm.Parameter("Positive electrode porosity"), "positive electrode"
)
b_e_s = self.param.s.b_e
b_e_w = self.param.p.b_e
# Interfacial reactions
j0_w = self.param.p.prim.j0(c_e_w, c_s_surf_w, T)
U_w = self.param.p.prim.U
ne_w = self.param.p.prim.ne
# Particle diffusion parameters
D_w = self.param.p.prim.D
c_w_init = pybamm.surf(self.param.p.prim.c_init)
# Electrode equation parameters
eps_s_w = pybamm.Parameter("Positive electrode active material volume fraction")
b_s_w = self.param.p.b_s
sigma_w = self.param.p.sigma
# Other parameters (for outputs)
c_w_max = self.param.p.prim.c_max
L_w = self.param.p.L
eps = pybamm.concatenation(eps_s, eps_w)
tor = pybamm.concatenation(eps_s**b_e_s, eps_w**b_e_w)
F_RT = self.param.F / (self.param.R * T)
RT_F = self.param.R * T / self.param.F
sto_surf_w = c_s_surf_w / c_w_max
j_w = (
2
* j0_w
* pybamm.sinh(ne_w / 2 * F_RT * (phi_s_w - phi_e_w - U_w(sto_surf_w, T)))
)
R_w = self.param.p.prim.R_typ
a_w = 3 * eps_s_w / R_w
a_j_w = a_w * j_w
a_j_s = pybamm.PrimaryBroadcast(0, "separator")
a_j = pybamm.concatenation(a_j_s, a_j_w)
# Yang params
j_p_spm = i_cell / (L_w * a_w) # TODO: delete me, here for debugging
l_s_p = R_w
D_s = pybamm.surf(D_w(c_s_w, T))
a_d_p = self.my_params["Positive electrode Yang fitting parameter"]
small_perturbation = 1e-200
t_dif_p = l_s_p ** 2 / D_s
t_cut_p = t_dif_p / (a_d_p ** 2)
regime_p = (pybamm.t < t_cut_p)
l_dif_p = (a_d_p * (D_s * (pybamm.t + small_perturbation))**0.5) * regime_p + l_s_p * (1 - regime_p)
if model_type == "averaged":
yang_corr_term = - (l_dif_p / (2 * l_s_p) - l_dif_p**2 / (6 * l_s_p**2)) * (j_p_spm * l_s_p / (D_s * self.param.F))
elif model_type == "full":
yang_corr_term = - (l_dif_p / (2 * l_s_p) - l_dif_p**2 / (6 * l_s_p**2)) * (j_w * l_s_p / (D_s * self.param.F))
self.algebraic[c_s_surf_w] = c_s_surf_w - (c_s_w + yang_corr_term)
######################
# State of Charge
######################
I = self.param.current_with_time
# The `rhs` dictionary contains differential equations, with the key being the
# variable in the d/dt
self.rhs[Q] = I / 3600
# Initial conditions must be provided for the ODEs
self.initial_conditions[Q] = pybamm.Scalar(0)
######################
# Particles
######################
# The div and grad operators will be converted to the appropriate matrix
# multiplication at the discretisation stage
self.rhs[c_s_w] = - a_j_w / (self.param.F * eps_s_w)
# Boundary conditions must be provided for equations with spatial
# derivatives
self.boundary_conditions[c_s_w] = {
"left": (pybamm.Scalar(0), "Neumann"),
#"right": (-j_w / pybamm.surf(D_w(c_s_w, T)) / self.param.F, "Neumann"),
"right": (pybamm.Scalar(0), "Neumann"),
}
self.initial_conditions[c_s_w] = c_w_init
self.boundary_conditions[c_s_surf_w] = {
"left": (pybamm.Scalar(0), "Neumann"),
"right": (pybamm.Scalar(0), "Neumann"),
}
self.initial_conditions[c_s_surf_w] = c_w_init
# Events specify points at which a solution should terminate
self.events += [
pybamm.Event(
"Minimum positive particle surface concentration",
pybamm.min(sto_surf_w) - 0.01,
),
pybamm.Event(
"Maximum positive particle surface concentration",
(1 - 0.01) - pybamm.max(sto_surf_w),
),
]
######################
# Current in the solid
######################
sigma_eff_w = sigma_w(T) * eps_s_w**b_s_w
i_s_w = -sigma_eff_w * pybamm.grad(phi_s_w)
self.boundary_conditions[phi_s_w] = {
"left": (pybamm.Scalar(0), "Neumann"),
"right": (
i_cell / pybamm.boundary_value(-sigma_eff_w, "right"),
"Neumann",
),
}
# multiply by Lx**2 to improve conditioning
self.algebraic[phi_s_w] = self.param.L_x**2 * (pybamm.div(i_s_w) + a_j_w)
# Initial conditions must also be provided for algebraic equations, as an
# initial guess for a root-finding algorithm which calculates consistent
# initial conditions
self.initial_conditions[phi_s_w] = self.param.p.prim.U_init
######################
# Electrolyte concentration
######################
N_e = -tor * self.param.D_e(c_e, T) * pybamm.grad(c_e)
self.rhs[c_e] = (1 / eps) * (
-pybamm.div(N_e) + (1 - self.param.t_plus(c_e, T)) * a_j / self.param.F
)
dce_dx = (
-(1 - self.param.t_plus(c_e, T))
* i_cell
/ (tor * self.param.F * self.param.D_e(c_e, T))
)
self.boundary_conditions[c_e] = {
"left": (pybamm.boundary_value(dce_dx, "left"), "Neumann"),
"right": (pybamm.Scalar(0), "Neumann"),
}
self.initial_conditions[c_e] = self.param.c_e_init
self.events.append(
pybamm.Event(
"Zero electrolyte concentration cut-off", pybamm.min(c_e) - 0.002
)
)
######################
# Current in the electrolyte
######################
i_e = (self.param.kappa_e(c_e, T) * tor) * (
self.param.chiRT_over_Fc(c_e, T) * pybamm.grad(c_e) - pybamm.grad(phi_e)
)
# multiply by Lx**2 to improve conditioning
self.algebraic[phi_e] = self.param.L_x**2 * (pybamm.div(i_e) - a_j)
# reference potential
L_Li = self.param.n.L
sigma_Li = self.param.n.sigma
j_Li = self.param.j0_Li_metal(pybamm.boundary_value(c_e, "left"), c_w_max, T)
eta_Li = 2 * RT_F * pybamm.arcsinh(i_cell / (2 * j_Li))
phi_s_cn = 0
delta_phi = eta_Li
delta_phis_Li = L_Li * i_cell / sigma_Li(T)
ref_potential = phi_s_cn - delta_phis_Li - delta_phi
self.boundary_conditions[phi_e] = {
"left": (ref_potential, "Dirichlet"),
"right": (pybamm.Scalar(0), "Neumann"),
}
self.initial_conditions[phi_e] = -self.param.n.prim.U_init
######################
# (Some) variables
######################
vdrop_cell = pybamm.boundary_value(phi_s_w, "right") - ref_potential
vdrop_Li = -eta_Li - delta_phis_Li
voltage = vdrop_cell + vdrop_Li
num_cells = pybamm.Parameter(
"Number of cells connected in series to make a battery"
)
c_e_total = pybamm.x_average(eps * c_e)
c_s_surf_w_av = pybamm.x_average(c_s_surf_w)
c_s_rav = pybamm.r_average(c_s_w)
c_s_vol_av = pybamm.x_average(eps_s_w * c_s_rav)
# Cut-off voltage
self.events.append(
pybamm.Event("Minimum voltage [V]", voltage - self.param.voltage_low_cut)
)
self.events.append(
pybamm.Event("Maximum voltage [V]", self.param.voltage_high_cut - voltage)
)
# Cut-off open-circuit voltage (for event switch with casadi 'fast with events'
# mode)
tol = 0.1
self.events.append(
pybamm.Event(
"Minimum voltage switch",
voltage - (self.param.voltage_low_cut - tol),
pybamm.EventType.SWITCH,
)
)
self.events.append(
pybamm.Event(
"Maximum voltage switch",
voltage - (self.param.voltage_high_cut + tol),
pybamm.EventType.SWITCH,
)
)
self.variables = {
"Time [s]": pybamm.t,
"Discharge capacity [A.h]": Q,
"Positive particle surface concentration [mol.m-3]": c_s_surf_w,
"X-averaged positive particle surface concentration "
"[mol.m-3]": c_s_surf_w_av,
"Positive particle concentration [mol.m-3]": c_s_w,
"Total lithium in positive electrode [mol]": c_s_vol_av
* L_w
* self.param.A_cc,
"Electrolyte concentration [mol.m-3]": c_e,
"Separator electrolyte concentration [mol.m-3]": c_e_s,
"Positive electrolyte concentration [mol.m-3]": c_e_w,
"Total lithium in electrolyte [mol]": c_e_total
* self.param.L_x
* self.param.A_cc,
"Current [A]": I,
"Current variable [A]": I, # for compatibility with pybamm.Experiment
"Current density [A.m-2]": i_cell,
"Positive electrode potential [V]": phi_s_w,
"Positive electrode open-circuit potential [V]": U_w(sto_surf_w, T),
"Electrolyte potential [V]": phi_e,
"Separator electrolyte potential [V]": phi_e_s,
"Positive electrolyte potential [V]": phi_e_w,
"Voltage drop in the cell [V]": vdrop_cell,
"Negative electrode exchange current density [A.m-2]": j_Li,
"Negative electrode reaction overpotential [V]": eta_Li,
"Negative electrode potential drop [V]": delta_phis_Li,
"Voltage [V]": voltage,
"Battery voltage [V]": voltage * num_cells,
"Instantaneous power [W.m-2]": i_cell * voltage,
"Yang correction term": pybamm.x_average(yang_corr_term),
}
parameter_values = pybamm.ParameterValues("Xu2019")
my_params = {
"Positive electrode Yang fitting parameter": 1.0,
}
models = {}
models["Yang"] = BasicYangModel(my_params=my_params, model_type="full")
models["DFN"] = pybamm.lithium_ion.DFN(options={"working electrode": "positive"})
sims = []
t_eval = [0, 6000]
for key, model in models.items():
solver = pybamm.IDAKLUSolver(atol=1e-8, rtol=1e-6)
#solver = pybamm.CasadiSolver()
sim = pybamm.Simulation(model, parameter_values=parameter_values, solver=solver)
sim.solve(t_eval=t_eval)
sims.append(sim)
pybamm.dynamic_plot(sims)
# Plot Yang's correction term
x = sims[0].solution["Time [s]"].entries / 3600
y = sims[0].solution["Yang correction term"].entries
plt.plot(x, y)
plt.xlabel("Time [h]")
plt.ylabel("Yang's correction term (surface conc - average conc)")
plt.show()
Relevant log output
Hi @isaacbasil, there are some known issues related to the time-dependent heaviside function (see https://github.com/pybamm-team/PyBaMM/pull/4994). As a workaround, can you please try evaluating the numeric value for t_cut_p and adding it to the t_eval vector?
Hi @MarcBerliner, thanks for your reply. I just tried that but unfortunately the same error is produced. The relevant code block now looks like this:
#======= Yang workaround ==========
l_s_p = parameter_values["Positive particle radius [m]"]
a_d_p = 1.0
D_s_p = parameter_values["Positive electrode diffusivity [m2.s-1]"]
t_dif_p = l_s_p ** 2 / D_s_p
t_cut_p = t_dif_p / (a_d_p ** 2)
# =================================
sims = []
t_eval = [0, 6000]
for key, model in models.items():
solver = pybamm.IDAKLUSolver()
#solver = pybamm.CasadiSolver()
sim = pybamm.Simulation(model, parameter_values=parameter_values, solver=solver)
if key == "Yang":
t_eval.insert(1, t_cut_p)
t_eval = np.array(t_eval)
sim.solve(t_eval=t_eval)
sims.append(sim)
@isaacbasil this workaround fixed it for me. Try replacing the t_cut_p line with
regime_p = (2 * pybamm.t < 2 * t_cut_p)
Pybamm tries to be fancy if it detects a heaviside function with a pybamm.t as one of the arguments. The multiple of 2 on the lhs and rhs is arbitrary, but it prevents this special casing, and it works as expected now.