Unable to solve the model after Discretisation

[karkulali28]

I am trying to solve single particle model using PyBamm and when I simulate the model, I get the following error message. Any idea on how to solve the issue.

KeyError Traceback (most recent call last)
Cell In [189], line 1
----> 1 disc.process_model(model)

File ~\AppData\Local\Programs\Python\Python310\lib\site-packages\pybamm\discretisations\discretisation.py:206, in Discretisation.process_model(self, model, inplace, check_model, remove_independent_variables_from_rhs)
203 model_disc.bcs = self.bcs
205 pybamm.logger.verbose("Discretise initial conditions for {}".format(model.name))
--> 206 ics, concat_ics = self.process_initial_conditions(model)
207 model_disc.initial_conditions = ics
208 model_disc.concatenated_initial_conditions = concat_ics

File ~\AppData\Local\Programs\Python\Python310\lib\site-packages\pybamm\discretisations\discretisation.py:425, in Discretisation.process_initial_conditions(self, model)
419 processed_initial_conditions = self.process_dict(
420 model.initial_conditions, ics=True
421 )
423 # Concatenate initial conditions into a single vector
424 # check that all initial conditions are set
--> 425 processed_concatenated_initial_conditions = self._concatenate_in_order(
426 processed_initial_conditions, check_complete=True
427 )
429 return processed_initial_conditions, processed_concatenated_initial_conditions

File ~\AppData\Local\Programs\Python\Python310\lib\site-packages\pybamm\discretisations\discretisation.py:984, in Discretisation._concatenate_in_order(self, var_eqn_dict, check_complete, sparse)
982 else:
983 unpacked_variables.append(symbol)
--> 984 slices.append(self.y_slices[symbol][0])
986 if check_complete:
987 # Check keys from the given var_eqn_dict against self.y_slices
988 unpacked_variables_set = set(unpacked_variables)

KeyError: Variable(-0x7786521cf4ca09d6, Negative electrode potential [V], children=[], domains={'primary': ['negative electrode']})

[brosaplanella]

Can you share a minimum working example to reproduce this error?

[karkulali28]

I have added the code below for your reference.. Do let me know if you are able to run it..

[karkulali28]

import numpy as np
import pybamm
import matplotlib.pyplot as plt

Define the parameters

Q = pybamm.Variable("Discharge capacity [A.h]")
c_s_n = pybamm.Variable("Negative particle concentration [mol.m-3]", domain="negative particle")
c_s_p = pybamm.Variable("Positive particle concentration [mol.m-3]", domain="positive particle")

D = pybamm.Parameter("Diffusion coefficient [m2.s-1]")
I_app = pybamm.InputParameter("Applied current [A]")

F = pybamm.Parameter("Faraday constant [C.mol-1]")
R = pybamm.Parameter("Molar gas constant [J.mol-1.K-1]")
T = pybamm.Parameter("Temperature [K]")

a = pybamm.Parameter("Surface area per unit volume [m-1]")
L_n = pybamm.Parameter("Negative electrode thickness [m]")
L_p = pybamm.Parameter("Positive electrode thickness [m]")
L_s = pybamm.Parameter("Separator thickness [m]")
R_n = pybamm.Parameter("Negative Particle Radius [m]")
R_p = pybamm.Parameter("Positive Particle Radius [m]")

#j = pybamm.InputParameter("Interfacial current density [A.m-2]")
c0_n = pybamm.Parameter("Initial concentration in negative electrode [mol.m-3]")
c0_p = pybamm.Parameter("Initial concentration in positive electrode [mol.m-3]")

sigma_p = pybamm.Parameter("Positive electrode conductivity [S.m-1]")
sigma_n = pybamm.Parameter("Negative electrode conductivity [S.m-1]")
kappa = pybamm.Parameter("Electrolyte conductivity [S.m-1]")

phi_n = pybamm.Variable("Negative electrode potential [V]",domain ="negative electrode")
phi_p = pybamm.Variable("Positive electrode potential [V]",domain ="positive electrode")

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_s, phi_e_p)

c_surf_p = pybamm.surf(c_s_p) # Positive surface concentration
c_surf_n = pybamm.surf(c_s_n) # Negative surface concentration

input_n = {"Negative particle surface concentration [mol.m-3]": c_surf_n}
input_p = {"Positive particle surface concentration [mol.m-3]": c_surf_p}

j0_n = pybamm.FunctionParameter("Negative electrode exchange-current density [A.m-2]",input_n)
j0_p = pybamm.FunctionParameter("Positive electrode exchange-current density [A.m-2]",input_p)

U_n = pybamm.FunctionParameter("Negative electrode OCP [V]", input_n)
U_p = pybamm.FunctionParameter("Positive electrode OCP [V]", input_p)

