Solver Error when ionic conductivity < 0.05 S/m in custom model

[AtishaySinghai]

I'm encountering a SolverError when running a custom PyBaMM model with ionic conductivity set below 0.05 S/m. The solver fails to find a consistent state during experiment initialization.
Here, is the code. Can anyone help

import pybamm
class BaseModel(pybamm.BaseBatteryModel):
"""
Overwrites default parameters from Base Model with default parameters for
lithium-ion models

Parameters
----------
options : dict-like, optional
    A dictionary of options to be passed to the model. If this is a dict (and not
    a subtype of dict), it will be processed by :class:`pybamm.BatteryModelOptions`
    to ensure that the options are valid. If this is a subtype of dict, it is
    assumed that the options have already been processed and are valid. This allows
    for the use of custom options classes. The default options are given by
    :class:`pybamm.BatteryModelOptions`.
name : str, optional
    The name of the model. The default is "Unnamed battery model".
build : bool, optional
    Whether to build the model on instantiation. Default is True. Setting this
    option to False allows users to change any number of the submodels before
    building the complete model (submodels cannot be changed after the model is
    built).
"""

def __init__(self, options=None, name="Unnamed lithium-ion model", build=False):
    super().__init__(options, name)
    self.param = pybamm.LithiumIonParameters(self.options)

    self.set_standard_output_variables()

def set_submodels(self, build):
    self.set_external_circuit_submodel()
    self.set_porosity_submodel()
    self.set_interface_utilisation_submodel()
    self.set_crack_submodel()
    self.set_active_material_submodel()
    self.set_transport_efficiency_submodels()
    self.set_convection_submodel()
    self.set_open_circuit_potential_submodel()
    self.set_intercalation_kinetics_submodel()
    self.set_particle_submodel()
    self.set_solid_submodel()
    self.set_electrolyte_concentration_submodel()
    self.set_electrolyte_potential_submodel()
    self.set_thermal_submodel()
    self.set_surface_temperature_submodel()
    self.set_current_collector_submodel()
    self.set_sei_submodel()
    self.set_sei_on_cracks_submodel()
    self.set_lithium_plating_submodel()
    self.set_li_metal_counter_electrode_submodels()
    self.set_total_interface_submodel()

    if build:
        self.build_model()

@property
def default_parameter_values(self):
    if self.options.whole_cell_domains == [
        "negative electrode",
        "separator",
        "positive electrode",
    ]:
        return pybamm.ParameterValues("Marquis2019")
    else:
        return pybamm.ParameterValues("Xu2019")

@property
def default_quick_plot_variables(self):
    if self.options.whole_cell_domains == ["separator", "positive electrode"]:
        return [
            "Electrolyte concentration [mol.m-3]",
            "Positive particle surface concentration [mol.m-3]",
            "Current [A]",
            "Electrolyte potential [V]",
            "Positive electrode potential [V]",
            "Voltage [V]",
        ]
    else:
        return [
            "Negative particle surface concentration [mol.m-3]",
            "Electrolyte concentration [mol.m-3]",
            "Positive particle surface concentration [mol.m-3]",
            "Current [A]",
            "Negative electrode potential [V]",
            "Electrolyte potential [V]",
            "Positive electrode potential [V]",
            "Voltage [V]",
        ]

def set_standard_output_variables(self):
    super().set_standard_output_variables()

    # Particle concentration position
    var = pybamm.standard_spatial_vars
    if self.options.electrode_types["negative"] == "porous":
        self.variables.update({"r_n [m]": var.r_n})
    if self.options.electrode_types["positive"] == "porous":
        self.variables.update({"r_p [m]": var.r_p})

