From e6743a5ec865b167ccf853f9b57f154f12f5aea0 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 6 Aug 2026 23:47:45 +0100 Subject: [PATCH 1/7] Add GroupedDFN model --- pybop/models/lithium_ion/__init__.py | 1 + pybop/models/lithium_ion/grouped_dfn.py | 730 ++++++++++++++++++ .../integration/models/test_grouped_models.py | 2 + tests/unit/test_models.py | 10 +- 4 files changed, 741 insertions(+), 2 deletions(-) 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 1745f9b42..25d20968a 100644 --- a/pybop/models/lithium_ion/__init__.py +++ b/pybop/models/lithium_ion/__init__.py @@ -3,4 +3,5 @@ # from .grouped_spm import GroupedSPM from .grouped_spme import GroupedSPMe +from .grouped_dfn import GroupedDFN from .cell_temperature import CellTemperature diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py new file mode 100644 index 000000000..42381c205 --- /dev/null +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -0,0 +1,730 @@ +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.alternative_functions import FunctionalDiffusionTime +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.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]") + + 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]"] + 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) + ) + + # 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 + ) + ) + j0_n = param.evaluate( + param["Negative electrode exchange-current density [A.m-2]"]( + ce0, c_max_n / 2, c_max_n, T + ) + ) + + # 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 / 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: + 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 / 3600 + 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 = 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 + + 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.h]": Q_meas, + "Reference electrolyte capacity [A.h]": 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 diff --git a/tests/integration/models/test_grouped_models.py b/tests/integration/models/test_grouped_models.py index deda9f635..8946831ee 100644 --- a/tests/integration/models/test_grouped_models.py +++ b/tests/integration/models/test_grouped_models.py @@ -32,6 +32,8 @@ 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 28db1fae7..6199cd6a1 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -6,12 +6,16 @@ import pybop -GROUPED_MODEL = pybop.lithium_ion.GroupedSPM | pybop.lithium_ion.GroupedSPMe +GROUPED_MODEL = ( + pybop.lithium_ion.GroupedSPM + | pybop.lithium_ion.GroupedSPMe + | pybop.lithium_ion.GroupedDFN +) class TestModels: """ - A class to test pybop created models. + A class to test pybop models. """ pytestmark = pytest.mark.unit @@ -24,6 +28,8 @@ 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"}), pybop.li_half_cell.WeppnerHuggins(), pybop.li_half_cell.SPDiffusion(), ], From a5f75f5daffec9ddafd203233ec8b2a4783567a1 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Thu, 6 Aug 2026 23:52:22 +0100 Subject: [PATCH 2/7] Allow electrolyte diffusivity function --- .../comparison_examples/grouped_SPMe.py | 5 +-- pybop/models/lithium_ion/grouped_dfn.py | 38 ++++++++++++++----- pybop/models/lithium_ion/grouped_spme.py | 33 +++++++++++----- 3 files changed, 54 insertions(+), 22 deletions(-) diff --git a/examples/scripts/comparison_examples/grouped_SPMe.py b/examples/scripts/comparison_examples/grouped_SPMe.py index 426e6cb82..3155357bc 100644 --- a/examples/scripts/comparison_examples/grouped_SPMe.py +++ b/examples/scripts/comparison_examples/grouped_SPMe.py @@ -18,12 +18,9 @@ # Use the Chen2020 parameters parameter_values = pybamm.ParameterValues("Chen2020") -# Fix the electrolyte diffusivity and conductivity +# Fix the electrolyte conductivity ce0 = parameter_values["Initial concentration in electrolyte [mol.m-3]"] T = parameter_values["Ambient temperature [K]"] -parameter_values["Electrolyte diffusivity [m2.s-1]"] = parameter_values[ - "Electrolyte diffusivity [m2.s-1]" -](ce0, T) parameter_values["Electrolyte conductivity [S.m-1]"] = parameter_values[ "Electrolyte conductivity [S.m-1]" ](ce0, T) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 42381c205..6c328b5cb 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -161,7 +161,6 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): 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]") ###################### @@ -283,13 +282,24 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): # Electrolyte ###################### self.rhs[sto_e_n] = ( - pybamm.div(pybamm.grad(sto_e_n) * beta_n / tau_e + (1 - t_plus) * i_e_n) + pybamm.div( + pybamm.grad(sto_e_n) + * beta_n + / self.tau_e(sto_e_n, T, "negative electrode") + + (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 + pybamm.grad(sto_e_sep) / self.tau_e(sto_e_sep, T, "separator") + - 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) + pybamm.div( + pybamm.grad(sto_e_p) + * beta_p + / self.tau_e(sto_e_p, T, "positive electrode") + + (1 - t_plus) * i_e_p + ) ) / zeta_p self.boundary_conditions[sto_e_n] = { @@ -451,6 +461,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 tau_e(self, sto_e, T, domain): + """ + Dimensional electolyte diffusion time scale [s]. + """ + Domain = domain.capitalize() + inputs = {f"{Domain} electrolyte stoichiometry": sto_e, "Temperature [K]": T} + return FunctionParameter("Electrolyte diffusion time scale [s]", inputs) + def j(self, sto_surf, sto_e, eta_RT_F, domain): """ Dimensionless exchange rate. @@ -473,9 +491,6 @@ def default_parameter_values(self) -> pybamm.ParameterValues: 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 @@ -612,7 +627,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Separator and electrolyte properties ce0 = param["Initial concentration in electrolyte [mol.m-3]"] - De = param["Electrolyte diffusivity [m2.s-1]"] # (ce0, T) + De = param["Electrolyte diffusivity [m2.s-1]"] L_s = param["Separator thickness [m]"] epsilon_sep = param["Separator porosity"] b_sep = param["Separator Bruggeman coefficient (electrolyte)"] @@ -677,7 +692,12 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal 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) + try: + tau_e = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) + except TypeError: + De_prefactor = epsilon_sep * L**2 / epsilon_sep**b_sep + tau_e = FunctionalDiffusionTime(De_prefactor, De, ce0) + beta_p = epsilon_p**b_p / epsilon_sep**b_sep beta_n = epsilon_n**b_n / epsilon_sep**b_sep diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index e9c027110..9df8ca7f8 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -170,7 +170,6 @@ def __init__( zeta_n = Parameter("Negative electrode relative porosity") zeta_p = Parameter("Positive electrode relative porosity") - tau_e = Parameter("Electrolyte diffusion time scale [s]") beta_n = Parameter("Negative electrode relative transport efficiency") beta_p = Parameter("Positive electrode relative transport efficiency") @@ -278,16 +277,22 @@ def __init__( ###################### self.rhs[sto_e_n] = ( pybamm.div( - pybamm.grad(sto_e_n) * beta_n / tau_e - (t_plus * I / Q_e) * x_n / l_n + pybamm.grad(sto_e_n) + * beta_n + / self.tau_e(sto_e_n, T, "negative electrode") + - (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 - t_plus * I / Q_e + pybamm.grad(sto_e_sep) / self.tau_e(sto_e_sep, T, "separator") + - t_plus * I / Q_e ) self.rhs[sto_e_p] = ( pybamm.div( - pybamm.grad(sto_e_p) * beta_p / tau_e + pybamm.grad(sto_e_p) + * beta_p + / self.tau_e(sto_e_p, T, "positive electrode") - (t_plus * I / Q_e) * (1 - x_p) / l_p ) + (3 / Q_e) * Q_th_p * j_p / l_p @@ -444,6 +449,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 tau_e(self, sto_e, T, domain): + """ + Dimensional electolyte diffusion time scale [s]. + """ + Domain = domain.capitalize() + inputs = {f"{Domain} electrolyte stoichiometry": sto_e, "Temperature [K]": T} + return FunctionParameter("Electrolyte diffusion time scale [s]", inputs) + def j(self, sto_surf, sto_e, eta_RT_F, domain): """ Dimensionless exchange rate. @@ -466,9 +479,6 @@ def default_parameter_values(self) -> ParameterValues: 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 @@ -605,7 +615,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Separator and electrolyte properties ce0 = param["Initial concentration in electrolyte [mol.m-3]"] - De = param["Electrolyte diffusivity [m2.s-1]"] # (ce0, T) + De = param["Electrolyte diffusivity [m2.s-1]"] L_s = param["Separator thickness [m]"] epsilon_sep = param["Separator porosity"] b_sep = param["Separator Bruggeman coefficient (electrolyte)"] @@ -671,7 +681,12 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal 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) + try: + tau_e = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) + except TypeError: + De_prefactor = epsilon_sep * L**2 / epsilon_sep**b_sep + tau_e = FunctionalDiffusionTime(De_prefactor, De, ce0) + beta_p = epsilon_p**b_p / epsilon_sep**b_sep beta_n = epsilon_n**b_n / epsilon_sep**b_sep From 1f4b91898122ce8f19d2a59ffa16c263137d0b73 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Tue, 14 Jul 2026 16:31:37 +0100 Subject: [PATCH 3/7] Fix DFN electrolyte flux --- pybop/models/lithium_ion/grouped_dfn.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 6c328b5cb..a30a1df7e 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -204,7 +204,7 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): 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 + # Electrolyte currents [s-1] 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 @@ -286,8 +286,9 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): pybamm.grad(sto_e_n) * beta_n / self.tau_e(sto_e_n, T, "negative electrode") - + (1 - t_plus) * i_e_n + - t_plus * i_e_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) / self.tau_e(sto_e_sep, T, "separator") @@ -298,8 +299,9 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): pybamm.grad(sto_e_p) * beta_p / self.tau_e(sto_e_p, T, "positive electrode") - + (1 - t_plus) * i_e_p + - t_plus * i_e_p ) + + (3 / Q_e) * Q_th_p * j_p / l_p ) / zeta_p self.boundary_conditions[sto_e_n] = { From 62cb66cb6e6718ed41553773965041dc9591e84d Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Wed, 15 Jul 2026 18:04:45 +0100 Subject: [PATCH 4/7] Remove interface approximations --- pybop/models/lithium_ion/grouped_dfn.py | 75 ++++++++++---------- pybop/models/lithium_ion/grouped_spme.py | 89 +++++++++++------------- 2 files changed, 75 insertions(+), 89 deletions(-) diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index a30a1df7e..1af29b4fe 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -1,5 +1,6 @@ import pybamm from pybamm import ( + ConcatenationVariable, Event, FunctionParameter, Parameter, @@ -95,6 +96,7 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): "Positive electrode electrolyte stoichiometry", domain="positive electrode", ) + sto_e = ConcatenationVariable(sto_e_n, sto_e_sep, sto_e_p) # 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 @@ -281,45 +283,43 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): ###################### # Electrolyte ###################### - self.rhs[sto_e_n] = ( - pybamm.div( - pybamm.grad(sto_e_n) - * beta_n - / self.tau_e(sto_e_n, T, "negative electrode") - - t_plus * i_e_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) / self.tau_e(sto_e_sep, T, "separator") - - t_plus * I / Q_e + beta = pybamm.concatenation( + PrimaryBroadcast(beta_n, "negative electrode"), + PrimaryBroadcast(Scalar(1), "separator"), + PrimaryBroadcast(beta_p, "positive electrode"), ) - self.rhs[sto_e_p] = ( - pybamm.div( - pybamm.grad(sto_e_p) - * beta_p - / self.tau_e(sto_e_p, T, "positive electrode") - - t_plus * i_e_p - ) - + (3 / Q_e) * Q_th_p * j_p / l_p - ) / zeta_p + tau_e = pybamm.concatenation( + self.tau_e(sto_e_n, T, "negative electrode"), + self.tau_e(sto_e_sep, T, "separator"), + self.tau_e(sto_e_p, T, "positive electrode"), + ) + i_e = pybamm.concatenation( + i_e_n, + PrimaryBroadcast(I / Q_e, "separator"), + i_e_p, + ) + N_e = -pybamm.grad(sto_e) * beta / tau_e + t_plus * i_e + + j_e_n = (3 / Q_e) * Q_th_n * j_n / l_n + j_e_sep = PrimaryBroadcast(Scalar(0), "separator") + j_e_p = (3 / Q_e) * Q_th_p * j_p / l_p + j_e = pybamm.concatenation(j_e_n, j_e_sep, j_e_p) + + zeta = pybamm.concatenation( + PrimaryBroadcast(zeta_n, "negative electrode"), + PrimaryBroadcast(Scalar(1), "separator"), + PrimaryBroadcast(zeta_p, "positive electrode"), + ) + + # The concatenated stoichiometry (sto_e) and flux (N_e) are continuous across interfaces + self.rhs[sto_e] = (-pybamm.div(N_e) + j_e) / zeta - self.boundary_conditions[sto_e_n] = { + self.boundary_conditions[sto_e] = { "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) + self.initial_conditions[sto_e] = Scalar(1) # Electrolyte overpotential eta_e = (2 * (1 - t_plus) * RT_F) * ( @@ -416,18 +416,15 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): "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 - ), + "Electrolyte stoichiometry": sto_e, "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 - ), + "Electrolyte scaled current density [s-1]": i_e, + "Electrolyte flux [s-1]": N_e, "Time [s]": pybamm_t, "Time [h]": pybamm_t / 3600, "Current [A]": I, diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 9df8ca7f8..80c4d46bb 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -1,5 +1,6 @@ import pybamm from pybamm import ( + ConcatenationVariable, Event, FunctionParameter, Parameter, @@ -90,6 +91,7 @@ def __init__( "Positive electrode electrolyte stoichiometry", domain="positive electrode", ) + sto_e = ConcatenationVariable(sto_e_n, sto_e_sep, sto_e_p) # Spatial variables x_n = SpatialVariable("x_n", domain=["negative electrode"]) @@ -122,16 +124,8 @@ def __init__( ), # 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, + "Minimum electrolyte stoichiometry", + pybamm.min(sto_e) - 0, ), ] @@ -275,45 +269,43 @@ def __init__( ###################### # Electrolyte ###################### - self.rhs[sto_e_n] = ( - pybamm.div( - pybamm.grad(sto_e_n) - * beta_n - / self.tau_e(sto_e_n, T, "negative electrode") - - (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) / self.tau_e(sto_e_sep, T, "separator") - - t_plus * I / Q_e + beta = pybamm.concatenation( + PrimaryBroadcast(beta_n, "negative electrode"), + PrimaryBroadcast(Scalar(1), "separator"), + PrimaryBroadcast(beta_p, "positive electrode"), ) - self.rhs[sto_e_p] = ( - pybamm.div( - pybamm.grad(sto_e_p) - * beta_p - / self.tau_e(sto_e_p, T, "positive electrode") - - (t_plus * I / Q_e) * (1 - x_p) / l_p - ) - + (3 / Q_e) * Q_th_p * j_p / l_p - ) / zeta_p + tau_e = pybamm.concatenation( + self.tau_e(sto_e_n, T, "negative electrode"), + self.tau_e(sto_e_sep, T, "separator"), + self.tau_e(sto_e_p, T, "positive electrode"), + ) + i_e = pybamm.concatenation( + (I / Q_e) * x_n / l_n, + PrimaryBroadcast(I / Q_e, "separator"), + (I / Q_e) * (1 - x_p) / l_p, + ) + N_e = -pybamm.grad(sto_e) * beta / tau_e + t_plus * i_e - self.boundary_conditions[sto_e_n] = { + j_e_n = (3 / Q_e) * Q_th_n * j_n / l_n + j_e_sep = PrimaryBroadcast(Scalar(0), "separator") + j_e_p = (3 / Q_e) * Q_th_p * j_p / l_p + j_e = pybamm.concatenation(j_e_n, j_e_sep, j_e_p) + + zeta = pybamm.concatenation( + PrimaryBroadcast(zeta_n, "negative electrode"), + PrimaryBroadcast(Scalar(1), "separator"), + PrimaryBroadcast(zeta_p, "positive electrode"), + ) + + # The concatenated stoichiometry (sto_e) and flux (N_e) are continuous across interfaces + self.rhs[sto_e] = (-pybamm.div(N_e) + j_e) / zeta + + self.boundary_conditions[sto_e] = { "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) + self.initial_conditions[sto_e] = Scalar(1) ###################### # Cell voltage @@ -394,9 +386,7 @@ def __init__( "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 - ), + "Electrolyte stoichiometry": sto_e, "Positive particle stoichiometry": sto_p, "Positive particle surface stoichiometry": PrimaryBroadcast( sto_p_surf, "positive electrode" @@ -410,10 +400,9 @@ def __init__( - pybamm.boundary_value(eta_p, "right"), "Electrolyte potential [V]": -v_s_n - (2 * RT_F * (1 - t_plus)) - * ( - pybamm.boundary_value(pybamm.log(sto_e_n), "left") - - pybamm.log(pybamm.concatenation(sto_e_n, sto_e_sep, sto_e_p)) - ), + * (pybamm.boundary_value(pybamm.log(sto_e_n), "left") - pybamm.log(sto_e)), + "Electrolyte scaled current density [s-1]": i_e, + "Electrolyte flux [s-1]": N_e, "Time [s]": pybamm_t, "Time [h]": pybamm_t / 3600, "Current [A]": I, From 4da53bd45c90d6e7f846adab812cb66d25d30714 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 7 Aug 2026 17:19:01 +0100 Subject: [PATCH 5/7] Enable flux concatenations --- pyproject.toml | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 9e7540668..65c7d087a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,7 +29,7 @@ classifiers = [ # versions in the tests in nightly_dependency_tests.yml as appropriate requires-python = ">=3.10, <3.15" dependencies = [ - "pybamm[plot]>=26.3.0", + "pybamm @ git+https://github.com/NicolaCourtier/PyBaMM@flux-concatenation", "numpy>=1.26", "scipy>=1.12", "pints>=0.6.0", @@ -78,6 +78,10 @@ dev = [ "ruff", ] +# Temporarily allow dependencies to be direct references +[tool.hatch.metadata] +allow-direct-references = true + [tool.hatch.build.targets.sdist] exclude = [ "papers/", From 81277fbb63c577ddddd79036f36026447ea89ae9 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Fri, 7 Aug 2026 17:19:13 +0100 Subject: [PATCH 6/7] Use flux boundary conditions --- pybop/models/li_half_cell/sp_diffusion.py | 7 +-- pybop/models/lithium_ion/grouped_dfn.py | 55 +++++++++-------------- pybop/models/lithium_ion/grouped_spm.py | 24 ++++------ pybop/models/lithium_ion/grouped_spme.py | 51 +++++++++------------ 4 files changed, 55 insertions(+), 82 deletions(-) diff --git a/pybop/models/li_half_cell/sp_diffusion.py b/pybop/models/li_half_cell/sp_diffusion.py index 366144e4c..49106c919 100644 --- a/pybop/models/li_half_cell/sp_diffusion.py +++ b/pybop/models/li_half_cell/sp_diffusion.py @@ -95,13 +95,14 @@ 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_p] = pybamm.div(pybamm.grad(sto_p) / self.tau_d(sto_p)) + N_s_p = -pybamm.grad(sto_p) / self.tau_d(sto_p) + self.rhs[sto_p] = -pybamm.div(N_s_p) # Boundary conditions must be provided for equations with spatial derivatives j_p = -I / (3 * Q_th_p) self.boundary_conditions[sto_p] = { - "left": (Scalar(0), "Neumann"), - "right": (-self.tau_d(sto_p_surf) * j_p, "Neumann"), + "left": (Scalar(0), ("Flux", N_s_p)), + "right": (j_p, ("Flux", N_s_p)), } self.initial_conditions[sto_p] = sto_p_init diff --git a/pybop/models/lithium_ion/grouped_dfn.py b/pybop/models/lithium_ion/grouped_dfn.py index 1af29b4fe..481b28ad5 100644 --- a/pybop/models/lithium_ion/grouped_dfn.py +++ b/pybop/models/lithium_ion/grouped_dfn.py @@ -6,6 +6,7 @@ Parameter, ParameterValues, PrimaryBroadcast, + PrimaryBroadcastToEdges, Scalar, Variable, ) @@ -235,24 +236,12 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): 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", - ), + "left": (Scalar(0), ("Flux", i_e_n)), + "right": (I / Q_e, ("Flux", i_e_n)), } 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"), + "left": (I / Q_e, ("Flux", i_e_p)), + "right": (Scalar(0), ("Flux", i_e_p)), } ###################### @@ -260,21 +249,19 @@ 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, T, "negative") - ) - self.rhs[sto_p] = pybamm.div( - pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") - ) + N_s_n = -pybamm.grad(sto_n) / self.tau_d(sto_n, T, "negative") + N_s_p = -pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") + self.rhs[sto_n] = -pybamm.div(N_s_n) + self.rhs[sto_p] = -pybamm.div(N_s_p) # 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"), + "left": (Scalar(0), ("Flux", N_s_n)), + "right": (j_n, ("Flux", N_s_n)), } self.boundary_conditions[sto_p] = { - "left": (Scalar(0), "Neumann"), - "right": (-self.tau_d(sto_p_surf, T, "positive") * j_p, "Neumann"), + "left": (Scalar(0), ("Flux", N_s_p)), + "right": (j_p, ("Flux", N_s_p)), } self.initial_conditions[sto_n] = sto_n_init @@ -295,16 +282,16 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): ) i_e = pybamm.concatenation( i_e_n, - PrimaryBroadcast(I / Q_e, "separator"), + PrimaryBroadcastToEdges(I / Q_e, "separator"), i_e_p, ) N_e = -pybamm.grad(sto_e) * beta / tau_e + t_plus * i_e - j_e_n = (3 / Q_e) * Q_th_n * j_n / l_n - j_e_sep = PrimaryBroadcast(Scalar(0), "separator") - j_e_p = (3 / Q_e) * Q_th_p * j_p / l_p - j_e = pybamm.concatenation(j_e_n, j_e_sep, j_e_p) - + j_e = pybamm.concatenation( + (3 / Q_e) * Q_th_n * j_n / l_n, + PrimaryBroadcast(Scalar(0), "separator"), + (3 / Q_e) * Q_th_p * j_p / l_p, + ) zeta = pybamm.concatenation( PrimaryBroadcast(zeta_n, "negative electrode"), PrimaryBroadcast(Scalar(1), "separator"), @@ -315,8 +302,8 @@ def __init__(self, name="Grouped Doyle Fuller Newman Model", **model_kwargs): self.rhs[sto_e] = (-pybamm.div(N_e) + j_e) / zeta self.boundary_conditions[sto_e] = { - "left": (Scalar(0), "Neumann"), - "right": (Scalar(0), "Neumann"), + "left": (Scalar(0), ("Flux", N_e)), + "right": (Scalar(0), ("Flux", N_e)), } self.initial_conditions[sto_e] = Scalar(1) diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 251854db0..de9912e29 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -192,27 +192,19 @@ 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, T, "negative") - ) - self.rhs[sto_p] = pybamm.div( - pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") - ) + N_s_n = -pybamm.grad(sto_n) / self.tau_d(sto_n, T, "negative") + N_s_p = -pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") + self.rhs[sto_n] = -pybamm.div(N_s_n) + self.rhs[sto_p] = -pybamm.div(N_s_p) # 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") * pybamm.x_average(j_n), - "Neumann", - ), + "left": (Scalar(0), ("Flux", N_s_n)), + "right": (pybamm.x_average(j_n), ("Flux", N_s_n)), } self.boundary_conditions[sto_p] = { - "left": (Scalar(0), "Neumann"), - "right": ( - -self.tau_d(sto_p_surf, T, "positive") * pybamm.x_average(j_p), - "Neumann", - ), + "left": (Scalar(0), ("Flux", N_s_p)), + "right": (pybamm.x_average(j_p), ("Flux", N_s_p)), } self.initial_conditions[sto_n] = sto_n_init diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index 80c4d46bb..74ccc60b3 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -6,8 +6,9 @@ Parameter, ParameterValues, PrimaryBroadcast, + PrimaryBroadcastToEdges, Scalar, - SpatialVariable, + SpatialVariableEdge, Variable, ) from pybamm import t as pybamm_t @@ -94,8 +95,8 @@ def __init__( sto_e = ConcatenationVariable(sto_e_n, sto_e_sep, sto_e_p) # Spatial variables - x_n = SpatialVariable("x_n", domain=["negative electrode"]) - x_p = SpatialVariable("x_p", domain=["positive electrode"]) + x_n_edge = SpatialVariableEdge("x_n", domain=["negative electrode"]) + x_p_edge = SpatialVariableEdge("x_p", 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 @@ -240,27 +241,19 @@ 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, T, "negative") - ) - self.rhs[sto_p] = pybamm.div( - pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") - ) + N_s_n = -pybamm.grad(sto_n) / self.tau_d(sto_n, T, "negative") + N_s_p = -pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") + self.rhs[sto_n] = -pybamm.div(N_s_n) + self.rhs[sto_p] = -pybamm.div(N_s_p) # 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") * pybamm.x_average(j_n), - "Neumann", - ), + "left": (Scalar(0), ("Flux", N_s_n)), + "right": (pybamm.x_average(j_n), ("Flux", N_s_n)), } self.boundary_conditions[sto_p] = { - "left": (Scalar(0), "Neumann"), - "right": ( - -self.tau_d(sto_p_surf, T, "positive") * pybamm.x_average(j_p), - "Neumann", - ), + "left": (Scalar(0), ("Flux", N_s_p)), + "right": (pybamm.x_average(j_p), ("Flux", N_s_p)), } self.initial_conditions[sto_n] = sto_n_init @@ -280,17 +273,17 @@ def __init__( self.tau_e(sto_e_p, T, "positive electrode"), ) i_e = pybamm.concatenation( - (I / Q_e) * x_n / l_n, - PrimaryBroadcast(I / Q_e, "separator"), - (I / Q_e) * (1 - x_p) / l_p, + (I / Q_e) * x_n_edge / l_n, + PrimaryBroadcastToEdges(I / Q_e, "separator"), + (I / Q_e) * (1 - x_p_edge) / l_p, ) N_e = -pybamm.grad(sto_e) * beta / tau_e + t_plus * i_e - j_e_n = (3 / Q_e) * Q_th_n * j_n / l_n - j_e_sep = PrimaryBroadcast(Scalar(0), "separator") - j_e_p = (3 / Q_e) * Q_th_p * j_p / l_p - j_e = pybamm.concatenation(j_e_n, j_e_sep, j_e_p) - + j_e = pybamm.concatenation( + (3 / Q_e) * Q_th_n * j_n / l_n, + PrimaryBroadcast(Scalar(0), "separator"), + (3 / Q_e) * Q_th_p * j_p / l_p, + ) zeta = pybamm.concatenation( PrimaryBroadcast(zeta_n, "negative electrode"), PrimaryBroadcast(Scalar(1), "separator"), @@ -301,8 +294,8 @@ def __init__( self.rhs[sto_e] = (-pybamm.div(N_e) + j_e) / zeta self.boundary_conditions[sto_e] = { - "left": (Scalar(0), "Neumann"), - "right": (Scalar(0), "Neumann"), + "left": (Scalar(0), ("Flux", N_e)), + "right": (Scalar(0), ("Flux", N_e)), } self.initial_conditions[sto_e] = Scalar(1) From cc19ffdef57749b445dbda97c75adca034f8e930 Mon Sep 17 00:00:00 2001 From: NicolaCourtier <45851982+NicolaCourtier@users.noreply.github.com> Date: Wed, 19 Aug 2026 12:47:36 +0100 Subject: [PATCH 7/7] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 641fdcd6c..6ebe3c6e0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#977](https://github.com/pybop-team/PyBOP/pull/977) - Updates the `lithium_ion` models to use flux boundary conditions and adds the `GroupedDFN` model. - [#973](https://github.com/pybop-team/PyBOP/pull/973) - Adds example scripts for EIS parameterisation. - [#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.