j_n = 2 * j0_n * pybamm.sinh((F / 2 / R / T) * (phi_n - phi_e_n - U_n))
j_s = pybamm.PrimaryBroadcast(0, "separator")
j_p = 2 * j0_p * pybamm.sinh((F / 2 / R / T) * (phi_p - phi_e_p - U_p))
j = pybamm.concatenation(j_s,j_p)

model = pybamm.BaseModel()

charge conservation equations

i_n = -sigma_n * pybamm.grad(phi_n)
i_e = -kappa * pybamm.grad(phi_e)
model.algebraic = {
phi_n: pybamm.div(i_n) + a * j_n,
phi_e: pybamm.div(i_e) - a * j,
}

#Govering Equations - Negative electrode
N_s_n = -D * pybamm.grad(c_s_n) # flux
dcndt = -pybamm.div(N_s_n)

charge conservation equations

i_p = -sigma_p * pybamm.grad(phi_p)
i_e = -kappa * pybamm.grad(phi_e)
model.algebraic = {
phi_p: pybamm.div(i_p) + a * j_p,
phi_e: pybamm.div(i_e) - a * j,
}

#Govering Equations - Positive electrode
N_s_p = -D * pybamm.grad(c_s_p) # flux
dcpdt = -pybamm.div(N_s_p)

model.rhs = {c_s_n: dcndt,c_s_p: dcpdt}
#Boundary Conditions - Negative electrode

lbc_n = pybamm.Scalar(0)
rbc_n = -I_app / (L_n * a)

#Boundary Conditions - Positive electrode

lbc_p = pybamm.Scalar(0)
rbc_p = I_app / (L_p * a)

model.boundary_conditions = {phi_n: {"left": (pybamm.Scalar(0), "Neumann"), "right": (-I_app / a / sigma_n, "Neumann")}, phi_e: {"left": (pybamm.Scalar(0), "Dirichlet"), "right": (pybamm.Scalar(0), "Neumann")}, phi_p: {"left": (pybamm.Scalar(0), "Neumann"), "right": (I_app / a / sigma_p, "Neumann")}, c_s_n: {"left": (lbc_n, "Neumann"), "right": (rbc_n, "Neumann")},
c_s_p: {"left": (lbc_p, "Neumann"), "right": (rbc_p, "Neumann")}}

input_p = {"Initial concentration [mol.m-3]": c0_p}
U_init_p = pybamm.FunctionParameter("Positive electrode OCP [V]", input_p)

input_n = {"Initial concentration [mol.m-3]": c0_n}
U_init_n = pybamm.FunctionParameter("Positive electrode OCP [V]", input_n)

model.initial_conditions = { phi_n: U_init_n, phi_p: U_init_p, c_s_n: c0_n, c_s_p: c0_p}

model.variables = {
"Negative particle concentration [mol.m-3]": c_s_n,
"Positive particle concentration [mol.m-3]": c_s_p,
"Negative particle surface concentration [mol.m-3]": pybamm.surf(c_s_n),
"Positive particle surface concentration [mol.m-3]": pybamm.surf(c_s_p),
"Positive electrode potential [V]": phi_p,
"Negative electrode potential [V]": phi_n,
"Electrolyte potential [V]": phi_e,
"Positive electrode interfacial current density [A.m-2]": j_p,
"Negative electrode interfacial current density [A.m-2]": j_n,
"Positive electrode OCP [V]":pybamm.boundary_value(U_p, "right"),
"Negative electrode OCP [V]":pybamm.boundary_value(U_n, "right"),
"Voltage [V]": pybamm.boundary_value(phi_p, "right"),
}

from pybamm import tanh

both functions will depend on the maximum concentration

c_max_p = pybamm.Parameter("Maximum concentration in positive electrode [mol.m-3]")
c_max_n = pybamm.Parameter("Maximum concentration in negative electrode [mol.m-3]")

def exchange_current_density_negative(c_surf_n):
k = 6 * 10 ** (-7) # reaction rate [(A/m2)(m3/mol)**1.5]
c_e = 1000 # (constant) electrolyte concentration [mol.m-3]

return k * c_e** 0.5 * c_surf_n ** 0.5 * (c_max_n - c_surf_n) ** 0.5

def exchange_current_density_positive(c_surf_p):
k = 6 * 10 ** (-7) # reaction rate [(A/m2)(m3/mol)**1.5]
c_e = 1000 # (constant) electrolyte concentration [mol.m-3]

return k * c_e** 0.5 * c_surf_p ** 0.5 * (c_max_p - c_surf_p) ** 0.5