def set_degradation_variables(self):
    """Sets variables that quantify degradation (LAM, LLI, etc)"""

    domains = [d for d in self.options.whole_cell_domains if d != "separator"]
    for domain in domains:
        Domain = domain.capitalize()
        self.variables[f"Total lithium in {domain} [mol]"] = sum(
            self.variables[f"Total lithium in {phase} phase in {domain} [mol]"]
            for phase in self.options.phases[domain.split()[0]]
        )

        # LAM
        Q_k = self.variables[f"{Domain} capacity [A.h]"]
        domain_param = getattr(self.param, domain[0])  # param.n or param.p
        LAM_k = (1 - Q_k / domain_param.Q_init) * 100
        self.variables.update(
            {
                f"LAM_{domain[0]}e [%]": LAM_k,
                f"Loss of active material in {domain} [%]": LAM_k,
            }
        )

    # LLI
    n_Li_e = self.variables["Total lithium in electrolyte [mol]"]
    n_Li_particles = sum(
        self.variables[f"Total lithium in {domain} [mol]"] for domain in domains
    )
    n_Li = n_Li_particles + n_Li_e

    # LLI is usually defined based only on the percentage lithium lost from
    # particles
    LLI = (1 - n_Li_particles / self.param.n_Li_particles_init) * 100
    LLI_tot = (1 - n_Li / self.param.n_Li_init) * 100

    self.variables.update(
        {
            "LLI [%]": LLI,
            "Loss of lithium inventory [%]": LLI,
            "Loss of lithium inventory, including electrolyte [%]": LLI_tot,
            # Total lithium
            "Total lithium [mol]": n_Li,
            "Total lithium in particles [mol]": n_Li_particles,
            "Total lithium capacity [A.h]": n_Li * self.param.F / 3600,
            "Total lithium capacity in particles [A.h]": n_Li_particles
            * self.param.F
            / 3600,
            # Lithium lost
            "Total lithium lost [mol]": self.param.n_Li_init - n_Li,
            "Total lithium lost from particles [mol]": self.param.n_Li_particles_init
            - n_Li_particles,
            "Total lithium lost from electrolyte [mol]": self.param.n_Li_e_init
            - n_Li_e,
        }
    )

    # Lithium lost to side reactions
    # Different way of measuring LLI but should give same value
    n_Li_lost_neg_sei = self.variables["Loss of lithium to negative SEI [mol]"]
    n_Li_lost_pos_sei = self.variables["Loss of lithium to positive SEI [mol]"]
    n_Li_lost_reactions = n_Li_lost_neg_sei + n_Li_lost_pos_sei
    for domain in domains:
        dom = domain.split()[0].lower()
        n_Li_lost_sei_cracks = self.variables[
            f"Loss of lithium to {dom} SEI on cracks [mol]"
        ]
        n_Li_lost_pl = self.variables[
            f"Loss of lithium to {dom} lithium plating [mol]"
        ]
        n_Li_lost_reactions += n_Li_lost_sei_cracks + n_Li_lost_pl

    self.variables.update(
        {
            "Total lithium lost to side reactions [mol]": n_Li_lost_reactions,
            "Total capacity lost to side reactions [A.h]": n_Li_lost_reactions
            * self.param.F
            / 3600,
        }
    )

def set_default_summary_variables(self):
    """
    Sets the default summary variables.
    """
    summary_variables = [
        "Time [s]",
        "Time [h]",
        "Throughput capacity [A.h]",
        "Throughput energy [W.h]",
        # LAM, LLI
        "Loss of lithium inventory [%]",
        "Loss of lithium inventory, including electrolyte [%]",
        # Total lithium
        "Total lithium [mol]",
        "Total lithium in electrolyte [mol]",
        "Total lithium in particles [mol]",
        # Lithium lost
        "Total lithium lost [mol]",
        "Total lithium lost from particles [mol]",
        "Total lithium lost from electrolyte [mol]",
        "Loss of lithium to negative SEI [mol]",
        "Loss of capacity to negative SEI [A.h]",
        "Loss of lithium to positive SEI [mol]",
        "Loss of capacity to positive SEI [A.h]",
        "Total lithium lost to side reactions [mol]",
        "Total capacity lost to side reactions [A.h]",
        # Resistance
        "Local ECM resistance [Ohm]",
    ]

    if self.options.electrode_types["negative"] == "porous":
        summary_variables += [
            "Negative electrode capacity [A.h]",
            "Loss of active material in negative electrode [%]",
            "Total lithium in negative electrode [mol]",
            "Loss of lithium to negative lithium plating [mol]",
            "Loss of capacity to negative lithium plating [A.h]",
            "Loss of lithium to negative SEI on cracks [mol]",
            "Loss of capacity to negative SEI on cracks [A.h]",
        ]
    if self.options.electrode_types["positive"] == "porous":
        summary_variables += [
            "Positive electrode capacity [A.h]",
            "Loss of active material in positive electrode [%]",
            "Total lithium in positive electrode [mol]",
            "Loss of lithium to positive lithium plating [mol]",
            "Loss of capacity to positive lithium plating [A.h]",
            "Loss of lithium to positive SEI on cracks [mol]",
            "Loss of capacity to positive SEI on cracks [A.h]",
        ]

    self.summary_variables = summary_variables

