Merge pull request #808 from pybamm-team/issue-759-non-smooth-dae

Issue 759 non smooth dae

docs/source/models/base_models/event.rst

@@ -0,0 +1,10 @@
+Event
+=====
+
+.. autoclass:: pybamm.Event
+    :members:
+
+.. autoclass:: pybamm.EventType
+    :members:
+
+

docs/source/models/base_models/index.rst

@@ -6,3 +6,4 @@ Base Models
 
   base_model
   base_battery_model
+  event

pybamm/__init__.py

@@ -148,6 +148,8 @@ def version(formatted=False):
 #
 from .models.base_model import BaseModel
 from .models import standard_variables
+from .models.event import Event
+from .models.event import EventType
 
 # Battery models
 from .models.full_battery_models.base_battery_model import BaseBatteryModel

pybamm/discretisations/discretisation.py

@@ -173,11 +173,16 @@ def process_model(self, model, inplace=True, check_model=True):
         model_disc.algebraic, model_disc.concatenated_algebraic = alg, concat_alg
 
         # Process events
-        processed_events = {}
+        processed_events = []
         pybamm.logger.info("Discretise events for {}".format(model.name))
-        for event, equation in model.events.items():
-            pybamm.logger.debug("Discretise event '{}'".format(event))
-            processed_events[event] = self.process_symbol(equation)
+        for event in model.events:
+            pybamm.logger.debug("Discretise event '{}'".format(event.name))
+            processed_event = pybamm.Event(
+                event.name,
+                self.process_symbol(event.expression),
+                event.event_type
+            )
+            processed_events.append(processed_event)
         model_disc.events = processed_events
 
         # Create mass matrix

pybamm/expression_tree/binary_operators.py

@@ -644,6 +644,16 @@ def inner(left, right):
 class Heaviside(BinaryOperator):
     """A node in the expression tree representing a heaviside step function.
 
+    Adding this operation to the rhs or algebraic equations in a model can often cause a
+    discontinuity in the solution. For the specific cases listed below, this will be
+    automatically handled by the solver. In the general case, you can explicitly tell
+    the solver of discontinuities by adding a :class:`Event` object with
+    :class:`EventType` DISCONTINUITY to the model's list of events.
+
+    In the case where the Heaviside function is of the form `pybamm.t < x`, `pybamm.t <=
+    x`, `x < pybamm.t`, or `x <= pybamm.t`, where `x` is any constant equation, this
+    DISCONTINUITY event will automatically be added by the solver.
+
     **Extends:** :class:`BinaryOperator`
     """
 

pybamm/models/base_model.py

@@ -46,9 +46,10 @@ class BaseModel(object):
     variables: dict
         A dictionary that maps strings to expressions that represent
         the useful variables
-    events: list
-        A list of events that should cause the solver to terminate (e.g. concentration
-        goes negative)
+    events: list of :class:`pybamm.Event`
+        A list of events. Each event can either cause the solver to terminate
+        (e.g. concentration goes negative), or be used to inform the solver of the
+        existance of a discontinuity (e.g. discontinuity in the input current)
     concatenated_rhs : :class:`pybamm.Concatenation`
         After discretisation, contains the expressions representing the rhs equations
         concatenated into a single expression
@@ -105,7 +106,7 @@ def __init__(self, name="Unnamed model"):
         self._initial_conditions = {}
         self._boundary_conditions = {}
         self._variables = pybamm.FuzzyDict()
-        self._events = {}
+        self._events = []
         self._concatenated_rhs = None
         self._concatenated_algebraic = None
         self._concatenated_initial_conditions = None
@@ -337,7 +338,7 @@ def update(self, *submodels):
                 self._boundary_conditions, submodel.boundary_conditions
             )
             self.variables.update(submodel.variables)  # keys are strings so no check
-            self._events.update(submodel.events)
+            self._events += submodel.events
 
     def check_and_combine_dict(self, dict1, dict2):
         # check that the key ids are distinct

pybamm/models/event.py

@@ -0,0 +1,66 @@
+from enum import Enum
+
+
+class EventType(Enum):
+    """
+    Defines the type of event, see :class:`pybamm.Event`
+
+    TERMINATION indicates an event that will terminate the solver, the expression should
+    return 0 when the event is triggered
+
+    DISCONTINUITY indicates an expected discontinuity in the solution, the expression
+    should return the time that the discontinuity occurs. The solver will integrate up
+    to the discontinuity and then restart just after the discontinuity.
+
+    """
+    TERMINATION = 0
+    DISCONTINUITY = 1
+
+
+class Event:
+    """
+
+    Defines an event for use within a pybamm model
+
+    Attributes
+    ----------
+
+    name: str
+        A string giving the name of the event
+    event_type: :class:`pybamm.EventType`
+        An enum defining the type of event
+    expression: :class:`pybamm.Symbol`
+        An expression that defines when the event occurs
+
+
+    """
+
+    def __init__(self, name, expression, event_type=EventType.TERMINATION):
+        self._name = name
+        self._expression = expression
+        self._event_type = event_type
+
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+        """
+        Acts as a drop-in replacement for :func:`pybamm.Symbol.evaluate`
+        """
+        return self._expression.evaluate(t, y, u, known_evals)
+
+    def __str__(self):
+        return self._name
+
+    @property
+    def name(self):
+        return self._name
+
+    @property
+    def expression(self):
+        return self._expression
+
+    @expression.setter
+    def expression(self, value):
+        self._expression = value
+
+    @property
+    def event_type(self):
+        return self._event_type