def open_circuit_potential_negative(c_surf_n):
stretch = 1.062
sto = stretch * c_surf_n / c_max_n

u_eq_n = (
    2.16216
    + 0.07645 * tanh(30.834 - 54.4806 * sto)
    + 2.1581 * tanh(52.294 - 50.294 * sto)
    - 0.14169 * tanh(11.0923 - 19.8543 * sto)
    + 0.2051 * tanh(1.4684 - 5.4888 * sto)
    + 0.2531 * tanh((-sto + 0.56478) / 0.1316)
    - 0.02167 * tanh((sto - 0.525) / 0.006)
)
return u_eq_n   

def open_circuit_potential_positive(c_surf_p):
stretch = 1.062
sto = stretch * c_surf_p / c_max_p

u_eq_p = (
    6.16216
    + 0.07645 * tanh(30.834 - 54.4806 * sto)
    + 2.1581 * tanh(52.294 - 50.294 * sto)
    - 0.14169 * tanh(11.0923 - 19.8543 * sto)
    + 0.2051 * tanh(1.4684 - 5.4888 * sto)
    + 0.2531 * tanh((-sto + 0.56478) / 0.1316)
    - 0.02167 * tanh((sto - 0.525) / 0.006)
)
return u_eq_p

param = pybamm.ParameterValues(
{
"Initial concentration in negative electrode [mol.m-3]": 19986.609595075,
"Initial concentration in positive electrode [mol.m-3]": 30730.7554385565,
#"Initial concentration in positive electrode [mol.m-3]": 25370,
#"Initial concentration in negative electrode [mol.m-3]": 25370,

"Negative electrode thickness [m]": 100e-6,
"Positive electrode thickness [m]": 100e-6,
"Separator thickness [m]": 25e-6,

"Applied current [A]": 0.680616,
"Diffusion coefficient [m2.s-1]": 1e-13,
"Surface area per unit volume [m-1]": 0.15e6,

"Negative Particle Radius [m]": 1e-06,
"Positive Particle Radius [m]": 1e-06,

"Positive electrode conductivity [S.m-1]": 10,
"Negative electrode conductivity [S.m-1]": 10,
"Electrolyte conductivity [S.m-1]": 1,

"Molar gas constant [J.mol-1.K-1]": 8.314,
"Temperature [K]": 298.15,
"Faraday constant [C.mol-1]": 96485,

"Maximum concentration in positive electrode [mol.m-3]": 51217,
"Maximum concentration in negative electrode [mol.m-3]": 51217, 

"Positive electrode exchange-current density [A.m-2]": exchange_current_density_positive,
"Negative electrode exchange-current density [A.m-2]": exchange_current_density_negative,

"Positive electrode OCP [V]": open_circuit_potential_positive,
"Negative electrode OCP [V]": open_circuit_potential_negative,

})

r_n = pybamm.SpatialVariable("r_n", domain=["negative particle"], auxiliary_domains={
"secondary": "negative electrode"}, coord_sys="spherical polar")
r_p = pybamm.SpatialVariable("r_p",auxiliary_domains={"secondary": "positive electrode"}, domain=["positive particle"], coord_sys="spherical polar")
x_s = pybamm.SpatialVariable("x_s", domain=["separator"], coord_sys="cartesian")
x_n = pybamm.SpatialVariable("x_n", domain=["negative electrode"], coord_sys="cartesian")
x_p = pybamm.SpatialVariable("x_p", domain=["positive electrode"], coord_sys="cartesian")

geometry = {

"negative electrode": {x_n: {"min": 0, "max": L_n}},
"negative particle": {r_n: {"min": 0, "max": R_n}},
"separator": {x_s: {"min": -L_s, "max": 0}}, 
"positive electrode": {x_p: {"min": 0, "max": L_p}},
"positive particle": {r_p: {"min": 0, "max": R_p}},

}
param.process_model(model)
param.process_geometry(geometry)

submesh_types = {
"negative electrode": pybamm.Uniform1DSubMesh,
"negative particle": pybamm.Uniform1DSubMesh,
"separator": pybamm.Uniform1DSubMesh,
"positive electrode": pybamm.Uniform1DSubMesh,
"positive particle": pybamm.Uniform1DSubMesh,
}
var_pts = {x_n: 10, x_p: 20,x_s: 10, r_n: 20, r_p: 20}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
spatial_methods = {
"negative electrode": pybamm.FiniteVolume(),
"negative particle": pybamm.FiniteVolume(),
"separator": pybamm.FiniteVolume(),
"positive electrode": pybamm.FiniteVolume(),
"positive particle": pybamm.FiniteVolume(),
}

disc = pybamm.Discretisation(mesh,spatial_methods)
disc.process_model(model)