def set_open_circuit_potential_submodel(self):
    for domain in ["negative", "positive"]:
        if self.options.electrode_types[domain] == "porous":
            reaction = "lithium-ion main"
        elif self.options.electrode_types[domain] == "planar":
            reaction = "lithium metal plating"
        domain_options = getattr(self.options, domain)
        for phase in self.options.phases[domain]:
            ocp_option = getattr(domain_options, phase)["open-circuit potential"]
            ocp_submodels = pybamm.open_circuit_potential
            if ocp_option == "single":
                ocp_model = ocp_submodels.SingleOpenCircuitPotential
            elif ocp_option == "current sigmoid":
                ocp_model = ocp_submodels.CurrentSigmoidOpenCircuitPotential
            elif ocp_option == "Wycisk":
                pybamm.citations.register("Wycisk2022")
                ocp_model = ocp_submodels.WyciskOpenCircuitPotential
            elif ocp_option == "MSMR":
                ocp_model = ocp_submodels.MSMROpenCircuitPotential
            self.submodels[f"{domain} {phase} open-circuit potential"] = ocp_model(
                self.param, domain, reaction, self.options, phase
            )

def set_sei_submodel(self):
    for domain in ["negative", "positive"]:
        if self.options.electrode_types[domain] == "planar":
            reaction_loc = "interface"
        elif self.options["x-average side reactions"] == "true":
            reaction_loc = "x-average"
        else:
            reaction_loc = "full electrode"
        phases = self.options.phases[domain]
        for phase in phases:
            sei_option = getattr(getattr(self.options, domain), phase)["SEI"]
            if sei_option == "none":
                submodel = pybamm.sei.NoSEI(self.param, domain, self.options, phase)
            elif sei_option == "constant":
                submodel = pybamm.sei.ConstantSEI(
                    self.param, domain, self.options, phase
                )
            else:
                submodel = pybamm.sei.SEIGrowth(
                    self.param,
                    domain,
                    reaction_loc,
                    self.options,
                    phase,
                    cracks=False,
                )
            self.submodels[f"{domain} {phase} sei"] = submodel
        if len(phases) > 1:
            self.submodels[f"{domain} total sei"] = pybamm.sei.TotalSEI(
                self.param, domain, self.options
            )

def set_sei_on_cracks_submodel(self):
    # Do not set "sei on cracks" submodel for a planar electrode. For porous
    # electrodes, "sei on cracks" submodel must be set, even if it is zero
    for domain in self.options.whole_cell_domains:
        if domain != "separator":
            domain = domain.split()[0].lower()
            sei_option = getattr(self.options, domain)["SEI"]
            sei_on_cracks_option = getattr(self.options, domain)["SEI on cracks"]
            phases = self.options.phases[domain]
            for phase in phases:
                if (
                    sei_option in ["none", "constant"]
                    or sei_on_cracks_option == "false"
                ):
                    submodel = pybamm.sei.NoSEI(
                        self.param, domain, self.options, phase, cracks=True
                    )
                else:
                    if self.options["x-average side reactions"] == "true":
                        reaction_loc = "x-average"
                    else:
                        reaction_loc = "full electrode"
                    submodel = pybamm.sei.SEIGrowth(
                        self.param,
                        domain,
                        reaction_loc,
                        self.options,
                        phase,
                        cracks=True,
                    )
                self.submodels[f"{domain} {phase} sei on cracks"] = submodel
            if len(phases) > 1:
                self.submodels[f"{domain} total sei on cracks"] = (
                    pybamm.sei.TotalSEI(
                        self.param, domain, self.options, cracks=True
                    )
                )

def set_lithium_plating_submodel(self):
    # Do not set "lithium plating" submodel for a planar electrode. For porous
    # electrodes, "lithium plating" submodel must be set, even if it is zero
    for domain in self.options.whole_cell_domains:
        if domain != "separator":
            domain = domain.split()[0].lower()
            phases = self.options.phases[domain]
            for phase in phases:
                lithium_plating_opt = getattr(getattr(self.options, domain), phase)[
                    "lithium plating"
                ]
                if lithium_plating_opt == "none":
                    submodel = pybamm.lithium_plating.NoPlating(
                        self.param, domain, self.options, phase
                    )
                else:
                    x_average = self.options["x-average side reactions"] == "true"
                    submodel = pybamm.lithium_plating.Plating(
                        self.param, domain, x_average, self.options, phase
                    )
                self.submodels[f"{domain} {phase} lithium plating"] = submodel
            if len(phases) > 1:
                self.submodels[f"{domain} total lithium plating"] = (
                    pybamm.lithium_plating.TotalLithiumPlating(
                        self.param, domain, self.options
                    )
                )

