Problem with initialisation of a custom model

[isaacbasil]

Hi, I am trying to code Yang's model, as described in this paper.

I am running into a strange problem, in which the model does not initialise. The surface concentration is estimated from the average concentration using the following equations,

  l_s_p = pybamm.Parameter("Positive particle radius [m]")
  a_d_p = pybamm.Parameter("Positive electrode Yang fitting parameter")
  t_dif_p = l_s_p ** 2 / self.D_s_p
  t_cut_p = t_dif_p / (a_d_p ** 2)
  regime_p = (pybamm.t < t_cut_p)
  l_dif_p = (a_d_p * (self.D_s_p * pybamm.t)**0.5) * regime_p + l_s_p * (1 - regime_p)

  self.c_s_p_surf = self.c_s_p - (l_dif_p / (2 * l_s_p) - l_dif_p**2 / (6 * l_s_p)) * (self.j_p_bar * l_s_p / (self.D_s_p * self.F))

which leads to the following error,

pybamm.expression_tree.exceptions.SolverError: Could not find consistent states: Could not find acceptable solution: solver terminated successfully, but maximum solution error (165.894) above tolerance (1e-06)

This is apparently caused by the presence of pybamm.t in the equation, because if I instead use,
l_dif_p = (a_d_p * (self.D_s_p)**0.5) * regime_p + l_s_p * (1 - regime_p)
then the model solves without any problems.

Does anyone have any ideas for a workaround?

Thanks,

Isaac

MWE below.

import pybamm
from matplotlib import pyplot as plt

def u0_nmc_Chen2020(sto):
    u_eq = (
            -0.8090 * sto
            + 4.4875
            - 0.0428 * pybamm.tanh(18.5138 * (sto - 0.5542))
            - 17.7326 * pybamm.tanh(15.7890 * (sto - 0.3117))
            + 17.5842 * pybamm.tanh(15.9308 * (sto - 0.3120))
    )
    return u_eq

def j_simplified(c_s, c_s_max, phi_s):
    '''a simplified Butler Volmer'''
    T = 298.15
    R = 8.314
    j0 = 100.0
    x = c_s / c_s_max
    eta = phi_s - u0_nmc_Chen2020(x)
    return j0 * (pybamm.exp(0.5 * F / (R * T) * eta)
                    - pybamm.exp(- 0.5 * F / (R * T) * eta))

model = pybamm.BaseModel()

x_p = pybamm.SpatialVariable("x", domain="positive electrode", coord_sys="cartesian")

L_p = pybamm.Scalar(40e-6)

geometry = {"positive electrode": {x_p: {"min": 0, "max": L_p }}}

var_pts = {x_p: 40}

submesh_types = {"positive electrode": pybamm.Uniform1DSubMesh, "positive particle": pybamm.Uniform1DSubMesh}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
spatial_methods = {"positive electrode": pybamm.FiniteVolume(),"positive particle": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)

c_s_p = pybamm.Variable("c_s_p", domain="positive electrode")
phi_s_p = pybamm.Variable("phi_s_p", domain="positive electrode")


D_s_p = pybamm.Scalar(1e-14) # AM diffusivity
sigma_p = pybamm.Scalar(3.0)
vf_am_p = pybamm.Scalar(0.6) # volume fraction of AM
i_app = pybamm.Scalar(-30.0) 
a = pybamm.PrimaryBroadcast(300000, "positive electrode") # specific surface area
c_s_p_max = pybamm.Scalar(50400) 
c_s_0 = pybamm.Scalar(25000) 

j_p_spm = i_app / (L_p * a)


# ----- Defining parameters for Yang's model --------
l_s_p = pybamm.Scalar(4e-6) # particle length scale
a_d_p = pybamm.Scalar(1.0) # Yang's fitting parameter
F = pybamm.Scalar(96485) # Faraday
t_dif_p = l_s_p ** 2 / D_s_p
t_cut_p = t_dif_p / (a_d_p ** 2)
regime_p = (pybamm.t < t_cut_p)

# This is the crucial line... 
l_dif_p = (a_d_p * (D_s_p * pybamm.t)**0.5) * regime_p + l_s_p * (1 - regime_p)

# if the commented line below is instead used, the model solves.
#l_dif_p = (a_d_p * (D_s_p)**0.5) * regime_p + l_s_p * (1 - regime_p)


c_s_p_surf = c_s_p - (l_dif_p / (2 * l_s_p) - l_dif_p**2 / (6 * l_s_p)) * (j_p_spm * l_s_p / (D_s_p * F))

j_p = j_simplified(c_s_p_surf, c_s_p_max, phi_s_p)


a_j_p = a * j_p

dc_s_dt = - 1/(vf_am_p * F) * a_j_p 
model.rhs[c_s_p] = dc_s_dt

i_p = - sigma_p * pybamm.grad(phi_s_p) 
model.algebraic[phi_s_p] = L_p ** 2 * (- pybamm.div(i_p) - a_j_p)

grad_phi_rhs = i_app / sigma_p
model.boundary_conditions = {
    c_s_p: {"left": (0, "Neumann"), "right": (0, "Neumann")},
    phi_s_p: {"left": (pybamm.Scalar(0), "Neumann"), "right": (grad_phi_rhs, "Neumann"),} 
    }


model.initial_conditions = {c_s_p: c_s_0, phi_s_p:u0_nmc_Chen2020((c_s_0/c_s_p_max))}

model.variables = {"c_s_p": c_s_p, "phi_s_p": phi_s_p, "Voltage [V]": pybamm.boundary_value(phi_s_p, "right"), "Time [s]": pybamm.t}

disc.process_model(model)
sim = pybamm.Simulation(model)
solution = sim.solve([0, 1000])

plt.plot(solution["Time [s]"].entries, solution["Voltage [V]"].entries)
plt.xlabel("Time [s]")
plt.ylabel("Voltage [V]")
plt.show()

[MarcBerliner]

Hi @isaacbasil, can you please provide a minimum working example?

[isaacbasil]

Hi @MarcBerliner , I've just added one. It's not supposed to a realistic battery model, it just demonstrates the error caused by pybamm.t during the initialisation. I've highlighted the key line in the MWE.

[MarcBerliner]

It's not immediately clear to me what the issue is. But, if I add a tiny perturbation to the time in the offending line, like

l_dif_p = (a_d_p * (D_s_p * (pybamm.t + 1e-308))**0.5) * regime_p + l_s_p * (1 - regime_p)

then it works for me.

[isaacbasil]

Thanks @MarcBerliner that works for me too.