#1082 basic 'fast with events' algorithm working but needs fine tuning

examples/scripts/DFN.py

@@ -8,12 +8,14 @@
 pybamm.set_logging_level("INFO")
 
 # load model
-model = pybamm.lithium_ion.DFN()
+model = pybamm.lithium_ion.DFN()  # {"operating mode": "power"})
 # create geometry
 geometry = model.default_geometry
 
 # load parameter values and process model and geometry
 param = model.default_parameter_values
+# param.update({"Power function [W]": 3.5}, check_already_exists=False)
+param["Current function [A]"] /= 10
 param.process_geometry(geometry)
 param.process_model(model)
 
@@ -27,21 +29,23 @@
 disc.process_model(model)
 
 # solve model
-t_eval = np.linspace(0, 5000, 100)
-solver = pybamm.CasadiSolver(mode="fast", atol=1e-6, rtol=1e-3)
-solution = solver.solve(model, t_eval)
+t_eval = np.linspace(0, 5000 * 10, 100)
+solver = pybamm.CasadiSolver(mode="safe", atol=1e-6, rtol=1e-6)
+solution1 = solver.solve(model, t_eval)
+solver = pybamm.CasadiSolver(mode="fast with events", atol=1e-6, rtol=1e-6)
+solution2 = solver.solve(model, t_eval)
 
 # plot
 plot = pybamm.QuickPlot(
-    solution,
+    [solution1, solution2],
     [
-        "Negative particle concentration [mol.m-3]",
+        # "Negative particle concentration [mol.m-3]",
         "Electrolyte concentration [mol.m-3]",
-        "Positive particle concentration [mol.m-3]",
+        # "Positive particle concentration [mol.m-3]",
         "Current [A]",
         "Negative electrode potential [V]",
-        "Electrolyte potential [V]",
-        "Positive electrode potential [V]",
+        # "Electrolyte potential [V]",
+        # "Positive electrode potential [V]",
         "Terminal voltage [V]",
     ],
     time_unit="seconds",

examples/scripts/experimental_protocols/cccv.py

@@ -4,7 +4,7 @@
 import pybamm
 import matplotlib.pyplot as plt
 
-pybamm.set_logging_level("VERBOSE")
+pybamm.set_logging_level("INFO")
 experiment = pybamm.Experiment(
     [
         (
@@ -19,24 +19,24 @@
 )
 model = pybamm.lithium_ion.SPM()
 sim = pybamm.Simulation(
-    model, experiment=experiment, solver=pybamm.CasadiSolver(dt_max=1e4)
+    model, experiment=experiment, solver=pybamm.CasadiSolver("fast with events")
 )
 sim.solve()
 
 # Plot voltages from the discharge segments only
-fig, ax = plt.subplots()
-for i in range(3):
-    # Extract sub solutions
-    sol = sim.solution.cycles[i]
-    # Extract variables
-    t = sol["Time [h]"].entries
-    V = sol["Terminal voltage [V]"].entries
-    # Plot
-    ax.plot(t - t[0], V, label="Discharge {}".format(i + 1))
-    ax.set_xlabel("Time [h]")
-    ax.set_ylabel("Voltage [V]")
-    ax.set_xlim([0, 10])
-ax.legend(loc="lower left")
+# fig, ax = plt.subplots()
+# for i in range(3):
+#     # Extract sub solutions
+#     sol = sim.solution.cycles[i]
+#     # Extract variables
+#     t = sol["Time [h]"].entries
+#     V = sol["Terminal voltage [V]"].entries
+#     # Plot
+#     ax.plot(t - t[0], V, label="Discharge {}".format(i + 1))
+#     ax.set_xlabel("Time [h]")
+#     ax.set_ylabel("Voltage [V]")
+#     ax.set_xlim([0, 10])
+# ax.legend(loc="lower left")
 
 # Save time, voltage, current, discharge capacity, temperature, and electrolyte
 # concentration to csv and matlab formats

pybamm/input/parameters/lithium-ion/negative_electrodes/graphite_UMBL_Mohtat2020/graphite_ocp_PeymanMPM.py

@@ -25,10 +25,7 @@ def graphite_ocp_PeymanMPM(sto):
     return u_eq
 
 
-# if __name__ == "__main__": # pragma: no cover
-#     import matplotlib.pyplot as plt
-#     import numpy as np
-
-#     x = np.linspace(0, 1)
-#     plt.plot(x, graphite_ocp_PeymanMPM(x))
-#     plt.show()
+if __name__ == "__main__":  # pragma: no cover
+    x = pybamm.linspace(1e-10, 1 - 1e-10, 1000)
+    # pybamm.plot(x, graphite_ocp_PeymanMPM(x))
+    pybamm.plot(x, -1e-8 * pybamm.log(x / (1 - x)))

pybamm/input/parameters/lithium-ion/positive_electrodes/NMC_UMBL_Mohtat2020/NMC_ocp_PeymanMPM.py

@@ -30,9 +30,6 @@ def NMC_ocp_PeymanMPM(sto):
     return u_eq
 
 
-# if __name__ == "__main__": # pragma: no cover
-#     import matplotlib.pyplot as plt
-
-#     x = np.linspace(0, 1)
-#     plt.plot(x, NMC_ocp_PeymanMPM(x))
-#     plt.show()
+# if __name__ == "__main__":  # pragma: no cover
+#     x = pybamm.linspace(0, 1)
+#     pybamm.plot(x, NMC_ocp_PeymanMPM(x))

pybamm/models/full_battery_models/base_battery_model.py

@@ -899,13 +899,13 @@ def set_voltage_variables(self):
         )
 
         # Cut-off voltage
-        # self.events.append(
-        #     pybamm.Event(
-        #         "Minimum voltage",
-        #         V - self.param.voltage_low_cut,
-        #         pybamm.EventType.TERMINATION,
-        #     )
-        # )
+        self.events.append(
+            pybamm.Event(
+                "Minimum voltage",
+                V - self.param.voltage_low_cut,
+                pybamm.EventType.TERMINATION,
+            )
+        )
         self.events.append(
             pybamm.Event(
                 "Maximum voltage",

pybamm/models/submodels/interface/base_interface.py

@@ -72,6 +72,11 @@ def _get_exchange_current_density(self, variables):
                 c_s_surf = c_s_surf.orphans[0]
                 c_e = c_e.orphans[0]
                 T = T.orphans[0]
+
+            tol = 1e-8
+            c_e = pybamm.maximum(tol, c_e)
+            c_s_surf = pybamm.maximum(tol, pybamm.minimum(c_s_surf, 1 - tol))
+
             if self.domain == "Negative":
                 j0 = self.param.j0_n(c_e, c_s_surf, T) / self.param.C_r_n
             elif self.domain == "Positive":

pybamm/parameters/lithium_ion_parameters.py

@@ -367,12 +367,20 @@ def U_n_dimensional(self, sto, T):
         """Dimensional open-circuit potential in the negative electrode [V]"""
         inputs = {"Negative particle stoichiometry": sto}
         u_ref = pybamm.FunctionParameter("Negative 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
+        # see #1435
+        u_ref -= 1e-6 * pybamm.log(sto / (1 - sto))
         return u_ref + (T - self.T_ref) * self.dUdT_n_dimensional(sto)
 
     def U_p_dimensional(self, sto, T):
         """Dimensional open-circuit potential in the positive electrode [V]"""
         inputs = {"Positive particle stoichiometry": sto}
         u_ref = pybamm.FunctionParameter("Positive 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
+        # see #1435
+        u_ref -= 1e-6 * pybamm.log(sto / (1 - sto))
         return u_ref + (T - self.T_ref) * self.dUdT_p_dimensional(sto)
 
     def dUdT_n_dimensional(self, sto):
@@ -880,13 +888,11 @@ def R_p(self, x):
         return self.R_p_dimensional(x_dim) / self.R_p_typ
 
     def c_n_init(self, x):
-        """Dimensionless initial concentration as a function of dimensionless position x
-        """
+        """Dimensionless initial concentration as a function of dimensionless position x"""
         return self.c_n_init_dimensional(x) / self.c_n_max
 
     def c_p_init(self, x):
-        """Dimensionless initial concentration as a function of dimensionless position x
-        """
+        """Dimensionless initial concentration as a function of dimensionless position x"""
         return self.c_p_init_dimensional(x) / self.c_p_max
 
     def rho(self, T):

pybamm/simulation.py

@@ -172,12 +172,14 @@ def set_up_experiment(self, model, experiment):
         for op, events in zip(experiment.operating_conditions, experiment.events):
             if op[1] in ["A", "C"]:
                 # Update inputs for constant current
+                capacity = self._parameter_values["Nominal cell capacity [A.h]"]
                 if op[1] == "A":
                     I = op[0]
+                    Crate = I / capacity
                 else:
                     # Scale C-rate with capacity to obtain current
-                    capacity = self._parameter_values["Nominal cell capacity [A.h]"]
-                    I = op[0] * capacity
+                    Crate = op[0]
+                    I = Crate * capacity
                 operating_inputs = {
                     "Current switch": 1,
                     "Voltage switch": 0,
@@ -238,8 +240,13 @@ def set_up_experiment(self, model, experiment):
             # Add time to the experiment times
             dt = op[2]
             if dt is None:
-                # max simulation time: 1 week
-                dt = 7 * 24 * 3600
+                if op[1] in ["A", "C"]:
+                    # Current control: max simulation time: 3 * max simulation time
+                    # based on C-rate
+                    dt = 3 / abs(Crate)
+                else:
+                    # max simulation time: 1 week
+                    dt = 7 * 24 * 3600
             self._experiment_times.append(dt)
 
         # Set up model for experiment

pybamm/solvers/base_solver.py

@@ -378,15 +378,17 @@ def report(string):
                     # events' mode of the casadi solver
                     # see #1082
                     k = 20
-                    if (
-                        event.name in ["Minimum voltage", "Maximum voltage"]
-                        or "[experiment]" in event.name
-                    ):
+                    if "voltage" in event.name.lower():
                         init_sign = float(
                             np.sign(event_eval(0, model.y0, inputs_stacked))
                         )
-                        event_sigmoid = pybamm.sigmoid(
-                            0, init_sign * event.expression, k
+                        # We create a sigmoid for each event which will multiply the
+                        # rhs. Doing * 2 - 1 ensures that when the event is crossed,
+                        # the sigmoid is zero. Hence the rhs is zero and the solution
+                        # stays constant for the rest of the simulation period
+                        # We can then cut off the part after the event was crossed
+                        event_sigmoid = (
+                            pybamm.sigmoid(0, init_sign * event.expression, k) * 2 - 1
                         )
                         event_casadi = process(
                             event_sigmoid, "event", use_jacobian=False