Merge branch 'develop' into issue-1429-simplifications

.gitignore

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

CHANGELOG.md

@@ -1,11 +1,16 @@
 # [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
 
-- Perform more automatic simplifications of the expression tree ([#1449](https://github.com/pybamm-team/PyBaMM/pull/1449))
-- Reduce time taken to hash a sparse `Matrix` object ([#1449](https://github.com/pybamm-team/PyBaMM/pull/1449))
+-   Improved how the CasADi solver's "safe" mode finds events ([#1450](https://github.com/pybamm-team/PyBaMM/pull/1450))
+-   Perform more automatic simplifications of the expression tree ([#1449](https://github.com/pybamm-team/PyBaMM/pull/1449))
+-   Reduce time taken to hash a sparse `Matrix` object ([#1449](https://github.com/pybamm-team/PyBaMM/pull/1449))
 
-# [v0.4.0](https://github.com/pybamm-team/PyBaMM/tree/v0.4.0) - 2021-03-27
+# [v0.4.0](https://github.com/pybamm-team/PyBaMM/tree/v0.4.0) - 2021-03-28
 
 This release introduces:
 - several new models, including reversible and irreversible plating submodels, submodels for loss of active material, Yang et al.'s (2017) coupled SEI/plating/pore clogging model, and the Newman-Tobias model

docs/conf.py

@@ -22,7 +22,7 @@
 # -- Project information -----------------------------------------------------
 
 project = "PyBaMM"
-copyright = "2020, The PyBaMM Team"
+copyright = "2021, The PyBaMM Team"
 author = "The PyBaMM Team"
 
 # The short X.Y version
@@ -66,8 +66,8 @@
 # source_suffix = ['.rst', '.md']
 source_suffix = ".rst"
 
-# The main toctree document.
-main_doc = "index"
+# The master toctree document.
+master_doc = "index"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -154,16 +154,14 @@
 # Grouping the document tree into LaTeX files. List of tuples
 # (source start file, target name, title,
 #  author, documentclass [howto, manual, or own class]).
-latex_documents = [
-    (main_doc, "PyBaMM.tex", "PyBaMM Documentation", "Valentin Sulzer", "manual")
-]
+latex_documents = [(master_doc, "PyBaMM.tex", "PyBaMM Documentation", author, "manual")]
 
 
 # -- Options for manual page output ------------------------------------------
 
 # One entry per manual page. List of tuples
 # (source start file, name, description, authors, manual section).
-man_pages = [(main_doc, "pybamm", "PyBaMM Documentation", [author], 1)]
+man_pages = [(master_doc, "pybamm", "PyBaMM Documentation", [author], 1)]
 
 
 # -- Options for Texinfo output ----------------------------------------------
@@ -173,7 +171,7 @@
 #  dir menu entry, description, category)
 texinfo_documents = [
     (
-        main_doc,
+        master_doc,
         "PyBaMM",
         "PyBaMM Documentation",
         author,

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()