def set_total_interface_submodel(self):
    self.submodels["total interface"] = pybamm.interface.TotalInterfacialCurrent(
        self.param, "lithium-ion", self.options
    )

def set_crack_submodel(self):
    for domain in self.options.whole_cell_domains:
        if domain != "separator":
            domain = domain.split()[0].lower()
            phases = self.options.phases[domain]
            for phase in phases:
                crack = getattr(self.options, domain)["particle mechanics"]
                if crack == "none":
                    self.submodels[f"{domain} {phase}particle mechanics"] = (
                        pybamm.particle_mechanics.NoMechanics(
                            self.param, domain, options=self.options, phase=phase
                        )
                    )
                elif crack == "swelling only":
                    self.submodels[f"{domain} {phase}particle mechanics"] = (
                        pybamm.particle_mechanics.SwellingOnly(
                            self.param, domain, options=self.options, phase=phase
                        )
                    )
                elif crack == "swelling and cracking":
                    self.submodels[f"{domain} {phase}particle mechanics"] = (
                        pybamm.particle_mechanics.CrackPropagation(
                            self.param,
                            domain,
                            self.x_average,
                            options=self.options,
                            phase=phase,
                        )
                    )

def set_active_material_submodel(self):
    for domain in ["negative", "positive"]:
        if self.options.electrode_types[domain] == "porous":
            lam = getattr(self.options, domain)["loss of active material"]
            phases = self.options.phases[domain]
            for phase in phases:
                if lam == "none":
                    submod = pybamm.active_material.Constant(
                        self.param, domain, self.options, phase
                    )
                else:
                    submod = pybamm.active_material.LossActiveMaterial(
                        self.param, domain, self.options, self.x_average, phase
                    )
                self.submodels[f"{domain} {phase} active material"] = submod

            # Submodel for the total active material, summing up each phase
            if len(phases) > 1:
                self.submodels[f"{domain} total active material"] = (
                    pybamm.active_material.Total(self.param, domain, self.options)
                )

def set_porosity_submodel(self):
    if (
        self.options["SEI porosity change"] == "false"
        and self.options["lithium plating porosity change"] == "false"
    ):
        self.submodels["porosity"] = pybamm.porosity.Constant(
            self.param, self.options
        )
    elif (
        self.options["SEI porosity change"] == "true"
        or self.options["lithium plating porosity change"] == "true"
    ):
        x_average = self.options["x-average side reactions"] == "true"
        self.submodels["porosity"] = pybamm.porosity.ReactionDriven(
            self.param, self.options, x_average
        )

def set_li_metal_counter_electrode_submodels(self):
    for domain in ["negative", "positive"]:
        if self.options.electrode_types[domain] == "porous":
            continue
        if (
            self.options["SEI"] in ["none", "constant"]
            and self.options["intercalation kinetics"] == "symmetric Butler-Volmer"
            and self.options["surface form"] == "false"
        ):
            # only symmetric Butler-Volmer can be inverted
            self.submodels[f"{domain} electrode potential"] = (
                pybamm.electrode.ohm.LithiumMetalExplicit(
                    self.param, domain, self.options
                )
            )
            self.submodels[f"{domain} electrode interface"] = (
                pybamm.kinetics.InverseButlerVolmer(
                    self.param, domain, "lithium metal plating", self.options
                )
            )  # assuming symmetric reaction for now so we can take the inverse
            self.submodels[f"{domain} electrode interface current"] = (
                pybamm.kinetics.CurrentForInverseButlerVolmerLithiumMetal(
                    self.param, domain, "lithium metal plating", self.options
                )
            )
        else:
            self.submodels[f"{domain} electrode potential"] = (
                pybamm.electrode.ohm.LithiumMetalSurfaceForm(
                    self.param, domain, self.options
                )
            )
            neg_intercalation_kinetics = self.get_intercalation_kinetics(domain)
            self.submodels[f"{domain} electrode interface"] = (
                neg_intercalation_kinetics(
                    self.param,
                    domain,
                    "lithium metal plating",
                    self.options,
                    "primary",
                )
            )

def set_convection_submodel(self):
    self.submodels["transverse convection"] = (
        pybamm.convection.transverse.NoConvection(self.param, self.options)
    )
    self.submodels["through-cell convection"] = (
        pybamm.convection.through_cell.NoConvection(self.param, self.options)
    )

