Fixed merge conflict in changelog

CHANGELOG.md

@@ -16,10 +16,12 @@
 
 ## Bug fixes
 
+- Fixed a bug where the "basic" lithium-ion models gave incorrect results when using nonlinear particle diffusivity ([#3207](https://github.com/pybamm-team/PyBaMM/pull/3207))
+- Particle size distributions now work with SPMe and NewmanTobias models ([#3207](https://github.com/pybamm-team/PyBaMM/pull/3207))
 - Negative half-cell models now work, but the procedure to run them has changed ([#3198](https://github.com/pybamm-team/PyBaMM/pull/3198))
 - Fix to simulate c_rate steps with drive cycles ([#3186](https://github.com/pybamm-team/PyBaMM/pull/3186))
 - Parameters in `Prada2013` have been updated to better match those given in the paper, which is a 2.3 Ah cell, instead of the mix-and-match with the 1.1 Ah cell from Lain2019.
-- Error generated when invalid parameter values are passed.
+- Error generated when invalid parameter values are passed. ([#3132](https://github.com/pybamm-team/PyBaMM/pull/3132))
 - Thevenin() model is now constructed with standard variables: `Time [s], Time [min], Time [h]` ([#3143](https://github.com/pybamm-team/PyBaMM/pull/3143))
 - Fix SEI Example Notebook ([#3166](https://github.com/pybamm-team/PyBaMM/pull/3166))
 

examples/scripts/cycling_ageing.py

@@ -1,7 +1,7 @@
 import pybamm as pb
 
 pb.set_logging_level("NOTICE")
-model = pb.lithium_ion.DFN(
+model = pb.lithium_ion.SPM(
     {
         "SEI": "ec reaction limited",
         "SEI film resistance": "distributed",

noxfile.py

@@ -114,13 +114,13 @@ def set_dev(session):
     session.install("cmake")
     if sys.platform == "linux" or sys.platform == "darwin":
         session.run(
-        "echo",
-        "export",
-        f"LD_LIBRARY_PATH={PYBAMM_ENV['LD_LIBRARY_PATH']}",
-        ">>",
-        f"{envbindir}/activate",
-        external=True,  # silence warning about echo being an external command
-    )
+            "echo",
+            "export",
+            f"LD_LIBRARY_PATH={PYBAMM_ENV['LD_LIBRARY_PATH']}",
+            ">>",
+            f"{envbindir}/activate",
+            external=True,  # silence warning about echo being an external command
+        )
 
 
 @nox.session(name="tests")

pybamm/models/full_battery_models/lithium_ion/basic_dfn.py

@@ -170,14 +170,14 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (
-                -j_n / (param.F * param.n.prim.D(c_s_surf_n, T)),
+                -j_n / (param.F * pybamm.surf(param.n.prim.D(c_s_n, T))),
                 "Neumann",
             ),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (
-                -j_p / (param.F * param.p.prim.D(c_s_surf_p, T)),
+                -j_p / (param.F * pybamm.surf(param.p.prim.D(c_s_p, T))),
                 "Neumann",
             ),
         }

pybamm/models/full_battery_models/lithium_ion/basic_dfn_composite.py

@@ -219,15 +219,24 @@ def __init__(self, name="Composite graphite/silicon Doyle-Fuller-Newman model"):
         # Boundary conditions must be provided for equations with spatial derivatives
         self.boundary_conditions[c_s_n_p1] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-j_n_p1 / param.F / param.n.prim.D(c_s_surf_n_p1, T), "Neumann"),
+            "right": (
+                -j_n_p1 / param.F / pybamm.surf(param.n.prim.D(c_s_n_p1, T)),
+                "Neumann",
+            ),
         }
         self.boundary_conditions[c_s_n_p2] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-j_n_p2 / param.F / param.n.sec.D(c_s_surf_n_p2, T), "Neumann"),
+            "right": (
+                -j_n_p2 / param.F / pybamm.surf(param.n.sec.D(c_s_n_p2, T)),
+                "Neumann",
+            ),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-j_p / param.F / param.p.prim.D(c_s_surf_p, T), "Neumann"),
+            "right": (
+                -j_p / param.F / pybamm.surf(param.p.prim.D(c_s_p, T)),
+                "Neumann",
+            ),
         }
         self.initial_conditions[c_s_n_p1] = param.n.prim.c_init
         self.initial_conditions[c_s_n_p2] = param.n.sec.c_init

pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py

@@ -153,7 +153,7 @@ def __init__(self, options=None, name="Doyle-Fuller-Newman half cell model"):
         # derivatives
         self.boundary_conditions[c_s_w] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-j_w / D_w(c_s_surf_w, T) / param.F, "Neumann"),
+            "right": (-j_w / pybamm.surf(D_w(c_s_w, T)) / param.F, "Neumann"),
         }
         self.initial_conditions[c_s_w] = c_w_init
 

pybamm/models/full_battery_models/lithium_ion/basic_spm.py

@@ -86,14 +86,14 @@ def __init__(self, name="Single Particle Model"):
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (
-                -j_n / param.F / param.n.prim.D(c_s_surf_n, T),
+                -j_n / (param.F * pybamm.surf(param.n.prim.D(c_s_n, T))),
                 "Neumann",
             ),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (
-                -j_p / param.F / param.p.prim.D(c_s_surf_p, T),
+                -j_p / (param.F * pybamm.surf(param.p.prim.D(c_s_p, T))),
                 "Neumann",
             ),
         }

pybamm/models/full_battery_models/lithium_ion/electrode_soh.py

@@ -268,12 +268,8 @@ def solve(self, inputs):
         sol_dict = {key: sol[key].data[0] for key in sol.all_models[0].variables.keys()}
 
         # Calculate theoretical energy
-        x_0 = sol_dict["x_0"]
-        y_0 = sol_dict["y_0"]
-        x_100 = sol_dict["x_100"]
-        y_100 = sol_dict["y_100"]
         energy = pybamm.lithium_ion.electrode_soh.theoretical_energy_integral(
-            self.parameter_values, x_100, x_0, y_100, y_0
+            self.parameter_values, sol_dict
         )
         sol_dict.update({"Maximum theoretical energy [W.h]": energy})
         return sol_dict
@@ -598,7 +594,7 @@ def get_min_max_stoichiometries(
     return esoh_solver.get_min_max_stoichiometries()
 
 
-def theoretical_energy_integral(parameter_values, n_i, n_f, p_i, p_f, points=100):
+def theoretical_energy_integral(parameter_values, inputs, points=100):
     """
     Calculate maximum energy possible from a cell given OCV, initial soc, and final soc
     given voltage limits, open-circuit potentials, etc defined by parameter_values
@@ -618,20 +614,25 @@ def theoretical_energy_integral(parameter_values, n_i, n_f, p_i, p_f, points=100
     E
         The total energy of the cell in Wh
     """
-    n_vals = np.linspace(n_i, n_f, num=points)
-    p_vals = np.linspace(p_i, p_f, num=points)
+    x_0 = inputs["x_0"]
+    y_0 = inputs["y_0"]
+    x_100 = inputs["x_100"]
+    y_100 = inputs["y_100"]
+    Q_p = inputs["Q_p"]
+    x_vals = np.linspace(x_100, x_0, num=points)
+    y_vals = np.linspace(y_100, y_0, num=points)
     # Calculate OCV at each stoichiometry
     param = pybamm.LithiumIonParameters()
     T = param.T_amb(0)
-    Vs = np.empty(n_vals.shape)
-    for i in range(n_vals.size):
+    Vs = np.empty(x_vals.shape)
+    for i in range(x_vals.size):
         Vs[i] = (
-            parameter_values.evaluate(param.p.prim.U(p_vals[i], T)).item()
-            - parameter_values.evaluate(param.n.prim.U(n_vals[i], T)).item()
+            parameter_values.evaluate(param.p.prim.U(y_vals[i], T)).item()
+            - parameter_values.evaluate(param.n.prim.U(x_vals[i], T)).item()
         )
     # Calculate dQ
-    Q_p = parameter_values.evaluate(param.p.prim.Q_init) * (p_f - p_i)
-    dQ = Q_p / (points - 1)
+    Q = Q_p * (y_0 - y_100)
+    dQ = Q / (points - 1)
     # Integrate and convert to W-h
     E = np.trapz(Vs, dx=dQ)
     return E
@@ -663,7 +664,10 @@ def calculate_theoretical_energy(
     # Get initial and final stoichiometric values.
     x_100, y_100 = get_initial_stoichiometries(initial_soc, parameter_values)
     x_0, y_0 = get_initial_stoichiometries(final_soc, parameter_values)
+    Q_p = parameter_values.evaluate(pybamm.LithiumIonParameters().p.prim.Q_init)
     E = theoretical_energy_integral(
-        parameter_values, x_100, x_0, y_100, y_0, points=points
+        parameter_values,
+        {"x_100": x_100, "x_0": x_0, "y_100": y_100, "y_0": y_0, "Q_p": Q_p},
+        points=points,
     )
     return E

pybamm/models/submodels/interface/kinetics/base_kinetics.py

@@ -57,6 +57,7 @@ def get_coupled_variables(self, variables):
         reaction_name = self.reaction_name
         phase_name = self.phase_name
 
+        # Get surface potential difference
         if self.reaction == "lithium metal plating":  # li metal electrode (half-cell)
             delta_phi = variables[
                 "Lithium metal interface surface potential difference [V]"
@@ -79,7 +80,8 @@ def get_coupled_variables(self, variables):
 
         # Get exchange-current density
         j0 = self._get_exchange_current_density(variables)
-        # Get open-circuit potential variables and reaction overpotential
+
+        # Get open-circuit potential
         if (
             domain_options["particle size"] == "distribution"
             and self.options.electrode_types[domain] == "porous"
@@ -92,9 +94,14 @@ def get_coupled_variables(self, variables):
             ocp = variables[
                 f"{Domain} electrode {reaction_name}open-circuit potential [V]"
             ]
-        # If ocp was broadcast, take only the orphan.
-        if isinstance(ocp, pybamm.Broadcast):
+        # If ocp was broadcast, and delta_phi's secondary domain is "current collector",
+        # then take only the orphan.
+        if isinstance(ocp, pybamm.Broadcast) and delta_phi.domains["secondary"] == [
+            "current collector"
+        ]:
             ocp = ocp.orphans[0]
+
+        # Get reaction overpotential
         eta_r = delta_phi - ocp
 
         # Get average interfacial current density
@@ -130,8 +137,15 @@ def get_coupled_variables(self, variables):
             eta_sei = pybamm.Scalar(0)
         eta_r += eta_sei
 
+        # Broadcast j0 to match eta_r's domain, if necessary
+        if j0.secondary_domain == ["current collector"] and eta_r.secondary_domain == [
+            f"{domain} electrode"
+        ]:
+            j0 = pybamm.SecondaryBroadcast(j0, [f"{domain} electrode"])
+
         # Get number of electrons in reaction
         ne = self._get_number_of_electrons_in_reaction()
+
         # Get kinetics. Note: T and u must have the same domain as j0 and eta_r
         if self.options.electrode_types[domain] == "planar":
             T = variables["X-averaged cell temperature [K]"]

tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_basic_models.py

@@ -10,17 +10,18 @@
 
 class TestCompareBasicModels(TestCase):
     def test_compare_dfns(self):
+        parameter_values = pybamm.ParameterValues("Ecker2015")
         basic_dfn = pybamm.lithium_ion.BasicDFN()
         dfn = pybamm.lithium_ion.DFN()
 
         # Solve basic DFN
-        basic_sim = pybamm.Simulation(basic_dfn)
+        basic_sim = pybamm.Simulation(basic_dfn, parameter_values=parameter_values)
         t_eval = np.linspace(0, 3600)
         basic_sim.solve(t_eval)
         basic_sol = basic_sim.solution
 
         # Solve main DFN
-        sim = pybamm.Simulation(dfn)
+        sim = pybamm.Simulation(dfn, parameter_values=parameter_values)
         t_eval = np.linspace(0, 3600)
         sim.solve(t_eval)
         sol = sim.solution
@@ -64,17 +65,18 @@ def test_compare_dfns_composite(self):
             )
 
     def test_compare_spms(self):
+        parameter_values = pybamm.ParameterValues("Ecker2015")
         basic_spm = pybamm.lithium_ion.BasicSPM()
         spm = pybamm.lithium_ion.SPM()
 
         # Solve basic SPM
-        basic_sim = pybamm.Simulation(basic_spm)
+        basic_sim = pybamm.Simulation(basic_spm, parameter_values=parameter_values)
         t_eval = np.linspace(0, 3600)
         basic_sim.solve(t_eval)
         basic_sol = basic_sim.solution
 
         # Solve main SPM
-        sim = pybamm.Simulation(spm)
+        sim = pybamm.Simulation(spm, parameter_values=parameter_values)
         t_eval = np.linspace(0, 3600)
         sim.solve(t_eval)
         sol = sim.solution
@@ -84,7 +86,7 @@ def test_compare_spms(self):
         # Compare variables
         for name in basic_spm.variables:
             np.testing.assert_allclose(
-                basic_sol[name].entries, sol[name].entries, rtol=1e-5
+                basic_sol[name].entries, sol[name].entries, rtol=1e-4
             )
 
 

tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py

@@ -367,3 +367,7 @@ def test_well_posed_current_sigmoid_exchange_current(self):
     def test_well_posed_current_sigmoid_diffusivity(self):
         options = {"diffusivity": "current sigmoid"}
         self.check_well_posedness(options)
+
+    def test_well_posed_psd(self):
+        options = {"particle size": "distribution", "surface form": "algebraic"}
+        self.check_well_posedness(options)

tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py

@@ -36,12 +36,12 @@ def test_known_solution(self):
                 self.assertAlmostEqual(sol[key], sol_split[key].data[0], places=5)
             else:
                 # theoretical_energy is not present in sol_split
-                x_0 = sol_split["x_0"].data[0]
-                y_0 = sol_split["y_0"].data[0]
-                x_100 = sol_split["x_100"].data[0]
-                y_100 = sol_split["y_100"].data[0]
+                inputs = {
+                    k: sol_split[k].data[0]
+                    for k in ["x_0", "y_0", "x_100", "y_100", "Q_p"]
+                }
                 energy = pybamm.lithium_ion.electrode_soh.theoretical_energy_integral(
-                    parameter_values, x_100, x_0, y_100, y_0
+                    parameter_values, inputs
                 )
                 self.assertAlmostEqual(sol[key], energy, places=5)