Merge pull request #1450 from pybamm-team/issue-1082-events

Issue 1082 events

.gitignore

@@ -44,6 +44,7 @@ out/
 config.py
 matplotlibrc
 *.pickle
+*.sav
 
 # ideas
 ideas/

CHANGELOG.md

@@ -1,3 +1,13 @@
+# [Unreleased](https://github.com/pybamm-team/PyBaMM)
+
+## Features
+
+-   Added "fast with events" mode for the CasADi solver, which solves a model and finds events more efficiently than "safe" mode. As of PR #1450 this feature is still being tested and "safe" mode remains the default ([#1450](https://github.com/pybamm-team/PyBaMM/pull/1450))
+
+## Optimizations
+
+-   Improved how the CasADi solver's "safe" mode finds events ([#1450](https://github.com/pybamm-team/PyBaMM/pull/1450))
+
 # [v0.4.0](https://github.com/pybamm-team/PyBaMM/tree/v0.4.0) - 2021-03-28
 
 This release introduces:

examples/scripts/DFN.py

@@ -9,7 +9,6 @@
 
 # load model
 model = pybamm.lithium_ion.DFN()
-
 # create geometry
 geometry = model.default_geometry
 

examples/scripts/experimental_protocols/cccv.py

@@ -4,21 +4,24 @@
 import pybamm
 import matplotlib.pyplot as plt
 
-pybamm.set_logging_level("INFO")
+pybamm.set_logging_level("NOTICE")
 experiment = pybamm.Experiment(
     [
         (
             "Discharge at C/5 for 10 hours or until 3.3 V",
             "Rest for 1 hour",
             "Charge at 1 A until 4.1 V",
-            "Hold at 4.1 V until 50 mA",
+            "Hold at 4.1 V until 10 mA",
             "Rest for 1 hour",
         ),
     ]
-    * 3,
+    * 3
 )
 model = pybamm.lithium_ion.DFN()
-sim = pybamm.Simulation(model, experiment=experiment, solver=pybamm.CasadiSolver())
+
+sim = pybamm.Simulation(
+    model, experiment=experiment, solver=pybamm.CasadiSolver("safe")
+)
 sim.solve()
 
 # Plot voltages from the discharge segments only

pybamm/input/parameters/lithium-ion/experiments/1C_discharge_from_full_Marquis2019/parameters.csv

@@ -15,7 +15,7 @@ Total heat transfer coefficient [W.m-2.K-1],10,,
 Number of electrodes connected in parallel to make a cell,1,,
 Number of cells connected in series to make a battery,1,,
 Lower voltage cut-off [V],3.105,,
-Upper voltage cut-off [V],4.7,,
+Upper voltage cut-off [V],4.2,,
 ,,,
 # Initial conditions
 Initial concentration in negative electrode [mol.m-3],19986.609595075,Scott Moura FastDFN,

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/event.py

@@ -17,6 +17,7 @@ class EventType(Enum):
     TERMINATION = 0
     DISCONTINUITY = 1
     INTERPOLANT_EXTRAPOLATION = 2
+    SWITCH = 3
 
 
 class Event:

pybamm/models/full_battery_models/base_battery_model.py

@@ -916,6 +916,26 @@ def set_voltage_variables(self):
             )
         )
 
+        # Cut-off open-circuit voltage (for event switch with casadi 'fast with events'
+        # mode)
+        # A tolerance of 1 is sufficiently small since the dimensionless voltage is
+        # scaled with the thermal voltage (0.025V) and hence has a range of around 60
+        tol = 1
+        self.events.append(
+            pybamm.Event(
+                "Minimum voltage switch",
+                V - (self.param.voltage_low_cut - tol),
+                pybamm.EventType.SWITCH,
+            )
+        )
+        self.events.append(
+            pybamm.Event(
+                "Maximum voltage switch",
+                V - (self.param.voltage_high_cut + tol),
+                pybamm.EventType.SWITCH,
+            )
+        )
+
         # Power
         I_dim = self.variables["Current [A]"]
         self.variables.update({"Terminal power [W]": I_dim * V_dim})

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/models/submodels/particle/base_particle.py

@@ -143,23 +143,3 @@ def _get_standard_flux_variables(self, N_s, N_s_xav):
         }
 
         return variables
-
-    def set_events(self, variables):
-        c_s_surf = variables[self.domain + " particle surface concentration"]
-        tol = 1e-4
-
-        self.events.append(
-            pybamm.Event(
-                "Minumum " + self.domain.lower() + " particle surface concentration",
-                pybamm.min(c_s_surf) - tol,
-                pybamm.EventType.TERMINATION,
-            )
-        )
-
-        self.events.append(
-            pybamm.Event(
-                "Maximum " + self.domain.lower() + " particle surface concentration",
-                (1 - tol) - pybamm.max(c_s_surf),
-                pybamm.EventType.TERMINATION,
-            )
-        )

pybamm/parameters/lithium_ion_parameters.py

@@ -381,12 +381,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):
@@ -918,13 +926,13 @@ def R_p(self, x):
 
     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
 

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) * 3600  # seconds
+                else:
+                    # max simulation time: 1 day
+                    dt = 24 * 3600  # seconds
             self._experiment_times.append(dt)
 
         # Set up model for experiment
@@ -415,13 +422,14 @@ def set_up_model_for_experiment_new(self, model):
                         )
                     )
 
-                # Remove upper and lower voltage cut-offs that are *not* part of the
+                # Keep the min and max voltages as safeguards but add some tolerances
+                # so that they are not triggered before the voltage limits in the
                 # experiment