pybamm/models/full_battery_models/base_battery_model.py

@@ -767,8 +767,16 @@ def set_voltage_variables(self):
 
         # Cut-off voltage
         voltage = self.variables["Terminal voltage"]
-        self.events["Minimum voltage"] = voltage - self.param.voltage_low_cut
-        self.events["Maximum voltage"] = voltage - self.param.voltage_high_cut
+        self.events.append(pybamm.Event(
+            "Minimum voltage",
+            voltage - self.param.voltage_low_cut,
+            pybamm.EventType.TERMINATION
+        ))
+        self.events.append(pybamm.Event(
+            "Maximum voltage",
+            voltage - self.param.voltage_high_cut,
+            pybamm.EventType.TERMINATION
+        ))
 
         # Power
         I_dim = self.variables["Current [A]"]

pybamm/models/full_battery_models/lithium_ion/basic_dfn.py

@@ -180,20 +180,16 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         self.initial_conditions[c_s_n] = param.c_n_init
         self.initial_conditions[c_s_p] = param.c_p_init
         # Events specify points at which a solution should terminate
-        self.events.update(
-            {
-                "Minimum negative particle surface concentration": (
-                    pybamm.min(c_s_surf_n) - 0.01
-                ),
-                "Maximum negative particle surface concentration": (1 - 0.01)
-                - pybamm.max(c_s_surf_n),
-                "Minimum positive particle surface concentration": (
-                    pybamm.min(c_s_surf_p) - 0.01
-                ),
-                "Maximum positive particle surface concentration": (1 - 0.01)
-                - pybamm.max(c_s_surf_p),
-            }
-        )
+        self.events += [
+            pybamm.Event("Minimum negative particle surface concentration",
+                         pybamm.min(c_s_surf_n) - 0.01),
+            pybamm.Event("Maximum negative particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_n)),
+            pybamm.Event("Minimum positive particle surface concentration",
+                         pybamm.min(c_s_surf_p) - 0.01),
+            pybamm.Event("Maximum positive particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_p)),
+        ]
         ######################
         # Current in the solid
         ######################
@@ -245,7 +241,8 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
             "right": (pybamm.Scalar(0), "Neumann"),
         }
         self.initial_conditions[c_e] = param.c_e_init
-        self.events["Zero electrolyte concentration cut-off"] = pybamm.min(c_e) - 0.002
+        self.events.append(pybamm.Event("Zero electrolyte concentration cut-off",
+                                        pybamm.min(c_e) - 0.002))
 
         ######################
         # (Some) variables
@@ -263,8 +260,10 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
             "Positive electrode potential": phi_s_p,
             "Terminal voltage": voltage,
         }
-        self.events["Minimum voltage"] = voltage - param.voltage_low_cut
-        self.events["Maximum voltage"] = voltage - param.voltage_high_cut
+        self.events += [
+            pybamm.Event("Minimum voltage", voltage - param.voltage_low_cut),
+            pybamm.Event("Maximum voltage", voltage - param.voltage_high_cut),
+        ]
 
     @property
     def default_geometry(self):

pybamm/models/full_battery_models/lithium_ion/basic_spm.py

@@ -97,20 +97,16 @@ def __init__(self, name="Single Particle Model"):
         c_s_surf_n = pybamm.surf(c_s_n)
         c_s_surf_p = pybamm.surf(c_s_p)
         # Events specify points at which a solution should terminate
-        self.events.update(
-            {
-                "Minimum negative particle surface concentration": (
-                    pybamm.min(c_s_surf_n) - 0.01
-                ),
-                "Maximum negative particle surface concentration": (1 - 0.01)
-                - pybamm.max(c_s_surf_n),
-                "Minimum positive particle surface concentration": (
-                    pybamm.min(c_s_surf_p) - 0.01
-                ),
-                "Maximum positive particle surface concentration": (1 - 0.01)
-                - pybamm.max(c_s_surf_p),
-            }
-        )
+        self.events += [
+            pybamm.Event("Minimum negative particle surface concentration",
+                         pybamm.min(c_s_surf_n) - 0.01),
+            pybamm.Event("Maximum negative particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_n)),
+            pybamm.Event("Minimum positive particle surface concentration",
+                         pybamm.min(c_s_surf_p) - 0.01),
+            pybamm.Event("Maximum positive particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_p)),
+        ]
 
         # Note that the SPM does not have any algebraic equations, so the `algebraic`
         # dictionary remains empty
@@ -164,8 +160,10 @@ def __init__(self, name="Single Particle Model"):
             ),
             "Terminal voltage": V,
         }