def insert_reference_electrode(self, position=None):
    """
    Insert a reference electrode to measure the electrolyte potential at a given
    position in space. Adds model variables for the electrolyte potential at the
    reference electrode and for the potential difference between the electrode
    potentials measured at the electrode/current collector interface and the
    reference electrode. Only implemented for 1D models (i.e. where the
    'dimensionality' option is 0).

    Parameters
    ----------
    position : :class:`pybamm.Symbol`, optional
        The position in space at which to measure the electrolyte potential. If
        None, defaults to the mid-point of the separator.
    """
    if self.options["dimensionality"] != 0:
        raise NotImplementedError(
            "Reference electrode can only be inserted for models where "
            "'dimensionality' is 0. For other models, please add a reference "
            "electrode manually."
        )

    if position is None:
        position = self.param.n.L + self.param.s.L / 2

    phi_e_ref = pybamm.EvaluateAt(
        self.variables["Electrolyte potential [V]"], position
    )
    phi_p = pybamm.boundary_value(
        self.variables["Positive electrode potential [V]"], "right"
    )
    variables = {
        "Positive electrode 3E potential [V]": phi_p - phi_e_ref,
        "Reference electrode potential [V]": phi_e_ref,
    }

    if self.options["working electrode"] == "both":
        phi_n = pybamm.boundary_value(
            self.variables["Negative electrode potential [V]"], "left"
        )
        variables.update(
            {
                "Negative electrode 3E potential [V]": phi_n - phi_e_ref,
            }
        )
    self.variables.update(variables)

class BasicDFN(BaseModel):
"""Doyle-Fuller-Newman (DFN) model of a lithium-ion battery, from
:footcite:t:Marquis2019.

This class differs from the :class:`pybamm.lithium_ion.DFN` model class in that it
shows the whole model in a single class. This comes at the cost of flexibility in
comparing different physical effects, and in general the main DFN class should be
used instead.

Parameters
----------
name : str, optional
    The name of the model.

"""