-                new_model.events = [
-                    event
-                    for event in new_model.events
-                    if event.name not in ["Minimum voltage", "Maximum voltage"]
-                ]
+                for event in new_model.events:
+                    if event.name == "Minimum voltage":
+                        event._expression += 1
+                    elif event.name == "Maximum voltage":
+                        event._expression -= 1
 
                 # Update parameter values
                 new_parameter_values = self.parameter_values.copy()

pybamm/solvers/base_solver.py

@@ -131,7 +131,7 @@ def set_up(self, model, inputs=None, t_eval=None):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
-        inputs_dict : dict, optional
+        inputs : dict, optional
             Any input parameters to pass to the model when solving
         t_eval : numeric type, optional
             The times (in seconds) at which to compute the solution
@@ -355,40 +355,68 @@ def report(string):
         algebraic, algebraic_eval, jac_algebraic = process(
             model.concatenated_algebraic, "algebraic"
         )
-        terminate_events_eval = [
-            process(event.expression, "event", use_jacobian=False)[1]
-            for event in model.events
-            if event.event_type == pybamm.EventType.TERMINATION
-        ]
 
-        interpolant_extrapolation_events_eval = [
-            process(event.expression, "event", use_jacobian=False)[1]
-            for event in model.events
-            if event.event_type == pybamm.EventType.INTERPOLANT_EXTRAPOLATION
-        ]
+        # Calculate initial conditions
+        model.y0 = init_eval(inputs)
 
-        # discontinuity events are evaluated before the solver is called, so don't need
-        # to process them
-        discontinuity_events_eval = [
-            event
-            for event in model.events
-            if event.event_type == pybamm.EventType.DISCONTINUITY
-        ]
+        casadi_terminate_events = []
+        terminate_events_eval = []
+        interpolant_extrapolation_events_eval = []
+        discontinuity_events_eval = []
+        for n, event in enumerate(model.events):
+            if event.event_type == pybamm.EventType.DISCONTINUITY:
+                # discontinuity events are evaluated before the solver is called,
+                # so don't need to process them
+                discontinuity_events_eval.append(event)
+            elif event.event_type == pybamm.EventType.SWITCH:
+                if (
+                    isinstance(self, pybamm.CasadiSolver)
+                    and self.mode == "fast with events"
+                    and model.algebraic != {}
+                ):
+                    # Save some events to casadi_terminate_events for the 'fast with
+                    # events' mode of the casadi solver
+                    # We only need to do this if the model is a DAE model
+                    # see #1082
+                    k = 20
+                    init_sign = float(
+                        np.sign(event.evaluate(0, model.y0.full(), inputs=inputs))
+                    )
+                    # 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, f"event_{n}", use_jacobian=False
+                    )[0]
+                    # use the actual casadi object as this will go into the rhs
+                    casadi_terminate_events.append(event_casadi)
+            else:
+                # use the function call
+                event_eval = process(
+                    event.expression, f"event_{n}", use_jacobian=False
+                )[1]
+                if event.event_type == pybamm.EventType.TERMINATION:
+                    terminate_events_eval.append(event_eval)
+                elif event.event_type == pybamm.EventType.INTERPOLANT_EXTRAPOLATION:
+                    interpolant_extrapolation_events_eval.append(event_eval)
 
         # Add the solver attributes
         model.init_eval = init_eval
         model.rhs_eval = rhs_eval
         model.algebraic_eval = algebraic_eval
         model.jac_algebraic_eval = jac_algebraic
+        model.casadi_terminate_events = casadi_terminate_events
         model.terminate_events_eval = terminate_events_eval
         model.discontinuity_events_eval = discontinuity_events_eval
         model.interpolant_extrapolation_events_eval = (
             interpolant_extrapolation_events_eval
         )
 
-        # Calculate initial conditions
-        model.y0 = init_eval(inputs)
-
         # Save CasADi functions for the CasADi solver
         # Note: when we pass to casadi the ode part of the problem must be in explicit
         # form so we pre-multiply by the inverse of the mass matrix
@@ -1039,6 +1067,7 @@ def get_termination_reason(self, solution, events):
                 "the solver successfully reached the end of the integration interval",
             )
         elif solution.termination == "event":
+            pybamm.logger.debug("Start post-processing events")
             # Get final event value
             final_event_values = {}
 
@@ -1052,6 +1081,14 @@ def get_termination_reason(self, solution, events):
                         )
                     )
             termination_event = min(final_event_values, key=final_event_values.get)
+
+            # Check that it's actually an event
+            if abs(final_event_values[termination_event]) > 0.1:  # pragma: no cover
+                # Hard to test this
+                raise pybamm.SolverError(
+                    "Could not determine which event was triggered "
+                    "(possibly due to NaNs)"
+                )
             # Add the event to the solution object
             solution.termination = "event: {}".format(termination_event)
             # Update t, y and inputs to include event time and state
@@ -1071,6 +1108,7 @@ def get_termination_reason(self, solution, events):
                 event_sol.integration_time = 0
                 solution = solution + event_sol
 
+            pybamm.logger.debug("Finish post-processing events")
             return solution, solution.termination
         elif solution.termination == "success":
             return solution, solution.termination

pybamm/solvers/casadi_algebraic_solver.py

@@ -123,7 +123,7 @@ def _integrate(self, model, t_eval, inputs_dict=None):
                             extrap_event_names.append(event.name[12:])
 
                     raise pybamm.SolverError(
-                        "CasADI solver failed because the following interpolation "
+                        "CasADi solver failed because the following interpolation "
                         "bounds were exceeded at the initial conditions: {}. "
                         "You may need to provide additional interpolation points "
                         "outside these bounds.".format(extrap_event_names)