-        self.events["Minimum voltage"] = V - param.voltage_low_cut
-        self.events["Maximum voltage"] = V - param.voltage_high_cut
+        self.events += [
+            pybamm.Event("Minimum voltage", V - param.voltage_low_cut),
+            pybamm.Event("Maximum voltage", V - param.voltage_high_cut),
+        ]
 
     @property
     def default_geometry(self):

pybamm/models/submodels/base_submodel.py

@@ -40,10 +40,10 @@ class BaseSubModel:
     variables: dict
         A dictionary that maps strings to expressions that represent
         the useful variables
-    events: dict
-        A dictionary of events that should cause the solver to terminate (e.g.
-        concentration goes negative). The keys are strings and the values are
-        symbols.
+    events: list
+        A list of events. Each event can either cause the solver to terminate
+        (e.g. concentration goes negative), or be used to inform the solver of the
+        existance of a discontinuity (e.g. discontinuity in the input current)
     """
 
     def __init__(
@@ -63,7 +63,7 @@ def __init__(
         self.boundary_conditions = {}
         self.initial_conditions = {}
         self.variables = {}
-        self.events = {}
+        self.events = []
 
         self.domain = domain
         self.set_domain_for_broadcast()

pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py

@@ -105,4 +105,8 @@ def _get_standard_flux_variables(self, N_e):
 
     def set_events(self, variables):
         c_e = variables["Electrolyte concentration"]
-        self.events["Zero electrolyte concentration cut-off"] = pybamm.min(c_e) - 0.002
+        self.events.append(pybamm.Event(
+            "Zero electrolyte concentration cut-off",
+            pybamm.min(c_e) - 0.002,
+            pybamm.EventType.TERMINATION
+        ))

pybamm/models/submodels/particle/base_particle.py

@@ -91,10 +91,15 @@ def set_events(self, variables):
         c_s_surf = variables[self.domain + " particle surface concentration"]
         tol = 0.01
 
-        self.events[
-            "Minimum " + self.domain.lower() + " particle surface concentration"
-        ] = (pybamm.min(c_s_surf) - tol)
+        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
+        ))
 
-        self.events[
-            "Maximum " + self.domain.lower() + " particle surface concentration"
-        ] = (1 - tol) - pybamm.max(c_s_surf)

pybamm/models/submodels/porosity/base_porosity.py

@@ -96,7 +96,25 @@ def _get_standard_porosity_change_variables(self, deps_dt, set_leading_order=Fal
     def set_events(self, variables):
         eps_n = variables["Negative electrode porosity"]
         eps_p = variables["Positive electrode porosity"]
-        self.events["Zero negative electrode porosity cut-off"] = pybamm.min(eps_n)
-        self.events["Max negative electrode porosity cut-off"] = pybamm.max(eps_n) - 1
-        self.events["Zero positive electrode porosity cut-off"] = pybamm.min(eps_p)
-        self.events["Max positive electrode porosity cut-off"] = pybamm.max(eps_p) - 1
+        self.events.append(pybamm.Event(
+            "Zero negative electrode porosity cut-off",
+            pybamm.min(eps_n),
+            pybamm.EventType.TERMINATION
+        ))
+        self.events.append(pybamm.Event(
+            "Max negative electrode porosity cut-off",
+            pybamm.max(eps_n) - 1,
+            pybamm.EventType.TERMINATION
+        ))
+
+        self.events.append(pybamm.Event(
+            "Zero positive electrode porosity cut-off",
+            pybamm.min(eps_p),
+            pybamm.EventType.TERMINATION
+        ))
+
+        self.events.append(pybamm.Event(
+            "Max positive electrode porosity cut-off",
+            pybamm.max(eps_p) - 1,
+            pybamm.EventType.TERMINATION
+        ))

pybamm/parameters/parameter_values.py

@@ -383,9 +383,11 @@ def process_model(self, unprocessed_model, inplace=True):
                 "Processing parameters for {!r} (variables)".format(variable)
             )
             model.variables[variable] = self.process_symbol(equation)
-        for event, equation in model.events.items():
-            pybamm.logger.debug("Processing parameters for event '{}''".format(event))
-            model.events[event] = self.process_symbol(equation)
+
+        for event in model.events:
+            pybamm.logger.debug("Processing parameters for event'{}''"
+                                .format(event.name))
+            event.expression = self.process_symbol(event.expression)
 
         pybamm.logger.info("Finish setting parameters for {}".format(model.name))