def __init__(self, name="Doyle-Fuller-Newman model"):
    super().__init__(name=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.

    ######################
    # 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_n = pybamm.Variable(
        "Negative electrolyte concentration [mol.m-3]",
        domain="negative electrode",
    )
    c_e_s = pybamm.Variable(
        "Separator electrolyte concentration [mol.m-3]",
        domain="separator",
    )
    c_e_p = pybamm.Variable(
        "Positive electrolyte concentration [mol.m-3]",
        domain="positive electrode",
    )
    # Concatenations combine several variables into a single variable, to simplify
    # implementing equations that hold over several domains
    c_e = pybamm.concatenation(c_e_n, c_e_s, c_e_p)

    # Electrolyte potential
    phi_e_n = pybamm.Variable(
        "Negative electrolyte potential [V]",
        domain="negative electrode",
    )
    phi_e_s = pybamm.Variable(
        "Separator electrolyte potential [V]",
        domain="separator",
    )
    phi_e_p = pybamm.Variable(
        "Positive electrolyte potential [V]",
        domain="positive electrode",
    )
    phi_e = pybamm.concatenation(phi_e_n, phi_e_s, phi_e_p)

    # Electrode potential
    phi_s_n = pybamm.Variable(
        "Negative electrode potential [V]", domain="negative electrode"
    )
    phi_s_p = pybamm.Variable(
        "Positive electrode potential [V]",
        domain="positive electrode",
    )
    # Particle concentrations are variables on the particle domain, but also vary in
    # the x-direction (electrode domain) and so must be provided with auxiliary
    # domains
    c_s_n = pybamm.Variable(
        "Negative particle concentration [mol.m-3]",
        domain="negative particle",
        auxiliary_domains={"secondary": "negative electrode"},
    )
    c_s_p = pybamm.Variable(
        "Positive particle concentration [mol.m-3]",
        domain="positive particle",
        auxiliary_domains={"secondary": "positive electrode"},
    )

    # Constant temperature
    T = self.param.T_init

    ######################
    # Other set-up
    ######################

    # Current density
    i_cell = self.param.current_density_with_time

    # Porosity
    # Primary broadcasts are used to broadcast scalar quantities across a domain
    # into a vector of the right shape, for multiplying with other vectors
    eps_n = pybamm.PrimaryBroadcast(
        pybamm.Parameter("Negative electrode porosity"), "negative electrode"
    )
    eps_s = pybamm.PrimaryBroadcast(
        pybamm.Parameter("Separator porosity"), "separator"
    )
    eps_p = pybamm.PrimaryBroadcast(
        pybamm.Parameter("Positive electrode porosity"), "positive electrode"
    )
    eps = pybamm.concatenation(eps_n, eps_s, eps_p)

    # Active material volume fraction (eps + eps_s + eps_inactive = 1)
    eps_s_n = pybamm.Parameter("Negative electrode active material volume fraction")
    eps_s_p = pybamm.Parameter("Positive electrode active material volume fraction")

    # transport_efficiency
    tor = pybamm.concatenation(
        eps_n**self.param.n.b_e, eps_s**self.param.s.b_e, eps_p**self.param.p.b_e)
    
    a_n = 3 * eps_s_n/ self.param.n.prim.R_typ
    a_p = 3 * eps_s_p / self.param.p.prim.R_typ
    #a_p =10e4
    #a_n=10e4

    # Interfacial reactions
    # Surf takes the surface value of a variable, i.e. its boundary value on the
    # right side. This is also accessible via `boundary_value(x, "right")`, with
    # "left" providing the boundary value of the left side
    c_s_surf_n = pybamm.surf(c_s_n)
    sto_surf_n = c_s_surf_n / self.param.n.prim.c_max
    j0_n = self.param.n.prim.j0(c_e_n, c_s_surf_n, T)
    eta_n = phi_s_n - phi_e_n - self.param.n.prim.U(sto_surf_n, T )
    Feta_RT_n = self.param.F * eta_n / (self.param.R * T)
    j_n = 2 * j0_n * pybamm.sinh(self.param.n.prim.ne / 2 * Feta_RT_n)

    c_s_surf_p = pybamm.surf(c_s_p)
    sto_surf_p = c_s_surf_p / self.param.p.prim.c_max
    j0_p = self.param.p.prim.j0(c_e_p, c_s_surf_p, T)
    eta_p = phi_s_p - phi_e_p - self.param.p.prim.U(sto_surf_p, T)
    Feta_RT_p = self.param.F * eta_p / (self.param.R * T)
    j_s = pybamm.PrimaryBroadcast(0, "separator")
    j_p = 2 * j0_p * pybamm.sinh(self.param.p.prim.ne / 2 * Feta_RT_p)

    a_j_n = a_n * j_n
    a_j_p = a_p * j_p
    a_j = pybamm.concatenation(a_j_n, j_s, a_j_p)

    ######################
    # 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
    N_s_n = -self.param.n.prim.D(c_s_n, T) * pybamm.grad(c_s_n)
    N_s_p = -self.param.p.prim.D(c_s_p, T) * pybamm.grad(c_s_p)
    self.rhs[c_s_n] = -pybamm.div(N_s_n)
    self.rhs[c_s_p] = -pybamm.div(N_s_p)
    # Boundary conditions must be provided for equations with spatial derivatives
    self.boundary_conditions[c_s_n] = {
        "left": (pybamm.Scalar(0), "Neumann"),
        "right": (
            -j_n / (self.param.F * pybamm.surf(self.param.n.prim.D(c_s_n, T))),
            "Neumann",
        ),
    }
    self.boundary_conditions[c_s_p] = {
        "left": (pybamm.Scalar(0), "Neumann"),
        "right": (
            -j_p / (self.param.F * pybamm.surf(self.param.p.prim.D(c_s_p, T))),
            "Neumann",
        ),
    }
    self.initial_conditions[c_s_n] = self.param.n.prim.c_init
    self.initial_conditions[c_s_p] = self.param.p.prim.c_init
    ######################
    # Current in the solid
    ######################
    sigma_eff_n = 1
    i_s_n = -sigma_eff_n * pybamm.grad(phi_s_n)
    sigma_eff_p = self.param.p.sigma(T) * eps_s_p**self.param.p.b_s
    i_s_p = -sigma_eff_p * pybamm.grad(phi_s_p)
    # The `algebraic` dictionary contains differential equations, with the key being
    # the main scalar variable of interest in the equation
    # multiply by Lx**2 to improve conditioning
    self.algebraic[phi_s_n] = self.param.L_x**2 * (pybamm.div(i_s_n) + a_j_n)
    self.algebraic[phi_s_p] = self.param.L_x**2 * (pybamm.div(i_s_p) + a_j_p)
    self.boundary_conditions[phi_s_n] = {
        "left": (pybamm.Scalar(0), "Dirichlet"),
        "right": (pybamm.Scalar(0), "Neumann"),
    }
    self.boundary_conditions[phi_s_p] = {
        "left": (pybamm.Scalar(0), "Neumann"),
        "right": (i_cell / pybamm.boundary_value(-sigma_eff_p, "right"), "Neumann"),
    }
    # 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_n] = pybamm.Scalar(0)
    self.initial_conditions[phi_s_p] = self.param.ocv_init

    ######################
    # Current in the electrolyte
    ######################
    i_e = (0.001) * (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)
    self.boundary_conditions[phi_e] = {
        "left": (pybamm.Scalar(0), "Neumann"),
        "right": (pybamm.Scalar(0), "Neumann"),
    }
    self.initial_conditions[phi_e] = -self.param.n.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
    )
    self.boundary_conditions[c_e] = {
        "left": (pybamm.Scalar(0), "Neumann"),
        "right": (pybamm.Scalar(0), "Neumann"),
    }
    self.initial_conditions[c_e] = self.param.c_e_init

    ######################
    # (Some) variables
    ######################
    voltage = pybamm.boundary_value(phi_s_p, "right")
    num_cells = pybamm.Parameter(
        "Number of cells connected in series to make a battery"
    )
    # The `variables` dictionary contains all variables that might be useful for
    # visualising the solution of the model
    self.variables = {
        "Negative particle concentration [mol.m-3]": c_s_n,
        "Negative particle surface concentration [mol.m-3]": c_s_surf_n,
        "Electrolyte concentration [mol.m-3]": c_e,
        "Negative electrolyte concentration [mol.m-3]": c_e_n,
        "Separator electrolyte concentration [mol.m-3]": c_e_s,
        "Positive electrolyte concentration [mol.m-3]": c_e_p,
        "Positive particle concentration [mol.m-3]": c_s_p,
        "Positive particle surface concentration [mol.m-3]": c_s_surf_p,
        "Current [A]": I,
        "Current variable [A]": I,  # for compatibility with pybamm.Experiment
        "Negative electrode potential [V]": phi_s_n,
        "Electrolyte potential [V]": phi_e,
        "Negative electrolyte potential [V]": phi_e_n,
        "Separator electrolyte potential [V]": phi_e_s,
        "Positive electrolyte potential [V]": phi_e_p,
        "Positive electrode potential [V]": phi_s_p,
        "Voltage [V]": voltage,
        "Battery voltage [V]": voltage * num_cells,
        "Time [s]": pybamm.t,
        "Discharge capacity [A.h]": Q,
    }
    # Events specify points at which a solution should terminate
    self.events += [
        pybamm.Event("Minimum voltage [V]", voltage - self.param.voltage_low_cut),
        pybamm.Event("Maximum voltage [V]", self.param.voltage_high_cut - voltage),
    ]

