From 5a6310598c2c170525e7889ddc2c710709eb07e4 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Tue, 6 Jan 2026 20:31:47 +0000 Subject: [PATCH 01/37] Add GroupedSPMe voltage components --- pybop/models/lithium_ion/grouped_spme.py | 47 ++++++++++++++++++++++++ tests/unit/test_models.py | 8 ++++ 2 files changed, 55 insertions(+) diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index c39616105..c32062194 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -327,6 +327,51 @@ def __init__( Event("Maximum voltage [V]", self.param.voltage_high_cut - V), ] + ###################### + # Voltage components + ###################### + # Include the following variables to enable plotting via PyBaMM's plot_voltage_components + ocp_n_bulk = self.U( + pybamm.x_average(zeta_n * pybamm.r_average(sto_n)) + / pybamm.x_average(zeta_n), + "negative", + ) + ocp_p_bulk = self.U( + pybamm.x_average(zeta_p * pybamm.r_average(sto_p)) + / pybamm.x_average(zeta_p), + "positive", + ) + voltage_components = { + "Battery open-circuit voltage [V]": ocp_p_bulk - ocp_n_bulk, + "Battery particle concentration overpotential [V]": ( + (pybamm.x_average(self.U(sto_p_surf, "positive")) - ocp_p_bulk) + - (pybamm.x_average(self.U(sto_n_surf, "negative")) - ocp_n_bulk) + ), + "X-averaged battery reaction overpotential [V]": pybamm.x_average(eta_p) + - pybamm.x_average(eta_n), + "X-averaged battery concentration overpotential [V]": eta_e, + "X-averaged battery electrolyte ohmic losses [V]": Scalar(0), + "X-averaged battery solid phase ohmic losses [V]": Scalar(0), + "Contact overpotential [V]": R0 * I, # includes Ohmic losses in this model + # and split by electrode + "Negative electrode bulk open-circuit potential [V]": ocp_n_bulk, + "Positive electrode bulk open-circuit potential [V]": ocp_p_bulk, + "Negative particle concentration overpotential [V]": pybamm.x_average( + self.U(sto_n_surf, "negative") + ) + - ocp_n_bulk, + "Positive particle concentration overpotential [V]": pybamm.x_average( + self.U(sto_p_surf, "positive") + ) + - ocp_p_bulk, + "X-averaged negative electrode reaction overpotential [V]" + "": pybamm.x_average(eta_n), + "X-averaged positive electrode reaction overpotential [V]" + "": pybamm.x_average(eta_p), + "X-averaged battery negative solid phase ohmic losses [V]": Scalar(0), + "X-averaged battery positive solid phase ohmic losses [V]": Scalar(0), + } + ###################### # (Some) variables ###################### @@ -367,6 +412,7 @@ def __init__( - pybamm.log(pybamm.concatenation(sto_e_n, sto_e_sep, sto_e_p)) ), "Time [s]": pybamm_t, + "Time [h]": pybamm_t / 3600, "Current [A]": I, "Current variable [A]": I, # for compatibility with pybamm.Experiment "Discharge capacity [A.h]": Q, @@ -374,6 +420,7 @@ def __init__( "Voltage [V]": V, "Battery voltage [V]": V, "Open-circuit voltage [V]": U_p - U_n, + **voltage_components, } def U(self, sto, domain): diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index 097dff64b..3eeaa9fbf 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -1,6 +1,8 @@ import numpy as np import pybamm import pytest +from matplotlib.axes._axes import Axes +from matplotlib.figure import Figure import pybop @@ -45,3 +47,9 @@ def test_model_simulation(self, model): fig = solution.plot() assert isinstance(fig, pybamm.QuickPlot) + + if isinstance(model, pybop.lithium_ion.GroupedSPMe): + for split in [False, True]: + fig, ax = solution.plot_voltage_components(split_by_electrode=split) + assert isinstance(fig, Figure) + assert isinstance(ax, Axes) From 97d82e843ba00a2fe42cdd9d52163db2e75c7358 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Wed, 7 Jan 2026 13:54:46 +0000 Subject: [PATCH 02/37] Change parameter grouping --- pybop/models/lithium_ion/grouped_spme.py | 43 +++++++++++------------- 1 file changed, 20 insertions(+), 23 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index c32062194..803797466 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -56,7 +56,7 @@ def __init__( doi = {10.1149/1945-7111/add41b} } """ - ) + ) # Note that the electrode electrolyte timescales have been replaced by relative transport efficiencies ###################### # Variables @@ -155,9 +155,9 @@ def __init__( zeta_n = Parameter("Negative electrode relative porosity") zeta_p = Parameter("Positive electrode relative porosity") - tau_e_n = Parameter("Negative electrode electrolyte diffusion time scale [s]") - tau_e_sep = Parameter("Separator electrolyte diffusion time scale [s]") - tau_e_p = Parameter("Positive electrode electrolyte diffusion time scale [s]") + tau_e = Parameter("Electrolyte diffusion time scale [s]") + beta_n = Parameter("Negative electrode relative transport efficiency") + beta_p = Parameter("Positive electrode relative transport efficiency") ###################### # Input current (positive on discharge) @@ -273,35 +273,32 @@ def __init__( # Electrolyte ###################### self.rhs[sto_e_n] = ( - pybamm.div(pybamm.grad(sto_e_n) / tau_e_n - (t_plus * I / Q_e) * x_n / l_n) + pybamm.div( + pybamm.grad(sto_e_n) * beta_n / tau_e - (t_plus * I / Q_e) * x_n / l_n + ) + (3 / Q_e) * Q_th_n * j_n / l_n ) / zeta_n self.rhs[sto_e_sep] = pybamm.div( - pybamm.grad(sto_e_sep) / tau_e_sep - t_plus * I / Q_e + pybamm.grad(sto_e_sep) / tau_e - t_plus * I / Q_e ) self.rhs[sto_e_p] = ( pybamm.div( - pybamm.grad(sto_e_p) / tau_e_p - (t_plus * I / Q_e) * (1 - x_p) / l_p + pybamm.grad(sto_e_p) * beta_p / tau_e + - (t_plus * I / Q_e) * (1 - x_p) / l_p ) + (3 / Q_e) * Q_th_p * j_p / l_p ) / zeta_p self.boundary_conditions[sto_e_n] = { "left": (Scalar(0), "Neumann"), - "right": ( - tau_e_n * pybamm.boundary_gradient(sto_e_sep, "left") / tau_e_sep, - "Neumann", - ), + "right": (pybamm.boundary_gradient(sto_e_sep, "left") / beta_n, "Neumann"), } self.boundary_conditions[sto_e_sep] = { "left": (pybamm.boundary_value(sto_e_n, "right"), "Dirichlet"), "right": (pybamm.boundary_value(sto_e_p, "left"), "Dirichlet"), } self.boundary_conditions[sto_e_p] = { - "left": ( - tau_e_p * pybamm.boundary_gradient(sto_e_sep, "right") / tau_e_sep, - "Neumann", - ), + "left": (pybamm.boundary_gradient(sto_e_sep, "right") / beta_p, "Neumann"), "right": (Scalar(0), "Neumann"), } @@ -610,7 +607,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal epsilon_sep = param["Separator porosity"] b_sep = param["Separator Bruggeman coefficient (electrolyte)"] t_plus = param["Cation transference number"] - sigma_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) + kappa_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) # Compute the cell area and thickness A = param["Electrode height [m]"] * param["Electrode width [m]"] @@ -621,7 +618,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal L_p / (3 * epsilon_p**b_p) + L_s / (epsilon_sep**b_sep) + L_n / (3 * epsilon_n**b_n) - ) / (sigma_e * A) + ) / (kappa_e * A) Rs = (L_p / sigma_p + L_n / sigma_n) / (3 * A) R0 = Re + Rs + param["Contact resistance [Ohm]"] @@ -652,9 +649,9 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal tau_d_p = R_p**2 / D_p tau_d_n = R_n**2 / D_n - tau_e_p = epsilon_sep * L**2 / (epsilon_p**b_p * De) - tau_e_n = epsilon_sep * L**2 / (epsilon_n**b_n * De) - tau_e_sep = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) + tau_e = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) + beta_p = epsilon_p**b_p / epsilon_sep**b_sep + beta_n = epsilon_n**b_n / epsilon_sep**b_sep tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) tau_ct_n = F * R_n / (m_n * np.sqrt(ce0)) @@ -684,9 +681,9 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Negative electrode relative porosity": zeta_n, "Positive particle diffusion time scale [s]": tau_d_p, "Negative particle diffusion time scale [s]": tau_d_n, - "Positive electrode electrolyte diffusion time scale [s]": tau_e_p, - "Negative electrode electrolyte diffusion time scale [s]": tau_e_n, - "Separator electrolyte diffusion time scale [s]": tau_e_sep, + "Electrolyte diffusion time scale [s]": tau_e, + "Positive electrode relative transport efficiency": beta_p, + "Negative electrode relative transport efficiency": beta_n, "Positive electrode charge transfer time scale [s]": tau_ct_p, "Negative electrode charge transfer time scale [s]": tau_ct_n, "Positive electrode capacitance [F]": C_p, From 5e236e2ca4ef3975ef9ed54d55533e703efb8682 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Wed, 7 Jan 2026 14:59:03 +0000 Subject: [PATCH 03/37] Add minimum electrolyte events --- pybop/models/lithium_ion/grouped_spme.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 803797466..84cfb875b 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -115,6 +115,19 @@ def __init__( "Maximum positive particle surface stoichiometry", (1 - 0.01) - pybamm.max(sto_p_surf), ), + # model does not capture electrolyte depletion, use the DFN instead + Event( + "Minimum negative electrode electrolyte stoichiometry", + pybamm.min(sto_e_n) - 0, + ), + Event( + "Minimum separator electrolyte stoichiometry", + pybamm.min(sto_e_sep) - 0, + ), + Event( + "Minimum positive electrode electrolyte stoichiometry", + pybamm.min(sto_e_p) - 0, + ), ] ###################### From 9a56fa7def8502f5d1330f735672cd528fde3f89 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Wed, 7 Jan 2026 13:55:06 +0000 Subject: [PATCH 04/37] Add GroupedDFN --- pybop/models/lithium_ion/__init__.py | 1 + pybop/models/lithium_ion/grouped_dfn.py | 672 ++++++++++++++++++ .../integration/models/test_grouped_models.py | 1 + tests/unit/test_models.py | 1 + 4 files changed, 675 insertions(+) create mode 100644 pybop/models/lithium_ion/grouped_dfn.py diff --git a/pybop/models/lithium_ion/__init__.py b/pybop/models/lithium_ion/__init__.py index 0a1deb13f..8b7f39631 100644 --- a/pybop/models/lithium_ion/__init__.py +++ b/pybop/models/lithium_ion/__init__.py @@ -4,4 +4,5 @@ from .sp_diffusion import SPDiffusion from .grouped_spm import GroupedSPM from .grouped_spme import GroupedSPMe +from .grouped_dfn import GroupedDFN from .weppner_huggins import WeppnerHuggins diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py new file mode 100644 index 000000000..a90fb0705 --- /dev/null +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -0,0 +1,672 @@ +import numpy as np +import pybamm +from pybamm import ( + Event, + FunctionParameter, + Parameter, + ParameterValues, + PrimaryBroadcast, + Scalar, + Variable, +) +from pybamm import lithium_ion as pybamm_lithium_ion +from pybamm import t as pybamm_t +from pybamm.models.full_battery_models.lithium_ion.electrode_soh import ( + get_min_max_stoichiometries, +) + + +class GroupedDFN(pybamm_lithium_ion.BaseModel): + """ + A grouped parameter version of the Doyle Fuller Newman (DFN) model. + + Parameters + ---------- + name : str, optional + The name of the model. + **model_kwargs : optional + Valid PyBaMM model option keys and their values, for example: + options : dict, optional + A dictionary of options to customise the behaviour of the PyBaMM model. + build : bool, optional + If True, the model is built upon creation (default: False). + """ + + def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): + super().__init__(name=name, **model_kwargs) + + pybamm.citations.register( + """ + @article{Hallemans2025, + title = {{Physics-Based Battery Model Parametrisation from Impedance Data}}, + author = {Hallemans, Noël and Courtier, Nicola E. and Please, Colin P. and Planden, Brady and Dhoot, Rishit and Timms, Robert and Chapman, S. Jon and Howey, David and Duncan, Stephen R.}, + journal = {Journal of the Electrochemical Society}, + volume = {172}, + number = {6}, + pages = {060507}, + year = {2025}, + publisher = {The Electrochemical Society}, + doi = {10.1149/1945-7111/add41b} + } + """ + ) + + ###################### + # Variables + ###################### + # Variables that depend on time only are created without a domain + Q = Variable("Discharge capacity [A.h]") + Qt = Variable("Throughput capacity [A.h]") + + # Variables that vary spatially are created with a domain + sto_n = Variable( + "Negative particle stoichiometry", + domain="negative particle", + auxiliary_domains={"secondary": "negative electrode"}, + ) + sto_p = Variable( + "Positive particle stoichiometry", + domain="positive particle", + auxiliary_domains={"secondary": "positive electrode"}, + ) + sto_e_n = Variable( + "Negative electrode electrolyte stoichiometry", + domain="negative electrode", + ) + sto_e_sep = Variable( + "Separator electrolyte stoichiometry", + domain="separator", + ) + sto_e_p = Variable( + "Positive electrode electrolyte stoichiometry", + domain="positive electrode", + ) + + # 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 + sto_n_surf = pybamm.surf(sto_n) + sto_p_surf = pybamm.surf(sto_p) + + # Events specify points at which a solution should terminate + self.events += [ + Event( + "Minimum negative particle surface stoichiometry", + pybamm.min(sto_n_surf) - 0.01, + ), + Event( + "Maximum negative particle surface stoichiometry", + (1 - 0.01) - pybamm.max(sto_n_surf), + ), + Event( + "Minimum positive particle surface stoichiometry", + pybamm.min(sto_p_surf) - 0.01, + ), + Event( + "Maximum positive particle surface stoichiometry", + (1 - 0.01) - pybamm.max(sto_p_surf), + ), + ] + + ###################### + # Parameters + ###################### + # Parameters are purely symbolic at this stage, and will be set by the + # `ParameterValues` class when the model is processed. + + F = self.param.F # Faraday constant + Rg = self.param.R # Universal gas constant + T = self.param.T_init # Temperature + RT_F = Rg * T / F # Thermal voltage + + soc_init = Parameter("Initial SoC") + x_0 = Parameter("Minimum negative stoichiometry") + x_100 = Parameter("Maximum negative stoichiometry") + y_100 = Parameter("Minimum positive stoichiometry") + y_0 = Parameter("Maximum positive stoichiometry") + + # Grouped parameters + Q_th_p = Parameter("Measured cell capacity [A.s]") / (y_0 - y_100) + Q_th_n = Parameter("Measured cell capacity [A.s]") / (x_100 - x_0) + Q_e = Parameter("Reference electrolyte capacity [A.s]") + + tau_d_p = Parameter("Positive particle diffusion time scale [s]") + tau_d_n = Parameter("Negative particle diffusion time scale [s]") + + tau_ct_p = Parameter("Positive electrode charge transfer time scale [s]") + tau_ct_n = Parameter("Negative electrode charge transfer time scale [s]") + + l_p = Parameter("Positive electrode relative thickness") + l_n = Parameter("Negative electrode relative thickness") + + t_plus = Parameter("Cation transference number") + + R0 = Parameter("Series resistance [Ohm]") + + zeta_n = Parameter("Negative electrode relative porosity") + zeta_p = Parameter("Positive electrode relative porosity") + + tau_e_n = Parameter("Negative electrode electrolyte diffusion time scale [s]") + tau_e_sep = Parameter("Separator electrolyte diffusion time scale [s]") + tau_e_p = Parameter("Positive electrode electrolyte diffusion time scale [s]") + + varsigma_e = Parameter("Reference electrolyte scaled conductivity [V-1.s-1]") + varsigma_e_n = varsigma_e * tau_e_sep / tau_e_n + varsigma_e_p = varsigma_e * tau_e_sep / tau_e_p + + ###################### + # Input current (positive on discharge) + ###################### + I = self.param.current_with_time + + ###################### + # State of Charge + ###################### + # The `rhs` dictionary contains differential equations, with the key being the + # variable in the d/dt + self.rhs[Q] = I / 3600 + self.rhs[Qt] = abs(I) / 3600 + # Initial conditions must be provided for the ODEs + self.initial_conditions[Q] = Scalar(0) + self.initial_conditions[Qt] = Scalar(0) + + ###################### + # Potentials + ###################### + U_n = self.U(sto_n_surf, "negative") + U_p = self.U(sto_p_surf, "positive") + + sto_n_init = x_0 + (x_100 - x_0) * soc_init + sto_p_init = y_0 + (y_100 - y_0) * soc_init + U_n_init = self.U(sto_n_init, "negative") + U_p_init = self.U(sto_p_init, "positive") + + ###################### + # Exchange current + ###################### + # Primary broadcasts are used to broadcast scalar quantities across a domain + # into a vector of the right shape, for multiplying with other vectors + alpha = 0.5 # cathodic transfer coefficient + j0_n = ( + sto_n_surf**alpha * (sto_e_n * (1 - sto_n_surf)) ** (1 - alpha) / tau_ct_n + ) + j0_p = ( + sto_p_surf**alpha * (sto_e_p * (1 - sto_p_surf)) ** (1 - alpha) / tau_ct_p + ) + + ###################### + # Double layer + ###################### + # Additional variables + v_s_n = Variable( + "Negative particle surface voltage [V]", + domain="negative electrode", + ) + v_s_p = Variable( + "Positive particle surface voltage [V]", + domain="positive electrode", + ) + + # Additional parameters + C_p = Parameter("Positive electrode capacitance [F]") + C_n = Parameter("Negative electrode capacitance [F]") + + # Overpotentials + eta_n = v_s_n - U_n + eta_p = v_s_p - U_p + + # Exchange current + j_n = j0_n * ( + pybamm.exp((1 - alpha) * eta_n / RT_F) - pybamm.exp(-alpha * eta_n / RT_F) + ) + j_p = j0_p * ( + pybamm.exp((1 - alpha) * eta_p / RT_F) - pybamm.exp(-alpha * eta_p / RT_F) + ) + + # Electrolyte current + i_e_n = varsigma_e_n * ( + pybamm.grad(v_s_n) + + (2 * RT_F * (1 - t_plus)) * pybamm.grad(sto_e_n) / sto_e_n + ) + i_e_p = varsigma_e_p * ( + pybamm.grad(v_s_p) + + (2 * RT_F * (1 - t_plus)) * pybamm.grad(sto_e_p) / sto_e_p + ) + + # Electrode surface potentials + self.rhs[v_s_n] = (l_n * Q_e * pybamm.div(i_e_n) - 3 * Q_th_n * j_n) / C_n + self.rhs[v_s_p] = (l_p * Q_e * pybamm.div(i_e_p) - 3 * Q_th_p * j_p) / C_p + + self.initial_conditions[v_s_n] = U_n_init + self.initial_conditions[v_s_p] = U_p_init + + self.boundary_conditions[v_s_n] = { + "left": (Scalar(0), "Neumann"), + "right": ( + I / (Q_e * varsigma_e_n) + - (2 * RT_F * (1 - t_plus)) + * pybamm.boundary_gradient(sto_e_n, "right") + / pybamm.boundary_value(sto_e_n, "right"), + "Neumann", + ), + } + self.boundary_conditions[v_s_p] = { + "left": ( + I / (Q_e * varsigma_e_p) + - (2 * RT_F * (1 - t_plus)) + * pybamm.boundary_gradient(sto_e_p, "left") + / pybamm.boundary_value(sto_e_p, "left"), + "Neumann", + ), + "right": (Scalar(0), "Neumann"), + } + + ###################### + # Particles + ###################### + # The div and grad operators will be converted to the appropriate matrix + # multiplication at the discretisation stage + self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / tau_d_n) + self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / tau_d_p) + + # Boundary conditions must be provided for equations with spatial derivatives + self.boundary_conditions[sto_n] = { + "left": (Scalar(0), "Neumann"), + "right": (-tau_d_n * j_n, "Neumann"), + } + self.boundary_conditions[sto_p] = { + "left": (Scalar(0), "Neumann"), + "right": (-tau_d_p * j_p, "Neumann"), + } + + self.initial_conditions[sto_n] = sto_n_init + self.initial_conditions[sto_p] = sto_p_init + + ###################### + # Electrolyte + ###################### + self.rhs[sto_e_n] = ( + pybamm.div(pybamm.grad(sto_e_n) / tau_e_n + (1 - t_plus) * i_e_n) + ) / zeta_n + self.rhs[sto_e_sep] = pybamm.div( + pybamm.grad(sto_e_sep) / tau_e_sep - t_plus * I / Q_e + ) + self.rhs[sto_e_p] = ( + pybamm.div(pybamm.grad(sto_e_p) / tau_e_p + (1 - t_plus) * i_e_p) + ) / zeta_p + + self.boundary_conditions[sto_e_n] = { + "left": (Scalar(0), "Neumann"), + "right": ( + tau_e_n * pybamm.boundary_gradient(sto_e_sep, "left") / tau_e_sep, + "Neumann", + ), + } + self.boundary_conditions[sto_e_sep] = { + "left": (pybamm.boundary_value(sto_e_n, "right"), "Dirichlet"), + "right": (pybamm.boundary_value(sto_e_p, "left"), "Dirichlet"), + } + self.boundary_conditions[sto_e_p] = { + "left": ( + tau_e_p * pybamm.boundary_gradient(sto_e_sep, "right") / tau_e_sep, + "Neumann", + ), + "right": (Scalar(0), "Neumann"), + } + + self.initial_conditions[sto_e_n] = PrimaryBroadcast( + Scalar(1), "negative electrode" + ) + self.initial_conditions[sto_e_sep] = PrimaryBroadcast(Scalar(1), "separator") + self.initial_conditions[sto_e_p] = PrimaryBroadcast( + Scalar(1), "positive electrode" + ) + + # Electrolyte overpotential + eta_e = ( + (2 * (1 - t_plus) * RT_F) + * ( + pybamm.log(pybamm.boundary_value(sto_e_sep, "right")) + - pybamm.log(pybamm.boundary_value(sto_e_sep, "left")) + ) + - (1 - l_p - l_n) * I / (varsigma_e * Q_e) + - ( + pybamm.boundary_value(v_s_p, "right") + - pybamm.boundary_value(v_s_p, "left") + ) + - ( + pybamm.boundary_value(v_s_n, "right") + - pybamm.boundary_value(v_s_n, "left") + ) + ) + + ###################### + # Cell voltage + ###################### + V = ( + pybamm.boundary_value(v_s_p, "right") + - pybamm.boundary_value(v_s_n, "left") + + eta_e + - R0 * I + ) + + # Save the initial OCV + self.param.ocv_init = U_p_init - U_n_init + + # Events specify points at which a solution should terminate + self.events += [ + Event("Minimum voltage [V]", V - self.param.voltage_low_cut), + Event("Maximum voltage [V]", self.param.voltage_high_cut - V), + ] + + ###################### + # (Some) variables + ###################### + # The `variables` dictionary contains all variables that might be useful for + # visualising the solution of the model + self.variables = { + "Negative particle stoichiometry": sto_n, + "Negative particle surface stoichiometry": sto_n_surf, + "Negative particle surface voltage [V]": v_s_n, + "Negative electrode potential [V]": eta_n + - pybamm.boundary_value(eta_n, "left"), + "Negative electrode electrolyte stoichiometry": sto_e_n, + "Separator electrolyte stoichiometry": sto_e_sep, + "Positive electrode electrolyte stoichiometry": sto_e_p, + "Electrolyte stoichiometry": pybamm.concatenation( + sto_e_n, sto_e_sep, sto_e_p + ), + "Positive particle stoichiometry": sto_p, + "Positive particle surface stoichiometry": sto_p_surf, + "Positive particle surface voltage [V]": v_s_p, + "Positive electrode potential [V]": V + + eta_p + - pybamm.boundary_value(eta_p, "right"), + "Electrolyte scaled current density [s-1]": pybamm.concatenation( + i_e_n, PrimaryBroadcast(I / Q_e, "separator"), i_e_p + ), + "Time [s]": pybamm_t, + "Current [A]": I, + "Current variable [A]": I, # for compatibility with pybamm.Experiment + "Discharge capacity [A.h]": Q, + "Throughput capacity [A.h]": Qt, + "Voltage [V]": V, + "Battery voltage [V]": V, + "Open-circuit voltage [V]": pybamm.boundary_value(U_p, "right") + - pybamm.boundary_value(U_n, "left"), + } + + def U(self, sto, domain): + """ + Dimensional open-circuit potential [V], calculated as U(x) = U_ref(x). + Credit: PyBaMM + """ + # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later + # will ensure that ocp goes to +- infinity if sto goes into that region + # anyway + Domain = domain.capitalize() + tol = pybamm.settings.tolerances["U__c_s"] + sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol) + inputs = {f"{Domain} particle surface stoichiometry": sto} + u_ref = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) + + # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 + # this will not affect the OCP for most values of sto + out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1)) + + if domain == "negative": + out.print_name = r"U_\mathrm{n}(c^\mathrm{surf}_\mathrm{s,n})" + elif domain == "positive": + out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" + return out + + def build_model(self): + """ + Build model variables and equations + Credit: PyBaMM + """ + self._build_model() + + self._built = True + pybamm.logger.info(f"Finish building {self.name}") + + @property + def default_parameter_values(self) -> pybamm.ParameterValues: + param = ParameterValues("Chen2020") + ce0 = param["Initial concentration in electrolyte [mol.m-3]"] + T = param["Ambient temperature [K]"] + param["Electrolyte conductivity [S.m-1]"] = param[ + "Electrolyte conductivity [S.m-1]" + ](ce0, T) + param["Electrolyte diffusivity [m2.s-1]"] = param[ + "Electrolyte diffusivity [m2.s-1]" + ](ce0, T) + return self.create_grouped_parameters(param) + + @property + def default_quick_plot_variables(self): + return [ + "Negative particle surface stoichiometry", + "Electrolyte stoichiometry", + "Positive particle surface stoichiometry", + "Current [A]", + { + "Negative electrode potential [V]", + "Negative particle surface voltage [V]", + }, + {"Electrolyte scaled current density [s-1]"}, + { + "Positive electrode potential [V]", + "Positive particle surface voltage [V]", + }, + {"Open-circuit voltage [V]", "Voltage [V]"}, + ] + + @property + def default_var_pts(self): + x_n = pybamm.SpatialVariable( + "x_n", + domain=["negative electrode"], + coord_sys="cartesian", + ) + x_s = pybamm.SpatialVariable( + "x_s", + domain=["separator"], + coord_sys="cartesian", + ) + x_p = pybamm.SpatialVariable( + "x_p", + domain=["positive electrode"], + coord_sys="cartesian", + ) + + # Add particle domains + r_n = pybamm.SpatialVariable( + "r_n", + domain=["negative particle"], + auxiliary_domains={"secondary": "negative electrode"}, + coord_sys="spherical polar", + ) + r_p = pybamm.SpatialVariable( + "r_p", + domain=["positive particle"], + auxiliary_domains={"secondary": "positive electrode"}, + coord_sys="spherical polar", + ) + + return {x_n: 20, x_s: 20, x_p: 20, r_n: 20, r_p: 20} + + @property + def default_geometry(self): + l_p = Parameter("Positive electrode relative thickness") + l_n = Parameter("Negative electrode relative thickness") + + return { + "negative electrode": {"x_n": {"min": 0, "max": l_n}}, + "separator": {"x_s": {"min": l_n, "max": 1 - l_p}}, + "positive electrode": {"x_p": {"min": 1 - l_p, "max": 1}}, + "negative particle": {"r_n": {"min": 0, "max": 1}}, + "positive particle": {"r_p": {"min": 0, "max": 1}}, + } + + @property + def default_submesh_types(self): + return { + "negative electrode": pybamm.Uniform1DSubMesh, + "separator": pybamm.Uniform1DSubMesh, + "positive electrode": pybamm.Uniform1DSubMesh, + "negative particle": pybamm.Uniform1DSubMesh, + "positive particle": pybamm.Uniform1DSubMesh, + } + + @property + def default_spatial_methods(self): + return { + "negative electrode": pybamm.FiniteVolume(), + "separator": pybamm.FiniteVolume(), + "positive electrode": pybamm.FiniteVolume(), + "negative particle": pybamm.FiniteVolume(), + "positive particle": pybamm.FiniteVolume(), + } + + @staticmethod + def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterValues: + """ + Create a parameter set for the Grouped Doyle Fuller Newman Model from a + PyBaMM lithium-ion ParameterValues object. + + Parameters + ---------- + parameter_values : pybamm.ParameterValues + Parameters and their corresponding values. + + Returns + ------- + parameter_values : pybamm.ParameterValues + A new set of parameters and their values. + """ + param = parameter_values + + # Unpack physical parameters + F = param["Faraday constant [C.mol-1]"] + alpha_p = param["Positive electrode active material volume fraction"] + alpha_n = param["Negative electrode active material volume fraction"] + c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] + c_max_n = param["Maximum concentration in negative electrode [mol.m-3]"] + L_p = param["Positive electrode thickness [m]"] + L_n = param["Negative electrode thickness [m]"] + epsilon_p = param["Positive electrode porosity"] + epsilon_n = param["Negative electrode porosity"] + R_p = param["Positive particle radius [m]"] + R_n = param["Negative particle radius [m]"] + D_p = param["Positive particle diffusivity [m2.s-1]"] + D_n = param["Negative particle diffusivity [m2.s-1]"] + b_p = param["Positive electrode Bruggeman coefficient (electrolyte)"] + b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] + Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] + Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] + m_p = 3.42e-6 # (A/m2)(m3/mol)**1.5 + m_n = 6.48e-7 # (A/m2)(m3/mol)**1.5 + sigma_p = ( + param["Positive electrode conductivity [S.m-1]"] + * alpha_p ** param["Positive electrode Bruggeman coefficient (electrode)"] + ) + sigma_n = ( + param["Negative electrode conductivity [S.m-1]"] + * alpha_n ** param["Negative electrode Bruggeman coefficient (electrode)"] + ) + + # Separator and electrolyte properties + ce0 = param["Initial concentration in electrolyte [mol.m-3]"] + De = param["Electrolyte diffusivity [m2.s-1]"] # (ce0, T) + L_s = param["Separator thickness [m]"] + epsilon_sep = param["Separator porosity"] + b_sep = param["Separator Bruggeman coefficient (electrolyte)"] + t_plus = param["Cation transference number"] + sigma_e = ( + param["Electrolyte conductivity [S.m-1]"] # (ce0, T) + * (epsilon_sep**b_sep) + ) + + # Compute the cell area and thickness + A = param["Electrode height [m]"] * param["Electrode width [m]"] + L = L_p + L_n + L_s + + # Compute the series resistance + Rs = (L_p / sigma_p + L_n / sigma_n) / (3 * A) + R0 = Rs + param["Contact resistance [Ohm]"] + + # Compute the stoichiometry limits and initial SOC + x_0, x_100, y_100, y_0 = get_min_max_stoichiometries(param) + sto_p_init = ( + param["Initial concentration in positive electrode [mol.m-3]"] / c_max_p + ) + soc_init = (sto_p_init - y_0) / (y_100 - y_0) + + # Compute the capacity within the stoichiometry limits + Q_th_p = F * alpha_p * c_max_p * L_p * A + Q_th_n = F * alpha_n * c_max_n * L_n * A + Q_meas_p = (y_0 - y_100) * Q_th_p + Q_meas_n = (x_100 - x_0) * Q_th_n + if abs(Q_meas_n / Q_meas_p - 1) > 1e-6: + raise ValueError( + "The measured capacity should be the same for both electrodes." + ) + + # Grouped parameters + Q_meas = (Q_meas_n + Q_meas_p) / 2 + Q_e = F * epsilon_sep * ce0 * L * A + varsigma_e = sigma_e / (F * epsilon_sep * ce0 * L**2) + + zeta_p = epsilon_p / epsilon_sep + zeta_n = epsilon_n / epsilon_sep + + tau_d_p = R_p**2 / D_p + tau_d_n = R_n**2 / D_n + + tau_e_p = epsilon_sep * L**2 / (epsilon_p**b_p * De) + tau_e_n = epsilon_sep * L**2 / (epsilon_n**b_n * De) + tau_e_sep = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) + + tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) + tau_ct_n = F * R_n / (m_n * np.sqrt(ce0)) + + C_p = 3 * alpha_p * Cdl_p * L_p * A / R_p + C_n = 3 * alpha_n * Cdl_n * L_n * A / R_n + + l_p = L_p / L + l_n = L_n / L + + parameter_dictionary = { + "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], + "Current function [A]": param["Current function [A]"], + "Initial temperature [K]": param["Ambient temperature [K]"], + "Initial SoC": soc_init, + "Minimum negative stoichiometry": x_0, + "Maximum negative stoichiometry": x_100, + "Minimum positive stoichiometry": y_100, + "Maximum positive stoichiometry": y_0, + "Lower voltage cut-off [V]": param["Lower voltage cut-off [V]"], + "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], + "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], + "Negative electrode OCP [V]": param["Negative electrode OCP [V]"], + "Measured cell capacity [A.s]": Q_meas, + "Reference electrolyte capacity [A.s]": Q_e, + "Reference electrolyte scaled conductivity [V-1.s-1]": varsigma_e, + "Positive electrode relative porosity": zeta_p, + "Negative electrode relative porosity": zeta_n, + "Positive particle diffusion time scale [s]": tau_d_p, + "Negative particle diffusion time scale [s]": tau_d_n, + "Positive electrode electrolyte diffusion time scale [s]": tau_e_p, + "Negative electrode electrolyte diffusion time scale [s]": tau_e_n, + "Separator electrolyte diffusion time scale [s]": tau_e_sep, + "Positive electrode charge transfer time scale [s]": tau_ct_p, + "Negative electrode charge transfer time scale [s]": tau_ct_n, + "Positive electrode capacitance [F]": C_p, + "Negative electrode capacitance [F]": C_n, + "Cation transference number": t_plus, + "Positive electrode relative thickness": l_p, + "Negative electrode relative thickness": l_n, + "Series resistance [Ohm]": R0, + } + return ParameterValues(values=parameter_dictionary) diff --git a/tests/integration/models/test_grouped_models.py b/tests/integration/models/test_grouped_models.py index ea57f6020..860d0b97e 100644 --- a/tests/integration/models/test_grouped_models.py +++ b/tests/integration/models/test_grouped_models.py @@ -32,6 +32,7 @@ class TestGroupedModels: pybop.lithium_ion.GroupedSPM(options={"surface form": "differential"}), pybop.lithium_ion.GroupedSPMe(), pybop.lithium_ion.GroupedSPMe(options={"surface form": "differential"}), + pybop.lithium_ion.GroupedDFN(), ], scope="module", ) diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index 3eeaa9fbf..2d3381cec 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -23,6 +23,7 @@ class TestModels: pybop.lithium_ion.GroupedSPM(options={"surface form": "differential"}), pybop.lithium_ion.GroupedSPMe(), pybop.lithium_ion.GroupedSPMe(options={"surface form": "differential"}), + pybop.lithium_ion.GroupedDFN(), ], scope="module", ) From 7d6e401f6132955c0b753f10e30378db9a5768f3 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Wed, 7 Jan 2026 13:57:11 +0000 Subject: [PATCH 05/37] Add voltage components --- pybop/models/lithium_ion/grouped_dfn.py | 136 +++++++++++++++--------- tests/unit/test_models.py | 4 +- 2 files changed, 90 insertions(+), 50 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index a90fb0705..0b796a1f2 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -49,7 +49,8 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): doi = {10.1149/1945-7111/add41b} } """ - ) + ) # Note that the electrode electrolyte timescales have been replaced by relative transport efficiencies + # and there is an additional variable (i_e) and associated grouped parameter (gamma_e) compared to the SPMe ###################### # Variables @@ -146,13 +147,11 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): zeta_n = Parameter("Negative electrode relative porosity") zeta_p = Parameter("Positive electrode relative porosity") - tau_e_n = Parameter("Negative electrode electrolyte diffusion time scale [s]") - tau_e_sep = Parameter("Separator electrolyte diffusion time scale [s]") - tau_e_p = Parameter("Positive electrode electrolyte diffusion time scale [s]") + beta_n = Parameter("Negative electrode relative transport efficiency") + beta_p = Parameter("Positive electrode relative transport efficiency") - varsigma_e = Parameter("Reference electrolyte scaled conductivity [V-1.s-1]") - varsigma_e_n = varsigma_e * tau_e_sep / tau_e_n - varsigma_e_p = varsigma_e * tau_e_sep / tau_e_p + tau_e = Parameter("Electrolyte diffusion time scale [s]") + gamma_e = Parameter("Reference electrolyte scaled conductivity [V-1.s-1]") ###################### # Input current (positive on discharge) @@ -224,11 +223,11 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): ) # Electrolyte current - i_e_n = varsigma_e_n * ( + i_e_n = (beta_n * gamma_e) * ( pybamm.grad(v_s_n) + (2 * RT_F * (1 - t_plus)) * pybamm.grad(sto_e_n) / sto_e_n ) - i_e_p = varsigma_e_p * ( + i_e_p = (beta_p * gamma_e) * ( pybamm.grad(v_s_p) + (2 * RT_F * (1 - t_plus)) * pybamm.grad(sto_e_p) / sto_e_p ) @@ -243,19 +242,19 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): self.boundary_conditions[v_s_n] = { "left": (Scalar(0), "Neumann"), "right": ( - I / (Q_e * varsigma_e_n) + I / (beta_n * gamma_e * Q_e) - (2 * RT_F * (1 - t_plus)) - * pybamm.boundary_gradient(sto_e_n, "right") - / pybamm.boundary_value(sto_e_n, "right"), + * (pybamm.boundary_gradient(sto_e_sep, "left") / beta_n) + / pybamm.boundary_value(sto_e_sep, "left"), "Neumann", ), } self.boundary_conditions[v_s_p] = { "left": ( - I / (Q_e * varsigma_e_p) + I / (beta_p * gamma_e * Q_e) - (2 * RT_F * (1 - t_plus)) - * pybamm.boundary_gradient(sto_e_p, "left") - / pybamm.boundary_value(sto_e_p, "left"), + * (pybamm.boundary_gradient(sto_e_sep, "right") / beta_p) + / pybamm.boundary_value(sto_e_sep, "right"), "Neumann", ), "right": (Scalar(0), "Neumann"), @@ -286,31 +285,25 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): # Electrolyte ###################### self.rhs[sto_e_n] = ( - pybamm.div(pybamm.grad(sto_e_n) / tau_e_n + (1 - t_plus) * i_e_n) + pybamm.div(pybamm.grad(sto_e_n) * beta_n / tau_e + (1 - t_plus) * i_e_n) ) / zeta_n self.rhs[sto_e_sep] = pybamm.div( - pybamm.grad(sto_e_sep) / tau_e_sep - t_plus * I / Q_e + pybamm.grad(sto_e_sep) / tau_e - t_plus * I / Q_e ) self.rhs[sto_e_p] = ( - pybamm.div(pybamm.grad(sto_e_p) / tau_e_p + (1 - t_plus) * i_e_p) + pybamm.div(pybamm.grad(sto_e_p) * beta_p / tau_e + (1 - t_plus) * i_e_p) ) / zeta_p self.boundary_conditions[sto_e_n] = { "left": (Scalar(0), "Neumann"), - "right": ( - tau_e_n * pybamm.boundary_gradient(sto_e_sep, "left") / tau_e_sep, - "Neumann", - ), + "right": (pybamm.boundary_gradient(sto_e_sep, "left") / beta_n, "Neumann"), } self.boundary_conditions[sto_e_sep] = { "left": (pybamm.boundary_value(sto_e_n, "right"), "Dirichlet"), "right": (pybamm.boundary_value(sto_e_p, "left"), "Dirichlet"), } self.boundary_conditions[sto_e_p] = { - "left": ( - tau_e_p * pybamm.boundary_gradient(sto_e_sep, "right") / tau_e_sep, - "Neumann", - ), + "left": (pybamm.boundary_gradient(sto_e_sep, "right") / beta_p, "Neumann"), "right": (Scalar(0), "Neumann"), } @@ -323,32 +316,32 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): ) # Electrolyte overpotential - eta_e = ( + eta_e = (2 * (1 - t_plus) * RT_F) * ( + pybamm.x_average(pybamm.log(sto_e_p)) + - pybamm.x_average(pybamm.log(sto_e_n)) + ) + + # Electrolyte Ohmic losses + DPhi_e = ( (2 * (1 - t_plus) * RT_F) * ( - pybamm.log(pybamm.boundary_value(sto_e_sep, "right")) - - pybamm.log(pybamm.boundary_value(sto_e_sep, "left")) + pybamm.log(pybamm.boundary_value(sto_e_p, "left")) + - pybamm.log(pybamm.boundary_value(sto_e_n, "right")) ) - - (1 - l_p - l_n) * I / (varsigma_e * Q_e) + - eta_e + - (1 - l_p - l_n) * I / (gamma_e * Q_e) - ( - pybamm.boundary_value(v_s_p, "right") + pybamm.x_average(v_s_p) + - pybamm.x_average(v_s_n) - pybamm.boundary_value(v_s_p, "left") - ) - - ( - pybamm.boundary_value(v_s_n, "right") - - pybamm.boundary_value(v_s_n, "left") + + pybamm.boundary_value(v_s_n, "right") ) ) ###################### # Cell voltage ###################### - V = ( - pybamm.boundary_value(v_s_p, "right") - - pybamm.boundary_value(v_s_n, "left") - + eta_e - - R0 * I - ) + V = pybamm.x_average(v_s_p) - pybamm.x_average(v_s_n) + eta_e + DPhi_e - R0 * I # Save the initial OCV self.param.ocv_init = U_p_init - U_n_init @@ -359,6 +352,49 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): Event("Maximum voltage [V]", self.param.voltage_high_cut - V), ] + ###################### + # Voltage components + ###################### + # Include the following variables to enable plotting via PyBaMM's plot_voltage_components + ocp_n_bulk = self.U( + pybamm.x_average(zeta_n * pybamm.r_average(sto_n)) + / pybamm.x_average(zeta_n), + "negative", + ) + ocp_p_bulk = self.U( + pybamm.x_average(zeta_p * pybamm.r_average(sto_p)) + / pybamm.x_average(zeta_p), + "positive", + ) + voltage_components = { + "Battery open-circuit voltage [V]": ocp_p_bulk - ocp_n_bulk, + "Battery particle concentration overpotential [V]": ( + (pybamm.x_average(self.U(sto_p_surf, "positive")) - ocp_p_bulk) + - (pybamm.x_average(self.U(sto_n_surf, "negative")) - ocp_n_bulk) + ), + "X-averaged battery reaction overpotential [V]": pybamm.x_average(eta_p) + - pybamm.x_average(eta_n), + "X-averaged battery concentration overpotential [V]": eta_e, + "X-averaged battery electrolyte ohmic losses [V]": DPhi_e, + "X-averaged battery solid phase ohmic losses [V]": Scalar(0), + "Contact overpotential [V]": R0 * I, # includes solid phase Ohmic losses + # and split by electrode + "Negative electrode bulk open-circuit potential [V]": ocp_n_bulk, + "Positive electrode bulk open-circuit potential [V]": ocp_p_bulk, + "Negative particle concentration overpotential [V]": ( + pybamm.x_average(self.U(sto_n_surf, "negative")) - ocp_n_bulk + ), + "Positive particle concentration overpotential [V]": ( + pybamm.x_average(self.U(sto_p_surf, "positive")) - ocp_p_bulk + ), + "X-averaged negative electrode reaction overpotential [V]" + "": pybamm.x_average(eta_n), + "X-averaged positive electrode reaction overpotential [V]" + "": pybamm.x_average(eta_p), + "X-averaged battery negative solid phase ohmic losses [V]": Scalar(0), + "X-averaged battery positive solid phase ohmic losses [V]": Scalar(0), + } + ###################### # (Some) variables ###################### @@ -386,6 +422,7 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): i_e_n, PrimaryBroadcast(I / Q_e, "separator"), i_e_p ), "Time [s]": pybamm_t, + "Time [h]": pybamm_t / 3600, "Current [A]": I, "Current variable [A]": I, # for compatibility with pybamm.Experiment "Discharge capacity [A.h]": Q, @@ -394,6 +431,7 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): "Battery voltage [V]": V, "Open-circuit voltage [V]": pybamm.boundary_value(U_p, "right") - pybamm.boundary_value(U_n, "left"), + **voltage_components, } def U(self, sto, domain): @@ -616,7 +654,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Grouped parameters Q_meas = (Q_meas_n + Q_meas_p) / 2 Q_e = F * epsilon_sep * ce0 * L * A - varsigma_e = sigma_e / (F * epsilon_sep * ce0 * L**2) + gamma_e = sigma_e / (F * epsilon_sep * ce0 * L**2) zeta_p = epsilon_p / epsilon_sep zeta_n = epsilon_n / epsilon_sep @@ -624,9 +662,9 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal tau_d_p = R_p**2 / D_p tau_d_n = R_n**2 / D_n - tau_e_p = epsilon_sep * L**2 / (epsilon_p**b_p * De) - tau_e_n = epsilon_sep * L**2 / (epsilon_n**b_n * De) - tau_e_sep = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) + tau_e = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) + beta_p = epsilon_p**b_p / epsilon_sep**b_sep + beta_n = epsilon_n**b_n / epsilon_sep**b_sep tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) tau_ct_n = F * R_n / (m_n * np.sqrt(ce0)) @@ -652,14 +690,14 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Negative electrode OCP [V]": param["Negative electrode OCP [V]"], "Measured cell capacity [A.s]": Q_meas, "Reference electrolyte capacity [A.s]": Q_e, - "Reference electrolyte scaled conductivity [V-1.s-1]": varsigma_e, + "Reference electrolyte scaled conductivity [V-1.s-1]": gamma_e, + "Positive electrode relative transport efficiency": beta_p, + "Negative electrode relative transport efficiency": beta_n, "Positive electrode relative porosity": zeta_p, "Negative electrode relative porosity": zeta_n, "Positive particle diffusion time scale [s]": tau_d_p, "Negative particle diffusion time scale [s]": tau_d_n, - "Positive electrode electrolyte diffusion time scale [s]": tau_e_p, - "Negative electrode electrolyte diffusion time scale [s]": tau_e_n, - "Separator electrolyte diffusion time scale [s]": tau_e_sep, + "Electrolyte diffusion time scale [s]": tau_e, "Positive electrode charge transfer time scale [s]": tau_ct_p, "Negative electrode charge transfer time scale [s]": tau_ct_n, "Positive electrode capacitance [F]": C_p, diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index 2d3381cec..6c3dfb886 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -49,7 +49,9 @@ def test_model_simulation(self, model): fig = solution.plot() assert isinstance(fig, pybamm.QuickPlot) - if isinstance(model, pybop.lithium_ion.GroupedSPMe): + if isinstance( + model, pybop.lithium_ion.GroupedSPMe | pybop.lithium_ion.GroupedDFN + ): for split in [False, True]: fig, ax = solution.plot_voltage_components(split_by_electrode=split) assert isinstance(fig, Figure) From f92085da9e9e2996e6dc62c0e3f9aa5bd311edb9 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 8 Jan 2026 00:04:00 +0000 Subject: [PATCH 06/37] Add set_initial_state --- pybop/models/lithium_ion/grouped_dfn.py | 54 ++++++++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 0b796a1f2..5a0c62f77 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -707,4 +707,56 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Negative electrode relative thickness": l_n, "Series resistance [Ohm]": R0, } - return ParameterValues(values=parameter_dictionary) + parameter_values = ParameterValues(values=parameter_dictionary) + parameter_values._set_initial_state = set_initial_state # noqa: SLF001 + return parameter_values + + +def set_initial_state( + initial_value, + parameter_values, + direction=None, + param=None, + inplace=True, + options=None, + inputs=None, + tol=1e-6, +): + """ + Set the value of the initial state of charge. + + Parameters + ---------- + initial_value : float + Target initial value. + If float, interpreted as SOC, must be between 0 and 1. + If string e.g. "4 V", interpreted as voltage, must be between V_min and V_max. + parameter_values : :class:`pybamm.ParameterValues` + Parameters and their corresponding values. + param : :class:`pybamm.LithiumIonParameters`, optional + The symbolic parameter set to use for the simulation. + If not provided, the default parameter set will be used. + inplace: bool, optional + If True, replace the parameters values in place. Otherwise, return a new set of + parameter values. Default is True. + options : dict-like, optional + A dictionary of options to be passed to the model, see + :class:`pybamm.BatteryModelOptions`. + inputs : dict, optional + A dictionary of input parameters to pass to the model when solving. + tol : float, optional + The tolerance for the solver used to compute the initial stoichiometries. + A lower value results in higher precision but may increase computation time. + Default is 1e-6. + """ + parameter_values = parameter_values if inplace else parameter_values.copy() + + if isinstance(initial_value, int | float): + if not 0 <= initial_value <= 1: + raise ValueError("Initial SOC should be between 0 and 1") + parameter_values["Initial SoC"] = initial_value + + else: + raise ValueError("Initial value must be a float between 0 and 1.") + + return parameter_values From b19d6fd5ae07113529dc58b24a98502951cb9893 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 8 Jan 2026 17:12:35 +0000 Subject: [PATCH 07/37] Add SPM voltage components --- pybop/models/lithium_ion/grouped_spm.py | 43 +++++++++++++++++++++++-- 1 file changed, 41 insertions(+), 2 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index f1771799b..22c4a4bf9 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -237,6 +237,43 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): Event("Maximum voltage [V]", self.param.voltage_high_cut - V), ] + ###################### + # Voltage components + ###################### + # Include the following variables to enable plotting via PyBaMM's plot_voltage_components + ocp_n_bulk = self.U(pybamm.x_average(pybamm.r_average(sto_n)), "negative") + ocp_p_bulk = self.U(pybamm.x_average(pybamm.r_average(sto_p)), "positive") + voltage_components = { + "Battery open-circuit voltage [V]": ocp_p_bulk - ocp_n_bulk, + "Battery particle concentration overpotential [V]": ( + (pybamm.x_average(self.U(sto_p_surf, "positive")) - ocp_p_bulk) + - (pybamm.x_average(self.U(sto_n_surf, "negative")) - ocp_n_bulk) + ), + "X-averaged battery reaction overpotential [V]": pybamm.x_average(eta_p) + - pybamm.x_average(eta_n), + "X-averaged battery concentration overpotential [V]": Scalar(0), + "X-averaged battery electrolyte ohmic losses [V]": Scalar(0), + "X-averaged battery solid phase ohmic losses [V]": Scalar(0), + "Contact overpotential [V]": R0 * I, # includes Ohmic losses in this model + # and split by electrode + "Negative electrode bulk open-circuit potential [V]": ocp_n_bulk, + "Positive electrode bulk open-circuit potential [V]": ocp_p_bulk, + "Negative particle concentration overpotential [V]": pybamm.x_average( + self.U(sto_n_surf, "negative") + ) + - ocp_n_bulk, + "Positive particle concentration overpotential [V]": pybamm.x_average( + self.U(sto_p_surf, "positive") + ) + - ocp_p_bulk, + "X-averaged negative electrode reaction overpotential [V]" + "": pybamm.x_average(eta_n), + "X-averaged positive electrode reaction overpotential [V]" + "": pybamm.x_average(eta_p), + "X-averaged battery negative solid phase ohmic losses [V]": Scalar(0), + "X-averaged battery positive solid phase ohmic losses [V]": Scalar(0), + } + ###################### # (Some) variables ###################### @@ -265,6 +302,7 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): + eta_p - pybamm.boundary_value(eta_p, "right"), "Time [s]": pybamm_t, + "Time [h]": pybamm_t / 3600, "Current [A]": I, "Current variable [A]": I, # for compatibility with pybamm.Experiment "Discharge capacity [A.h]": Q, @@ -272,6 +310,7 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): "Voltage [V]": V, "Battery voltage [V]": V, "Open-circuit voltage [V]": U_p - U_n, + **voltage_components, } def U(self, sto, domain): @@ -454,7 +493,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal L_s = param["Separator thickness [m]"] epsilon_sep = param["Separator porosity"] b_sep = param["Separator Bruggeman coefficient (electrolyte)"] - sigma_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) + kappa_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) # Compute the cell area and thickness A = param["Electrode height [m]"] * param["Electrode width [m]"] @@ -465,7 +504,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal L_p / (3 * epsilon_p**b_p) + L_s / (epsilon_sep**b_sep) + L_n / (3 * epsilon_n**b_n) - ) / (sigma_e * A) + ) / (kappa_e * A) Rs = (L_p / sigma_p + L_n / sigma_n) / (3 * A) R0 = Re + Rs + param["Contact resistance [Ohm]"] From 5df89f9384d7279da307b2c99ef8c9b1ce33356d Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Wed, 18 Feb 2026 15:14:35 +0000 Subject: [PATCH 08/37] Finish merging --- pybop/models/lithium_ion/grouped_dfn.py | 69 ++----------------------- 1 file changed, 5 insertions(+), 64 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 5a0c62f77..a4bd24086 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -9,14 +9,15 @@ Scalar, Variable, ) -from pybamm import lithium_ion as pybamm_lithium_ion from pybamm import t as pybamm_t from pybamm.models.full_battery_models.lithium_ion.electrode_soh import ( get_min_max_stoichiometries, ) +from pybop.models.lithium_ion.base_model import BaseGroupedModel -class GroupedDFN(pybamm_lithium_ion.BaseModel): + +class GroupedDFN(BaseGroupedModel): """ A grouped parameter version of the Doyle Fuller Newman (DFN) model. @@ -458,16 +459,6 @@ def U(self, sto, domain): out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" return out - def build_model(self): - """ - Build model variables and equations - Credit: PyBaMM - """ - self._build_model() - - self._built = True - pybamm.logger.info(f"Finish building {self.name}") - @property def default_parameter_values(self) -> pybamm.ParameterValues: param = ParameterValues("Chen2020") @@ -586,7 +577,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal param = parameter_values # Unpack physical parameters - F = param["Faraday constant [C.mol-1]"] + F = pybamm.constants.F.value alpha_p = param["Positive electrode active material volume fraction"] alpha_n = param["Negative electrode active material volume fraction"] c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] @@ -708,55 +699,5 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Series resistance [Ohm]": R0, } parameter_values = ParameterValues(values=parameter_dictionary) - parameter_values._set_initial_state = set_initial_state # noqa: SLF001 + parameter_values._set_initial_state = GroupedDFN.set_initial_state # noqa: SLF001 return parameter_values - - -def set_initial_state( - initial_value, - parameter_values, - direction=None, - param=None, - inplace=True, - options=None, - inputs=None, - tol=1e-6, -): - """ - Set the value of the initial state of charge. - - Parameters - ---------- - initial_value : float - Target initial value. - If float, interpreted as SOC, must be between 0 and 1. - If string e.g. "4 V", interpreted as voltage, must be between V_min and V_max. - parameter_values : :class:`pybamm.ParameterValues` - Parameters and their corresponding values. - param : :class:`pybamm.LithiumIonParameters`, optional - The symbolic parameter set to use for the simulation. - If not provided, the default parameter set will be used. - inplace: bool, optional - If True, replace the parameters values in place. Otherwise, return a new set of - parameter values. Default is True. - options : dict-like, optional - A dictionary of options to be passed to the model, see - :class:`pybamm.BatteryModelOptions`. - inputs : dict, optional - A dictionary of input parameters to pass to the model when solving. - tol : float, optional - The tolerance for the solver used to compute the initial stoichiometries. - A lower value results in higher precision but may increase computation time. - Default is 1e-6. - """ - parameter_values = parameter_values if inplace else parameter_values.copy() - - if isinstance(initial_value, int | float): - if not 0 <= initial_value <= 1: - raise ValueError("Initial SOC should be between 0 and 1") - parameter_values["Initial SoC"] = initial_value - - else: - raise ValueError("Initial value must be a float between 0 and 1.") - - return parameter_values From 2030fa320225578013a9a1a26993a6bc1e81e8a4 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Tue, 17 Mar 2026 23:16:32 +0000 Subject: [PATCH 09/37] Enable sto-dependent diffusion timescales --- pybop/models/lithium_ion/grouped_dfn.py | 19 +++++++++++------- pybop/models/lithium_ion/grouped_spm.py | 25 +++++++++++++++++------- pybop/models/lithium_ion/grouped_spme.py | 25 +++++++++++++++++------- 3 files changed, 48 insertions(+), 21 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index a4bd24086..bd2006f2a 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -132,9 +132,6 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): Q_th_n = Parameter("Measured cell capacity [A.s]") / (x_100 - x_0) Q_e = Parameter("Reference electrolyte capacity [A.s]") - tau_d_p = Parameter("Positive particle diffusion time scale [s]") - tau_d_n = Parameter("Negative particle diffusion time scale [s]") - tau_ct_p = Parameter("Positive electrode charge transfer time scale [s]") tau_ct_n = Parameter("Negative electrode charge transfer time scale [s]") @@ -266,17 +263,17 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): ###################### # The div and grad operators will be converted to the appropriate matrix # multiplication at the discretisation stage - self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / tau_d_n) - self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / tau_d_p) + self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / self.tau_d(sto_n, "negative")) + self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / self.tau_d(sto_p, "positive")) # Boundary conditions must be provided for equations with spatial derivatives self.boundary_conditions[sto_n] = { "left": (Scalar(0), "Neumann"), - "right": (-tau_d_n * j_n, "Neumann"), + "right": (-self.tau_d(sto_n_surf, "negative") * j_n, "Neumann"), } self.boundary_conditions[sto_p] = { "left": (Scalar(0), "Neumann"), - "right": (-tau_d_p * j_p, "Neumann"), + "right": (-self.tau_d(sto_p_surf, "positive") * j_p, "Neumann"), } self.initial_conditions[sto_n] = sto_n_init @@ -459,6 +456,14 @@ def U(self, sto, domain): out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" return out + def tau_d(self, sto, domain): + """ + Dimensional diffusion time scale [s]. + """ + Domain = domain.capitalize() + inputs = {f"{Domain} particle surface stoichiometry": sto} + return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) + @property def default_parameter_values(self) -> pybamm.ParameterValues: param = ParameterValues("Chen2020") diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 9108dc1db..51a00c094 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -119,9 +119,6 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): Q_th_p = Parameter("Measured cell capacity [A.s]") / (y_0 - y_100) Q_th_n = Parameter("Measured cell capacity [A.s]") / (x_100 - x_0) - tau_d_p = Parameter("Positive particle diffusion time scale [s]") - tau_d_n = Parameter("Negative particle diffusion time scale [s]") - tau_ct_p = Parameter("Positive electrode charge transfer time scale [s]") tau_ct_n = Parameter("Negative electrode charge transfer time scale [s]") @@ -208,17 +205,23 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): ###################### # The div and grad operators will be converted to the appropriate matrix # multiplication at the discretisation stage - self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / tau_d_n) - self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / tau_d_p) + self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / self.tau_d(sto_n, "negative")) + self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / self.tau_d(sto_p, "positive")) # Boundary conditions must be provided for equations with spatial derivatives self.boundary_conditions[sto_n] = { "left": (Scalar(0), "Neumann"), - "right": (-tau_d_n * pybamm.x_average(j_n), "Neumann"), + "right": ( + -self.tau_d(sto_n_surf, "negative") * pybamm.x_average(j_n), + "Neumann", + ), } self.boundary_conditions[sto_p] = { "left": (Scalar(0), "Neumann"), - "right": (-tau_d_p * pybamm.x_average(j_p), "Neumann"), + "right": ( + -self.tau_d(sto_p_surf, "positive") * pybamm.x_average(j_p), + "Neumann", + ), } self.initial_conditions[sto_n] = sto_n_init @@ -338,6 +341,14 @@ def U(self, sto, domain): out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" return out + def tau_d(self, sto, domain): + """ + Dimensional diffusion time scale [s]. + """ + Domain = domain.capitalize() + inputs = {f"{Domain} particle surface stoichiometry": sto} + return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) + @property def default_parameter_values(self) -> ParameterValues: param = ParameterValues("Chen2020") diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index aace961b0..35c44b755 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -153,9 +153,6 @@ def __init__( Q_th_n = Parameter("Measured cell capacity [A.s]") / (x_100 - x_0) Q_e = Parameter("Reference electrolyte capacity [A.s]") - tau_d_p = Parameter("Positive particle diffusion time scale [s]") - tau_d_n = Parameter("Negative particle diffusion time scale [s]") - tau_ct_p = Parameter("Positive electrode charge transfer time scale [s]") tau_ct_n = Parameter("Negative electrode charge transfer time scale [s]") @@ -267,17 +264,23 @@ def __init__( ###################### # The div and grad operators will be converted to the appropriate matrix # multiplication at the discretisation stage - self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / tau_d_n) - self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / tau_d_p) + self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / self.tau_d(sto_n, "negative")) + self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / self.tau_d(sto_p, "positive")) # Boundary conditions must be provided for equations with spatial derivatives self.boundary_conditions[sto_n] = { "left": (Scalar(0), "Neumann"), - "right": (-tau_d_n * pybamm.x_average(j_n), "Neumann"), + "right": ( + -self.tau_d(sto_n_surf, "negative") * pybamm.x_average(j_n), + "Neumann", + ), } self.boundary_conditions[sto_p] = { "left": (Scalar(0), "Neumann"), - "right": (-tau_d_p * pybamm.x_average(j_p), "Neumann"), + "right": ( + -self.tau_d(sto_p_surf, "positive") * pybamm.x_average(j_p), + "Neumann", + ), } self.initial_conditions[sto_n] = sto_n_init @@ -458,6 +461,14 @@ def U(self, sto, domain): out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" return out + def tau_d(self, sto, domain): + """ + Dimensional diffusion time scale [s]. + """ + Domain = domain.capitalize() + inputs = {f"{Domain} particle surface stoichiometry": sto} + return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) + @property def default_parameter_values(self) -> ParameterValues: param = ParameterValues("Chen2020") From ee636b569b5bc2e0344d08167102fa75143f6c8e Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Mon, 30 Mar 2026 13:14:19 +0100 Subject: [PATCH 10/37] Pass ambient temperature --- pybop/models/lithium_ion/grouped_dfn.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index bd2006f2a..15e43647c 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -674,6 +674,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_dictionary = { "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], "Current function [A]": param["Current function [A]"], + "Ambient temperature [K]": param["Ambient temperature [K]"], "Initial temperature [K]": param["Ambient temperature [K]"], "Initial SoC": soc_init, "Minimum negative stoichiometry": x_0, From 7ae8baff271dbabda77c3ddd42f779487e77e17c Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 15 May 2026 11:31:20 +0100 Subject: [PATCH 11/37] Update event tolerance --- pybop/models/lithium_ion/grouped_dfn.py | 9 +++++---- pybop/models/lithium_ion/grouped_spm.py | 9 +++++---- pybop/models/lithium_ion/grouped_spme.py | 9 +++++---- pybop/models/lithium_ion/sp_diffusion.py | 5 +++-- 4 files changed, 18 insertions(+), 14 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 15e43647c..cb5f82125 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -91,22 +91,23 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): sto_p_surf = pybamm.surf(sto_p) # Events specify points at which a solution should terminate + tol = pybamm.settings.tolerances["U__c_s"] self.events += [ Event( "Minimum negative particle surface stoichiometry", - pybamm.min(sto_n_surf) - 0.01, + pybamm.min(sto_n_surf) - tol, ), Event( "Maximum negative particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_n_surf), + (1 - tol) - pybamm.max(sto_n_surf), ), Event( "Minimum positive particle surface stoichiometry", - pybamm.min(sto_p_surf) - 0.01, + pybamm.min(sto_p_surf) - tol, ), Event( "Maximum positive particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_p_surf), + (1 - tol) - pybamm.max(sto_p_surf), ), ] diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 37c3c5d63..0908ce3b4 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -79,22 +79,23 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): sto_p_surf = pybamm.surf(sto_p) # Events specify points at which a solution should terminate + tol = pybamm.settings.tolerances["U__c_s"] self.events += [ Event( "Minimum negative particle surface stoichiometry", - pybamm.min(sto_n_surf) - 0.01, + pybamm.min(sto_n_surf) - tol, ), Event( "Maximum negative particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_n_surf), + (1 - tol) - pybamm.max(sto_n_surf), ), Event( "Minimum positive particle surface stoichiometry", - pybamm.min(sto_p_surf) - 0.01, + pybamm.min(sto_p_surf) - tol, ), Event( "Maximum positive particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_p_surf), + (1 - tol) - pybamm.max(sto_p_surf), ), ] diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 00396bbaf..8505abdde 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -99,22 +99,23 @@ def __init__( sto_p_surf = pybamm.surf(sto_p) # Events specify points at which a solution should terminate + tol = pybamm.settings.tolerances["U__c_s"] self.events += [ Event( "Minimum negative particle surface stoichiometry", - pybamm.min(sto_n_surf) - 0.01, + pybamm.min(sto_n_surf) - tol, ), Event( "Maximum negative particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_n_surf), + (1 - tol) - pybamm.max(sto_n_surf), ), Event( "Minimum positive particle surface stoichiometry", - pybamm.min(sto_p_surf) - 0.01, + pybamm.min(sto_p_surf) - tol, ), Event( "Maximum positive particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_p_surf), + (1 - tol) - pybamm.max(sto_p_surf), ), # model does not capture electrolyte depletion, use the DFN instead Event( diff --git a/pybop/models/lithium_ion/sp_diffusion.py b/pybop/models/lithium_ion/sp_diffusion.py index 6798aec17..4af588fbc 100644 --- a/pybop/models/lithium_ion/sp_diffusion.py +++ b/pybop/models/lithium_ion/sp_diffusion.py @@ -48,14 +48,15 @@ def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): sto_surf = pybamm.surf(sto) # Events specify points at which a solution should terminate + tol = pybamm.settings.tolerances["U__c_s"] self.events += [ Event( "Minimum particle surface stoichiometry", - pybamm.min(sto_surf) - 0.01, + pybamm.min(sto_surf) - tol, ), Event( "Maximum particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_surf), + (1 - tol) - pybamm.max(sto_surf), ), ] From b8d1a7a9797c93831958ca53c603dd2fd2291ac1 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 15 May 2026 11:31:35 +0100 Subject: [PATCH 12/37] Simplify open-circuit voltage --- pybop/models/lithium_ion/cell_temperature.py | 13 ++----------- pybop/models/lithium_ion/grouped_dfn.py | 13 ++----------- pybop/models/lithium_ion/grouped_spm.py | 13 ++----------- pybop/models/lithium_ion/grouped_spme.py | 13 ++----------- pybop/models/lithium_ion/sp_diffusion.py | 13 ++----------- 5 files changed, 10 insertions(+), 55 deletions(-) diff --git a/pybop/models/lithium_ion/cell_temperature.py b/pybop/models/lithium_ion/cell_temperature.py index 3fbb32f31..e2317bb48 100644 --- a/pybop/models/lithium_ion/cell_temperature.py +++ b/pybop/models/lithium_ion/cell_temperature.py @@ -145,21 +145,12 @@ def __init__(self, name="Cell Temperature Model", **model_kwargs): def U(self, sto, domain): """ - Dimensional open-circuit potential [V], calculated as U(x) = U_ref(x). + Dimensional open-circuit potential [V]. Credit: PyBaMM """ - # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later - # will ensure that ocp goes to +- infinity if sto goes into that region - # anyway Domain = domain.capitalize() - tol = pybamm.settings.tolerances["U__c_s"] - sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol) inputs = {f"{Domain} particle surface stoichiometry": sto} - u_ref = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) - - # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 - # this will not affect the OCP for most values of sto - out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1)) + out = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) if domain == "negative": out.print_name = r"U_\mathrm{n}(c^\mathrm{surf}_\mathrm{s,n})" diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index cb5f82125..505b587c0 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -435,21 +435,12 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): def U(self, sto, domain): """ - Dimensional open-circuit potential [V], calculated as U(x) = U_ref(x). + Dimensional open-circuit potential [V]. Credit: PyBaMM """ - # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later - # will ensure that ocp goes to +- infinity if sto goes into that region - # anyway Domain = domain.capitalize() - tol = pybamm.settings.tolerances["U__c_s"] - sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol) inputs = {f"{Domain} particle surface stoichiometry": sto} - u_ref = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) - - # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 - # this will not affect the OCP for most values of sto - out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1)) + out = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) if domain == "negative": out.print_name = r"U_\mathrm{n}(c^\mathrm{surf}_\mathrm{s,n})" diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 0908ce3b4..581defe71 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -321,21 +321,12 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): def U(self, sto, domain): """ - Dimensional open-circuit potential [V], calculated as U(x) = U_ref(x). + Dimensional open-circuit potential [V]. Credit: PyBaMM """ - # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later - # will ensure that ocp goes to +- infinity if sto goes into that region - # anyway Domain = domain.capitalize() - tol = pybamm.settings.tolerances["U__c_s"] - sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol) inputs = {f"{Domain} particle surface stoichiometry": sto} - u_ref = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) - - # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 - # this will not affect the OCP for most values of sto - out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1)) + out = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) if domain == "negative": out.print_name = r"U_\mathrm{n}(c^\mathrm{surf}_\mathrm{s,n})" diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 8505abdde..e19d6cf71 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -441,21 +441,12 @@ def __init__( def U(self, sto, domain): """ - Dimensional open-circuit potential [V], calculated as U(x) = U_ref(x). + Dimensional open-circuit potential [V]. Credit: PyBaMM """ - # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later - # will ensure that ocp goes to +- infinity if sto goes into that region - # anyway Domain = domain.capitalize() - tol = pybamm.settings.tolerances["U__c_s"] - sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol) inputs = {f"{Domain} particle surface stoichiometry": sto} - u_ref = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) - - # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 - # this will not affect the OCP for most values of sto - out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1)) + out = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) if domain == "negative": out.print_name = r"U_\mathrm{n}(c^\mathrm{surf}_\mathrm{s,n})" diff --git a/pybop/models/lithium_ion/sp_diffusion.py b/pybop/models/lithium_ion/sp_diffusion.py index 4af588fbc..1ba955034 100644 --- a/pybop/models/lithium_ion/sp_diffusion.py +++ b/pybop/models/lithium_ion/sp_diffusion.py @@ -129,20 +129,11 @@ def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): def U(self, sto): """ - Dimensional open-circuit potential [V], calculated as U(x) = U_ref(x). + Dimensional open-circuit potential [V]. Credit: PyBaMM """ - # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later - # will ensure that ocp goes to +- infinity if sto goes into that region - # anyway - tol = pybamm.settings.tolerances["U__c_s"] - sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol) inputs = {"Particle surface stoichiometry": sto} - u_ref = FunctionParameter("Electrode OCP [V]", inputs) - - # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 - # this will not affect the OCP for most values of sto - out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1)) + out = FunctionParameter("Electrode OCP [V]", inputs) out.print_name = r"U(c^\mathrm{surf}_\mathrm{s})" return out From faa5fa05528e44de04788f05df11ea9a16db5814 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 28 May 2026 13:26:00 +0100 Subject: [PATCH 13/37] Enable alpha not 0.5 --- pybop/models/lithium_ion/grouped_dfn.py | 49 ++++++++------- pybop/models/lithium_ion/grouped_spm.py | 52 +++++++--------- pybop/models/lithium_ion/grouped_spme.py | 60 +++++++++---------- .../integration/models/test_grouped_models.py | 1 + tests/unit/test_models.py | 1 + 5 files changed, 81 insertions(+), 82 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 505b587c0..70b2c2384 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -36,6 +36,9 @@ class GroupedDFN(BaseGroupedModel): def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): super().__init__(name=name, **model_kwargs) + # Unpack model options + include_double_layer = self.options["surface form"] == "differential" + pybamm.citations.register( """ @article{Hallemans2025, @@ -61,6 +64,15 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): Qt = Variable("Throughput capacity [A.h]") # Variables that vary spatially are created with a domain + v_s_n = Variable( + "Negative particle surface voltage [V]", + domain="negative electrode", + ) + v_s_p = Variable( + "Positive particle surface voltage [V]", + domain="positive electrode", + ) + sto_n = Variable( "Negative particle stoichiometry", domain="negative particle", @@ -185,6 +197,8 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): # Primary broadcasts are used to broadcast scalar quantities across a domain # into a vector of the right shape, for multiplying with other vectors alpha = 0.5 # cathodic transfer coefficient + + # Reference exchange current j0_n = ( sto_n_surf**alpha * (sto_e_n * (1 - sto_n_surf)) ** (1 - alpha) / tau_ct_n ) @@ -192,23 +206,6 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): sto_p_surf**alpha * (sto_e_p * (1 - sto_p_surf)) ** (1 - alpha) / tau_ct_p ) - ###################### - # Double layer - ###################### - # Additional variables - v_s_n = Variable( - "Negative particle surface voltage [V]", - domain="negative electrode", - ) - v_s_p = Variable( - "Positive particle surface voltage [V]", - domain="positive electrode", - ) - - # Additional parameters - C_p = Parameter("Positive electrode capacitance [F]") - C_n = Parameter("Negative electrode capacitance [F]") - # Overpotentials eta_n = v_s_n - U_n eta_p = v_s_p - U_p @@ -231,9 +228,21 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): + (2 * RT_F * (1 - t_plus)) * pybamm.grad(sto_e_p) / sto_e_p ) - # Electrode surface potentials - self.rhs[v_s_n] = (l_n * Q_e * pybamm.div(i_e_n) - 3 * Q_th_n * j_n) / C_n - self.rhs[v_s_p] = (l_p * Q_e * pybamm.div(i_e_p) - 3 * Q_th_p * j_p) / C_p + ###################### + # Double layer + ###################### + if include_double_layer: + # Additional parameters + C_p = Parameter("Positive electrode capacitance [F]") + C_n = Parameter("Negative electrode capacitance [F]") + + # Electrode surface potentials + self.rhs[v_s_n] = (l_n * Q_e * pybamm.div(i_e_n) - 3 * Q_th_n * j_n) / C_n + self.rhs[v_s_p] = (l_p * Q_e * pybamm.div(i_e_p) - 3 * Q_th_p * j_p) / C_p + else: + # Electrode surface potentials + self.algebraic[v_s_n] = l_n * Q_e * pybamm.div(i_e_n) - 3 * Q_th_n * j_n + self.algebraic[v_s_p] = l_p * Q_e * pybamm.div(i_e_p) - 3 * Q_th_p * j_p self.initial_conditions[v_s_n] = U_n_init self.initial_conditions[v_s_p] = U_p_init diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 581defe71..64fc1513e 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -62,6 +62,9 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): Q = Variable("Discharge capacity [A.h]") Qt = Variable("Throughput capacity [A.h]") + v_s_n = Variable("Negative particle surface voltage variable [V]") + v_s_p = Variable("Positive particle surface voltage variable [V]") + # Variables that vary spatially are created with a domain sto_n = Variable( "Negative particle stoichiometry", @@ -158,48 +161,39 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): # Primary broadcasts are used to broadcast scalar quantities across a domain # into a vector of the right shape, for multiplying with other vectors alpha = 0.5 # cathodic transfer coefficient + + # Reference exchange current j0_n = sto_n_surf**alpha * (1 - sto_n_surf) ** (1 - alpha) / tau_ct_n j0_p = sto_p_surf**alpha * (1 - sto_p_surf) ** (1 - alpha) / tau_ct_p - if not include_double_layer: - # Assuming alpha = 0.5 - j_n = PrimaryBroadcast(I / (3 * Q_th_n), "negative electrode") - j_p = PrimaryBroadcast(-I / (3 * Q_th_p), "positive electrode") - eta_n = 2 * RT_F * pybamm.arcsinh(j_n / (2 * j0_n)) - eta_p = 2 * RT_F * pybamm.arcsinh(j_p / (2 * j0_p)) - v_s_n = pybamm.x_average(eta_n + U_n) - v_s_p = pybamm.x_average(eta_p + U_p) + + # Overpotentials + eta_n = PrimaryBroadcast(v_s_n - U_n, "negative electrode") + eta_p = PrimaryBroadcast(v_s_p - U_p, "positive electrode") + + # Exchange current + j_n = j0_n * ( + pybamm.exp((1 - alpha) * eta_n / RT_F) - pybamm.exp(-alpha * eta_n / RT_F) + ) + j_p = j0_p * ( + pybamm.exp((1 - alpha) * eta_p / RT_F) - pybamm.exp(-alpha * eta_p / RT_F) + ) ###################### # Double layer ###################### if include_double_layer: - # Additional variables - v_s_n = Variable("Negative particle surface voltage variable [V]") - v_s_p = Variable("Positive particle surface voltage variable [V]") - # Additional parameters C_p = Parameter("Positive electrode capacitance [F]") C_n = Parameter("Negative electrode capacitance [F]") - # Overpotentials - eta_n = PrimaryBroadcast(v_s_n - U_n, "negative electrode") - eta_p = PrimaryBroadcast(v_s_p - U_p, "positive electrode") - - # Exchange current - j_n = j0_n * ( - pybamm.exp((1 - alpha) * eta_n / RT_F) - - pybamm.exp(-alpha * eta_n / RT_F) - ) - j_p = j0_p * ( - pybamm.exp((1 - alpha) * eta_p / RT_F) - - pybamm.exp(-alpha * eta_p / RT_F) - ) - - # Electrode surface potentials self.rhs[v_s_n] = 1 / C_n * (I - 3 * Q_th_n * pybamm.x_average(j_n)) self.rhs[v_s_p] = 1 / C_p * (-I - 3 * Q_th_p * pybamm.x_average(j_p)) - self.initial_conditions[v_s_n] = U_n_init - self.initial_conditions[v_s_p] = U_p_init + else: + self.algebraic[v_s_n] = I - 3 * Q_th_n * pybamm.x_average(j_n) + self.algebraic[v_s_p] = -I - 3 * Q_th_p * pybamm.x_average(j_p) + + self.initial_conditions[v_s_n] = U_n_init + self.initial_conditions[v_s_p] = U_p_init ###################### # Particles diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index e19d6cf71..46554ddf8 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -66,6 +66,9 @@ def __init__( Q = Variable("Discharge capacity [A.h]") Qt = Variable("Throughput capacity [A.h]") + v_s_n = Variable("Negative particle surface voltage variable [V]") + v_s_p = Variable("Positive particle surface voltage variable [V]") + # Variables that vary spatially are created with a domain sto_n = Variable( "Negative particle stoichiometry", @@ -209,56 +212,47 @@ def __init__( # Primary broadcasts are used to broadcast scalar quantities across a domain # into a vector of the right shape, for multiplying with other vectors alpha = 0.5 # cathodic transfer coefficient + + # Reference exchange current j0_n = ( sto_n_surf**alpha * (sto_e_n * (1 - sto_n_surf)) ** (1 - alpha) / tau_ct_n ) j0_p = ( sto_p_surf**alpha * (sto_e_p * (1 - sto_p_surf)) ** (1 - alpha) / tau_ct_p ) - if not include_double_layer: - # Assuming alpha = 0.5 - j_n = PrimaryBroadcast(I / (3 * Q_th_n), "negative electrode") - j_p = PrimaryBroadcast(-I / (3 * Q_th_p), "positive electrode") - eta_n = 2 * RT_F * pybamm.arcsinh(j_n / (2 * j0_n)) - eta_p = 2 * RT_F * pybamm.arcsinh(j_p / (2 * j0_p)) - v_s_n = pybamm.x_average(eta_n + U_n) - v_s_p = pybamm.x_average(eta_p + U_p) + + # Overpotentials + eta_n = (v_s_n - U_n) + (2 * RT_F * (1 - t_plus)) * ( + pybamm.x_average(pybamm.log(sto_e_n)) - pybamm.log(sto_e_n) + ) + eta_p = (v_s_p - U_p) + (2 * RT_F * (1 - t_plus)) * ( + pybamm.x_average(pybamm.log(sto_e_p)) - pybamm.log(sto_e_p) + ) + + # Exchange current + j_n = j0_n * ( + pybamm.exp((1 - alpha) * eta_n / RT_F) - pybamm.exp(-alpha * eta_n / RT_F) + ) + j_p = j0_p * ( + pybamm.exp((1 - alpha) * eta_p / RT_F) - pybamm.exp(-alpha * eta_p / RT_F) + ) ###################### # Double layer ###################### if include_double_layer: - # Additional variables - v_s_n = Variable("Negative particle surface voltage variable [V]") - v_s_p = Variable("Positive particle surface voltage variable [V]") - # Additional parameters C_p = Parameter("Positive electrode capacitance [F]") C_n = Parameter("Negative electrode capacitance [F]") - # Overpotentials - eta_n = (v_s_n - U_n) + (2 * RT_F * (1 - t_plus)) * ( - pybamm.x_average(pybamm.log(sto_e_n)) - pybamm.log(sto_e_n) - ) - eta_p = (v_s_p - U_p) + (2 * RT_F * (1 - t_plus)) * ( - pybamm.x_average(pybamm.log(sto_e_p)) - pybamm.log(sto_e_p) - ) - - # Exchange current - j_n = j0_n * ( - pybamm.exp((1 - alpha) * eta_n / RT_F) - - pybamm.exp(-alpha * eta_n / RT_F) - ) - j_p = j0_p * ( - pybamm.exp((1 - alpha) * eta_p / RT_F) - - pybamm.exp(-alpha * eta_p / RT_F) - ) - - # Electrode surface potentials self.rhs[v_s_n] = 1 / C_n * (I - 3 * Q_th_n * pybamm.x_average(j_n)) self.rhs[v_s_p] = 1 / C_p * (-I - 3 * Q_th_p * pybamm.x_average(j_p)) - self.initial_conditions[v_s_n] = U_n_init - self.initial_conditions[v_s_p] = U_p_init + else: + self.algebraic[v_s_n] = I - 3 * Q_th_n * pybamm.x_average(j_n) + self.algebraic[v_s_p] = -I - 3 * Q_th_p * pybamm.x_average(j_p) + + self.initial_conditions[v_s_n] = U_n_init + self.initial_conditions[v_s_p] = U_p_init ###################### # Particles diff --git a/tests/integration/models/test_grouped_models.py b/tests/integration/models/test_grouped_models.py index e8dcf0596..8946831ee 100644 --- a/tests/integration/models/test_grouped_models.py +++ b/tests/integration/models/test_grouped_models.py @@ -33,6 +33,7 @@ class TestGroupedModels: pybop.lithium_ion.GroupedSPMe(), pybop.lithium_ion.GroupedSPMe(options={"surface form": "differential"}), pybop.lithium_ion.GroupedDFN(), + pybop.lithium_ion.GroupedDFN(options={"surface form": "differential"}), ], scope="module", ) diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index 15ab0a97c..2f4652f1d 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -25,6 +25,7 @@ class TestModels: pybop.lithium_ion.GroupedSPMe(), pybop.lithium_ion.GroupedSPMe(options={"surface form": "differential"}), pybop.lithium_ion.GroupedDFN(), + pybop.lithium_ion.GroupedDFN(options={"surface form": "differential"}), ], scope="module", ) From 1447a0102cbd8ea47565c3b64d7abd7fbe33a44b Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 28 May 2026 14:36:20 +0100 Subject: [PATCH 14/37] Refactor Butler Volmer --- pybop/models/lithium_ion/grouped_dfn.py | 27 ++++++++++-------------- pybop/models/lithium_ion/grouped_spm.py | 23 ++++++++++---------- pybop/models/lithium_ion/grouped_spme.py | 27 ++++++++++-------------- 3 files changed, 33 insertions(+), 44 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 70b2c2384..cc3698ca3 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -196,27 +196,14 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): ###################### # Primary broadcasts are used to broadcast scalar quantities across a domain # into a vector of the right shape, for multiplying with other vectors - alpha = 0.5 # cathodic transfer coefficient - - # Reference exchange current - j0_n = ( - sto_n_surf**alpha * (sto_e_n * (1 - sto_n_surf)) ** (1 - alpha) / tau_ct_n - ) - j0_p = ( - sto_p_surf**alpha * (sto_e_p * (1 - sto_p_surf)) ** (1 - alpha) / tau_ct_p - ) # Overpotentials eta_n = v_s_n - U_n eta_p = v_s_p - U_p - # Exchange current - j_n = j0_n * ( - pybamm.exp((1 - alpha) * eta_n / RT_F) - pybamm.exp(-alpha * eta_n / RT_F) - ) - j_p = j0_p * ( - pybamm.exp((1 - alpha) * eta_p / RT_F) - pybamm.exp(-alpha * eta_p / RT_F) - ) + # Exchange rates + j_n = self.butler_volmer(sto_n_surf, sto_e_n, eta_n / RT_F) / tau_ct_n + j_p = self.butler_volmer(sto_p_surf, sto_e_p, eta_p / RT_F) / tau_ct_p # Electrolyte current i_e_n = (beta_n * gamma_e) * ( @@ -465,6 +452,14 @@ def tau_d(self, sto, domain): inputs = {f"{Domain} particle surface stoichiometry": sto} return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) + def butler_volmer(self, sto_surf, sto_e, eta_RT_F): + """ + Dimensionless Butler-Volmer exchange rate. + """ + alpha = 0.5 # cathodic transfer coefficient + j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) + @property def default_parameter_values(self) -> pybamm.ParameterValues: param = ParameterValues("Chen2020") diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 64fc1513e..dda80a0f0 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -160,23 +160,14 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): ###################### # Primary broadcasts are used to broadcast scalar quantities across a domain # into a vector of the right shape, for multiplying with other vectors - alpha = 0.5 # cathodic transfer coefficient - - # Reference exchange current - j0_n = sto_n_surf**alpha * (1 - sto_n_surf) ** (1 - alpha) / tau_ct_n - j0_p = sto_p_surf**alpha * (1 - sto_p_surf) ** (1 - alpha) / tau_ct_p # Overpotentials eta_n = PrimaryBroadcast(v_s_n - U_n, "negative electrode") eta_p = PrimaryBroadcast(v_s_p - U_p, "positive electrode") - # Exchange current - j_n = j0_n * ( - pybamm.exp((1 - alpha) * eta_n / RT_F) - pybamm.exp(-alpha * eta_n / RT_F) - ) - j_p = j0_p * ( - pybamm.exp((1 - alpha) * eta_p / RT_F) - pybamm.exp(-alpha * eta_p / RT_F) - ) + # Exchange rates + j_n = self.butler_volmer(sto_n_surf, eta_n / RT_F) / tau_ct_n + j_p = self.butler_volmer(sto_p_surf, eta_p / RT_F) / tau_ct_p ###################### # Double layer @@ -336,6 +327,14 @@ def tau_d(self, sto, domain): inputs = {f"{Domain} particle surface stoichiometry": sto} return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) + def butler_volmer(self, sto_surf, eta_RT_F): + """ + Dimensionless Butler-Volmer exchange rate. + """ + alpha = 0.5 # cathodic transfer coefficient + j0 = sto_surf**alpha * (1 - sto_surf) ** (1 - alpha) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) + @property def default_parameter_values(self) -> ParameterValues: param = ParameterValues("Chen2020") diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 46554ddf8..ff626f77c 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -211,15 +211,6 @@ def __init__( ###################### # Primary broadcasts are used to broadcast scalar quantities across a domain # into a vector of the right shape, for multiplying with other vectors - alpha = 0.5 # cathodic transfer coefficient - - # Reference exchange current - j0_n = ( - sto_n_surf**alpha * (sto_e_n * (1 - sto_n_surf)) ** (1 - alpha) / tau_ct_n - ) - j0_p = ( - sto_p_surf**alpha * (sto_e_p * (1 - sto_p_surf)) ** (1 - alpha) / tau_ct_p - ) # Overpotentials eta_n = (v_s_n - U_n) + (2 * RT_F * (1 - t_plus)) * ( @@ -229,13 +220,9 @@ def __init__( pybamm.x_average(pybamm.log(sto_e_p)) - pybamm.log(sto_e_p) ) - # Exchange current - j_n = j0_n * ( - pybamm.exp((1 - alpha) * eta_n / RT_F) - pybamm.exp(-alpha * eta_n / RT_F) - ) - j_p = j0_p * ( - pybamm.exp((1 - alpha) * eta_p / RT_F) - pybamm.exp(-alpha * eta_p / RT_F) - ) + # Exchange rates + j_n = self.butler_volmer(sto_n_surf, sto_e_n, eta_n / RT_F) / tau_ct_n + j_p = self.butler_volmer(sto_p_surf, sto_e_p, eta_p / RT_F) / tau_ct_p ###################### # Double layer @@ -456,6 +443,14 @@ def tau_d(self, sto, domain): inputs = {f"{Domain} particle surface stoichiometry": sto} return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) + def butler_volmer(self, sto_surf, sto_e, eta_RT_F): + """ + Dimensionless Butler-Volmer exchange rate. + """ + alpha = 0.5 # cathodic transfer coefficient + j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) + @property def default_parameter_values(self) -> ParameterValues: param = ParameterValues("Chen2020") From 21fb48fceaf7de8b599a40eaf00bf163684c0c7b Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 28 May 2026 15:00:24 +0100 Subject: [PATCH 15/37] Move exchange rate to parameter set --- pybop/models/lithium_ion/grouped_dfn.py | 31 ++++++++++++++++++------ pybop/models/lithium_ion/grouped_spm.py | 30 +++++++++++++++++------ pybop/models/lithium_ion/grouped_spme.py | 31 ++++++++++++++++++------ 3 files changed, 71 insertions(+), 21 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index cc3698ca3..da9b42890 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -202,8 +202,8 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): eta_p = v_s_p - U_p # Exchange rates - j_n = self.butler_volmer(sto_n_surf, sto_e_n, eta_n / RT_F) / tau_ct_n - j_p = self.butler_volmer(sto_p_surf, sto_e_p, eta_p / RT_F) / tau_ct_p + j_n = self.j(sto_n_surf, sto_e_n, eta_n / RT_F, "negative") / tau_ct_n + j_p = self.j(sto_p_surf, sto_e_p, eta_p / RT_F, "positive") / tau_ct_p # Electrolyte current i_e_n = (beta_n * gamma_e) * ( @@ -452,13 +452,19 @@ def tau_d(self, sto, domain): inputs = {f"{Domain} particle surface stoichiometry": sto} return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) - def butler_volmer(self, sto_surf, sto_e, eta_RT_F): + def j(self, sto_surf, sto_e, eta_RT_F, domain): """ - Dimensionless Butler-Volmer exchange rate. + Dimensionless exchange rate. """ - alpha = 0.5 # cathodic transfer coefficient - j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) - return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) + Domain = domain.capitalize() + inputs = { + f"{Domain} particle surface stoichiometry": sto_surf, + f"{Domain} electrode electrolyte stoichiometry": sto_e, + f"{Domain} electrode dimensionless overpotential": eta_RT_F, + } + return FunctionParameter( + f"{Domain} electrode dimensionless exchange rate", inputs + ) @property def default_parameter_values(self) -> pybamm.ParameterValues: @@ -691,6 +697,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Positive particle diffusion time scale [s]": tau_d_p, "Negative particle diffusion time scale [s]": tau_d_n, "Electrolyte diffusion time scale [s]": tau_e, + "Positive electrode dimensionless exchange rate": GroupedDFN.symmetric_butler_volmer, + "Negative electrode dimensionless exchange rate": GroupedDFN.symmetric_butler_volmer, "Positive electrode charge transfer time scale [s]": tau_ct_p, "Negative electrode charge transfer time scale [s]": tau_ct_n, "Positive electrode capacitance [F]": C_p, @@ -703,3 +711,12 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_values = ParameterValues(values=parameter_dictionary) parameter_values._set_initial_state = GroupedDFN.set_initial_state # noqa: SLF001 return parameter_values + + @staticmethod + def symmetric_butler_volmer(sto_surf, sto_e, eta_RT_F): + """ + Dimensionless Butler-Volmer exchange rate. + """ + alpha = 0.5 # cathodic transfer coefficient + j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index dda80a0f0..020558ca4 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -166,8 +166,8 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): eta_p = PrimaryBroadcast(v_s_p - U_p, "positive electrode") # Exchange rates - j_n = self.butler_volmer(sto_n_surf, eta_n / RT_F) / tau_ct_n - j_p = self.butler_volmer(sto_p_surf, eta_p / RT_F) / tau_ct_p + j_n = self.j(sto_n_surf, eta_n / RT_F, "negative") / tau_ct_n + j_p = self.j(sto_p_surf, eta_p / RT_F, "positive") / tau_ct_p ###################### # Double layer @@ -327,13 +327,18 @@ def tau_d(self, sto, domain): inputs = {f"{Domain} particle surface stoichiometry": sto} return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) - def butler_volmer(self, sto_surf, eta_RT_F): + def j(self, sto_surf, eta_RT_F, domain): """ - Dimensionless Butler-Volmer exchange rate. + Dimensionless exchange rate. """ - alpha = 0.5 # cathodic transfer coefficient - j0 = sto_surf**alpha * (1 - sto_surf) ** (1 - alpha) - return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) + Domain = domain.capitalize() + inputs = { + f"{Domain} particle surface stoichiometry": sto_surf, + f"{Domain} electrode dimensionless overpotential": eta_RT_F, + } + return FunctionParameter( + f"{Domain} electrode dimensionless exchange rate", inputs + ) @property def default_parameter_values(self) -> ParameterValues: @@ -545,6 +550,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Measured cell capacity [A.s]": Q_meas, "Positive particle diffusion time scale [s]": tau_d_p, "Negative particle diffusion time scale [s]": tau_d_n, + "Positive electrode dimensionless exchange rate": GroupedSPM.symmetric_butler_volmer, + "Negative electrode dimensionless exchange rate": GroupedSPM.symmetric_butler_volmer, "Positive electrode charge transfer time scale [s]": tau_ct_p, "Negative electrode charge transfer time scale [s]": tau_ct_n, "Positive electrode capacitance [F]": C_p, @@ -556,3 +563,12 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_values = ParameterValues(values=parameter_dictionary) parameter_values._set_initial_state = GroupedSPM.set_initial_state # noqa: SLF001 return parameter_values + + @staticmethod + def symmetric_butler_volmer(sto_surf, eta_RT_F): + """ + Dimensionless Butler-Volmer exchange rate. + """ + alpha = 0.5 # cathodic transfer coefficient + j0 = sto_surf**alpha * (1 - sto_surf) ** (1 - alpha) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index ff626f77c..d126efcd8 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -221,8 +221,8 @@ def __init__( ) # Exchange rates - j_n = self.butler_volmer(sto_n_surf, sto_e_n, eta_n / RT_F) / tau_ct_n - j_p = self.butler_volmer(sto_p_surf, sto_e_p, eta_p / RT_F) / tau_ct_p + j_n = self.j(sto_n_surf, sto_e_n, eta_n / RT_F, "negative") / tau_ct_n + j_p = self.j(sto_p_surf, sto_e_p, eta_p / RT_F, "positive") / tau_ct_p ###################### # Double layer @@ -443,13 +443,19 @@ def tau_d(self, sto, domain): inputs = {f"{Domain} particle surface stoichiometry": sto} return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) - def butler_volmer(self, sto_surf, sto_e, eta_RT_F): + def j(self, sto_surf, sto_e, eta_RT_F, domain): """ - Dimensionless Butler-Volmer exchange rate. + Dimensionless exchange rate. """ - alpha = 0.5 # cathodic transfer coefficient - j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) - return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) + Domain = domain.capitalize() + inputs = { + f"{Domain} particle surface stoichiometry": sto_surf, + f"{Domain} electrode electrolyte stoichiometry": sto_e, + f"{Domain} electrode dimensionless overpotential": eta_RT_F, + } + return FunctionParameter( + f"{Domain} electrode dimensionless exchange rate", inputs + ) @property def default_parameter_values(self) -> ParameterValues: @@ -682,6 +688,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Electrolyte diffusion time scale [s]": tau_e, "Positive electrode relative transport efficiency": beta_p, "Negative electrode relative transport efficiency": beta_n, + "Positive electrode dimensionless exchange rate": GroupedSPMe.symmetric_butler_volmer, + "Negative electrode dimensionless exchange rate": GroupedSPMe.symmetric_butler_volmer, "Positive electrode charge transfer time scale [s]": tau_ct_p, "Negative electrode charge transfer time scale [s]": tau_ct_n, "Positive electrode capacitance [F]": C_p, @@ -694,3 +702,12 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_values = ParameterValues(values=parameter_dictionary) parameter_values._set_initial_state = GroupedSPMe.set_initial_state # noqa: SLF001 return parameter_values + + @staticmethod + def symmetric_butler_volmer(sto_surf, sto_e, eta_RT_F): + """ + Dimensionless Butler-Volmer exchange rate. + """ + alpha = 0.5 # cathodic transfer coefficient + j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) From 27920e358836b124cdae3ce71845917290a84159 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Wed, 3 Jun 2026 23:09:47 +0100 Subject: [PATCH 16/37] Add Butler-Volmer variations --- pybop/models/lithium_ion/grouped_dfn.py | 52 +++++++++++++++++++++++- pybop/models/lithium_ion/grouped_spm.py | 48 +++++++++++++++++++++- pybop/models/lithium_ion/grouped_spme.py | 52 +++++++++++++++++++++++- 3 files changed, 149 insertions(+), 3 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index da9b42890..a12609ac1 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -717,6 +717,56 @@ def symmetric_butler_volmer(sto_surf, sto_e, eta_RT_F): """ Dimensionless Butler-Volmer exchange rate. """ - alpha = 0.5 # cathodic transfer coefficient + j0 = (sto_surf * sto_e * (1 - sto_surf)) ** 0.5 + return 2 * j0 * pybamm.sinh(0.5 * eta_RT_F) + + @staticmethod + def get_asymmetric_butler_volmer(domain: str): + """ + Get the asymmetric Butler-Volmer exchange rate function for the "positive" or + "negative" domain. + """ + Domain = domain.capitalize() + alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") + + return AsymmetricButlerVolmer(alpha) + + @staticmethod + def get_multiphase_butler_volmer(domain: str): + """ + Get the multiphase Butler-Volmer exchange rate function for the "positive" or + "negative" domain. + """ + Domain = domain.capitalize() + alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") + omega = pybamm.Parameter(f"{Domain} electrode charge transfer ideality factor") + + return MultiphaseButlerVolmer(alpha, omega) + + +""" Alternative, dimensionless exchange rate functions written as classes to allow pickling. """ + + +class AsymmetricButlerVolmer: + def __init__(self, alpha): + self.alpha = alpha # cathodic transfer coefficient + + def __call__(self, sto_surf, sto_e, eta_RT_F): + alpha = self.alpha j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) + + +class MultiphaseButlerVolmer: + def __init__(self, alpha, omega): + self.alpha = alpha + self.omega = omega + + def __call__(self, sto_surf, sto_e, eta_RT_F): + alpha, omega = self.alpha, self.omega + j0 = ( + sto_surf ** (alpha * omega) + * (1 - sto_surf) ** ((1 - alpha) * omega) + * sto_e ** (1 - alpha) + ) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 020558ca4..192f23eef 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -569,6 +569,52 @@ def symmetric_butler_volmer(sto_surf, eta_RT_F): """ Dimensionless Butler-Volmer exchange rate. """ - alpha = 0.5 # cathodic transfer coefficient + j0 = (sto_surf * (1 - sto_surf)) ** 0.5 + return 2 * j0 * pybamm.sinh(0.5 * eta_RT_F) + + @staticmethod + def get_asymmetric_butler_volmer(domain: str): + """ + Get the asymmetric Butler-Volmer exchange rate function for the "positive" or + "negative" domain. + """ + Domain = domain.capitalize() + alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") + + return AsymmetricButlerVolmer(alpha) + + @staticmethod + def get_multiphase_butler_volmer(domain: str): + """ + Get the multiphase Butler-Volmer exchange rate function for the "positive" or + "negative" domain. + """ + Domain = domain.capitalize() + alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") + omega = pybamm.Parameter(f"{Domain} electrode charge transfer ideality factor") + + return MultiphaseButlerVolmer(alpha, omega) + + +""" Alternative, dimensionless exchange rate functions written as classes to allow pickling. """ + + +class AsymmetricButlerVolmer: + def __init__(self, alpha): + self.alpha = alpha # cathodic transfer coefficient + + def __call__(self, sto_surf, eta_RT_F): + alpha = self.alpha j0 = sto_surf**alpha * (1 - sto_surf) ** (1 - alpha) return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) + + +class MultiphaseButlerVolmer: + def __init__(self, alpha, omega): + self.alpha = alpha + self.omega = omega + + def __call__(self, sto_surf, eta_RT_F): + alpha, omega = self.alpha, self.omega + j0 = sto_surf ** (alpha * omega) * (1 - sto_surf) ** ((1 - alpha) * omega) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index d126efcd8..2b553c7b4 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -708,6 +708,56 @@ def symmetric_butler_volmer(sto_surf, sto_e, eta_RT_F): """ Dimensionless Butler-Volmer exchange rate. """ - alpha = 0.5 # cathodic transfer coefficient + j0 = (sto_surf * sto_e * (1 - sto_surf)) ** 0.5 + return 2 * j0 * pybamm.sinh(0.5 * eta_RT_F) + + @staticmethod + def get_asymmetric_butler_volmer(domain: str): + """ + Get the asymmetric Butler-Volmer exchange rate function for the "positive" or + "negative" domain. + """ + Domain = domain.capitalize() + alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") + + return AsymmetricButlerVolmer(alpha) + + @staticmethod + def get_multiphase_butler_volmer(domain: str): + """ + Get the multiphase Butler-Volmer exchange rate function for the "positive" or + "negative" domain. + """ + Domain = domain.capitalize() + alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") + omega = pybamm.Parameter(f"{Domain} electrode charge transfer ideality factor") + + return MultiphaseButlerVolmer(alpha, omega) + + +""" Alternative, dimensionless exchange rate functions written as classes to allow pickling. """ + + +class AsymmetricButlerVolmer: + def __init__(self, alpha): + self.alpha = alpha # cathodic transfer coefficient + + def __call__(self, sto_surf, sto_e, eta_RT_F): + alpha = self.alpha j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) + + +class MultiphaseButlerVolmer: + def __init__(self, alpha, omega): + self.alpha = alpha + self.omega = omega + + def __call__(self, sto_surf, sto_e, eta_RT_F): + alpha, omega = self.alpha, self.omega + j0 = ( + sto_surf ** (alpha * omega) + * (1 - sto_surf) ** ((1 - alpha) * omega) + * sto_e ** (1 - alpha) + ) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) From 3abd7145953515a048e00d5d26af90fe40605676 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 4 Jun 2026 16:13:14 +0100 Subject: [PATCH 17/37] Squeeze array to get singular cost --- pybop/models/lithium_ion/base_model.py | 3 +-- pybop/models/lithium_ion/sp_diffusion.py | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/pybop/models/lithium_ion/base_model.py b/pybop/models/lithium_ion/base_model.py index 2f4733060..0cec8753a 100644 --- a/pybop/models/lithium_ion/base_model.py +++ b/pybop/models/lithium_ion/base_model.py @@ -104,8 +104,7 @@ def ocv_function(soc): "Negative electrode OCP [V]", {"Negative particle stoichiometry": sto_n}, ) - - return parameter_values.evaluate(U_p - U_n, inputs=inputs) + return parameter_values.evaluate(U_p - U_n, inputs=inputs).squeeze() inverse_ocv = InverseOCV(ocv_function) soc = inverse_ocv(V_init) diff --git a/pybop/models/lithium_ion/sp_diffusion.py b/pybop/models/lithium_ion/sp_diffusion.py index 1ba955034..0e3504072 100644 --- a/pybop/models/lithium_ion/sp_diffusion.py +++ b/pybop/models/lithium_ion/sp_diffusion.py @@ -283,8 +283,7 @@ def ocv_function(sto): U = FunctionParameter( "Electrode OCP [V]", {"Particle stoichiometry": sto} ) - - return parameter_values.evaluate(U, inputs=inputs) + return parameter_values.evaluate(U, inputs=inputs).squeeze() inverse_ocv = InverseOCV(ocv_function) sto = inverse_ocv(V_init) From 25a4eb462055cd05692610a8f6abb1b61472594f Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 6 Aug 2026 20:19:22 +0100 Subject: [PATCH 18/37] Test simulation of experiment --- pybop/models/lithium_ion/grouped_dfn.py | 2 +- pybop/models/lithium_ion/grouped_spm.py | 2 +- pybop/models/lithium_ion/grouped_spme.py | 2 +- tests/unit/test_models.py | 12 +++++++++++- 4 files changed, 14 insertions(+), 4 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index a12609ac1..4441e247e 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -362,6 +362,7 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): "positive", ) voltage_components = { + "Battery voltage [V]": V, "Battery open-circuit voltage [V]": ocp_p_bulk - ocp_n_bulk, "Battery particle concentration overpotential [V]": ( (pybamm.x_average(self.U(sto_p_surf, "positive")) - ocp_p_bulk) @@ -423,7 +424,6 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): "Discharge capacity [A.h]": Q, "Throughput capacity [A.h]": Qt, "Voltage [V]": V, - "Battery voltage [V]": V, "Open-circuit voltage [V]": pybamm.boundary_value(U_p, "right") - pybamm.boundary_value(U_n, "left"), **voltage_components, diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 192f23eef..6da6cdb71 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -234,6 +234,7 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): ocp_n_bulk = self.U(pybamm.x_average(pybamm.r_average(sto_n)), "negative") ocp_p_bulk = self.U(pybamm.x_average(pybamm.r_average(sto_p)), "positive") voltage_components = { + "Battery voltage [V]": V, "Battery open-circuit voltage [V]": ocp_p_bulk - ocp_n_bulk, "Battery particle concentration overpotential [V]": ( (pybamm.x_average(self.U(sto_p_surf, "positive")) - ocp_p_bulk) @@ -299,7 +300,6 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): "Throughput capacity [A.h]": Qt, "Voltage [V]": V, "Voltage expression [V]": V, # for compatibility with "voltage as a state" - "Battery voltage [V]": V, "Open-circuit voltage [V]": U_p - U_n, **voltage_components, } diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 2b553c7b4..9d114c7ef 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -338,6 +338,7 @@ def __init__( "positive", ) voltage_components = { + "Battery voltage [V]": V, "Battery open-circuit voltage [V]": ocp_p_bulk - ocp_n_bulk, "Battery particle concentration overpotential [V]": ( (pybamm.x_average(self.U(sto_p_surf, "positive")) - ocp_p_bulk) @@ -415,7 +416,6 @@ def __init__( "Throughput capacity [A.h]": Qt, "Voltage [V]": V, "Voltage expression [V]": V, # for compatibility with "voltage as a state" - "Battery voltage [V]": V, "Open-circuit voltage [V]": U_p - U_n, **voltage_components, } diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index 2f4652f1d..276db3f9e 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -52,13 +52,23 @@ def test_model_simulation(self, model): assert isinstance(fig, pybamm.QuickPlot) if isinstance( - model, pybop.lithium_ion.GroupedSPMe | pybop.lithium_ion.GroupedDFN + model, + pybop.lithium_ion.GroupedSPM + | pybop.lithium_ion.GroupedSPMe + | pybop.lithium_ion.GroupedDFN, ): for split in [False, True]: fig, ax = solution.plot_voltage_components(split_by_electrode=split) assert isinstance(fig, Figure) assert isinstance(ax, Axes) + if not isinstance(model, pybop.ExponentialDecayModel): + experiment = pybamm.Experiment(["Discharge at 1C for 1 minute"]) + solution = pybamm.Simulation(model, experiment=experiment).solve( + calc_esoh=False + ) + assert solution["Time [s]"].data[-1] == 60 + def test_set_initial_state(self, model): if isinstance(model, pybop.ExponentialDecayModel): pass # Only testing the battery models for now From 2b7d0ef5a5bdd00ba194db1ca4448d270b0864d5 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Wed, 10 Jun 2026 23:59:53 +0100 Subject: [PATCH 19/37] Get reaction rates from exchange current functions --- pybop/models/lithium_ion/grouped_dfn.py | 13 +++++++++---- pybop/models/lithium_ion/grouped_spm.py | 13 +++++++++---- pybop/models/lithium_ion/grouped_spme.py | 13 +++++++++---- 3 files changed, 27 insertions(+), 12 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 4441e247e..3cf994f9d 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -585,6 +585,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Unpack physical parameters F = pybamm.constants.F.value + T = param["Ambient temperature [K]"] alpha_p = param["Positive electrode active material volume fraction"] alpha_n = param["Negative electrode active material volume fraction"] c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] @@ -601,8 +602,12 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] - m_p = 3.42e-6 # (A/m2)(m3/mol)**1.5 - m_n = 6.48e-7 # (A/m2)(m3/mol)**1.5 + m_p = param["Positive electrode exchange-current density [A.m-2]"]( + 1, 1, 2, T + ) # (A/m2)(m3/mol)**1.5 + m_n = param["Negative electrode exchange-current density [A.m-2]"]( + 1, 1, 2, T + ) # (A/m2)(m3/mol)**1.5 sigma_p = ( param["Positive electrode conductivity [S.m-1]"] * alpha_p ** param["Positive electrode Bruggeman coefficient (electrode)"] @@ -676,8 +681,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_dictionary = { "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], "Current function [A]": param["Current function [A]"], - "Ambient temperature [K]": param["Ambient temperature [K]"], - "Initial temperature [K]": param["Ambient temperature [K]"], + "Ambient temperature [K]": T, + "Initial temperature [K]": T, "Initial SoC": soc_init, "Minimum negative stoichiometry": x_0, "Maximum negative stoichiometry": x_100, diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 6da6cdb71..ab8ad6658 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -454,6 +454,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Unpack physical parameters F = pybamm.constants.F.value + T = param["Ambient temperature [K]"] alpha_p = param["Positive electrode active material volume fraction"] alpha_n = param["Negative electrode active material volume fraction"] c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] @@ -470,8 +471,12 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] - m_p = 3.42e-6 # (A/m2)(m3/mol)**1.5 - m_n = 6.48e-7 # (A/m2)(m3/mol)**1.5 + m_p = param["Positive electrode exchange-current density [A.m-2]"]( + 1, 1, 2, T + ) # (A/m2)(m3/mol)**1.5 + m_n = param["Negative electrode exchange-current density [A.m-2]"]( + 1, 1, 2, T + ) # (A/m2)(m3/mol)**1.5 sigma_p = ( param["Positive electrode conductivity [S.m-1]"] * alpha_p ** param["Positive electrode Bruggeman coefficient (electrode)"] @@ -536,8 +541,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_dictionary = { "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], "Current function [A]": param["Current function [A]"], - "Ambient temperature [K]": param["Ambient temperature [K]"], - "Initial temperature [K]": param["Ambient temperature [K]"], + "Ambient temperature [K]": T, + "Initial temperature [K]": T, "Initial SoC": soc_init, "Minimum negative stoichiometry": x_0, "Maximum negative stoichiometry": x_100, diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 9d114c7ef..7ff96cd46 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -576,6 +576,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Unpack physical parameters F = pybamm.constants.F.value + T = param["Ambient temperature [K]"] alpha_p = param["Positive electrode active material volume fraction"] alpha_n = param["Negative electrode active material volume fraction"] c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] @@ -592,8 +593,12 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] - m_p = 3.42e-6 # (A/m2)(m3/mol)**1.5 - m_n = 6.48e-7 # (A/m2)(m3/mol)**1.5 + m_p = param["Positive electrode exchange-current density [A.m-2]"]( + 1, 1, 2, T + ) # (A/m2)(m3/mol)**1.5 + m_n = param["Negative electrode exchange-current density [A.m-2]"]( + 1, 1, 2, T + ) # (A/m2)(m3/mol)**1.5 sigma_p = ( param["Positive electrode conductivity [S.m-1]"] * alpha_p ** param["Positive electrode Bruggeman coefficient (electrode)"] @@ -668,8 +673,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_dictionary = { "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], "Current function [A]": param["Current function [A]"], - "Ambient temperature [K]": param["Ambient temperature [K]"], - "Initial temperature [K]": param["Ambient temperature [K]"], + "Ambient temperature [K]": T, + "Initial temperature [K]": T, "Initial SoC": soc_init, "Minimum negative stoichiometry": x_0, "Maximum negative stoichiometry": x_100, From 1bb821e2917452ce5b8409d2f2003cd6fd1b7dc7 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 12 Jun 2026 20:15:45 +0100 Subject: [PATCH 20/37] Update plot variables --- pybop/models/lithium_ion/grouped_dfn.py | 13 ++++-------- pybop/models/lithium_ion/grouped_spm.py | 25 ++++++++++++++---------- pybop/models/lithium_ion/grouped_spme.py | 19 ++++++++---------- 3 files changed, 27 insertions(+), 30 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 3cf994f9d..41f9da1c2 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -227,7 +227,6 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): self.rhs[v_s_n] = (l_n * Q_e * pybamm.div(i_e_n) - 3 * Q_th_n * j_n) / C_n self.rhs[v_s_p] = (l_p * Q_e * pybamm.div(i_e_p) - 3 * Q_th_p * j_p) / C_p else: - # Electrode surface potentials self.algebraic[v_s_n] = l_n * Q_e * pybamm.div(i_e_n) - 3 * Q_th_n * j_n self.algebraic[v_s_p] = l_p * Q_e * pybamm.div(i_e_p) - 3 * Q_th_p * j_p @@ -302,13 +301,9 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): "right": (Scalar(0), "Neumann"), } - self.initial_conditions[sto_e_n] = PrimaryBroadcast( - Scalar(1), "negative electrode" - ) - self.initial_conditions[sto_e_sep] = PrimaryBroadcast(Scalar(1), "separator") - self.initial_conditions[sto_e_p] = PrimaryBroadcast( - Scalar(1), "positive electrode" - ) + self.initial_conditions[sto_e_n] = Scalar(1) + self.initial_conditions[sto_e_sep] = Scalar(1) + self.initial_conditions[sto_e_p] = Scalar(1) # Electrolyte overpotential eta_e = (2 * (1 - t_plus) * RT_F) * ( @@ -490,7 +485,7 @@ def default_quick_plot_variables(self): "Negative electrode potential [V]", "Negative particle surface voltage [V]", }, - {"Electrolyte scaled current density [s-1]"}, + "Electrolyte scaled current density [s-1]", { "Positive electrode potential [V]", "Positive particle surface voltage [V]", diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index ab8ad6658..4f49f43ce 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -177,8 +177,9 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): C_p = Parameter("Positive electrode capacitance [F]") C_n = Parameter("Negative electrode capacitance [F]") - self.rhs[v_s_n] = 1 / C_n * (I - 3 * Q_th_n * pybamm.x_average(j_n)) - self.rhs[v_s_p] = 1 / C_p * (-I - 3 * Q_th_p * pybamm.x_average(j_p)) + # Electrode surface potentials + self.rhs[v_s_n] = (I - 3 * Q_th_n * pybamm.x_average(j_n)) / C_n + self.rhs[v_s_p] = (-I - 3 * Q_th_p * pybamm.x_average(j_p)) / C_p else: self.algebraic[v_s_n] = I - 3 * Q_th_n * pybamm.x_average(j_n) self.algebraic[v_s_p] = -I - 3 * Q_th_p * pybamm.x_average(j_p) @@ -191,21 +192,25 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): ###################### # The div and grad operators will be converted to the appropriate matrix # multiplication at the discretisation stage - self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / self.tau_d(sto_n, "negative")) - self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / self.tau_d(sto_p, "positive")) + self.rhs[sto_n] = pybamm.div( + pybamm.grad(sto_n) / self.tau_d(sto_n, T, "negative") + ) + self.rhs[sto_p] = pybamm.div( + pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") + ) # Boundary conditions must be provided for equations with spatial derivatives self.boundary_conditions[sto_n] = { "left": (Scalar(0), "Neumann"), "right": ( - -self.tau_d(sto_n_surf, "negative") * pybamm.x_average(j_n), + -self.tau_d(sto_n_surf, T, "negative") * pybamm.x_average(j_n), "Neumann", ), } self.boundary_conditions[sto_p] = { "left": (Scalar(0), "Neumann"), "right": ( - -self.tau_d(sto_p_surf, "positive") * pybamm.x_average(j_p), + -self.tau_d(sto_p_surf, T, "positive") * pybamm.x_average(j_p), "Neumann", ), } @@ -319,12 +324,12 @@ def U(self, sto, domain): out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" return out - def tau_d(self, sto, domain): + def tau_d(self, sto, T, domain): """ Dimensional diffusion time scale [s]. """ Domain = domain.capitalize() - inputs = {f"{Domain} particle surface stoichiometry": sto} + inputs = {f"{Domain} particle surface stoichiometry": sto, "Temperature [K]": T} return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) def j(self, sto_surf, eta_RT_F, domain): @@ -353,8 +358,8 @@ def default_parameter_values(self) -> ParameterValues: @property def default_quick_plot_variables(self): return [ - "Negative particle surface stoichiometry", - "Positive particle surface stoichiometry", + "Negative particle stoichiometry", + "Positive particle stoichiometry", "Current [A]", { "Negative electrode potential [V]", diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 7ff96cd46..58b092de2 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -232,8 +232,9 @@ def __init__( C_p = Parameter("Positive electrode capacitance [F]") C_n = Parameter("Negative electrode capacitance [F]") - self.rhs[v_s_n] = 1 / C_n * (I - 3 * Q_th_n * pybamm.x_average(j_n)) - self.rhs[v_s_p] = 1 / C_p * (-I - 3 * Q_th_p * pybamm.x_average(j_p)) + # Electrode surface potentials + self.rhs[v_s_n] = (I - 3 * Q_th_n * pybamm.x_average(j_n)) / C_n + self.rhs[v_s_p] = (-I - 3 * Q_th_p * pybamm.x_average(j_p)) / C_p else: self.algebraic[v_s_n] = I - 3 * Q_th_n * pybamm.x_average(j_n) self.algebraic[v_s_p] = -I - 3 * Q_th_p * pybamm.x_average(j_p) @@ -301,13 +302,9 @@ def __init__( "right": (Scalar(0), "Neumann"), } - self.initial_conditions[sto_e_n] = PrimaryBroadcast( - Scalar(1), "negative electrode" - ) - self.initial_conditions[sto_e_sep] = PrimaryBroadcast(Scalar(1), "separator") - self.initial_conditions[sto_e_p] = PrimaryBroadcast( - Scalar(1), "positive electrode" - ) + self.initial_conditions[sto_e_n] = Scalar(1) + self.initial_conditions[sto_e_sep] = Scalar(1) + self.initial_conditions[sto_e_p] = Scalar(1) ###################### # Cell voltage @@ -473,9 +470,9 @@ def default_parameter_values(self) -> ParameterValues: @property def default_quick_plot_variables(self): return [ - "Negative particle surface stoichiometry", + "Negative particle stoichiometry", "Electrolyte stoichiometry", - "Positive particle surface stoichiometry", + "Positive particle stoichiometry", "Current [A]", { "Negative electrode potential [V]", From c5786e7cee20045d9c5d28d037c39f4bfd2f348f Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Tue, 14 Jul 2026 12:49:57 +0100 Subject: [PATCH 21/37] Update test --- pybop/models/lithium_ion/utils.py | 2 +- tests/unit/test_models.py | 13 +++++++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/pybop/models/lithium_ion/utils.py b/pybop/models/lithium_ion/utils.py index 6448306b0..e5402429c 100644 --- a/pybop/models/lithium_ion/utils.py +++ b/pybop/models/lithium_ion/utils.py @@ -133,7 +133,7 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): for x in inputs: diff = np.abs(ocv_function(x["Root"]) - self.ocv_value) sol = Solution() - sol.set_solution_variable("Difference", data=np.asarray([diff])) + sol.set_solution_variable("Difference", data=np.atleast_1d(diff)) solutions.append(sol) return solutions diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index 276db3f9e..423b4f136 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -6,6 +6,12 @@ import pybop +GROUPED_MODEL = ( + pybop.lithium_ion.GroupedSPM + | pybop.lithium_ion.GroupedSPMe + | pybop.lithium_ion.GroupedDFN +) + class TestModels: """ @@ -51,12 +57,7 @@ def test_model_simulation(self, model): fig = solution.plot() assert isinstance(fig, pybamm.QuickPlot) - if isinstance( - model, - pybop.lithium_ion.GroupedSPM - | pybop.lithium_ion.GroupedSPMe - | pybop.lithium_ion.GroupedDFN, - ): + if isinstance(model, GROUPED_MODEL): for split in [False, True]: fig, ax = solution.plot_voltage_components(split_by_electrode=split) assert isinstance(fig, Figure) From 2d7590b071687b2e8754e6074a51008a88d28661 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Tue, 14 Jul 2026 16:30:28 +0100 Subject: [PATCH 22/37] Update parameter setting --- pybop/models/lithium_ion/base_model.py | 4 ++-- pybop/models/lithium_ion/grouped_dfn.py | 8 ++++---- pybop/models/lithium_ion/grouped_spm.py | 8 ++++---- pybop/models/lithium_ion/grouped_spme.py | 8 ++++---- pybop/models/lithium_ion/sp_diffusion.py | 14 +++++++++++++- 5 files changed, 27 insertions(+), 15 deletions(-) diff --git a/pybop/models/lithium_ion/base_model.py b/pybop/models/lithium_ion/base_model.py index 0cec8753a..f555d9302 100644 --- a/pybop/models/lithium_ion/base_model.py +++ b/pybop/models/lithium_ion/base_model.py @@ -77,10 +77,10 @@ def set_initial_state( if isinstance(initial_value, str) and initial_value.endswith("V"): V_init = float(initial_value[:-1]) V_min = parameter_values.evaluate( - pybamm.Parameter("Lower voltage cut-off [V]"), inputs=inputs + Parameter("Lower voltage cut-off [V]"), inputs=inputs ) V_max = parameter_values.evaluate( - pybamm.Parameter("Upper voltage cut-off [V]"), inputs=inputs + Parameter("Upper voltage cut-off [V]"), inputs=inputs ) if not V_min - tol <= V_init <= V_max + tol: diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 41f9da1c2..6a7d8469c 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -597,11 +597,11 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] - m_p = param["Positive electrode exchange-current density [A.m-2]"]( - 1, 1, 2, T + m_p = param.evaluate( + param["Positive electrode exchange-current density [A.m-2]"](1, 1, 2, T) ) # (A/m2)(m3/mol)**1.5 - m_n = param["Negative electrode exchange-current density [A.m-2]"]( - 1, 1, 2, T + m_n = param.evaluate( + param["Negative electrode exchange-current density [A.m-2]"](1, 1, 2, T) ) # (A/m2)(m3/mol)**1.5 sigma_p = ( param["Positive electrode conductivity [S.m-1]"] diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 4f49f43ce..9acfad3bf 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -476,11 +476,11 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] - m_p = param["Positive electrode exchange-current density [A.m-2]"]( - 1, 1, 2, T + m_p = param.evaluate( + param["Positive electrode exchange-current density [A.m-2]"](1, 1, 2, T) ) # (A/m2)(m3/mol)**1.5 - m_n = param["Negative electrode exchange-current density [A.m-2]"]( - 1, 1, 2, T + m_n = param.evaluate( + param["Negative electrode exchange-current density [A.m-2]"](1, 1, 2, T) ) # (A/m2)(m3/mol)**1.5 sigma_p = ( param["Positive electrode conductivity [S.m-1]"] diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 58b092de2..67f63ae11 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -590,11 +590,11 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] - m_p = param["Positive electrode exchange-current density [A.m-2]"]( - 1, 1, 2, T + m_p = param.evaluate( + param["Positive electrode exchange-current density [A.m-2]"](1, 1, 2, T) ) # (A/m2)(m3/mol)**1.5 - m_n = param["Negative electrode exchange-current density [A.m-2]"]( - 1, 1, 2, T + m_n = param.evaluate( + param["Negative electrode exchange-current density [A.m-2]"](1, 1, 2, T) ) # (A/m2)(m3/mol)**1.5 sigma_p = ( param["Positive electrode conductivity [S.m-1]"] diff --git a/pybop/models/lithium_ion/sp_diffusion.py b/pybop/models/lithium_ion/sp_diffusion.py index 0e3504072..2bef6e533 100644 --- a/pybop/models/lithium_ion/sp_diffusion.py +++ b/pybop/models/lithium_ion/sp_diffusion.py @@ -1,3 +1,4 @@ +import numpy as np import pybamm from pybamm import ( Event, @@ -206,11 +207,15 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Unpack physical parameters F = pybamm.constants.F.value + T = param["Ambient temperature [K]"] alpha = param["Positive electrode active material volume fraction"] c_max = param["Maximum concentration in positive electrode [mol.m-3]"] L = param["Positive electrode thickness [m]"] R = param["Positive particle radius [m]"] D = param["Positive particle diffusivity [m2.s-1]"] + m = param.evaluate( + param["Positive electrode exchange-current density [A.m-2]"](1, 1, 2, T) + ) # (A/m2)(m3/mol)**1.5 sto_init = ( param["Initial concentration in positive electrode [mol.m-3]"] / c_max ) @@ -223,6 +228,13 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal Q_th = F * alpha * c_max * L * A tau_d = R**2 / D + # Estimate the series resistance, neglecting conductivity losses + RT_F = pybamm.constants.R.value * param["Ambient temperature [K]"] / F + ce0 = param["Initial concentration in electrolyte [mol.m-3]"] + tau_ct = F * R / (m * np.sqrt(ce0)) + Rct_typ = (2 * RT_F * tau_ct) / (3 * Q_th) + R0 = Rct_typ + param["Contact resistance [Ohm]"] + parameter_dictionary = { "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], "Current function [A]": param["Current function [A]"], @@ -230,7 +242,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Electrode OCP [V]": ocp, "Theoretical electrode capacity [A.s]": Q_th, "Particle diffusion time scale [s]": tau_d, - "Series resistance [Ohm]": 1, + "Series resistance [Ohm]": R0, } parameter_values = ParameterValues(values=parameter_dictionary) parameter_values._set_initial_state = SPDiffusion.set_initial_state # noqa: SLF001 From 0e317690f2332f3728b7924fbb4c437a3381ed0c Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Mon, 20 Jul 2026 20:12:50 +0100 Subject: [PATCH 23/37] Use positive domain name --- .../battery_parameterisation/gitt_fitting.py | 6 +- .../comparison_examples/gitt_models.py | 7 +- pybop/applications/gitt_methods.py | 22 ++- pybop/models/lithium_ion/sp_diffusion.py | 184 ++++++++++++------ pybop/models/lithium_ion/weppner_huggins.py | 26 +-- tests/integration/models/test_gitt_models.py | 2 +- tests/integration/test_applications.py | 11 +- tests/unit/test_models.py | 11 +- 8 files changed, 169 insertions(+), 100 deletions(-) diff --git a/examples/scripts/battery_parameterisation/gitt_fitting.py b/examples/scripts/battery_parameterisation/gitt_fitting.py index 6abe016f2..d6ea887eb 100644 --- a/examples/scripts/battery_parameterisation/gitt_fitting.py +++ b/examples/scripts/battery_parameterisation/gitt_fitting.py @@ -60,7 +60,9 @@ gitt_parameter_data = gitt_fit() # Plot the functional parameters -pybop.plot.dataset(gitt_parameter_data, signal=["Particle diffusion time scale [s]"]) +pybop.plot.dataset( + gitt_parameter_data, signal=["Positive particle diffusion time scale [s]"] +) pybop.plot.dataset(gitt_parameter_data, signal=["Series resistance [Ohm]"]) # Run the identified model @@ -76,7 +78,7 @@ # Return to the original model and update the diffusivity value diffusivity = np.mean( parameter_values["Positive particle radius [m]"] ** 2 - / gitt_parameter_data["Particle diffusion time scale [s]"] + / gitt_parameter_data["Positive particle diffusion time scale [s]"] ) parameter_values.update({"Positive particle diffusivity [m2.s-1]": diffusivity}) diff --git a/examples/scripts/comparison_examples/gitt_models.py b/examples/scripts/comparison_examples/gitt_models.py index d5df9d907..27071dbe3 100644 --- a/examples/scripts/comparison_examples/gitt_models.py +++ b/examples/scripts/comparison_examples/gitt_models.py @@ -61,7 +61,7 @@ # Fitting parameters grouped_parameter_values.update( { - "Particle diffusion time scale [s]": diffusion_parameter, + "Positive particle diffusion time scale [s]": diffusion_parameter, "Reference voltage [V]": pybop.Parameter( initial_value=grouped_parameter_values["Reference voltage [V]"], ), @@ -77,7 +77,7 @@ # Fitting parameters grouped_parameter_values.update( { - "Particle diffusion time scale [s]": diffusion_parameter, + "Positive particle diffusion time scale [s]": diffusion_parameter, "Series resistance [Ohm]": pybop.Parameter( initial_value=grouped_parameter_values["Series resistance [Ohm]"], ), @@ -103,7 +103,8 @@ result = optim.run() print(result) print( - "Diffusion time [s]:", result.best_inputs["Particle diffusion time scale [s]"] + "Diffusion time [s]:", + result.best_inputs["Positive particle diffusion time scale [s]"], ) # Plot the timeseries output diff --git a/pybop/applications/gitt_methods.py b/pybop/applications/gitt_methods.py index 02151c4d1..e580392af 100644 --- a/pybop/applications/gitt_methods.py +++ b/pybop/applications/gitt_methods.py @@ -36,7 +36,9 @@ def __init__( ): self.parameter_values = parameter_values.copy() self.parameters = { - "Particle diffusion time scale [s]": pybop.Parameter(bounds=[0, np.inf]), + "Positive particle diffusion time scale [s]": pybop.Parameter( + bounds=[0, np.inf] + ), "Series resistance [Ohm]": pybop.Parameter(bounds=[0, np.inf]), } self.cost = cost or pybop.RootMeanSquaredError @@ -110,7 +112,9 @@ def __init__( self.optimiser_options = optimiser_options or self.optimiser.default_options() # Set up OCV root-finding function - self.inverse_ocp = pybop.InverseOCV(parameter_values["Electrode OCP [V]"]) + self.inverse_ocp = pybop.InverseOCV( + parameter_values["Positive electrode OCP [V]"] + ) # Initialise single pulse fitter self.pulse_fit = GITTPulseFit( @@ -148,7 +152,9 @@ def __call__(self) -> pybop.Dataset: # Log the result self.pulses.append(pulse_result) diffusion_time.append( - pulse_result.best_inputs["Particle diffusion time scale [s]"] + pulse_result.best_inputs[ + "Positive particle diffusion time scale [s]" + ] ) series_resistance.append( pulse_result.best_inputs["Series resistance [Ohm]"] @@ -169,14 +175,16 @@ def __call__(self) -> pybop.Dataset: self.parameter_data = pybop.Dataset( { "Stoichiometry": np.asarray(stoichiometry), - "Particle diffusion time scale [s]": np.asarray(diffusion_time), + "Positive particle diffusion time scale [s]": np.asarray( + diffusion_time + ), "Series resistance [Ohm]": np.asarray(series_resistance), cost_name: np.asarray(best_cost), } if len(stoichiometry) > 1 and stoichiometry[-1] > stoichiometry[0] else { "Stoichiometry": np.flipud(np.asarray(stoichiometry)), - "Particle diffusion time scale [s]": np.flipud( + "Positive particle diffusion time scale [s]": np.flipud( np.asarray(diffusion_time) ), "Series resistance [Ohm]": np.flipud(np.asarray(series_resistance)), @@ -187,8 +195,8 @@ def __call__(self) -> pybop.Dataset: # Compute mean values self.best_inputs = { - "Particle diffusion time scale [s]": np.mean( - self.parameter_data["Particle diffusion time scale [s]"] + "Positive particle diffusion time scale [s]": np.mean( + self.parameter_data["Positive particle diffusion time scale [s]"] ), "Series resistance [Ohm]": np.mean( self.parameter_data["Series resistance [Ohm]"] diff --git a/pybop/models/lithium_ion/sp_diffusion.py b/pybop/models/lithium_ion/sp_diffusion.py index 2bef6e533..e6b384ada 100644 --- a/pybop/models/lithium_ion/sp_diffusion.py +++ b/pybop/models/lithium_ion/sp_diffusion.py @@ -11,6 +11,9 @@ Variable, ) from pybamm import t as pybamm_t +from pybamm.models.full_battery_models.lithium_ion.electrode_soh_half_cell import ( + get_min_max_stoichiometries, +) from pybop.models.lithium_ion.base_model import BaseGroupedModel from pybop.models.lithium_ion.utils import InverseOCV @@ -45,19 +48,19 @@ def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): Qt = Variable("Throughput capacity [A.h]") # Variables that vary spatially are created with a domain - sto = Variable("Particle stoichiometry", domain="particle") - sto_surf = pybamm.surf(sto) + sto_p = Variable("Positive particle stoichiometry", domain="positive particle") + sto_p_surf = pybamm.surf(sto_p) # Events specify points at which a solution should terminate tol = pybamm.settings.tolerances["U__c_s"] self.events += [ Event( - "Minimum particle surface stoichiometry", - pybamm.min(sto_surf) - tol, + "Minimum positive particle surface stoichiometry", + pybamm.min(sto_p_surf) - tol, ), Event( - "Maximum particle surface stoichiometry", - (1 - tol) - pybamm.max(sto_surf), + "Maximum positive particle surface stoichiometry", + (1 - tol) - pybamm.max(sto_p_surf), ), ] @@ -67,9 +70,12 @@ def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): # Parameters are purely symbolic at this stage, and will be set by the # `ParameterValues` class when the model is processed. - Q_th = Parameter("Theoretical electrode capacity [A.s]") + soc_init = Parameter("Initial SoC") + y_100 = Parameter("Minimum positive stoichiometry") + y_0 = Parameter("Maximum positive stoichiometry") - sto_init = Parameter("Initial stoichiometry") + # Grouped parameters + Q_th_p = Parameter("Measured cell capacity [A.s]") / (y_0 - y_100) ###################### # Input current (positive on discharge) @@ -92,32 +98,44 @@ def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): ###################### # The div and grad operators will be converted to the appropriate matrix # multiplication at the discretisation stage - self.rhs[sto] = pybamm.div(pybamm.grad(sto) / self.tau_d(sto)) + self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / self.tau_d(sto_p)) # Boundary conditions must be provided for equations with spatial derivatives - j = -I / (3 * Q_th) - self.boundary_conditions[sto] = { + j_p = -I / (3 * Q_th_p) + self.boundary_conditions[sto_p] = { "left": (Scalar(0), "Neumann"), - "right": (-self.tau_d(sto_surf) * j, "Neumann"), + "right": (-self.tau_d(sto_p_surf) * j_p, "Neumann"), } - self.initial_conditions[sto] = sto_init + sto_p_init = y_0 + (y_100 - y_0) * soc_init + self.initial_conditions[sto_p] = sto_p_init ###################### # Cell voltage ###################### - U = self.U(sto_surf) - V = U - self.R0(sto_surf) * I + sto_p_average = sto_p_init + Q * 3600 / Q_th_p # pybamm.r_average(sto_p) + U = self.U(sto_p_surf, "positive") - self.U(sto_p_average, "negative") + V = U - self.R0(sto_p_surf) * I # Save the initial OCV - self.param.ocv_init = self.U(sto_init) + self.param.ocv_init = self.U(sto_p_init, "positive") - self.U( + sto_p_init, "negative" + ) + + # Events specify points at which a solution should terminate + self.events += [ + Event("Minimum voltage [V]", V - self.param.voltage_low_cut), + Event("Maximum voltage [V]", self.param.voltage_high_cut - V), + ] ###################### # (Some) variables ###################### self.variables = { - "Particle stoichiometry": sto, - "Particle surface stoichiometry": PrimaryBroadcast(sto_surf, "particle"), + "Positive particle stoichiometry": sto_p, + "Positive particle surface stoichiometry": PrimaryBroadcast( + sto_p_surf, "positive particle" + ), "Time [s]": pybamm_t, "Current [A]": I, "Current variable [A]": I, # for compatibility with pybamm.Experiment @@ -128,29 +146,36 @@ def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): "Open-circuit voltage [V]": U, } - def U(self, sto): + def U(self, sto, domain): """ Dimensional open-circuit potential [V]. Credit: PyBaMM """ - inputs = {"Particle surface stoichiometry": sto} - out = FunctionParameter("Electrode OCP [V]", inputs) + Domain = domain.capitalize() + if domain == "negative": + inputs = {"Average positive particle stoichiometry": sto} + else: + inputs = {f"{Domain} particle surface stoichiometry": sto} + out = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) - out.print_name = r"U(c^\mathrm{surf}_\mathrm{s})" + if domain == "negative": + out.print_name = r"U_\mathrm{n}(c^\mathrm{av}_\mathrm{s,p})" + elif domain == "positive": + out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" return out def tau_d(self, sto): """ - Diffusion time scale [s] for lithium in the particles. + Dimensional solid-state diffusion time scale [s]. """ - inputs = {"Particle surface stoichiometry": sto} - return FunctionParameter("Particle diffusion time scale [s]", inputs) + inputs = {"Positive particle surface stoichiometry": sto} + return FunctionParameter("Positive particle diffusion time scale [s]", inputs) def R0(self, sto): """ Series resistance [Ohm]. """ - inputs = {"Particle surface stoichiometry": sto} + inputs = {"Positive particle surface stoichiometry": sto} return FunctionParameter("Series resistance [Ohm]", inputs) @property @@ -161,29 +186,30 @@ def default_parameter_values(self) -> ParameterValues: @property def default_quick_plot_variables(self): return [ - "Particle stoichiometry", - "Particle surface stoichiometry", + "Positive particle stoichiometry", + "Positive particle surface stoichiometry", "Current [A]", {"Open-circuit voltage [V]", "Voltage [V]"}, ] @property def default_var_pts(self): - r = SpatialVariable("r", domain=["particle"], coord_sys="spherical polar") - return {r: 20} + r_p = SpatialVariable( + "r_p", domain=["positive particle"], coord_sys="spherical polar" + ) + return {r_p: 20} @property def default_geometry(self): - r = SpatialVariable("r", domain=["particle"], coord_sys="spherical polar") - return {"particle": {r: {"min": 0, "max": 1}}} + return {"positive particle": {"r_p": {"min": 0, "max": 1}}} @property def default_submesh_types(self): - return {"particle": pybamm.Uniform1DSubMesh} + return {"positive particle": pybamm.Uniform1DSubMesh} @property def default_spatial_methods(self): - return {"particle": pybamm.FiniteVolume()} + return {"positive particle": pybamm.FiniteVolume()} @staticmethod def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterValues: @@ -208,40 +234,52 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Unpack physical parameters F = pybamm.constants.F.value T = param["Ambient temperature [K]"] - alpha = param["Positive electrode active material volume fraction"] - c_max = param["Maximum concentration in positive electrode [mol.m-3]"] - L = param["Positive electrode thickness [m]"] - R = param["Positive particle radius [m]"] - D = param["Positive particle diffusivity [m2.s-1]"] - m = param.evaluate( + alpha_p = param["Positive electrode active material volume fraction"] + c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] + L_p = param["Positive electrode thickness [m]"] + R_p = param["Positive particle radius [m]"] + D_p = param["Positive particle diffusivity [m2.s-1]"] + m_p = param.evaluate( param["Positive electrode exchange-current density [A.m-2]"](1, 1, 2, T) ) # (A/m2)(m3/mol)**1.5 - sto_init = ( - param["Initial concentration in positive electrode [mol.m-3]"] / c_max - ) - ocp = param["Positive electrode OCP [V]"] # Compute the cell area A = param["Electrode height [m]"] * param["Electrode width [m]"] + # Compute the stoichiometry limits and initial SOC + d = get_min_max_stoichiometries(param) + y_0, y_100 = d["x_0"], d["x_100"] + sto_p_init = ( + param["Initial concentration in positive electrode [mol.m-3]"] / c_max_p + ) + soc_init = (sto_p_init - y_0) / (y_100 - y_0) + + # Compute the capacity within the stoichiometry limits + Q_th_p = F * alpha_p * c_max_p * L_p * A + Q_meas = (y_0 - y_100) * Q_th_p + # Grouped parameters - Q_th = F * alpha * c_max * L * A - tau_d = R**2 / D + tau_d_p = R_p**2 / D_p # Estimate the series resistance, neglecting conductivity losses RT_F = pybamm.constants.R.value * param["Ambient temperature [K]"] / F ce0 = param["Initial concentration in electrolyte [mol.m-3]"] - tau_ct = F * R / (m * np.sqrt(ce0)) - Rct_typ = (2 * RT_F * tau_ct) / (3 * Q_th) + tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) + Rct_typ = (2 * RT_F * tau_ct_p) / (3 * Q_th_p) R0 = Rct_typ + param["Contact resistance [Ohm]"] parameter_dictionary = { "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], "Current function [A]": param["Current function [A]"], - "Initial stoichiometry": sto_init, - "Electrode OCP [V]": ocp, - "Theoretical electrode capacity [A.s]": Q_th, - "Particle diffusion time scale [s]": tau_d, + "Initial SoC": soc_init, + "Minimum positive stoichiometry": y_100, + "Maximum positive stoichiometry": y_0, + "Lower voltage cut-off [V]": param["Lower voltage cut-off [V]"], + "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], + "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], + "Negative electrode OCP [V]": 0.0, + "Measured cell capacity [A.s]": Q_meas, + "Positive particle diffusion time scale [s]": tau_d_p, "Series resistance [Ohm]": R0, } parameter_values = ParameterValues(values=parameter_dictionary) @@ -290,25 +328,49 @@ def set_initial_state( if isinstance(initial_value, str) and initial_value.endswith("V"): V_init = float(initial_value[:-1]) + V_min = parameter_values.evaluate( + Parameter("Lower voltage cut-off [V]"), inputs=inputs + ) + V_max = parameter_values.evaluate( + Parameter("Upper voltage cut-off [V]"), inputs=inputs + ) + + if not V_min - tol <= V_init <= V_max + tol: + raise ValueError( + f"Initial voltage {V_init}V is outside the voltage limits ({V_min}, {V_max})." + ) - def ocv_function(sto): - U = FunctionParameter( - "Electrode OCP [V]", {"Particle stoichiometry": sto} + y_100 = parameter_values.evaluate( + Parameter("Minimum positive stoichiometry"), inputs=inputs + ) + y_0 = parameter_values.evaluate( + Parameter("Maximum positive stoichiometry"), inputs=inputs + ) + + def ocv_function(soc): + sto_p = y_0 - soc * (y_0 - y_100) + U_p = FunctionParameter( + "Positive electrode OCP [V]", + {"Positive particle stoichiometry": sto_p}, + ) + U_n = FunctionParameter( + "Negative electrode OCP [V]", + {"Positive particle stoichiometry": sto_p}, ) - return parameter_values.evaluate(U, inputs=inputs).squeeze() + return parameter_values.evaluate(U_p - U_n, inputs=inputs).squeeze() inverse_ocv = InverseOCV(ocv_function) - sto = inverse_ocv(V_init) + soc = inverse_ocv(V_init) elif isinstance(initial_value, int | float): - sto = initial_value + soc = initial_value else: raise ValueError("Initial value must be a float or a string ending in 'V'.") - if not 0 <= sto <= 1: - raise ValueError("Initial stoichiometry should be between 0 and 1.") + if not 0 <= soc <= 1: + raise ValueError("Initial SOC should be between 0 and 1.") - parameter_values["Initial stoichiometry"] = sto + parameter_values["Initial SoC"] = soc return parameter_values diff --git a/pybop/models/lithium_ion/weppner_huggins.py b/pybop/models/lithium_ion/weppner_huggins.py index 37b752a42..b57e206f8 100644 --- a/pybop/models/lithium_ion/weppner_huggins.py +++ b/pybop/models/lithium_ion/weppner_huggins.py @@ -48,12 +48,12 @@ def __init__(self, name="Weppner & Huggins model", **model_kwargs): # Parameters are purely symbolic at this stage, and will be set by the # `ParameterValues` class when the model is processed. - Q_th = Parameter("Theoretical electrode capacity [A.s]") + Q_th_p = Parameter("Theoretical electrode capacity [A.s]") U = Parameter("Reference voltage [V]") U_prime = Parameter("Derivative of the OCP wrt stoichiometry [V]") - tau_d = Parameter("Particle diffusion time scale [s]") + tau_d_p = Parameter("Positive particle diffusion time scale [s]") ###################### # Input current (positive on discharge) @@ -64,9 +64,9 @@ def __init__(self, name="Weppner & Huggins model", **model_kwargs): # Governing equations ###################### # Surface stoichiometry - sto_surf = 2 * I / (3 * Q_th) * (pybamm_t * tau_d / np.pi) ** 0.5 + sto_p_surf = 2 * I / (3 * Q_th_p) * (pybamm_t * tau_d_p / np.pi) ** 0.5 # Linearised voltage - V = U + U_prime * sto_surf + V = U + U_prime * sto_p_surf ###################### # (Some) variables @@ -136,25 +136,25 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Unpack physical parameters F = pybamm.constants.F.value - alpha = param["Positive electrode active material volume fraction"] - c_max = param["Maximum concentration in positive electrode [mol.m-3]"] - L = param["Positive electrode thickness [m]"] - R = param["Positive particle radius [m]"] - D = param["Positive particle diffusivity [m2.s-1]"] + alpha_p = param["Positive electrode active material volume fraction"] + c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] + L_p = param["Positive electrode thickness [m]"] + R_p = param["Positive particle radius [m]"] + D_p = param["Positive particle diffusivity [m2.s-1]"] # Compute the cell area A = param["Electrode height [m]"] * param["Electrode width [m]"] # Grouped parameters - Q_th = F * alpha * c_max * L * A - tau_d = R**2 / D + Q_th_p = F * alpha_p * c_max_p * L_p * A + tau_d_p = R_p**2 / D_p parameter_dictionary = { "Current function [A]": param["Current function [A]"], "Reference voltage [V]": 4, "Derivative of the OCP wrt stoichiometry [V]": -1, - "Theoretical electrode capacity [A.s]": Q_th, - "Particle diffusion time scale [s]": tau_d, + "Theoretical electrode capacity [A.s]": Q_th_p, + "Positive particle diffusion time scale [s]": tau_d_p, } parameter_values = ParameterValues(values=parameter_dictionary) parameter_values._set_initial_state = WeppnerHuggins.set_initial_state # noqa: SLF001 diff --git a/tests/integration/models/test_gitt_models.py b/tests/integration/models/test_gitt_models.py index 66f5d5f96..794cf9286 100644 --- a/tests/integration/models/test_gitt_models.py +++ b/tests/integration/models/test_gitt_models.py @@ -15,7 +15,7 @@ # Parameter configurations DIFFUSION_PARAMS = [ ("Theoretical electrode capacity [A.s]", 10), - ("Particle diffusion time scale [s]", 2000), + ("Positive particle diffusion time scale [s]", 2000), ] diff --git a/tests/integration/test_applications.py b/tests/integration/test_applications.py index 248305a41..e80fa320a 100644 --- a/tests/integration/test_applications.py +++ b/tests/integration/test_applications.py @@ -51,9 +51,10 @@ def test_interpolant(self, parameter_values, discharge_dataset): parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( parameter_values ) - parameter_values["Electrode OCP [V]"] = pybop.Interpolant( + parameter_values["Positive electrode OCP [V]"] = pybop.Interpolant( discharge_dataset["Stoichiometry"], discharge_dataset["Voltage [V]"] ) + parameter_values.set_initial_state(0.9) model = pybop.lithium_ion.SPDiffusion(build=True) t_eval = np.linspace(0, 10, 100) solution = pybamm.Simulation(model, parameter_values=parameter_values).solve( @@ -189,13 +190,13 @@ def test_gitt_pulse_fit( parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( half_cell_parameter_values ) - diffusion_time = parameter_values["Particle diffusion time scale [s]"] + diffusion_time = parameter_values["Positive particle diffusion time scale [s]"] gitt_fit = pybop.GITTPulseFit(parameter_values=parameter_values) gitt_result = gitt_fit(gitt_pulse=pulse_data) np.testing.assert_allclose( - gitt_result.best_inputs["Particle diffusion time scale [s]"], + gitt_result.best_inputs["Positive particle diffusion time scale [s]"], diffusion_time, rtol=5e-2, ) @@ -204,7 +205,7 @@ def test_gitt_fit(self, half_cell_model, half_cell_parameter_values, pulse_data) parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( half_cell_parameter_values ) - diffusion_time = parameter_values["Particle diffusion time scale [s]"] + diffusion_time = parameter_values["Positive particle diffusion time scale [s]"] with pytest.raises( ValueError, match="The initial current in the pulse dataset must be zero." @@ -224,7 +225,7 @@ def test_gitt_fit(self, half_cell_model, half_cell_parameter_values, pulse_data) gitt_parameter_data = gitt_fit() np.testing.assert_allclose( - gitt_parameter_data["Particle diffusion time scale [s]"], + gitt_parameter_data["Positive particle diffusion time scale [s]"], np.asarray([diffusion_time]), rtol=5e-2, ) diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index 423b4f136..4101adff8 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -83,17 +83,12 @@ def test_set_initial_state(self, model): param.set_initial_state(0.5) else: - if isinstance(model, pybop.lithium_ion.SPDiffusion): - initial_state = "Initial stoichiometry" - else: - initial_state = "Initial SoC" - param = model.default_parameter_values param.set_initial_state(0.5) - assert param[initial_state] == 0.5 + assert param["Initial SoC"] == 0.5 - param.set_initial_state("2.8 V") - assert 0 <= param[initial_state] <= 1 + param.set_initial_state("3.8 V") + assert 0 <= param["Initial SoC"] <= 1 with pytest.raises( ValueError, From a7cd2cd27b50f6280ec4a313d3ec91263eafe6c7 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 6 Aug 2026 20:49:08 +0100 Subject: [PATCH 24/37] Fix OCP averaging --- pybop/models/lithium_ion/grouped_dfn.py | 8 ++++---- pybop/models/lithium_ion/grouped_spm.py | 12 ++++++++++-- pybop/models/lithium_ion/grouped_spme.py | 8 ++++---- 3 files changed, 18 insertions(+), 10 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 6a7d8469c..078090b92 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -347,13 +347,13 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): ###################### # Include the following variables to enable plotting via PyBaMM's plot_voltage_components ocp_n_bulk = self.U( - pybamm.x_average(zeta_n * pybamm.r_average(sto_n)) - / pybamm.x_average(zeta_n), + pybamm.x_average(Q_th_n * pybamm.r_average(sto_n)) + / pybamm.x_average(Q_th_n), "negative", ) ocp_p_bulk = self.U( - pybamm.x_average(zeta_p * pybamm.r_average(sto_p)) - / pybamm.x_average(zeta_p), + pybamm.x_average(Q_th_p * pybamm.r_average(sto_p)) + / pybamm.x_average(Q_th_p), "positive", ) voltage_components = { diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 9acfad3bf..689d97ee0 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -236,8 +236,16 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): # Voltage components ###################### # Include the following variables to enable plotting via PyBaMM's plot_voltage_components - ocp_n_bulk = self.U(pybamm.x_average(pybamm.r_average(sto_n)), "negative") - ocp_p_bulk = self.U(pybamm.x_average(pybamm.r_average(sto_p)), "positive") + ocp_n_bulk = self.U( + pybamm.x_average(Q_th_n * pybamm.r_average(sto_n)) + / pybamm.x_average(Q_th_n), + "negative", + ) + ocp_p_bulk = self.U( + pybamm.x_average(Q_th_p * pybamm.r_average(sto_p)) + / pybamm.x_average(Q_th_p), + "positive", + ) voltage_components = { "Battery voltage [V]": V, "Battery open-circuit voltage [V]": ocp_p_bulk - ocp_n_bulk, diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 67f63ae11..a2c757400 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -325,13 +325,13 @@ def __init__( ###################### # Include the following variables to enable plotting via PyBaMM's plot_voltage_components ocp_n_bulk = self.U( - pybamm.x_average(zeta_n * pybamm.r_average(sto_n)) - / pybamm.x_average(zeta_n), + pybamm.x_average(Q_th_n * pybamm.r_average(sto_n)) + / pybamm.x_average(Q_th_n), "negative", ) ocp_p_bulk = self.U( - pybamm.x_average(zeta_p * pybamm.r_average(sto_p)) - / pybamm.x_average(zeta_p), + pybamm.x_average(Q_th_p * pybamm.r_average(sto_p)) + / pybamm.x_average(Q_th_p), "positive", ) voltage_components = { From f93e55445afeca96931071a4253af870fec741e5 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 6 Aug 2026 21:01:06 +0100 Subject: [PATCH 25/37] Allow diffusivity function of temperature --- pybop/models/lithium_ion/grouped_dfn.py | 18 +++++++++++------- pybop/models/lithium_ion/grouped_spme.py | 18 +++++++++++------- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 078090b92..476ecfc56 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -259,17 +259,21 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): ###################### # The div and grad operators will be converted to the appropriate matrix # multiplication at the discretisation stage - self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / self.tau_d(sto_n, "negative")) - self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / self.tau_d(sto_p, "positive")) + self.rhs[sto_n] = pybamm.div( + pybamm.grad(sto_n) / self.tau_d(sto_n, T, "negative") + ) + self.rhs[sto_p] = pybamm.div( + pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") + ) # Boundary conditions must be provided for equations with spatial derivatives self.boundary_conditions[sto_n] = { "left": (Scalar(0), "Neumann"), - "right": (-self.tau_d(sto_n_surf, "negative") * j_n, "Neumann"), + "right": (-self.tau_d(sto_n_surf, T, "negative") * j_n, "Neumann"), } self.boundary_conditions[sto_p] = { "left": (Scalar(0), "Neumann"), - "right": (-self.tau_d(sto_p_surf, "positive") * j_p, "Neumann"), + "right": (-self.tau_d(sto_p_surf, T, "positive") * j_p, "Neumann"), } self.initial_conditions[sto_n] = sto_n_init @@ -439,12 +443,12 @@ def U(self, sto, domain): out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" return out - def tau_d(self, sto, domain): + def tau_d(self, sto, T, domain): """ - Dimensional diffusion time scale [s]. + Dimensional solid-state diffusion time scale [s]. """ Domain = domain.capitalize() - inputs = {f"{Domain} particle surface stoichiometry": sto} + inputs = {f"{Domain} particle surface stoichiometry": sto, "Temperature [K]": T} return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) def j(self, sto_surf, sto_e, eta_RT_F, domain): diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index a2c757400..11f624d9e 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -247,21 +247,25 @@ def __init__( ###################### # The div and grad operators will be converted to the appropriate matrix # multiplication at the discretisation stage - self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / self.tau_d(sto_n, "negative")) - self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / self.tau_d(sto_p, "positive")) + self.rhs[sto_n] = pybamm.div( + pybamm.grad(sto_n) / self.tau_d(sto_n, T, "negative") + ) + self.rhs[sto_p] = pybamm.div( + pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") + ) # Boundary conditions must be provided for equations with spatial derivatives self.boundary_conditions[sto_n] = { "left": (Scalar(0), "Neumann"), "right": ( - -self.tau_d(sto_n_surf, "negative") * pybamm.x_average(j_n), + -self.tau_d(sto_n_surf, T, "negative") * pybamm.x_average(j_n), "Neumann", ), } self.boundary_conditions[sto_p] = { "left": (Scalar(0), "Neumann"), "right": ( - -self.tau_d(sto_p_surf, "positive") * pybamm.x_average(j_p), + -self.tau_d(sto_p_surf, T, "positive") * pybamm.x_average(j_p), "Neumann", ), } @@ -432,12 +436,12 @@ def U(self, sto, domain): out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" return out - def tau_d(self, sto, domain): + def tau_d(self, sto, T, domain): """ - Dimensional diffusion time scale [s]. + Dimensional solid-state diffusion time scale [s]. """ Domain = domain.capitalize() - inputs = {f"{Domain} particle surface stoichiometry": sto} + inputs = {f"{Domain} particle surface stoichiometry": sto, "Temperature [K]": T} return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) def j(self, sto_surf, sto_e, eta_RT_F, domain): From 11961765c18554fbd1aec587b78d4a5242adce3a Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 6 Aug 2026 21:12:32 +0100 Subject: [PATCH 26/37] Create alternative_functions.py --- .../lithium_ion/alternative_functions.py | 39 +++++++++++++++ pybop/models/lithium_ion/grouped_dfn.py | 44 ++++++----------- pybop/models/lithium_ion/grouped_spm.py | 49 +++++++------------ pybop/models/lithium_ion/grouped_spme.py | 44 ++++++----------- 4 files changed, 86 insertions(+), 90 deletions(-) create mode 100644 pybop/models/lithium_ion/alternative_functions.py diff --git a/pybop/models/lithium_ion/alternative_functions.py b/pybop/models/lithium_ion/alternative_functions.py new file mode 100644 index 000000000..e565b22fb --- /dev/null +++ b/pybop/models/lithium_ion/alternative_functions.py @@ -0,0 +1,39 @@ +import pybamm + +""" Alternative functions written as classes to allow pickling. """ + + +class FunctionalDiffusionTime: + def __init__(self, D_prefactor, D, c_scale): + self.D_prefactor = D_prefactor + self.D = D + self.c_scale = c_scale + + def __call__(self, sto, T): + D_prefactor, D, c_scale = self.D_prefactor, self.D, self.c_scale + return D_prefactor / D(sto * c_scale, T) + + +class AsymmetricButlerVolmer: + def __init__(self, alpha): + self.alpha = alpha # cathodic transfer coefficient + + def __call__(self, sto_surf, sto_e, eta_RT_F): + alpha = self.alpha + j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) + + +class MultiphaseButlerVolmer: + def __init__(self, alpha, omega): + self.alpha = alpha + self.omega = omega + + def __call__(self, sto_surf, sto_e, eta_RT_F): + alpha, omega = self.alpha, self.omega + j0 = ( + sto_surf ** (alpha * omega) + * (1 - sto_surf) ** ((1 - alpha) * omega) + * sto_e ** (1 - alpha) + ) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 476ecfc56..2869262b3 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -14,6 +14,11 @@ get_min_max_stoichiometries, ) +from pybop.models.lithium_ion.alternative_functions import ( + AsymmetricButlerVolmer, + FunctionalDiffusionTime, + MultiphaseButlerVolmer, +) from pybop.models.lithium_ion.base_model import BaseGroupedModel @@ -661,8 +666,15 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal zeta_p = epsilon_p / epsilon_sep zeta_n = epsilon_n / epsilon_sep - tau_d_p = R_p**2 / D_p - tau_d_n = R_n**2 / D_n + try: + tau_d_p = R_p**2 / D_p + except TypeError: + tau_d_p = FunctionalDiffusionTime(R_p**2, D_p, c_max_p) + + try: + tau_d_n = R_n**2 / D_n + except TypeError: + tau_d_n = FunctionalDiffusionTime(R_n**2, D_n, c_max_n) tau_e = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) beta_p = epsilon_p**b_p / epsilon_sep**b_sep @@ -746,31 +758,3 @@ def get_multiphase_butler_volmer(domain: str): omega = pybamm.Parameter(f"{Domain} electrode charge transfer ideality factor") return MultiphaseButlerVolmer(alpha, omega) - - -""" Alternative, dimensionless exchange rate functions written as classes to allow pickling. """ - - -class AsymmetricButlerVolmer: - def __init__(self, alpha): - self.alpha = alpha # cathodic transfer coefficient - - def __call__(self, sto_surf, sto_e, eta_RT_F): - alpha = self.alpha - j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) - return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) - - -class MultiphaseButlerVolmer: - def __init__(self, alpha, omega): - self.alpha = alpha - self.omega = omega - - def __call__(self, sto_surf, sto_e, eta_RT_F): - alpha, omega = self.alpha, self.omega - j0 = ( - sto_surf ** (alpha * omega) - * (1 - sto_surf) ** ((1 - alpha) * omega) - * sto_e ** (1 - alpha) - ) - return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 689d97ee0..f0980e180 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -14,6 +14,11 @@ get_min_max_stoichiometries, ) +from pybop.models.lithium_ion.alternative_functions import ( + AsymmetricButlerVolmer, + FunctionalDiffusionTime, + MultiphaseButlerVolmer, +) from pybop.models.lithium_ion.base_model import BaseGroupedModel @@ -166,8 +171,8 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): eta_p = PrimaryBroadcast(v_s_p - U_p, "positive electrode") # Exchange rates - j_n = self.j(sto_n_surf, eta_n / RT_F, "negative") / tau_ct_n - j_p = self.j(sto_p_surf, eta_p / RT_F, "positive") / tau_ct_p + j_n = self.j(sto_n_surf, 1.0, eta_n / RT_F, "negative") / tau_ct_n + j_p = self.j(sto_p_surf, 1.0, eta_p / RT_F, "positive") / tau_ct_p ###################### # Double layer @@ -340,13 +345,14 @@ def tau_d(self, sto, T, domain): inputs = {f"{Domain} particle surface stoichiometry": sto, "Temperature [K]": T} return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) - def j(self, sto_surf, eta_RT_F, domain): + def j(self, sto_surf, sto_e, eta_RT_F, domain): """ Dimensionless exchange rate. """ Domain = domain.capitalize() inputs = { f"{Domain} particle surface stoichiometry": sto_surf, + f"{Domain} electrode electrolyte stoichiometry": sto_e, f"{Domain} electrode dimensionless overpotential": eta_RT_F, } return FunctionParameter( @@ -539,8 +545,15 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Grouped parameters Q_meas = (Q_meas_n + Q_meas_p) / 2 - tau_d_p = R_p**2 / D_p - tau_d_n = R_n**2 / D_n + try: + tau_d_p = R_p**2 / D_p + except TypeError: + tau_d_p = FunctionalDiffusionTime(R_p**2, D_p, c_max_p) + + try: + tau_d_n = R_n**2 / D_n + except TypeError: + tau_d_n = FunctionalDiffusionTime(R_n**2, D_n, c_max_n) tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) tau_ct_n = F * R_n / (m_n * np.sqrt(ce0)) @@ -583,7 +596,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal return parameter_values @staticmethod - def symmetric_butler_volmer(sto_surf, eta_RT_F): + def symmetric_butler_volmer(sto_surf, sto_e, eta_RT_F): """ Dimensionless Butler-Volmer exchange rate. """ @@ -612,27 +625,3 @@ def get_multiphase_butler_volmer(domain: str): omega = pybamm.Parameter(f"{Domain} electrode charge transfer ideality factor") return MultiphaseButlerVolmer(alpha, omega) - - -""" Alternative, dimensionless exchange rate functions written as classes to allow pickling. """ - - -class AsymmetricButlerVolmer: - def __init__(self, alpha): - self.alpha = alpha # cathodic transfer coefficient - - def __call__(self, sto_surf, eta_RT_F): - alpha = self.alpha - j0 = sto_surf**alpha * (1 - sto_surf) ** (1 - alpha) - return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) - - -class MultiphaseButlerVolmer: - def __init__(self, alpha, omega): - self.alpha = alpha - self.omega = omega - - def __call__(self, sto_surf, eta_RT_F): - alpha, omega = self.alpha, self.omega - j0 = sto_surf ** (alpha * omega) * (1 - sto_surf) ** ((1 - alpha) * omega) - return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 11f624d9e..09e37ffb3 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -15,6 +15,11 @@ get_min_max_stoichiometries, ) +from pybop.models.lithium_ion.alternative_functions import ( + AsymmetricButlerVolmer, + FunctionalDiffusionTime, + MultiphaseButlerVolmer, +) from pybop.models.lithium_ion.base_model import BaseGroupedModel @@ -655,8 +660,15 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal zeta_p = epsilon_p / epsilon_sep zeta_n = epsilon_n / epsilon_sep - tau_d_p = R_p**2 / D_p - tau_d_n = R_n**2 / D_n + try: + tau_d_p = R_p**2 / D_p + except TypeError: + tau_d_p = FunctionalDiffusionTime(R_p**2, D_p, c_max_p) + + try: + tau_d_n = R_n**2 / D_n + except TypeError: + tau_d_n = FunctionalDiffusionTime(R_n**2, D_n, c_max_n) tau_e = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) beta_p = epsilon_p**b_p / epsilon_sep**b_sep @@ -739,31 +751,3 @@ def get_multiphase_butler_volmer(domain: str): omega = pybamm.Parameter(f"{Domain} electrode charge transfer ideality factor") return MultiphaseButlerVolmer(alpha, omega) - - -""" Alternative, dimensionless exchange rate functions written as classes to allow pickling. """ - - -class AsymmetricButlerVolmer: - def __init__(self, alpha): - self.alpha = alpha # cathodic transfer coefficient - - def __call__(self, sto_surf, sto_e, eta_RT_F): - alpha = self.alpha - j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) - return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) - - -class MultiphaseButlerVolmer: - def __init__(self, alpha, omega): - self.alpha = alpha - self.omega = omega - - def __call__(self, sto_surf, sto_e, eta_RT_F): - alpha, omega = self.alpha, self.omega - j0 = ( - sto_surf ** (alpha * omega) - * (1 - sto_surf) ** ((1 - alpha) * omega) - * sto_e ** (1 - alpha) - ) - return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) From 1a2678dfda2bbce5163aff1781fa5e792f3b3f80 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 6 Aug 2026 21:26:25 +0100 Subject: [PATCH 27/37] Delay DFN --- pybop/models/lithium_ion/__init__.py | 1 - pybop/models/lithium_ion/grouped_dfn.py | 760 ------------------ .../integration/models/test_grouped_models.py | 2 - tests/unit/test_models.py | 8 +- 4 files changed, 1 insertion(+), 770 deletions(-) delete mode 100644 pybop/models/lithium_ion/grouped_dfn.py diff --git a/pybop/models/lithium_ion/__init__.py b/pybop/models/lithium_ion/__init__.py index ff71b962d..7dcbbc032 100644 --- a/pybop/models/lithium_ion/__init__.py +++ b/pybop/models/lithium_ion/__init__.py @@ -4,6 +4,5 @@ from .sp_diffusion import SPDiffusion from .grouped_spm import GroupedSPM from .grouped_spme import GroupedSPMe -from .grouped_dfn import GroupedDFN from .weppner_huggins import WeppnerHuggins from .cell_temperature import CellTemperature diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py deleted file mode 100644 index 2869262b3..000000000 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ /dev/null @@ -1,760 +0,0 @@ -import numpy as np -import pybamm -from pybamm import ( - Event, - FunctionParameter, - Parameter, - ParameterValues, - PrimaryBroadcast, - Scalar, - Variable, -) -from pybamm import t as pybamm_t -from pybamm.models.full_battery_models.lithium_ion.electrode_soh import ( - get_min_max_stoichiometries, -) - -from pybop.models.lithium_ion.alternative_functions import ( - AsymmetricButlerVolmer, - FunctionalDiffusionTime, - MultiphaseButlerVolmer, -) -from pybop.models.lithium_ion.base_model import BaseGroupedModel - - -class GroupedDFN(BaseGroupedModel): - """ - A grouped parameter version of the Doyle Fuller Newman (DFN) model. - - Parameters - ---------- - name : str, optional - The name of the model. - **model_kwargs : optional - Valid PyBaMM model option keys and their values, for example: - options : dict, optional - A dictionary of options to customise the behaviour of the PyBaMM model. - build : bool, optional - If True, the model is built upon creation (default: False). - """ - - def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): - super().__init__(name=name, **model_kwargs) - - # Unpack model options - include_double_layer = self.options["surface form"] == "differential" - - pybamm.citations.register( - """ - @article{Hallemans2025, - title = {{Physics-Based Battery Model Parametrisation from Impedance Data}}, - author = {Hallemans, Noël and Courtier, Nicola E. and Please, Colin P. and Planden, Brady and Dhoot, Rishit and Timms, Robert and Chapman, S. Jon and Howey, David and Duncan, Stephen R.}, - journal = {Journal of the Electrochemical Society}, - volume = {172}, - number = {6}, - pages = {060507}, - year = {2025}, - publisher = {The Electrochemical Society}, - doi = {10.1149/1945-7111/add41b} - } - """ - ) # Note that the electrode electrolyte timescales have been replaced by relative transport efficiencies - # and there is an additional variable (i_e) and associated grouped parameter (gamma_e) compared to the SPMe - - ###################### - # Variables - ###################### - # Variables that depend on time only are created without a domain - Q = Variable("Discharge capacity [A.h]") - Qt = Variable("Throughput capacity [A.h]") - - # Variables that vary spatially are created with a domain - v_s_n = Variable( - "Negative particle surface voltage [V]", - domain="negative electrode", - ) - v_s_p = Variable( - "Positive particle surface voltage [V]", - domain="positive electrode", - ) - - sto_n = Variable( - "Negative particle stoichiometry", - domain="negative particle", - auxiliary_domains={"secondary": "negative electrode"}, - ) - sto_p = Variable( - "Positive particle stoichiometry", - domain="positive particle", - auxiliary_domains={"secondary": "positive electrode"}, - ) - sto_e_n = Variable( - "Negative electrode electrolyte stoichiometry", - domain="negative electrode", - ) - sto_e_sep = Variable( - "Separator electrolyte stoichiometry", - domain="separator", - ) - sto_e_p = Variable( - "Positive electrode electrolyte stoichiometry", - domain="positive electrode", - ) - - # 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 - sto_n_surf = pybamm.surf(sto_n) - sto_p_surf = pybamm.surf(sto_p) - - # Events specify points at which a solution should terminate - tol = pybamm.settings.tolerances["U__c_s"] - self.events += [ - Event( - "Minimum negative particle surface stoichiometry", - pybamm.min(sto_n_surf) - tol, - ), - Event( - "Maximum negative particle surface stoichiometry", - (1 - tol) - pybamm.max(sto_n_surf), - ), - Event( - "Minimum positive particle surface stoichiometry", - pybamm.min(sto_p_surf) - tol, - ), - Event( - "Maximum positive particle surface stoichiometry", - (1 - tol) - pybamm.max(sto_p_surf), - ), - ] - - ###################### - # Parameters - ###################### - # Parameters are purely symbolic at this stage, and will be set by the - # `ParameterValues` class when the model is processed. - - F = self.param.F # Faraday constant - Rg = self.param.R # Universal gas constant - T = self.param.T_init # Temperature - RT_F = Rg * T / F # Thermal voltage - - soc_init = Parameter("Initial SoC") - x_0 = Parameter("Minimum negative stoichiometry") - x_100 = Parameter("Maximum negative stoichiometry") - y_100 = Parameter("Minimum positive stoichiometry") - y_0 = Parameter("Maximum positive stoichiometry") - - # Grouped parameters - Q_th_p = Parameter("Measured cell capacity [A.s]") / (y_0 - y_100) - Q_th_n = Parameter("Measured cell capacity [A.s]") / (x_100 - x_0) - Q_e = Parameter("Reference electrolyte capacity [A.s]") - - tau_ct_p = Parameter("Positive electrode charge transfer time scale [s]") - tau_ct_n = Parameter("Negative electrode charge transfer time scale [s]") - - l_p = Parameter("Positive electrode relative thickness") - l_n = Parameter("Negative electrode relative thickness") - - t_plus = Parameter("Cation transference number") - - R0 = Parameter("Series resistance [Ohm]") - - zeta_n = Parameter("Negative electrode relative porosity") - zeta_p = Parameter("Positive electrode relative porosity") - - beta_n = Parameter("Negative electrode relative transport efficiency") - beta_p = Parameter("Positive electrode relative transport efficiency") - - tau_e = Parameter("Electrolyte diffusion time scale [s]") - gamma_e = Parameter("Reference electrolyte scaled conductivity [V-1.s-1]") - - ###################### - # Input current (positive on discharge) - ###################### - I = self.param.current_with_time - - ###################### - # State of Charge - ###################### - # The `rhs` dictionary contains differential equations, with the key being the - # variable in the d/dt - self.rhs[Q] = I / 3600 - self.rhs[Qt] = abs(I) / 3600 - # Initial conditions must be provided for the ODEs - self.initial_conditions[Q] = Scalar(0) - self.initial_conditions[Qt] = Scalar(0) - - ###################### - # Potentials - ###################### - U_n = self.U(sto_n_surf, "negative") - U_p = self.U(sto_p_surf, "positive") - - sto_n_init = x_0 + (x_100 - x_0) * soc_init - sto_p_init = y_0 + (y_100 - y_0) * soc_init - U_n_init = self.U(sto_n_init, "negative") - U_p_init = self.U(sto_p_init, "positive") - - ###################### - # Exchange current - ###################### - # Primary broadcasts are used to broadcast scalar quantities across a domain - # into a vector of the right shape, for multiplying with other vectors - - # Overpotentials - eta_n = v_s_n - U_n - eta_p = v_s_p - U_p - - # Exchange rates - j_n = self.j(sto_n_surf, sto_e_n, eta_n / RT_F, "negative") / tau_ct_n - j_p = self.j(sto_p_surf, sto_e_p, eta_p / RT_F, "positive") / tau_ct_p - - # Electrolyte current - i_e_n = (beta_n * gamma_e) * ( - pybamm.grad(v_s_n) - + (2 * RT_F * (1 - t_plus)) * pybamm.grad(sto_e_n) / sto_e_n - ) - i_e_p = (beta_p * gamma_e) * ( - pybamm.grad(v_s_p) - + (2 * RT_F * (1 - t_plus)) * pybamm.grad(sto_e_p) / sto_e_p - ) - - ###################### - # Double layer - ###################### - if include_double_layer: - # Additional parameters - C_p = Parameter("Positive electrode capacitance [F]") - C_n = Parameter("Negative electrode capacitance [F]") - - # Electrode surface potentials - self.rhs[v_s_n] = (l_n * Q_e * pybamm.div(i_e_n) - 3 * Q_th_n * j_n) / C_n - self.rhs[v_s_p] = (l_p * Q_e * pybamm.div(i_e_p) - 3 * Q_th_p * j_p) / C_p - else: - self.algebraic[v_s_n] = l_n * Q_e * pybamm.div(i_e_n) - 3 * Q_th_n * j_n - self.algebraic[v_s_p] = l_p * Q_e * pybamm.div(i_e_p) - 3 * Q_th_p * j_p - - self.initial_conditions[v_s_n] = U_n_init - self.initial_conditions[v_s_p] = U_p_init - - self.boundary_conditions[v_s_n] = { - "left": (Scalar(0), "Neumann"), - "right": ( - I / (beta_n * gamma_e * Q_e) - - (2 * RT_F * (1 - t_plus)) - * (pybamm.boundary_gradient(sto_e_sep, "left") / beta_n) - / pybamm.boundary_value(sto_e_sep, "left"), - "Neumann", - ), - } - self.boundary_conditions[v_s_p] = { - "left": ( - I / (beta_p * gamma_e * Q_e) - - (2 * RT_F * (1 - t_plus)) - * (pybamm.boundary_gradient(sto_e_sep, "right") / beta_p) - / pybamm.boundary_value(sto_e_sep, "right"), - "Neumann", - ), - "right": (Scalar(0), "Neumann"), - } - - ###################### - # Particles - ###################### - # The div and grad operators will be converted to the appropriate matrix - # multiplication at the discretisation stage - self.rhs[sto_n] = pybamm.div( - pybamm.grad(sto_n) / self.tau_d(sto_n, T, "negative") - ) - self.rhs[sto_p] = pybamm.div( - pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") - ) - - # Boundary conditions must be provided for equations with spatial derivatives - self.boundary_conditions[sto_n] = { - "left": (Scalar(0), "Neumann"), - "right": (-self.tau_d(sto_n_surf, T, "negative") * j_n, "Neumann"), - } - self.boundary_conditions[sto_p] = { - "left": (Scalar(0), "Neumann"), - "right": (-self.tau_d(sto_p_surf, T, "positive") * j_p, "Neumann"), - } - - self.initial_conditions[sto_n] = sto_n_init - self.initial_conditions[sto_p] = sto_p_init - - ###################### - # Electrolyte - ###################### - self.rhs[sto_e_n] = ( - pybamm.div(pybamm.grad(sto_e_n) * beta_n / tau_e + (1 - t_plus) * i_e_n) - ) / zeta_n - self.rhs[sto_e_sep] = pybamm.div( - pybamm.grad(sto_e_sep) / tau_e - t_plus * I / Q_e - ) - self.rhs[sto_e_p] = ( - pybamm.div(pybamm.grad(sto_e_p) * beta_p / tau_e + (1 - t_plus) * i_e_p) - ) / zeta_p - - self.boundary_conditions[sto_e_n] = { - "left": (Scalar(0), "Neumann"), - "right": (pybamm.boundary_gradient(sto_e_sep, "left") / beta_n, "Neumann"), - } - self.boundary_conditions[sto_e_sep] = { - "left": (pybamm.boundary_value(sto_e_n, "right"), "Dirichlet"), - "right": (pybamm.boundary_value(sto_e_p, "left"), "Dirichlet"), - } - self.boundary_conditions[sto_e_p] = { - "left": (pybamm.boundary_gradient(sto_e_sep, "right") / beta_p, "Neumann"), - "right": (Scalar(0), "Neumann"), - } - - self.initial_conditions[sto_e_n] = Scalar(1) - self.initial_conditions[sto_e_sep] = Scalar(1) - self.initial_conditions[sto_e_p] = Scalar(1) - - # Electrolyte overpotential - eta_e = (2 * (1 - t_plus) * RT_F) * ( - pybamm.x_average(pybamm.log(sto_e_p)) - - pybamm.x_average(pybamm.log(sto_e_n)) - ) - - # Electrolyte Ohmic losses - DPhi_e = ( - (2 * (1 - t_plus) * RT_F) - * ( - pybamm.log(pybamm.boundary_value(sto_e_p, "left")) - - pybamm.log(pybamm.boundary_value(sto_e_n, "right")) - ) - - eta_e - - (1 - l_p - l_n) * I / (gamma_e * Q_e) - - ( - pybamm.x_average(v_s_p) - - pybamm.x_average(v_s_n) - - pybamm.boundary_value(v_s_p, "left") - + pybamm.boundary_value(v_s_n, "right") - ) - ) - - ###################### - # Cell voltage - ###################### - V = pybamm.x_average(v_s_p) - pybamm.x_average(v_s_n) + eta_e + DPhi_e - R0 * I - - # Save the initial OCV - self.param.ocv_init = U_p_init - U_n_init - - # Events specify points at which a solution should terminate - self.events += [ - Event("Minimum voltage [V]", V - self.param.voltage_low_cut), - Event("Maximum voltage [V]", self.param.voltage_high_cut - V), - ] - - ###################### - # Voltage components - ###################### - # Include the following variables to enable plotting via PyBaMM's plot_voltage_components - ocp_n_bulk = self.U( - pybamm.x_average(Q_th_n * pybamm.r_average(sto_n)) - / pybamm.x_average(Q_th_n), - "negative", - ) - ocp_p_bulk = self.U( - pybamm.x_average(Q_th_p * pybamm.r_average(sto_p)) - / pybamm.x_average(Q_th_p), - "positive", - ) - voltage_components = { - "Battery voltage [V]": V, - "Battery open-circuit voltage [V]": ocp_p_bulk - ocp_n_bulk, - "Battery particle concentration overpotential [V]": ( - (pybamm.x_average(self.U(sto_p_surf, "positive")) - ocp_p_bulk) - - (pybamm.x_average(self.U(sto_n_surf, "negative")) - ocp_n_bulk) - ), - "X-averaged battery reaction overpotential [V]": pybamm.x_average(eta_p) - - pybamm.x_average(eta_n), - "X-averaged battery concentration overpotential [V]": eta_e, - "X-averaged battery electrolyte ohmic losses [V]": DPhi_e, - "X-averaged battery solid phase ohmic losses [V]": Scalar(0), - "Contact overpotential [V]": R0 * I, # includes solid phase Ohmic losses - # and split by electrode - "Negative electrode bulk open-circuit potential [V]": ocp_n_bulk, - "Positive electrode bulk open-circuit potential [V]": ocp_p_bulk, - "Negative particle concentration overpotential [V]": ( - pybamm.x_average(self.U(sto_n_surf, "negative")) - ocp_n_bulk - ), - "Positive particle concentration overpotential [V]": ( - pybamm.x_average(self.U(sto_p_surf, "positive")) - ocp_p_bulk - ), - "X-averaged negative electrode reaction overpotential [V]" - "": pybamm.x_average(eta_n), - "X-averaged positive electrode reaction overpotential [V]" - "": pybamm.x_average(eta_p), - "X-averaged battery negative solid phase ohmic losses [V]": Scalar(0), - "X-averaged battery positive solid phase ohmic losses [V]": Scalar(0), - } - - ###################### - # (Some) variables - ###################### - # The `variables` dictionary contains all variables that might be useful for - # visualising the solution of the model - self.variables = { - "Negative particle stoichiometry": sto_n, - "Negative particle surface stoichiometry": sto_n_surf, - "Negative particle surface voltage [V]": v_s_n, - "Negative electrode potential [V]": eta_n - - pybamm.boundary_value(eta_n, "left"), - "Negative electrode electrolyte stoichiometry": sto_e_n, - "Separator electrolyte stoichiometry": sto_e_sep, - "Positive electrode electrolyte stoichiometry": sto_e_p, - "Electrolyte stoichiometry": pybamm.concatenation( - sto_e_n, sto_e_sep, sto_e_p - ), - "Positive particle stoichiometry": sto_p, - "Positive particle surface stoichiometry": sto_p_surf, - "Positive particle surface voltage [V]": v_s_p, - "Positive electrode potential [V]": V - + eta_p - - pybamm.boundary_value(eta_p, "right"), - "Electrolyte scaled current density [s-1]": pybamm.concatenation( - i_e_n, PrimaryBroadcast(I / Q_e, "separator"), i_e_p - ), - "Time [s]": pybamm_t, - "Time [h]": pybamm_t / 3600, - "Current [A]": I, - "Current variable [A]": I, # for compatibility with pybamm.Experiment - "Discharge capacity [A.h]": Q, - "Throughput capacity [A.h]": Qt, - "Voltage [V]": V, - "Open-circuit voltage [V]": pybamm.boundary_value(U_p, "right") - - pybamm.boundary_value(U_n, "left"), - **voltage_components, - } - - def U(self, sto, domain): - """ - Dimensional open-circuit potential [V]. - Credit: PyBaMM - """ - Domain = domain.capitalize() - inputs = {f"{Domain} particle surface stoichiometry": sto} - out = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) - - if domain == "negative": - out.print_name = r"U_\mathrm{n}(c^\mathrm{surf}_\mathrm{s,n})" - elif domain == "positive": - out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" - return out - - def tau_d(self, sto, T, domain): - """ - Dimensional solid-state diffusion time scale [s]. - """ - Domain = domain.capitalize() - inputs = {f"{Domain} particle surface stoichiometry": sto, "Temperature [K]": T} - return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) - - def j(self, sto_surf, sto_e, eta_RT_F, domain): - """ - Dimensionless exchange rate. - """ - Domain = domain.capitalize() - inputs = { - f"{Domain} particle surface stoichiometry": sto_surf, - f"{Domain} electrode electrolyte stoichiometry": sto_e, - f"{Domain} electrode dimensionless overpotential": eta_RT_F, - } - return FunctionParameter( - f"{Domain} electrode dimensionless exchange rate", inputs - ) - - @property - def default_parameter_values(self) -> pybamm.ParameterValues: - param = ParameterValues("Chen2020") - ce0 = param["Initial concentration in electrolyte [mol.m-3]"] - T = param["Ambient temperature [K]"] - param["Electrolyte conductivity [S.m-1]"] = param[ - "Electrolyte conductivity [S.m-1]" - ](ce0, T) - param["Electrolyte diffusivity [m2.s-1]"] = param[ - "Electrolyte diffusivity [m2.s-1]" - ](ce0, T) - return self.create_grouped_parameters(param) - - @property - def default_quick_plot_variables(self): - return [ - "Negative particle surface stoichiometry", - "Electrolyte stoichiometry", - "Positive particle surface stoichiometry", - "Current [A]", - { - "Negative electrode potential [V]", - "Negative particle surface voltage [V]", - }, - "Electrolyte scaled current density [s-1]", - { - "Positive electrode potential [V]", - "Positive particle surface voltage [V]", - }, - {"Open-circuit voltage [V]", "Voltage [V]"}, - ] - - @property - def default_var_pts(self): - x_n = pybamm.SpatialVariable( - "x_n", - domain=["negative electrode"], - coord_sys="cartesian", - ) - x_s = pybamm.SpatialVariable( - "x_s", - domain=["separator"], - coord_sys="cartesian", - ) - x_p = pybamm.SpatialVariable( - "x_p", - domain=["positive electrode"], - coord_sys="cartesian", - ) - - # Add particle domains - r_n = pybamm.SpatialVariable( - "r_n", - domain=["negative particle"], - auxiliary_domains={"secondary": "negative electrode"}, - coord_sys="spherical polar", - ) - r_p = pybamm.SpatialVariable( - "r_p", - domain=["positive particle"], - auxiliary_domains={"secondary": "positive electrode"}, - coord_sys="spherical polar", - ) - - return {x_n: 20, x_s: 20, x_p: 20, r_n: 20, r_p: 20} - - @property - def default_geometry(self): - l_p = Parameter("Positive electrode relative thickness") - l_n = Parameter("Negative electrode relative thickness") - - return { - "negative electrode": {"x_n": {"min": 0, "max": l_n}}, - "separator": {"x_s": {"min": l_n, "max": 1 - l_p}}, - "positive electrode": {"x_p": {"min": 1 - l_p, "max": 1}}, - "negative particle": {"r_n": {"min": 0, "max": 1}}, - "positive particle": {"r_p": {"min": 0, "max": 1}}, - } - - @property - def default_submesh_types(self): - return { - "negative electrode": pybamm.Uniform1DSubMesh, - "separator": pybamm.Uniform1DSubMesh, - "positive electrode": pybamm.Uniform1DSubMesh, - "negative particle": pybamm.Uniform1DSubMesh, - "positive particle": pybamm.Uniform1DSubMesh, - } - - @property - def default_spatial_methods(self): - return { - "negative electrode": pybamm.FiniteVolume(), - "separator": pybamm.FiniteVolume(), - "positive electrode": pybamm.FiniteVolume(), - "negative particle": pybamm.FiniteVolume(), - "positive particle": pybamm.FiniteVolume(), - } - - @staticmethod - def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterValues: - """ - Create a parameter set for the Grouped Doyle Fuller Newman Model from a - PyBaMM lithium-ion ParameterValues object. - - Parameters - ---------- - parameter_values : pybamm.ParameterValues - Parameters and their corresponding values. - - Returns - ------- - parameter_values : pybamm.ParameterValues - A new set of parameters and their values. - """ - param = parameter_values - - # Unpack physical parameters - F = pybamm.constants.F.value - T = param["Ambient temperature [K]"] - alpha_p = param["Positive electrode active material volume fraction"] - alpha_n = param["Negative electrode active material volume fraction"] - c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] - c_max_n = param["Maximum concentration in negative electrode [mol.m-3]"] - L_p = param["Positive electrode thickness [m]"] - L_n = param["Negative electrode thickness [m]"] - epsilon_p = param["Positive electrode porosity"] - epsilon_n = param["Negative electrode porosity"] - R_p = param["Positive particle radius [m]"] - R_n = param["Negative particle radius [m]"] - D_p = param["Positive particle diffusivity [m2.s-1]"] - D_n = param["Negative particle diffusivity [m2.s-1]"] - b_p = param["Positive electrode Bruggeman coefficient (electrolyte)"] - b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] - Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] - Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] - m_p = param.evaluate( - param["Positive electrode exchange-current density [A.m-2]"](1, 1, 2, T) - ) # (A/m2)(m3/mol)**1.5 - m_n = param.evaluate( - param["Negative electrode exchange-current density [A.m-2]"](1, 1, 2, T) - ) # (A/m2)(m3/mol)**1.5 - sigma_p = ( - param["Positive electrode conductivity [S.m-1]"] - * alpha_p ** param["Positive electrode Bruggeman coefficient (electrode)"] - ) - sigma_n = ( - param["Negative electrode conductivity [S.m-1]"] - * alpha_n ** param["Negative electrode Bruggeman coefficient (electrode)"] - ) - - # Separator and electrolyte properties - ce0 = param["Initial concentration in electrolyte [mol.m-3]"] - De = param["Electrolyte diffusivity [m2.s-1]"] # (ce0, T) - L_s = param["Separator thickness [m]"] - epsilon_sep = param["Separator porosity"] - b_sep = param["Separator Bruggeman coefficient (electrolyte)"] - t_plus = param["Cation transference number"] - sigma_e = ( - param["Electrolyte conductivity [S.m-1]"] # (ce0, T) - * (epsilon_sep**b_sep) - ) - - # Compute the cell area and thickness - A = param["Electrode height [m]"] * param["Electrode width [m]"] - L = L_p + L_n + L_s - - # Compute the series resistance - Rs = (L_p / sigma_p + L_n / sigma_n) / (3 * A) - R0 = Rs + param["Contact resistance [Ohm]"] - - # Compute the stoichiometry limits and initial SOC - x_0, x_100, y_100, y_0 = get_min_max_stoichiometries(param) - sto_p_init = ( - param["Initial concentration in positive electrode [mol.m-3]"] / c_max_p - ) - soc_init = (sto_p_init - y_0) / (y_100 - y_0) - - # Compute the capacity within the stoichiometry limits - Q_th_p = F * alpha_p * c_max_p * L_p * A - Q_th_n = F * alpha_n * c_max_n * L_n * A - Q_meas_p = (y_0 - y_100) * Q_th_p - Q_meas_n = (x_100 - x_0) * Q_th_n - if abs(Q_meas_n / Q_meas_p - 1) > 1e-6: - raise ValueError( - "The measured capacity should be the same for both electrodes." - ) - - # Grouped parameters - Q_meas = (Q_meas_n + Q_meas_p) / 2 - Q_e = F * epsilon_sep * ce0 * L * A - gamma_e = sigma_e / (F * epsilon_sep * ce0 * L**2) - - zeta_p = epsilon_p / epsilon_sep - zeta_n = epsilon_n / epsilon_sep - - try: - tau_d_p = R_p**2 / D_p - except TypeError: - tau_d_p = FunctionalDiffusionTime(R_p**2, D_p, c_max_p) - - try: - tau_d_n = R_n**2 / D_n - except TypeError: - tau_d_n = FunctionalDiffusionTime(R_n**2, D_n, c_max_n) - - tau_e = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) - beta_p = epsilon_p**b_p / epsilon_sep**b_sep - beta_n = epsilon_n**b_n / epsilon_sep**b_sep - - tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) - tau_ct_n = F * R_n / (m_n * np.sqrt(ce0)) - - C_p = 3 * alpha_p * Cdl_p * L_p * A / R_p - C_n = 3 * alpha_n * Cdl_n * L_n * A / R_n - - l_p = L_p / L - l_n = L_n / L - - parameter_dictionary = { - "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], - "Current function [A]": param["Current function [A]"], - "Ambient temperature [K]": T, - "Initial temperature [K]": T, - "Initial SoC": soc_init, - "Minimum negative stoichiometry": x_0, - "Maximum negative stoichiometry": x_100, - "Minimum positive stoichiometry": y_100, - "Maximum positive stoichiometry": y_0, - "Lower voltage cut-off [V]": param["Lower voltage cut-off [V]"], - "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], - "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], - "Negative electrode OCP [V]": param["Negative electrode OCP [V]"], - "Measured cell capacity [A.s]": Q_meas, - "Reference electrolyte capacity [A.s]": Q_e, - "Reference electrolyte scaled conductivity [V-1.s-1]": gamma_e, - "Positive electrode relative transport efficiency": beta_p, - "Negative electrode relative transport efficiency": beta_n, - "Positive electrode relative porosity": zeta_p, - "Negative electrode relative porosity": zeta_n, - "Positive particle diffusion time scale [s]": tau_d_p, - "Negative particle diffusion time scale [s]": tau_d_n, - "Electrolyte diffusion time scale [s]": tau_e, - "Positive electrode dimensionless exchange rate": GroupedDFN.symmetric_butler_volmer, - "Negative electrode dimensionless exchange rate": GroupedDFN.symmetric_butler_volmer, - "Positive electrode charge transfer time scale [s]": tau_ct_p, - "Negative electrode charge transfer time scale [s]": tau_ct_n, - "Positive electrode capacitance [F]": C_p, - "Negative electrode capacitance [F]": C_n, - "Cation transference number": t_plus, - "Positive electrode relative thickness": l_p, - "Negative electrode relative thickness": l_n, - "Series resistance [Ohm]": R0, - } - parameter_values = ParameterValues(values=parameter_dictionary) - parameter_values._set_initial_state = GroupedDFN.set_initial_state # noqa: SLF001 - return parameter_values - - @staticmethod - def symmetric_butler_volmer(sto_surf, sto_e, eta_RT_F): - """ - Dimensionless Butler-Volmer exchange rate. - """ - j0 = (sto_surf * sto_e * (1 - sto_surf)) ** 0.5 - return 2 * j0 * pybamm.sinh(0.5 * eta_RT_F) - - @staticmethod - def get_asymmetric_butler_volmer(domain: str): - """ - Get the asymmetric Butler-Volmer exchange rate function for the "positive" or - "negative" domain. - """ - Domain = domain.capitalize() - alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") - - return AsymmetricButlerVolmer(alpha) - - @staticmethod - def get_multiphase_butler_volmer(domain: str): - """ - Get the multiphase Butler-Volmer exchange rate function for the "positive" or - "negative" domain. - """ - Domain = domain.capitalize() - alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") - omega = pybamm.Parameter(f"{Domain} electrode charge transfer ideality factor") - - return MultiphaseButlerVolmer(alpha, omega) diff --git a/tests/integration/models/test_grouped_models.py b/tests/integration/models/test_grouped_models.py index 8946831ee..deda9f635 100644 --- a/tests/integration/models/test_grouped_models.py +++ b/tests/integration/models/test_grouped_models.py @@ -32,8 +32,6 @@ class TestGroupedModels: pybop.lithium_ion.GroupedSPM(options={"surface form": "differential"}), pybop.lithium_ion.GroupedSPMe(), pybop.lithium_ion.GroupedSPMe(options={"surface form": "differential"}), - pybop.lithium_ion.GroupedDFN(), - pybop.lithium_ion.GroupedDFN(options={"surface form": "differential"}), ], scope="module", ) diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index 4101adff8..47eb1ed49 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -6,11 +6,7 @@ import pybop -GROUPED_MODEL = ( - pybop.lithium_ion.GroupedSPM - | pybop.lithium_ion.GroupedSPMe - | pybop.lithium_ion.GroupedDFN -) +GROUPED_MODEL = pybop.lithium_ion.GroupedSPM | pybop.lithium_ion.GroupedSPMe class TestModels: @@ -30,8 +26,6 @@ class TestModels: pybop.lithium_ion.GroupedSPM(options={"surface form": "differential"}), pybop.lithium_ion.GroupedSPMe(), pybop.lithium_ion.GroupedSPMe(options={"surface form": "differential"}), - pybop.lithium_ion.GroupedDFN(), - pybop.lithium_ion.GroupedDFN(options={"surface form": "differential"}), ], scope="module", ) From 5b8a72f2d11394e393cf2a2e42813a0145e0d559 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 6 Aug 2026 21:47:44 +0100 Subject: [PATCH 28/37] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index f683a96db..b75e88531 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#974](https://github.com/pybop-team/PyBOP/pull/974) - Adds voltage components to each grouped model as well as asymmetric and multiphase Butler-Volmer kinetics. - [#969](https://github.com/pybop-team/PyBOP/pull/969) - Updates synthetic data and adds example script for thermal parameterisation. - [#965](https://github.com/pybop-team/PyBOP/pull/965) - Adds synthetic data and example scripts for OCV parameterisation. - [#963](https://github.com/pybop-team/PyBOP/pull/963) - Adds an example for generating synthetic data from a specification and exporting it to a PyProBE-compatible parquet file. From 73233c3c0608f72e0b573f693e4877e6bd8fb640 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 7 Aug 2026 00:00:31 +0100 Subject: [PATCH 29/37] Update description --- pybop/models/lithium_ion/grouped_spm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index f0980e180..4fe9535d7 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -339,7 +339,7 @@ def U(self, sto, domain): def tau_d(self, sto, T, domain): """ - Dimensional diffusion time scale [s]. + Dimensional solid-state diffusion time scale [s]. """ Domain = domain.capitalize() inputs = {f"{Domain} particle surface stoichiometry": sto, "Temperature [K]": T} From 90e58d0f42a12add1fd0e26f417b0142e837b1ec Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 7 Aug 2026 12:12:36 +0100 Subject: [PATCH 30/37] Add unit test --- tests/unit/test_models.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index 47eb1ed49..e49dac0ad 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -99,6 +99,27 @@ def test_set_initial_state(self, model): ): param.set_initial_state("-1 V") + def test_alternative_functions(self, model): + if isinstance(model, GROUPED_MODEL): + parameter_values = model.default_parameter_values + + # Use alternative Butler-Volmer kinetics + parameter_values["Negative electrode charge transfer coefficient"] = 0.6 + parameter_values["Negative electrode dimensionless exchange rate"] = ( + model.get_asymmetric_butler_volmer("negative") + ) + parameter_values["Positive electrode charge transfer coefficient"] = 0.6 + parameter_values["Positive electrode charge transfer ideality factor"] = 0.6 + parameter_values["Positive electrode dimensionless exchange rate"] = ( + model.get_multiphase_butler_volmer("positive") + ) + + t_eval = np.linspace(0, 10, 11) + solution = pybamm.Simulation( + model, parameter_values=parameter_values + ).solve(t_eval=t_eval, t_interp=t_eval) + np.testing.assert_allclose(solution["Time [s]"].data, t_eval) + class TestModelUtils: """ From 227ec5ca3bbba28032681ef8530589e2dcae86fd Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 7 Aug 2026 16:02:22 +0100 Subject: [PATCH 31/37] Move Butler-Volmer functions to base model --- pybop/models/lithium_ion/base_model.py | 35 ++++++++++++++++++++++++ pybop/models/lithium_ion/grouped_spm.py | 33 ---------------------- pybop/models/lithium_ion/grouped_spme.py | 33 ---------------------- 3 files changed, 35 insertions(+), 66 deletions(-) diff --git a/pybop/models/lithium_ion/base_model.py b/pybop/models/lithium_ion/base_model.py index f555d9302..857939589 100644 --- a/pybop/models/lithium_ion/base_model.py +++ b/pybop/models/lithium_ion/base_model.py @@ -2,6 +2,10 @@ from pybamm import FunctionParameter, Parameter from pybamm import lithium_ion as pybamm_lithium_ion +from pybop.models.lithium_ion.alternative_functions import ( + AsymmetricButlerVolmer, + MultiphaseButlerVolmer, +) from pybop.models.lithium_ion.utils import InverseOCV @@ -121,3 +125,34 @@ def ocv_function(soc): parameter_values["Initial SoC"] = soc return parameter_values + + @staticmethod + def symmetric_butler_volmer(sto_surf, sto_e, eta_RT_F): + """ + Dimensionless Butler-Volmer exchange rate. + """ + j0 = (sto_surf * sto_e * (1 - sto_surf)) ** 0.5 + return 2 * j0 * pybamm.sinh(0.5 * eta_RT_F) + + @staticmethod + def get_asymmetric_butler_volmer(domain: str): + """ + Get the asymmetric Butler-Volmer exchange rate function for the "positive" or + "negative" domain. + """ + Domain = domain.capitalize() + alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") + + return AsymmetricButlerVolmer(alpha) + + @staticmethod + def get_multiphase_butler_volmer(domain: str): + """ + Get the multiphase Butler-Volmer exchange rate function for the "positive" or + "negative" domain. + """ + Domain = domain.capitalize() + alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") + omega = pybamm.Parameter(f"{Domain} electrode charge transfer ideality factor") + + return MultiphaseButlerVolmer(alpha, omega) diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 4fe9535d7..459845939 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -15,9 +15,7 @@ ) from pybop.models.lithium_ion.alternative_functions import ( - AsymmetricButlerVolmer, FunctionalDiffusionTime, - MultiphaseButlerVolmer, ) from pybop.models.lithium_ion.base_model import BaseGroupedModel @@ -594,34 +592,3 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_values = ParameterValues(values=parameter_dictionary) parameter_values._set_initial_state = GroupedSPM.set_initial_state # noqa: SLF001 return parameter_values - - @staticmethod - def symmetric_butler_volmer(sto_surf, sto_e, eta_RT_F): - """ - Dimensionless Butler-Volmer exchange rate. - """ - j0 = (sto_surf * (1 - sto_surf)) ** 0.5 - return 2 * j0 * pybamm.sinh(0.5 * eta_RT_F) - - @staticmethod - def get_asymmetric_butler_volmer(domain: str): - """ - Get the asymmetric Butler-Volmer exchange rate function for the "positive" or - "negative" domain. - """ - Domain = domain.capitalize() - alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") - - return AsymmetricButlerVolmer(alpha) - - @staticmethod - def get_multiphase_butler_volmer(domain: str): - """ - Get the multiphase Butler-Volmer exchange rate function for the "positive" or - "negative" domain. - """ - Domain = domain.capitalize() - alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") - omega = pybamm.Parameter(f"{Domain} electrode charge transfer ideality factor") - - return MultiphaseButlerVolmer(alpha, omega) diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 09e37ffb3..2898dbad8 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -16,9 +16,7 @@ ) from pybop.models.lithium_ion.alternative_functions import ( - AsymmetricButlerVolmer, FunctionalDiffusionTime, - MultiphaseButlerVolmer, ) from pybop.models.lithium_ion.base_model import BaseGroupedModel @@ -720,34 +718,3 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_values = ParameterValues(values=parameter_dictionary) parameter_values._set_initial_state = GroupedSPMe.set_initial_state # noqa: SLF001 return parameter_values - - @staticmethod - def symmetric_butler_volmer(sto_surf, sto_e, eta_RT_F): - """ - Dimensionless Butler-Volmer exchange rate. - """ - j0 = (sto_surf * sto_e * (1 - sto_surf)) ** 0.5 - return 2 * j0 * pybamm.sinh(0.5 * eta_RT_F) - - @staticmethod - def get_asymmetric_butler_volmer(domain: str): - """ - Get the asymmetric Butler-Volmer exchange rate function for the "positive" or - "negative" domain. - """ - Domain = domain.capitalize() - alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") - - return AsymmetricButlerVolmer(alpha) - - @staticmethod - def get_multiphase_butler_volmer(domain: str): - """ - Get the multiphase Butler-Volmer exchange rate function for the "positive" or - "negative" domain. - """ - Domain = domain.capitalize() - alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") - omega = pybamm.Parameter(f"{Domain} electrode charge transfer ideality factor") - - return MultiphaseButlerVolmer(alpha, omega) From ee1285a7fefb7ed3d2e111b775d5e9d49fe17771 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Sat, 8 Aug 2026 00:38:51 +0100 Subject: [PATCH 32/37] Fix reaction rate calculation --- pybop/models/lithium_ion/grouped_spm.py | 18 ++++++++++++------ pybop/models/lithium_ion/grouped_spme.py | 18 ++++++++++++------ pybop/models/lithium_ion/sp_diffusion.py | 10 +++++++--- 3 files changed, 31 insertions(+), 15 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 459845939..900ca6492 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -488,12 +488,6 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] - m_p = param.evaluate( - param["Positive electrode exchange-current density [A.m-2]"](1, 1, 2, T) - ) # (A/m2)(m3/mol)**1.5 - m_n = param.evaluate( - param["Negative electrode exchange-current density [A.m-2]"](1, 1, 2, T) - ) # (A/m2)(m3/mol)**1.5 sigma_p = ( param["Positive electrode conductivity [S.m-1]"] * alpha_p ** param["Positive electrode Bruggeman coefficient (electrode)"] @@ -510,6 +504,18 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_sep = param["Separator Bruggeman coefficient (electrolyte)"] kappa_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) + # Get reaction rates + m_p = param.evaluate( + param["Positive electrode exchange-current density [A.m-2]"]( + ce0, c_max_p / 2, c_max_p, T + ) + ) * (2 / (c_max_p * np.sqrt(ce0))) # (A/m2)(m3/mol)**1.5 + m_n = param.evaluate( + param["Negative electrode exchange-current density [A.m-2]"]( + ce0, c_max_n / 2, c_max_n, T + ) + ) * (2 / (c_max_n * np.sqrt(ce0))) # (A/m2)(m3/mol)**1.5 + # Compute the cell area and thickness A = param["Electrode height [m]"] * param["Electrode width [m]"] L = L_p + L_n + L_s diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 2898dbad8..1bec749a5 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -597,12 +597,6 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] - m_p = param.evaluate( - param["Positive electrode exchange-current density [A.m-2]"](1, 1, 2, T) - ) # (A/m2)(m3/mol)**1.5 - m_n = param.evaluate( - param["Negative electrode exchange-current density [A.m-2]"](1, 1, 2, T) - ) # (A/m2)(m3/mol)**1.5 sigma_p = ( param["Positive electrode conductivity [S.m-1]"] * alpha_p ** param["Positive electrode Bruggeman coefficient (electrode)"] @@ -621,6 +615,18 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal t_plus = param["Cation transference number"] kappa_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) + # Get reaction rates + m_p = param.evaluate( + param["Positive electrode exchange-current density [A.m-2]"]( + ce0, c_max_p / 2, c_max_p, T + ) + ) * (2 / (c_max_p * np.sqrt(ce0))) # (A/m2)(m3/mol)**1.5 + m_n = param.evaluate( + param["Negative electrode exchange-current density [A.m-2]"]( + ce0, c_max_n / 2, c_max_n, T + ) + ) * (2 / (c_max_n * np.sqrt(ce0))) # (A/m2)(m3/mol)**1.5 + # Compute the cell area and thickness A = param["Electrode height [m]"] * param["Electrode width [m]"] L = L_p + L_n + L_s diff --git a/pybop/models/lithium_ion/sp_diffusion.py b/pybop/models/lithium_ion/sp_diffusion.py index e6b384ada..dd3aa5c0f 100644 --- a/pybop/models/lithium_ion/sp_diffusion.py +++ b/pybop/models/lithium_ion/sp_diffusion.py @@ -239,9 +239,14 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal L_p = param["Positive electrode thickness [m]"] R_p = param["Positive particle radius [m]"] D_p = param["Positive particle diffusivity [m2.s-1]"] + + # Get reaction rate + ce0 = param["Initial concentration in electrolyte [mol.m-3]"] m_p = param.evaluate( - param["Positive electrode exchange-current density [A.m-2]"](1, 1, 2, T) - ) # (A/m2)(m3/mol)**1.5 + param["Positive electrode exchange-current density [A.m-2]"]( + ce0, c_max_p / 2, c_max_p, T + ) + ) * (2 / (c_max_p * np.sqrt(ce0))) # (A/m2)(m3/mol)**1.5 # Compute the cell area A = param["Electrode height [m]"] * param["Electrode width [m]"] @@ -263,7 +268,6 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Estimate the series resistance, neglecting conductivity losses RT_F = pybamm.constants.R.value * param["Ambient temperature [K]"] / F - ce0 = param["Initial concentration in electrolyte [mol.m-3]"] tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) Rct_typ = (2 * RT_F * tau_ct_p) / (3 * Q_th_p) R0 = Rct_typ + param["Contact resistance [Ohm]"] From eef65c0e78264a54490d6e2d94b2629e4e2fee0d Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Sat, 8 Aug 2026 00:52:43 +0100 Subject: [PATCH 33/37] Switch reaction rate to reference current --- pybop/models/lithium_ion/grouped_spm.py | 15 +++++++-------- pybop/models/lithium_ion/grouped_spme.py | 15 +++++++-------- pybop/models/lithium_ion/sp_diffusion.py | 9 ++++----- 3 files changed, 18 insertions(+), 21 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 900ca6492..282d6863b 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -1,4 +1,3 @@ -import numpy as np import pybamm from pybamm import ( Event, @@ -504,17 +503,17 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_sep = param["Separator Bruggeman coefficient (electrolyte)"] kappa_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) - # Get reaction rates - m_p = param.evaluate( + # Get reference exchange current density [A.m-2] + j0_p = param.evaluate( param["Positive electrode exchange-current density [A.m-2]"]( ce0, c_max_p / 2, c_max_p, T ) - ) * (2 / (c_max_p * np.sqrt(ce0))) # (A/m2)(m3/mol)**1.5 - m_n = param.evaluate( + ) + j0_n = param.evaluate( param["Negative electrode exchange-current density [A.m-2]"]( ce0, c_max_n / 2, c_max_n, T ) - ) * (2 / (c_max_n * np.sqrt(ce0))) # (A/m2)(m3/mol)**1.5 + ) # Compute the cell area and thickness A = param["Electrode height [m]"] * param["Electrode width [m]"] @@ -559,8 +558,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal except TypeError: tau_d_n = FunctionalDiffusionTime(R_n**2, D_n, c_max_n) - tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) - tau_ct_n = F * R_n / (m_n * np.sqrt(ce0)) + tau_ct_p = c_max_p * F * R_p / (2 * j0_p) + tau_ct_n = c_max_n * F * R_n / (2 * j0_n) C_p = 3 * alpha_p * Cdl_p * L_p * A / R_p C_n = 3 * alpha_n * Cdl_n * L_n * A / R_n diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 1bec749a5..aabe23f79 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -1,4 +1,3 @@ -import numpy as np import pybamm from pybamm import ( Event, @@ -615,17 +614,17 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal t_plus = param["Cation transference number"] kappa_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) - # Get reaction rates - m_p = param.evaluate( + # Get reference exchange current density [A.m-2] + j0_p = param.evaluate( param["Positive electrode exchange-current density [A.m-2]"]( ce0, c_max_p / 2, c_max_p, T ) - ) * (2 / (c_max_p * np.sqrt(ce0))) # (A/m2)(m3/mol)**1.5 - m_n = param.evaluate( + ) + j0_n = param.evaluate( param["Negative electrode exchange-current density [A.m-2]"]( ce0, c_max_n / 2, c_max_n, T ) - ) * (2 / (c_max_n * np.sqrt(ce0))) # (A/m2)(m3/mol)**1.5 + ) # Compute the cell area and thickness A = param["Electrode height [m]"] * param["Electrode width [m]"] @@ -678,8 +677,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal beta_p = epsilon_p**b_p / epsilon_sep**b_sep beta_n = epsilon_n**b_n / epsilon_sep**b_sep - tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) - tau_ct_n = F * R_n / (m_n * np.sqrt(ce0)) + tau_ct_p = c_max_p * F * R_p / (2 * j0_p) + tau_ct_n = c_max_n * F * R_n / (2 * j0_n) C_p = 3 * alpha_p * Cdl_p * L_p * A / R_p C_n = 3 * alpha_n * Cdl_n * L_n * A / R_n diff --git a/pybop/models/lithium_ion/sp_diffusion.py b/pybop/models/lithium_ion/sp_diffusion.py index dd3aa5c0f..7b699bcf1 100644 --- a/pybop/models/lithium_ion/sp_diffusion.py +++ b/pybop/models/lithium_ion/sp_diffusion.py @@ -1,4 +1,3 @@ -import numpy as np import pybamm from pybamm import ( Event, @@ -240,13 +239,13 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal R_p = param["Positive particle radius [m]"] D_p = param["Positive particle diffusivity [m2.s-1]"] - # Get reaction rate + # Get reference exchange current density [A.m-2] ce0 = param["Initial concentration in electrolyte [mol.m-3]"] - m_p = param.evaluate( + j0_p = param.evaluate( param["Positive electrode exchange-current density [A.m-2]"]( ce0, c_max_p / 2, c_max_p, T ) - ) * (2 / (c_max_p * np.sqrt(ce0))) # (A/m2)(m3/mol)**1.5 + ) # Compute the cell area A = param["Electrode height [m]"] * param["Electrode width [m]"] @@ -268,7 +267,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Estimate the series resistance, neglecting conductivity losses RT_F = pybamm.constants.R.value * param["Ambient temperature [K]"] / F - tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) + tau_ct_p = c_max_p * F * R_p / (2 * j0_p) Rct_typ = (2 * RT_F * tau_ct_p) / (3 * Q_th_p) R0 = Rct_typ + param["Contact resistance [Ohm]"] From abc079cfb7e2d1df25abdf7f5e86516d89f01254 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Mon, 17 Aug 2026 15:51:10 +0100 Subject: [PATCH 34/37] Rename prefactor to r2_scale --- pybop/models/lithium_ion/alternative_functions.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pybop/models/lithium_ion/alternative_functions.py b/pybop/models/lithium_ion/alternative_functions.py index e565b22fb..991e93d82 100644 --- a/pybop/models/lithium_ion/alternative_functions.py +++ b/pybop/models/lithium_ion/alternative_functions.py @@ -4,14 +4,14 @@ class FunctionalDiffusionTime: - def __init__(self, D_prefactor, D, c_scale): - self.D_prefactor = D_prefactor + def __init__(self, r2_scale, D, c_scale): + self.r2_scale = r2_scale self.D = D self.c_scale = c_scale def __call__(self, sto, T): - D_prefactor, D, c_scale = self.D_prefactor, self.D, self.c_scale - return D_prefactor / D(sto * c_scale, T) + r2_scale, D, c_scale = self.r2_scale, self.D, self.c_scale + return r2_scale / D(sto * c_scale, T) class AsymmetricButlerVolmer: From d028501db56cbff0de3f9623bd08c635479e350a Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Mon, 17 Aug 2026 15:51:21 +0100 Subject: [PATCH 35/37] Move alternative_functions.py --- pybop/models/{lithium_ion => }/alternative_functions.py | 0 pybop/models/lithium_ion/base_model.py | 2 +- pybop/models/lithium_ion/grouped_spm.py | 4 +--- pybop/models/lithium_ion/grouped_spme.py | 4 +--- 4 files changed, 3 insertions(+), 7 deletions(-) rename pybop/models/{lithium_ion => }/alternative_functions.py (100%) diff --git a/pybop/models/lithium_ion/alternative_functions.py b/pybop/models/alternative_functions.py similarity index 100% rename from pybop/models/lithium_ion/alternative_functions.py rename to pybop/models/alternative_functions.py diff --git a/pybop/models/lithium_ion/base_model.py b/pybop/models/lithium_ion/base_model.py index 857939589..d0af05b4e 100644 --- a/pybop/models/lithium_ion/base_model.py +++ b/pybop/models/lithium_ion/base_model.py @@ -2,7 +2,7 @@ from pybamm import FunctionParameter, Parameter from pybamm import lithium_ion as pybamm_lithium_ion -from pybop.models.lithium_ion.alternative_functions import ( +from pybop.models.alternative_functions import ( AsymmetricButlerVolmer, MultiphaseButlerVolmer, ) diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 282d6863b..5ba251979 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -13,9 +13,7 @@ get_min_max_stoichiometries, ) -from pybop.models.lithium_ion.alternative_functions import ( - FunctionalDiffusionTime, -) +from pybop.models.alternative_functions import FunctionalDiffusionTime from pybop.models.lithium_ion.base_model import BaseGroupedModel diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index aabe23f79..2f1cf035d 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -14,9 +14,7 @@ get_min_max_stoichiometries, ) -from pybop.models.lithium_ion.alternative_functions import ( - FunctionalDiffusionTime, -) +from pybop.models.alternative_functions import FunctionalDiffusionTime from pybop.models.lithium_ion.base_model import BaseGroupedModel From 70f13cf2220022557366fb65c96baf77cafde6ce Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Mon, 17 Aug 2026 17:16:51 +0100 Subject: [PATCH 36/37] Separate half-cell models --- .../battery_parameterisation/gitt_fitting.py | 4 +- .../battery_parameterisation/gitt_pulse.py | 2 +- .../comparison_examples/gitt_models.py | 18 +-- pybop/__init__.py | 1 + pybop/applications/gitt_methods.py | 2 +- pybop/models/li_half_cell/__init__.py | 5 + pybop/models/li_half_cell/base_model.py | 120 +++++++++++++++++ .../sp_diffusion.py | 123 ++---------------- .../weppner_huggins.py | 4 +- pybop/models/lithium_ion/__init__.py | 2 - tests/integration/models/test_gitt_models.py | 4 +- tests/integration/test_applications.py | 8 +- tests/unit/test_models.py | 8 +- 13 files changed, 165 insertions(+), 136 deletions(-) create mode 100644 pybop/models/li_half_cell/__init__.py create mode 100644 pybop/models/li_half_cell/base_model.py rename pybop/models/{lithium_ion => li_half_cell}/sp_diffusion.py (68%) rename pybop/models/{lithium_ion => li_half_cell}/weppner_huggins.py (98%) diff --git a/examples/scripts/battery_parameterisation/gitt_fitting.py b/examples/scripts/battery_parameterisation/gitt_fitting.py index d6ea887eb..8c067cec6 100644 --- a/examples/scripts/battery_parameterisation/gitt_fitting.py +++ b/examples/scripts/battery_parameterisation/gitt_fitting.py @@ -47,7 +47,7 @@ ) # Group the parameters -grouped_parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( +grouped_parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters( parameter_values ) @@ -66,7 +66,7 @@ pybop.plot.dataset(gitt_parameter_data, signal=["Series resistance [Ohm]"]) # Run the identified model -identified_model = pybop.lithium_ion.SPDiffusion(build=True) +identified_model = pybop.li_half_cell.SPDiffusion(build=True) grouped_parameter_values.update(gitt_fit.best_inputs) grouped_parameter_values["Current function [A]"] = pybamm.Interpolant( dataset["Time [s]"], dataset["Current [A]"], pybamm.t diff --git a/examples/scripts/battery_parameterisation/gitt_pulse.py b/examples/scripts/battery_parameterisation/gitt_pulse.py index 1ce651e42..be48cd90d 100644 --- a/examples/scripts/battery_parameterisation/gitt_pulse.py +++ b/examples/scripts/battery_parameterisation/gitt_pulse.py @@ -30,7 +30,7 @@ ) # Group the parameters -grouped_parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( +grouped_parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters( parameter_values ) diff --git a/examples/scripts/comparison_examples/gitt_models.py b/examples/scripts/comparison_examples/gitt_models.py index 27071dbe3..d6481139a 100644 --- a/examples/scripts/comparison_examples/gitt_models.py +++ b/examples/scripts/comparison_examples/gitt_models.py @@ -2,6 +2,7 @@ import pybamm import pybop +from pybop.models.li_half_cell import SPDiffusion, WeppnerHuggins # Define model and parameter values model_options = {"working electrode": "positive"} @@ -30,13 +31,13 @@ } ) -for model in [pybop.lithium_ion.WeppnerHuggins(), pybop.lithium_ion.SPDiffusion()]: +for model in [WeppnerHuggins(), SPDiffusion()]: # GITT target parameter diffusion_parameter = pybop.Parameter(pybop.Gaussian(5000, 1000)) - if isinstance(model, pybop.lithium_ion.WeppnerHuggins): + if isinstance(model, WeppnerHuggins): # Group parameter values - grouped_parameter_values = ( - pybop.lithium_ion.WeppnerHuggins.create_grouped_parameters(parameter_values) + grouped_parameter_values = WeppnerHuggins.create_grouped_parameters( + parameter_values ) # We can fit only the duration of the pulse @@ -70,8 +71,8 @@ else: # Group parameter values - grouped_parameter_values = ( - pybop.lithium_ion.SPDiffusion.create_grouped_parameters(parameter_values) + grouped_parameter_values = SPDiffusion.create_grouped_parameters( + parameter_values ) # Fitting parameters @@ -87,7 +88,7 @@ # Build the problem gitt_dataset = ( dataset.get_subset(pulse_index) - if isinstance(model, pybop.lithium_ion.WeppnerHuggins) + if isinstance(model, WeppnerHuggins) else dataset ) simulator = pybop.pybamm.Simulator( @@ -97,7 +98,8 @@ problem = pybop.Problem(simulator, cost) # Build the optimisation problem - optim = pybop.SciPyMinimize(problem) + options = pybop.SciPyMinimizeOptions(method="Nelder-Mead") + optim = pybop.SciPyMinimize(problem, options) # Run the optimisation problem result = optim.run() diff --git a/pybop/__init__.py b/pybop/__init__.py index 228d1e7f2..8d78e365f 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -56,6 +56,7 @@ # Model classes # from .models import lithium_ion +from .models import li_half_cell from .models._exponential_decay import ExponentialDecayModel from .models.lithium_ion.utils import Interpolant, InverseOCV diff --git a/pybop/applications/gitt_methods.py b/pybop/applications/gitt_methods.py index e580392af..ca6fa3499 100644 --- a/pybop/applications/gitt_methods.py +++ b/pybop/applications/gitt_methods.py @@ -46,7 +46,7 @@ def __init__( self.optimiser_options = optimiser_options or self.optimiser.default_options() # Create model - self.model = pybop.lithium_ion.SPDiffusion() + self.model = pybop.li_half_cell.SPDiffusion() self.problem = None def __call__( diff --git a/pybop/models/li_half_cell/__init__.py b/pybop/models/li_half_cell/__init__.py new file mode 100644 index 000000000..5ecba98cd --- /dev/null +++ b/pybop/models/li_half_cell/__init__.py @@ -0,0 +1,5 @@ +# +# Import lithium half cell models +# +from .sp_diffusion import SPDiffusion +from .weppner_huggins import WeppnerHuggins diff --git a/pybop/models/li_half_cell/base_model.py b/pybop/models/li_half_cell/base_model.py new file mode 100644 index 000000000..b5a2027ec --- /dev/null +++ b/pybop/models/li_half_cell/base_model.py @@ -0,0 +1,120 @@ +import pybamm +from pybamm import FunctionParameter, Parameter +from pybamm import lithium_ion as pybamm_lithium_ion + +from pybop.models.lithium_ion.utils import InverseOCV + + +class BaseHalfCellModel(pybamm_lithium_ion.BaseModel): + """ + A base model for PyBOP's grouped-parameter lithium-ion half-cell models. + + Parameters + ---------- + name : str, optional + The name of the model. + **model_kwargs : optional + Valid PyBaMM model option keys and their values, for example: + options : dict, optional + A dictionary of options to customise the behaviour of the PyBaMM model. + build : bool, optional + If True, the model is built upon creation (default: False). + """ + + def __init__(self, name="Base Model", **model_kwargs): + super().__init__(name=name, **model_kwargs) + + def build_model(self): + """ + Build model variables and equations + Credit: PyBaMM + """ + self._build_model() + + self._built = True + pybamm.logger.info(f"Finish building {self.name}") + + @staticmethod + def set_initial_state( + initial_value, + parameter_values, + direction=None, + param=None, + inplace=True, + options=None, + inputs=None, + tol=1e-6, + ): + """ + Set the value of the initial state of charge. + + Parameters + ---------- + initial_value : float + Target initial value. + If float, interpreted as SOC, must be between 0 and 1. + If string e.g. "4 V", interpreted as voltage, must be between V_min and V_max. + parameter_values : :class:`pybamm.ParameterValues` + Parameters and their corresponding values. + param : :class:`pybamm.LithiumIonParameters`, optional + The symbolic parameter set to use for the simulation. + If not provided, the default parameter set will be used. + inplace: bool, optional + If True, replace the parameters values in place. Otherwise, return a new set of + parameter values. Default is True. + options : dict-like, optional + A dictionary of options to be passed to the model, see + :class:`pybamm.BatteryModelOptions`. + inputs : dict, optional + A dictionary of input parameters to pass to the model when solving. + tol : float, optional + The tolerance for the solver used to compute the initial stoichiometries. + A lower value results in higher precision but may increase computation time. + Default is 1e-6. + """ + parameter_values = parameter_values if inplace else parameter_values.copy() + + if isinstance(initial_value, str) and initial_value.endswith("V"): + V_init = float(initial_value[:-1]) + V_min = parameter_values.evaluate( + Parameter("Lower voltage cut-off [V]"), inputs=inputs + ) + V_max = parameter_values.evaluate( + Parameter("Upper voltage cut-off [V]"), inputs=inputs + ) + + if not V_min - tol <= V_init <= V_max + tol: + raise ValueError( + f"Initial voltage {V_init}V is outside the voltage limits ({V_min}, {V_max})." + ) + + y_100 = parameter_values.evaluate( + Parameter("Minimum positive stoichiometry"), inputs=inputs + ) + y_0 = parameter_values.evaluate( + Parameter("Maximum positive stoichiometry"), inputs=inputs + ) + + def ocv_function(soc): + sto_p = y_0 - soc * (y_0 - y_100) + U_p = FunctionParameter( + "Positive electrode OCP [V]", + {"Positive particle stoichiometry": sto_p}, + ) + return parameter_values.evaluate(U_p, inputs=inputs).squeeze() + + inverse_ocv = InverseOCV(ocv_function) + soc = inverse_ocv(V_init) + + elif isinstance(initial_value, int | float): + soc = initial_value + + else: + raise ValueError("Initial value must be a float or a string ending in 'V'.") + + if not 0 <= soc <= 1: + raise ValueError("Initial SOC should be between 0 and 1.") + + parameter_values["Initial SoC"] = soc + + return parameter_values diff --git a/pybop/models/lithium_ion/sp_diffusion.py b/pybop/models/li_half_cell/sp_diffusion.py similarity index 68% rename from pybop/models/lithium_ion/sp_diffusion.py rename to pybop/models/li_half_cell/sp_diffusion.py index 7b699bcf1..d4bf49aa4 100644 --- a/pybop/models/lithium_ion/sp_diffusion.py +++ b/pybop/models/li_half_cell/sp_diffusion.py @@ -14,11 +14,10 @@ get_min_max_stoichiometries, ) -from pybop.models.lithium_ion.base_model import BaseGroupedModel -from pybop.models.lithium_ion.utils import InverseOCV +from pybop.models.li_half_cell.base_model import BaseHalfCellModel -class SPDiffusion(BaseGroupedModel): +class SPDiffusion(BaseHalfCellModel): """ Diffusion model for a single, spherical particle representing a half-cell for GITT. @@ -48,6 +47,10 @@ def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): # Variables that vary spatially are created with a domain sto_p = Variable("Positive particle stoichiometry", domain="positive particle") + + # 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 sto_p_surf = pybamm.surf(sto_p) # Events specify points at which a solution should terminate @@ -112,14 +115,11 @@ def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): ###################### # Cell voltage ###################### - sto_p_average = sto_p_init + Q * 3600 / Q_th_p # pybamm.r_average(sto_p) - U = self.U(sto_p_surf, "positive") - self.U(sto_p_average, "negative") + U = self.U(sto_p_surf) V = U - self.R0(sto_p_surf) * I # Save the initial OCV - self.param.ocv_init = self.U(sto_p_init, "positive") - self.U( - sto_p_init, "negative" - ) + self.param.ocv_init = self.U(sto_p_init) # Events specify points at which a solution should terminate self.events += [ @@ -145,22 +145,15 @@ def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): "Open-circuit voltage [V]": U, } - def U(self, sto, domain): + def U(self, sto): """ Dimensional open-circuit potential [V]. Credit: PyBaMM """ - Domain = domain.capitalize() - if domain == "negative": - inputs = {"Average positive particle stoichiometry": sto} - else: - inputs = {f"{Domain} particle surface stoichiometry": sto} - out = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) - - if domain == "negative": - out.print_name = r"U_\mathrm{n}(c^\mathrm{av}_\mathrm{s,p})" - elif domain == "positive": - out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" + inputs = {"Positive particle surface stoichiometry": sto} + out = FunctionParameter("Positive electrode OCP [V]", inputs) + + out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" return out def tau_d(self, sto): @@ -280,7 +273,6 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Lower voltage cut-off [V]": param["Lower voltage cut-off [V]"], "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], - "Negative electrode OCP [V]": 0.0, "Measured cell capacity [A.s]": Q_meas, "Positive particle diffusion time scale [s]": tau_d_p, "Series resistance [Ohm]": R0, @@ -288,92 +280,3 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_values = ParameterValues(values=parameter_dictionary) parameter_values._set_initial_state = SPDiffusion.set_initial_state # noqa: SLF001 return parameter_values - - @staticmethod - def set_initial_state( - initial_value, - parameter_values, - direction=None, - param=None, - inplace=True, - options=None, - inputs=None, - tol=1e-6, - ): - """ - Set the value of the initial state of charge. - - Parameters - ---------- - initial_value : float - Target initial value. - If float, interpreted as SOC, must be between 0 and 1. - If string e.g. "4 V", interpreted as voltage, must be between V_min and V_max. - parameter_values : :class:`pybamm.ParameterValues` - Parameters and their corresponding values. - param : :class:`pybamm.LithiumIonParameters`, optional - The symbolic parameter set to use for the simulation. - If not provided, the default parameter set will be used. - inplace: bool, optional - If True, replace the parameters values in place. Otherwise, return a new set of - parameter values. Default is True. - options : dict-like, optional - A dictionary of options to be passed to the model, see - :class:`pybamm.BatteryModelOptions`. - inputs : dict, optional - A dictionary of input parameters to pass to the model when solving. - tol : float, optional - The tolerance for the solver used to compute the initial stoichiometries. - A lower value results in higher precision but may increase computation time. - Default is 1e-6. - """ - parameter_values = parameter_values if inplace else parameter_values.copy() - - if isinstance(initial_value, str) and initial_value.endswith("V"): - V_init = float(initial_value[:-1]) - V_min = parameter_values.evaluate( - Parameter("Lower voltage cut-off [V]"), inputs=inputs - ) - V_max = parameter_values.evaluate( - Parameter("Upper voltage cut-off [V]"), inputs=inputs - ) - - if not V_min - tol <= V_init <= V_max + tol: - raise ValueError( - f"Initial voltage {V_init}V is outside the voltage limits ({V_min}, {V_max})." - ) - - y_100 = parameter_values.evaluate( - Parameter("Minimum positive stoichiometry"), inputs=inputs - ) - y_0 = parameter_values.evaluate( - Parameter("Maximum positive stoichiometry"), inputs=inputs - ) - - def ocv_function(soc): - sto_p = y_0 - soc * (y_0 - y_100) - U_p = FunctionParameter( - "Positive electrode OCP [V]", - {"Positive particle stoichiometry": sto_p}, - ) - U_n = FunctionParameter( - "Negative electrode OCP [V]", - {"Positive particle stoichiometry": sto_p}, - ) - return parameter_values.evaluate(U_p - U_n, inputs=inputs).squeeze() - - inverse_ocv = InverseOCV(ocv_function) - soc = inverse_ocv(V_init) - - elif isinstance(initial_value, int | float): - soc = initial_value - - else: - raise ValueError("Initial value must be a float or a string ending in 'V'.") - - if not 0 <= soc <= 1: - raise ValueError("Initial SOC should be between 0 and 1.") - - parameter_values["Initial SoC"] = soc - - return parameter_values diff --git a/pybop/models/lithium_ion/weppner_huggins.py b/pybop/models/li_half_cell/weppner_huggins.py similarity index 98% rename from pybop/models/lithium_ion/weppner_huggins.py rename to pybop/models/li_half_cell/weppner_huggins.py index b57e206f8..96c62b820 100644 --- a/pybop/models/lithium_ion/weppner_huggins.py +++ b/pybop/models/li_half_cell/weppner_huggins.py @@ -3,10 +3,10 @@ from pybamm import DummySolver, Parameter, ParameterValues from pybamm import t as pybamm_t -from pybop.models.lithium_ion.base_model import BaseGroupedModel +from pybop.models.li_half_cell.base_model import BaseHalfCellModel -class WeppnerHuggins(BaseGroupedModel): +class WeppnerHuggins(BaseHalfCellModel): """ Represents the Weppner & Huggins model to fit diffusion coefficients to GITT data. diff --git a/pybop/models/lithium_ion/__init__.py b/pybop/models/lithium_ion/__init__.py index 7dcbbc032..1745f9b42 100644 --- a/pybop/models/lithium_ion/__init__.py +++ b/pybop/models/lithium_ion/__init__.py @@ -1,8 +1,6 @@ # # Import lithium ion models # -from .sp_diffusion import SPDiffusion from .grouped_spm import GroupedSPM from .grouped_spme import GroupedSPMe -from .weppner_huggins import WeppnerHuggins from .cell_temperature import CellTemperature diff --git a/tests/integration/models/test_gitt_models.py b/tests/integration/models/test_gitt_models.py index 794cf9286..225e977f1 100644 --- a/tests/integration/models/test_gitt_models.py +++ b/tests/integration/models/test_gitt_models.py @@ -28,8 +28,8 @@ class TestGITTModels: @pytest.fixture( params=[ - pybop.lithium_ion.WeppnerHuggins(), - pybop.lithium_ion.SPDiffusion(), + pybop.li_half_cell.WeppnerHuggins(), + pybop.li_half_cell.SPDiffusion(), ], scope="module", ) diff --git a/tests/integration/test_applications.py b/tests/integration/test_applications.py index e80fa320a..f3173c9d6 100644 --- a/tests/integration/test_applications.py +++ b/tests/integration/test_applications.py @@ -48,14 +48,14 @@ def charge_dataset(self, parameter_values): ) def test_interpolant(self, parameter_values, discharge_dataset): - parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( + parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters( parameter_values ) parameter_values["Positive electrode OCP [V]"] = pybop.Interpolant( discharge_dataset["Stoichiometry"], discharge_dataset["Voltage [V]"] ) parameter_values.set_initial_state(0.9) - model = pybop.lithium_ion.SPDiffusion(build=True) + model = pybop.li_half_cell.SPDiffusion(build=True) t_eval = np.linspace(0, 10, 100) solution = pybamm.Simulation(model, parameter_values=parameter_values).solve( t_eval=t_eval, t_interp=t_eval @@ -187,7 +187,7 @@ def pulse_data(self, half_cell_model, half_cell_parameter_values): def test_gitt_pulse_fit( self, half_cell_model, half_cell_parameter_values, pulse_data ): - parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( + parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters( half_cell_parameter_values ) diffusion_time = parameter_values["Positive particle diffusion time scale [s]"] @@ -202,7 +202,7 @@ def test_gitt_pulse_fit( ) def test_gitt_fit(self, half_cell_model, half_cell_parameter_values, pulse_data): - parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( + parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters( half_cell_parameter_values ) diffusion_time = parameter_values["Positive particle diffusion time scale [s]"] diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index e49dac0ad..9aa28c3a0 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -20,8 +20,8 @@ class TestModels: params=[ pybop.ExponentialDecayModel(), pybop.lithium_ion.CellTemperature(), - pybop.lithium_ion.WeppnerHuggins(), - pybop.lithium_ion.SPDiffusion(), + pybop.li_half_cell.WeppnerHuggins(), + pybop.li_half_cell.SPDiffusion(), pybop.lithium_ion.GroupedSPM(), pybop.lithium_ion.GroupedSPM(options={"surface form": "differential"}), pybop.lithium_ion.GroupedSPMe(), @@ -68,7 +68,7 @@ def test_set_initial_state(self, model): if isinstance(model, pybop.ExponentialDecayModel): pass # Only testing the battery models for now - elif isinstance(model, pybop.lithium_ion.WeppnerHuggins): + elif isinstance(model, pybop.li_half_cell.WeppnerHuggins): param = model.default_parameter_values with pytest.raises( ValueError, @@ -93,7 +93,7 @@ def test_set_initial_state(self, model): with pytest.raises(ValueError, match="should be between 0 and 1."): param.set_initial_state(-1) - if not isinstance(model, pybop.lithium_ion.SPDiffusion): + if not isinstance(model, pybop.li_half_cell.SPDiffusion): with pytest.raises( ValueError, match=r"V is outside the voltage limits" ): From 9f8432279af5aff43a48fd544019d814320f8f2f Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Mon, 17 Aug 2026 20:07:33 +0100 Subject: [PATCH 37/37] Reset set_initial_state and use A.h --- .../lgm50_pulse_validation.ipynb | 2 +- .../comparison_examples/gitt_models.py | 2 +- .../dfn_parameterisation/3_cell_balancing.py | 2 +- pybop/models/li_half_cell/base_model.py | 22 ++++--------- pybop/models/li_half_cell/sp_diffusion.py | 33 +++++-------------- pybop/models/li_half_cell/weppner_huggins.py | 6 ++-- pybop/models/lithium_ion/base_model.py | 4 +-- pybop/models/lithium_ion/cell_temperature.py | 8 ++--- pybop/models/lithium_ion/grouped_spm.py | 10 +++--- pybop/models/lithium_ion/grouped_spme.py | 16 ++++----- tests/integration/models/test_gitt_models.py | 2 +- tests/unit/test_models.py | 14 +++++--- 12 files changed, 52 insertions(+), 69 deletions(-) diff --git a/examples/notebooks/battery_parameterisation/lgm50_pulse_validation.ipynb b/examples/notebooks/battery_parameterisation/lgm50_pulse_validation.ipynb index 9509c47cf..b091a7316 100644 --- a/examples/notebooks/battery_parameterisation/lgm50_pulse_validation.ipynb +++ b/examples/notebooks/battery_parameterisation/lgm50_pulse_validation.ipynb @@ -223,7 +223,7 @@ "id": "12", "metadata": {}, "source": [ - "We can construct the model and parameter values for a two-RC circuit model with initial values as listed below. Note, the \"Initial SOC\" is shifted slightly to better match the zero degree data." + "We can construct the model and parameter values for a two-RC circuit model with initial values as listed below. Note, the \"Initial SoC\" is shifted slightly to better match the zero degree data." ] }, { diff --git a/examples/scripts/comparison_examples/gitt_models.py b/examples/scripts/comparison_examples/gitt_models.py index d6481139a..913385a85 100644 --- a/examples/scripts/comparison_examples/gitt_models.py +++ b/examples/scripts/comparison_examples/gitt_models.py @@ -50,7 +50,7 @@ dataset["Discharge capacity [A.h]"][-1] - dataset["Discharge capacity [A.h]"][0] ) - * (grouped_parameter_values["Theoretical electrode capacity [A.s]"] / 3600) + * grouped_parameter_values["Theoretical electrode capacity [A.h]"] ) grouped_parameter_values.update( { diff --git a/examples/scripts/dfn_parameterisation/3_cell_balancing.py b/examples/scripts/dfn_parameterisation/3_cell_balancing.py index c67fd8407..acd05271c 100644 --- a/examples/scripts/dfn_parameterisation/3_cell_balancing.py +++ b/examples/scripts/dfn_parameterisation/3_cell_balancing.py @@ -354,7 +354,7 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): "Maximum negative stoichiometry": x_100, "Minimum positive stoichiometry": y_100, "Maximum positive stoichiometry": y_0, - "Measured cell capacity [A.s]": Q_soc / CE * 3600, + "Measured cell capacity [A.h]": Q_soc / CE, "Coulombic efficiency": CE, "Negative electrode pOCP [V]": negative_ocp_function, "Positive electrode pOCP [V]": positive_ocp_function, diff --git a/pybop/models/li_half_cell/base_model.py b/pybop/models/li_half_cell/base_model.py index b5a2027ec..d75de6089 100644 --- a/pybop/models/li_half_cell/base_model.py +++ b/pybop/models/li_half_cell/base_model.py @@ -52,7 +52,7 @@ def set_initial_state( ---------- initial_value : float Target initial value. - If float, interpreted as SOC, must be between 0 and 1. + If float, interpreted as stoichiometry, must be between 0 and 1. If string e.g. "4 V", interpreted as voltage, must be between V_min and V_max. parameter_values : :class:`pybamm.ParameterValues` Parameters and their corresponding values. @@ -88,15 +88,7 @@ def set_initial_state( f"Initial voltage {V_init}V is outside the voltage limits ({V_min}, {V_max})." ) - y_100 = parameter_values.evaluate( - Parameter("Minimum positive stoichiometry"), inputs=inputs - ) - y_0 = parameter_values.evaluate( - Parameter("Maximum positive stoichiometry"), inputs=inputs - ) - - def ocv_function(soc): - sto_p = y_0 - soc * (y_0 - y_100) + def ocv_function(sto_p): U_p = FunctionParameter( "Positive electrode OCP [V]", {"Positive particle stoichiometry": sto_p}, @@ -104,17 +96,17 @@ def ocv_function(soc): return parameter_values.evaluate(U_p, inputs=inputs).squeeze() inverse_ocv = InverseOCV(ocv_function) - soc = inverse_ocv(V_init) + sto_p = inverse_ocv(V_init) elif isinstance(initial_value, int | float): - soc = initial_value + sto_p = initial_value else: raise ValueError("Initial value must be a float or a string ending in 'V'.") - if not 0 <= soc <= 1: - raise ValueError("Initial SOC should be between 0 and 1.") + if not 0 <= sto_p <= 1: + raise ValueError("Initial stoichiometry should be between 0 and 1.") - parameter_values["Initial SoC"] = soc + parameter_values["Initial stoichiometry"] = sto_p return parameter_values diff --git a/pybop/models/li_half_cell/sp_diffusion.py b/pybop/models/li_half_cell/sp_diffusion.py index d4bf49aa4..366144e4c 100644 --- a/pybop/models/li_half_cell/sp_diffusion.py +++ b/pybop/models/li_half_cell/sp_diffusion.py @@ -10,9 +10,6 @@ Variable, ) from pybamm import t as pybamm_t -from pybamm.models.full_battery_models.lithium_ion.electrode_soh_half_cell import ( - get_min_max_stoichiometries, -) from pybop.models.li_half_cell.base_model import BaseHalfCellModel @@ -72,12 +69,10 @@ def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): # Parameters are purely symbolic at this stage, and will be set by the # `ParameterValues` class when the model is processed. - soc_init = Parameter("Initial SoC") - y_100 = Parameter("Minimum positive stoichiometry") - y_0 = Parameter("Maximum positive stoichiometry") - # Grouped parameters - Q_th_p = Parameter("Measured cell capacity [A.s]") / (y_0 - y_100) + Q_th_p = Parameter("Theoretical electrode capacity [A.h]") * 3600 + + sto_p_init = Parameter("Initial stoichiometry") ###################### # Input current (positive on discharge) @@ -109,7 +104,6 @@ def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): "right": (-self.tau_d(sto_p_surf) * j_p, "Neumann"), } - sto_p_init = y_0 + (y_100 - y_0) * soc_init self.initial_conditions[sto_p] = sto_p_init ###################### @@ -243,37 +237,28 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Compute the cell area A = param["Electrode height [m]"] * param["Electrode width [m]"] - # Compute the stoichiometry limits and initial SOC - d = get_min_max_stoichiometries(param) - y_0, y_100 = d["x_0"], d["x_100"] + # Compute the initial stoichiometry sto_p_init = ( param["Initial concentration in positive electrode [mol.m-3]"] / c_max_p ) - soc_init = (sto_p_init - y_0) / (y_100 - y_0) - - # Compute the capacity within the stoichiometry limits - Q_th_p = F * alpha_p * c_max_p * L_p * A - Q_meas = (y_0 - y_100) * Q_th_p # Grouped parameters + Q_th_p = F * alpha_p * c_max_p * L_p * A / 3600 tau_d_p = R_p**2 / D_p # Estimate the series resistance, neglecting conductivity losses - RT_F = pybamm.constants.R.value * param["Ambient temperature [K]"] / F - tau_ct_p = c_max_p * F * R_p / (2 * j0_p) - Rct_typ = (2 * RT_F * tau_ct_p) / (3 * Q_th_p) + RT_F = pybamm.constants.R.value * T / F + Rct_typ = (RT_F * R_p) / (3 * alpha_p * L_p * A * j0_p) R0 = Rct_typ + param["Contact resistance [Ohm]"] parameter_dictionary = { "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], "Current function [A]": param["Current function [A]"], - "Initial SoC": soc_init, - "Minimum positive stoichiometry": y_100, - "Maximum positive stoichiometry": y_0, + "Initial stoichiometry": sto_p_init, "Lower voltage cut-off [V]": param["Lower voltage cut-off [V]"], "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], - "Measured cell capacity [A.s]": Q_meas, + "Theoretical electrode capacity [A.h]": Q_th_p, "Positive particle diffusion time scale [s]": tau_d_p, "Series resistance [Ohm]": R0, } diff --git a/pybop/models/li_half_cell/weppner_huggins.py b/pybop/models/li_half_cell/weppner_huggins.py index 96c62b820..78f007f2a 100644 --- a/pybop/models/li_half_cell/weppner_huggins.py +++ b/pybop/models/li_half_cell/weppner_huggins.py @@ -48,7 +48,7 @@ def __init__(self, name="Weppner & Huggins model", **model_kwargs): # Parameters are purely symbolic at this stage, and will be set by the # `ParameterValues` class when the model is processed. - Q_th_p = Parameter("Theoretical electrode capacity [A.s]") + Q_th_p = Parameter("Theoretical electrode capacity [A.h]") * 3600 U = Parameter("Reference voltage [V]") U_prime = Parameter("Derivative of the OCP wrt stoichiometry [V]") @@ -146,14 +146,14 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal A = param["Electrode height [m]"] * param["Electrode width [m]"] # Grouped parameters - Q_th_p = F * alpha_p * c_max_p * L_p * A + Q_th_p = F * alpha_p * c_max_p * L_p * A / 3600 tau_d_p = R_p**2 / D_p parameter_dictionary = { "Current function [A]": param["Current function [A]"], "Reference voltage [V]": 4, "Derivative of the OCP wrt stoichiometry [V]": -1, - "Theoretical electrode capacity [A.s]": Q_th_p, + "Theoretical electrode capacity [A.h]": Q_th_p, "Positive particle diffusion time scale [s]": tau_d_p, } parameter_values = ParameterValues(values=parameter_dictionary) diff --git a/pybop/models/lithium_ion/base_model.py b/pybop/models/lithium_ion/base_model.py index d0af05b4e..9844b29aa 100644 --- a/pybop/models/lithium_ion/base_model.py +++ b/pybop/models/lithium_ion/base_model.py @@ -56,7 +56,7 @@ def set_initial_state( ---------- initial_value : float Target initial value. - If float, interpreted as SOC, must be between 0 and 1. + If float, interpreted as SoC, must be between 0 and 1. If string e.g. "4 V", interpreted as voltage, must be between V_min and V_max. parameter_values : :class:`pybamm.ParameterValues` Parameters and their corresponding values. @@ -120,7 +120,7 @@ def ocv_function(soc): raise ValueError("Initial value must be a float or a string ending in 'V'.") if not 0 <= soc <= 1: - raise ValueError("Initial SOC should be between 0 and 1.") + raise ValueError("Initial SoC should be between 0 and 1.") parameter_values["Initial SoC"] = soc diff --git a/pybop/models/lithium_ion/cell_temperature.py b/pybop/models/lithium_ion/cell_temperature.py index e2317bb48..3b5c40132 100644 --- a/pybop/models/lithium_ion/cell_temperature.py +++ b/pybop/models/lithium_ion/cell_temperature.py @@ -52,7 +52,7 @@ def __init__(self, name="Cell Temperature Model", **model_kwargs): # Parameters are purely symbolic at this stage, and will be set by the # `ParameterValues` class when the model is processed. - Q_meas = Parameter("Measured cell capacity [A.s]") + Q_meas = Parameter("Measured cell capacity [A.h]") * 3600 soc_init = Parameter("Initial SoC") x_0 = Parameter("Minimum negative stoichiometry") @@ -242,8 +242,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal soc_init = (sto_p_init - y_0) / (y_100 - y_0) # Compute the capacity within the stoichiometry limits - Q_th_p = F * alpha_p * c_max_p * L_p * A - Q_th_n = F * alpha_n * c_max_n * L_n * A + Q_th_p = F * alpha_p * c_max_p * L_p * A / 3600 + Q_th_n = F * alpha_n * c_max_n * L_n * A / 3600 Q_meas_p = (y_0 - y_100) * Q_th_p Q_meas_n = (x_100 - x_0) * Q_th_n if abs(Q_meas_n / Q_meas_p - 1) > 1e-6: @@ -270,7 +270,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], "Negative electrode OCP [V]": param["Negative electrode OCP [V]"], - "Measured cell capacity [A.s]": Q_meas, + "Measured cell capacity [A.h]": Q_meas, "Negative electrode OCP entropic change [V.K-1]": param[ "Negative electrode OCP entropic change [V.K-1]" ], diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 5ba251979..251854db0 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -120,8 +120,8 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): y_0 = Parameter("Maximum positive stoichiometry") # Grouped parameters - Q_th_p = Parameter("Measured cell capacity [A.s]") / (y_0 - y_100) - Q_th_n = Parameter("Measured cell capacity [A.s]") / (x_100 - x_0) + Q_th_p = Parameter("Measured cell capacity [A.h]") * 3600 / (y_0 - y_100) + Q_th_n = Parameter("Measured cell capacity [A.h]") * 3600 / (x_100 - x_0) tau_ct_p = Parameter("Positive electrode charge transfer time scale [s]") tau_ct_n = Parameter("Negative electrode charge transfer time scale [s]") @@ -534,8 +534,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal soc_init = (sto_p_init - y_0) / (y_100 - y_0) # Compute the capacity within the stoichiometry limits - Q_th_p = F * alpha_p * c_max_p * L_p * A - Q_th_n = F * alpha_n * c_max_n * L_n * A + Q_th_p = F * alpha_p * c_max_p * L_p * A / 3600 + Q_th_n = F * alpha_n * c_max_n * L_n * A / 3600 Q_meas_p = (y_0 - y_100) * Q_th_p Q_meas_n = (x_100 - x_0) * Q_th_n if abs(Q_meas_n / Q_meas_p - 1) > 1e-6: @@ -579,7 +579,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], "Negative electrode OCP [V]": param["Negative electrode OCP [V]"], - "Measured cell capacity [A.s]": Q_meas, + "Measured cell capacity [A.h]": Q_meas, "Positive particle diffusion time scale [s]": tau_d_p, "Negative particle diffusion time scale [s]": tau_d_n, "Positive electrode dimensionless exchange rate": GroupedSPM.symmetric_butler_volmer, diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 2f1cf035d..e9c027110 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -153,9 +153,9 @@ def __init__( y_0 = Parameter("Maximum positive stoichiometry") # Grouped parameters - Q_th_p = Parameter("Measured cell capacity [A.s]") / (y_0 - y_100) - Q_th_n = Parameter("Measured cell capacity [A.s]") / (x_100 - x_0) - Q_e = Parameter("Reference electrolyte capacity [A.s]") + Q_th_p = Parameter("Measured cell capacity [A.h]") * 3600 / (y_0 - y_100) + Q_th_n = Parameter("Measured cell capacity [A.h]") * 3600 / (x_100 - x_0) + Q_e = Parameter("Reference electrolyte capacity [A.h]") * 3600 tau_ct_p = Parameter("Positive electrode charge transfer time scale [s]") tau_ct_n = Parameter("Negative electrode charge transfer time scale [s]") @@ -645,8 +645,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal soc_init = (sto_p_init - y_0) / (y_100 - y_0) # Compute the capacity within the stoichiometry limits - Q_th_p = F * alpha_p * c_max_p * L_p * A - Q_th_n = F * alpha_n * c_max_n * L_n * A + Q_th_p = F * alpha_p * c_max_p * L_p * A / 3600 + Q_th_n = F * alpha_n * c_max_n * L_n * A / 3600 Q_meas_p = (y_0 - y_100) * Q_th_p Q_meas_n = (x_100 - x_0) * Q_th_n if abs(Q_meas_n / Q_meas_p - 1) > 1e-6: @@ -656,7 +656,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Grouped parameters Q_meas = (Q_meas_n + Q_meas_p) / 2 - Q_e = F * epsilon_sep * ce0 * L * A + Q_e = F * epsilon_sep * ce0 * L * A / 3600 zeta_p = epsilon_p / epsilon_sep zeta_n = epsilon_n / epsilon_sep @@ -698,8 +698,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], "Negative electrode OCP [V]": param["Negative electrode OCP [V]"], - "Measured cell capacity [A.s]": Q_meas, - "Reference electrolyte capacity [A.s]": Q_e, + "Measured cell capacity [A.h]": Q_meas, + "Reference electrolyte capacity [A.h]": Q_e, "Positive electrode relative porosity": zeta_p, "Negative electrode relative porosity": zeta_n, "Positive particle diffusion time scale [s]": tau_d_p, diff --git a/tests/integration/models/test_gitt_models.py b/tests/integration/models/test_gitt_models.py index 225e977f1..7c690c5a8 100644 --- a/tests/integration/models/test_gitt_models.py +++ b/tests/integration/models/test_gitt_models.py @@ -14,7 +14,7 @@ # Parameter configurations DIFFUSION_PARAMS = [ - ("Theoretical electrode capacity [A.s]", 10), + ("Theoretical electrode capacity [A.h]", 0.003), ("Positive particle diffusion time scale [s]", 2000), ] diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index 9aa28c3a0..28db1fae7 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -20,13 +20,14 @@ class TestModels: params=[ pybop.ExponentialDecayModel(), pybop.lithium_ion.CellTemperature(), - pybop.li_half_cell.WeppnerHuggins(), - pybop.li_half_cell.SPDiffusion(), pybop.lithium_ion.GroupedSPM(), pybop.lithium_ion.GroupedSPM(options={"surface form": "differential"}), pybop.lithium_ion.GroupedSPMe(), pybop.lithium_ion.GroupedSPMe(options={"surface form": "differential"}), + pybop.li_half_cell.WeppnerHuggins(), + pybop.li_half_cell.SPDiffusion(), ], + ids=lambda val: f"{type(val).__name__}", scope="module", ) def model(self, request): @@ -77,12 +78,17 @@ def test_set_initial_state(self, model): param.set_initial_state(0.5) else: + if isinstance(model, pybop.li_half_cell.SPDiffusion): + initial_state = "Initial stoichiometry" + else: + initial_state = "Initial SoC" + param = model.default_parameter_values param.set_initial_state(0.5) - assert param["Initial SoC"] == 0.5 + assert param[initial_state] == 0.5 param.set_initial_state("3.8 V") - assert 0 <= param["Initial SoC"] <= 1 + assert 0 <= param[initial_state] <= 1 with pytest.raises( ValueError,