model =BasicDFN()

param = pybamm.ParameterValues("Chen2020")

import numpy as np

def a(sto):
u_eq = (
13.4905 - 10.96038 * sto +
8.203617 * sto1.358699 -
3.10758e-6 * np.exp(127.1216 * sto - 114.2593) -
7.033556 * sto(-0.03362749)
)

return u_eq

def b(sto):

u_eq = 1

return u_eq

def c(c_e, c_s_surf, c_s_max, T):
m_ref = 6 * 10 ** (-7) # (A/m2)(m3/mol)**1.5 - includes ref concentrations
E_r = 39570
arrhenius = np.exp(E_r / pybamm.constants.R * (1 / 298.15 - 1 / T))

return 96500 * 2.57e-11* 1000**0.5 * c_s_surf**0.5 * (c_s_max - c_s_surf) ** 0.5

def d(c_e, c_s_surf, c_s_max, T):
m_ref = 6 * 10 ** (-7) # (A/m2)(m3/mol)**1.5 - includes ref concentrations
E_r = 39570
arrhenius = np.exp(E_r / pybamm.constants.R * (1 / 298.15 - 1 / T))

return 96500 * 2.57e-11 * 1000**0.5 * c_s_surf**0.5 * (c_s_max - c_s_surf) ** 0.5

param.update(
{
"Negative electrode OCP [V]": b,
"Negative electrode exchange-current density [A.m-2]": d,
"Positive electrode OCP [V]": a,
"Positive electrode exchange-current density [A.m-2]": c,
# cell
"Separator thickness [m]": 4e-05,
"Positive electrode thickness [m]": 7e-05,
"Positive current collector thickness [m]": 1.6e-05,
"Electrode height [m]": 0.02,
"Electrode width [m]": 1,
"Cell cooling surface area [m2]": 0.00531,
"Cell volume [m3]": 2.42e-5,
"Cell thermal expansion coefficient [m.K-1]": 1.1e-06,
"Negative current collector conductivity [S.m-1]": 58411000.0,
"Maximum concentration in negative electrode [mol.m-3]": 33133.0,
"Maximum concentration in positive electrode [mol.m-3]": 52500.0,
"Positive current collector conductivity [S.m-1]": 36914000.0,
"Negative current collector density [kg.m-3]": 8960.0,
"Positive current collector density [kg.m-3]": 2700.0,
"Negative current collector specific heat capacity [J.kg-1.K-1]": 385.0,
"Positive current collector specific heat capacity [J.kg-1.K-1]": 897.0,
"Negative current collector thermal conductivity [W.m-1.K-1]": 401.0,
"Positive current collector thermal conductivity [W.m-1.K-1]": 237.0,
"Nominal cell capacity [A.h]": 1,
"Current function [A]": 0.65,
"Contact resistance [Ohm]": 0,
# negative electrode
"Negative electrode conductivity [S.m-1]": 215.0,
"Negative particle diffusivity [m2.s-1]": 3.3e-14,
"Negative electrode porosity": 0.25,
"Negative electrode active material volume fraction": 0.8,
"Negative particle radius [m]": 5.86e-06,
"Negative electrode Bruggeman coefficient (electrolyte)": 1.5,
"Negative electrode Bruggeman coefficient (electrode)": 0,
"Negative electrode charge transfer coefficient": 0.5,
"Negative electrode double-layer capacity [F.m-2]": 0.2,
"Negative electrode density [kg.m-3]": 1657.0,
"Negative electrode specific heat capacity [J.kg-1.K-1]": 700.0,
"Negative electrode thermal conductivity [W.m-1.K-1]": 1.7,
"Negative electrode OCP entropic change [V.K-1]": 0.0,
# positive electrode
"Positive electrode conductivity [S.m-1]": 0.18,
"Positive particle diffusivity [m2.s-1]": 4e-15,
"Positive electrode porosity": 0.335,
"Positive electrode active material volume fraction": 0.6,
"Positive particle radius [m]": 5e-06,
"Positive electrode Bruggeman coefficient (electrolyte)": 1.5,
"Positive electrode Bruggeman coefficient (electrode)": 0,
"Positive electrode charge transfer coefficient": 0.5,
"Positive electrode double-layer capacity [F.m-2]": 0.2,
"Positive electrode density [kg.m-3]": 3262.0,
"Positive electrode specific heat capacity [J.kg-1.K-1]": 700.0,
"Positive electrode thermal conductivity [W.m-1.K-1]": 2.1,
"Positive electrode OCP entropic change [V.K-1]": 0.0,
# separator
"Separator porosity": 0.47,
"Separator Bruggeman coefficient (electrolyte)": 1.5,
"Separator density [kg.m-3]": 397.0,
"Separator specific heat capacity [J.kg-1.K-1]": 700.0,
"Separator thermal conductivity [W.m-1.K-1]": 0.16,
# electrolyte
"Initial concentration in electrolyte [mol.m-3]": 1000.0,
"Cation transference number": 0.5,
"Thermodynamic factor": 1.0, # experiment
"Reference temperature [K]": 298.15,
"Total heat transfer coefficient [W.m-2.K-1]": 10.0,
"Ambient temperature [K]": 298.15,
"Number of electrodes connected in parallel to make a cell": 1.0,
"Number of cells connected in series to make a battery": 1.0,
"Lower voltage cut-off [V]": 1,
"Upper voltage cut-off [V]": 4,
"Open-circuit voltage at 0% SOC [V]": 1,
"Open-circuit voltage at 100% SOC [V]": 4,
"Initial concentration in negative electrode [mol.m-3]": 29866.0,
"Initial concentration in positive electrode [mol.m-3]": 0.35*52500,
"Initial temperature [K]": 298.15,
"Lower voltage cut-off [V]": 1,
"Open-circuit voltage at 0% SOC [V]": 1
}
)

experiment = pybamm.Experiment(
[
("Discharge at 1.4 A until 2.2 V")
]
)

sim = pybamm.Simulation(model, experiment=experiment, parameter_values=param)

sim.solve()
sim.plot()