Repository: tinosulzer/PyBaMM

Issue:
Implement cut-off events in solvers
Solver needs to stop when certain events occur such as concentration going below zero or voltage goes below a threshold.
These events could be implemented in the model, e.g. 
```python3
self.events = [c_e, V-V_threshold]
```
and then the solver knows to stop if any of the values in `self.events` go negative.

**Motivation**
Solutions with concentration below zero are unphysical. In particular, for lead-acid, the open-circuit voltage is a logarithmic function of concentration, so we get a run-time warning if the concentration goes below zero, e.g. [here](https://travis-ci.org/tinosulzer/PyBaMM/jobs/500904077#L2387)

Discussion:

Patch:
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index dd090985..012628a3 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -70,6 +70,11 @@ class Discretisation(object):
         # model.initial_conditions and model.boundary_conditions
         model.variables = self.process_dict(model.variables)
 
+        # Process events
+        for idx, event in enumerate(model.events):
+            model.events[idx] = self.process_symbol(event)
+        model.concatenated_events = self.concatenate(*model.events)
+
         # Check that resulting model makes sense
         self.check_model(model)
 
@@ -79,26 +84,17 @@ class Discretisation(object):
         variables : iterable of :class:`pybamm.Variables`
             The variables for which to set slices
         """
-        # Unpack symbols in variables that are concatenations of variables
-        unpacked_variables = []
-        for symbol in variables:
-            if isinstance(symbol, pybamm.Concatenation):
-                unpacked_variables.extend([var for var in symbol.children])
-            else:
-                unpacked_variables.append(symbol)
-        # Set up y_slices
-        y_slices = {variable.id: None for variable in unpacked_variables}
+        y_slices = {variable.id: None for variable in variables}
         start = 0
         end = 0
-        # Iterate through unpacked variables, adding appropriate slices to y_slices
-        for variable in unpacked_variables:
+        for variable in variables:
             # If domain is empty then variable has size 1
             if variable.domain == []:
                 end += 1
             # Otherwise, add up the size of all the domains in variable.domain
             else:
                 for dom in variable.domain:
-                    end += self._spatial_methods[dom].mesh[dom].npts_for_broadcast
+                    end += self.mesh[dom].npts
             y_slices[variable.id] = slice(start, end)
             start = end
         self._y_slices = y_slices
@@ -342,15 +338,7 @@ class Discretisation(object):
             Discretised right-hand side equations
 
         """
-        # Unpack symbols in variables that are concatenations of variables
-        unpacked_variables = []
-        for symbol in var_eqn_dict.keys():
-            if isinstance(symbol, pybamm.Concatenation):
-                unpacked_variables.extend([var for var in symbol.children])
-            else:
-                unpacked_variables.append(symbol)
-        # Check keys from the given var_eqn_dict against self._y_slices
-        ids = {v.id for v in unpacked_variables}
+        ids = {v.id for v in var_eqn_dict.keys()}
         if ids != self._y_slices.keys():
             given_variable_names = [v.name for v in var_eqn_dict.keys()]
             raise pybamm.ModelError(
@@ -359,7 +347,7 @@ class Discretisation(object):
             )
 
         equations = list(var_eqn_dict.values())
-        slices = [self._y_slices[var.id] for var in unpacked_variables]
+        slices = [self._y_slices[var.id] for var in var_eqn_dict.keys()]
 
         # sort equations according to slices
         sorted_equations = [eq for _, eq in sorted(zip(slices, equations))]
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index ae33e20c..9bb84ac7 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -66,11 +66,9 @@ class NumpyConcatenation(pybamm.Symbol):
 
     def __init__(self, *children):
         children = list(children)
-        # Turn objects that evaluate to scalars to objects that evaluate to vectors,
-        # so that we can concatenate them
         for i, child in enumerate(children):
             if child.evaluates_to_number():
-                children[i] = pybamm.NumpyBroadcast(child, [], None)
+                children[i] = pybamm.Vector(np.array([child.evaluate()]))
         super().__init__("model concatenation", children, domain=[])
 
     def evaluate(self, t=None, y=None):
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 60d0df17..86693b96 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -35,6 +35,9 @@ 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)
 
     """
 
@@ -46,6 +49,7 @@ class BaseModel(object):
         self._initial_conditions_ydot = {}
         self._boundary_conditions = {}
         self._variables = {}
+        self._events = []
         self._concatenated_rhs = None
         self._concatenated_initial_conditions = None
 
@@ -83,7 +87,7 @@ class BaseModel(object):
             "negative particle": pybamm.FiniteVolume,
             "positive particle": pybamm.FiniteVolume,
         }
-        self.default_solver = pybamm.ScipySolver(method="RK45")
+        self.default_solver = pybamm.ScikitsOdeSolver()
 
     def _set_dict(self, dict, name):
         """
@@ -165,6 +169,14 @@ class BaseModel(object):
     def variables(self, variables):
         self._variables = variables
 
+    @property
+    def events(self):
+        return self._events
+
+    @events.setter
+    def events(self, events):
+        self._events = events
+
     @property
     def concatenated_rhs(self):
         return self._concatenated_rhs
@@ -232,20 +244,13 @@ class BaseModel(object):
         vars_in_algebraic_keys = set()
         vars_in_eqns = set()
         # Get all variables ids from rhs and algebraic keys and equations
-        # For equations we look through the whole expression tree.
-        # "Variables" can be Concatenations so we also have to look in the whole
-        # expression tree
         for var, eqn in self.rhs.items():
-            vars_in_rhs_keys.update(
-                [x.id for x in var.pre_order() if isinstance(x, pybamm.Variable)]
-            )
+            vars_in_rhs_keys.add(var.id)
             vars_in_eqns.update(
                 [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
             )
         for var, eqn in self.algebraic.items():
-            vars_in_algebraic_keys.update(
-                [x.id for x in var.pre_order() if isinstance(x, pybamm.Variable)]
-            )
+            vars_in_algebraic_keys.add(var.id)
             vars_in_eqns.update(
                 [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
             )
diff --git a/pybamm/models/lead_acid/loqs.py b/pybamm/models/lead_acid/loqs.py
index 74047e0f..6a259504 100644
--- a/pybamm/models/lead_acid/loqs.py
+++ b/pybamm/models/lead_acid/loqs.py
@@ -96,6 +96,9 @@ class LOQS(pybamm.BaseModel):
             "int(epsilon_times_c)dx": (ln * epsn + ls * epss + lp * epsp) * c,
         }
 
+        # Terminate if concentration goes below zero
+        self.events = [c]
+
         # Overwrite default parameter values
         self.default_parameter_values = pybamm.ParameterValues(
             "input/parameters/lead-acid/default.csv",
diff --git a/pybamm/models/li_ion/single_particle_model.py b/pybamm/models/li_ion/single_particle_model.py
index 830b0c12..268e9e91 100644
--- a/pybamm/models/li_ion/single_particle_model.py
+++ b/pybamm/models/li_ion/single_particle_model.py
@@ -4,6 +4,7 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 import pybamm
+import numpy as np
 
 
 class SPM(pybamm.BaseModel):
@@ -95,5 +96,5 @@ class SPM(pybamm.BaseModel):
             "V": V,
         }
 
-        # Overwrite default solver
-        self.default_solver = pybamm.ScipySolver(method="BDF")
+        # Cut-off if either concentration goes negative
+        self.events = [pybamm.Function(np.min, cn), pybamm.Function(np.min, cp)]
diff --git a/pybamm/models/reaction_diffusion.py b/pybamm/models/reaction_diffusion.py
index a0784cbb..43ed1c28 100644
--- a/pybamm/models/reaction_diffusion.py
+++ b/pybamm/models/reaction_diffusion.py
@@ -38,6 +38,3 @@ class ReactionDiffusionModel(pybamm.BaseModel):
 
         # Create own model from diffusion model
         self.update(diffusion_model)
-
-        # Overwrite default solver for faster solution
-        self.default_solver = pybamm.ScipySolver(method="BDF")
diff --git a/pybamm/models/simple_ode_model.py b/pybamm/models/simple_ode_model.py
index bd98142f..7a5f540b 100644
--- a/pybamm/models/simple_ode_model.py
+++ b/pybamm/models/simple_ode_model.py
@@ -54,6 +54,3 @@ class SimpleODEModel(pybamm.BaseModel):
             "b broadcasted": pybamm.Broadcast(b, whole_cell),
             "c broadcasted": pybamm.Broadcast(c, ["negative electrode", "separator"]),
         }
-
-        # Overwrite default solver for faster solution
-        # self.default_solver = pybamm.ScipySolver(method="BDF")
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 568f1b60..53d98c62 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -111,6 +111,9 @@ class ParameterValues(dict):
         for variable, equation in model.variables.items():
             model.variables[variable] = self.process_symbol(equation)
 
+        for idx, equation in enumerate(model.events):
+            model.events[idx] = self.process_symbol(equation)
+
     def process_geometry(self, geometry):
         """Assign parameter values to a geometry.
             Currently inplace, could be changed to return a new model.
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index deaf1146..0cd06db5 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -56,7 +56,7 @@ class DaeSolver(pybamm.BaseSolver):
 
         self.t, self.y = self.integrate(residuals, y0, ydot0, t_eval)
 
-    def integrate(self, residuals, y0, ydot0, t_eval):
+    def integrate(self, residuals, y0, ydot0, t_eval, events=None):
         """
         Solve a DAE model defined by residuals with initial conditions y0 and ydot0.
 
@@ -69,6 +69,9 @@ class DaeSolver(pybamm.BaseSolver):
             The initial conditions
         t_eval : numeric type
             The times at which to compute the solution
+        events : method, optional
+            A function that takes in t and y and returns conditions for the solver to
+            stop
 
         """
         raise NotImplementedError
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 04af0d30..374469d7 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -35,10 +35,12 @@ class OdeSolver(pybamm.BaseSolver):
         def dydt(t, y):
             return model.concatenated_rhs.evaluate(t, y)
 
+        events = [lambda t, y: event.evaluate(t, y) for event in model.events]
+
         y0 = model.concatenated_initial_conditions
-        self.t, self.y = self.integrate(dydt, y0, t_eval)
+        self.t, self.y = self.integrate(dydt, y0, t_eval, events=events)
 
-    def integrate(self, derivs, y0, t_eval):
+    def integrate(self, derivs, y0, t_eval, events=None):
         """
         Solve a model defined by dydt with initial conditions y0.
 
@@ -50,6 +52,9 @@ class OdeSolver(pybamm.BaseSolver):
             The initial conditions
         t_eval : numeric type
             The times at which to compute the solution
+        events : method, optional
+            A function that takes in t and y and returns conditions for the solver to
+            stop
 
         """
         raise NotImplementedError
diff --git a/pybamm/solvers/scikits_dae_solver.py b/pybamm/solvers/scikits_dae_solver.py
index bc4dc97b..6ceaab3b 100644
--- a/pybamm/solvers/scikits_dae_solver.py
+++ b/pybamm/solvers/scikits_dae_solver.py
@@ -7,6 +7,7 @@ import pybamm
 
 import numpy as np
 import importlib
+
 scikits_odes_spec = importlib.util.find_spec("scikits")
 if scikits_odes_spec is not None:
     scikits_odes_spec = importlib.util.find_spec("scikits.odes")
@@ -42,7 +43,7 @@ class ScikitsDaeSolver(pybamm.DaeSolver):
     def method(self, value):
         self._method = value
 
-    def integrate(self, residuals, y0, ydot0, t_eval):
+    def integrate(self, residuals, y0, ydot0, t_eval, events=None):
         """
         Solve a DAE model defined by residuals with initial conditions y0 and ydot_0.
 
@@ -55,16 +56,21 @@ class ScikitsDaeSolver(pybamm.DaeSolver):
             The initial conditions
         t_eval : numeric type
             The times at which to compute the solution
+        events : method, optional
+            A function that takes in t and y and returns conditions for the solver to
+            stop
 
         """
+
         def eqsres(t, y, ydot, return_residuals):
             return_residuals[:] = residuals(t, y, ydot)
 
-        extra_options = {
-            'old_api': False,
-            'rtol': self.tol,
-            'atol': self.tol,
-        }
+        def rootfn(t, y, ydot, return_root):
+            return_root[:] = [event(t, y) for event in events]
+
+        extra_options = {"old_api": False, "rtol": self.tol, "atol": self.tol}
+        if events:
+            extra_options.update({"rootfn": rootfn, "nr_rootfns": len(events)})
 
         dae_solver = scikits_odes.dae(self.method, eqsres, **extra_options)
         sol = dae_solver.solve(t_eval, y0, ydot0)
diff --git a/pybamm/solvers/scikits_ode_solver.py b/pybamm/solvers/scikits_ode_solver.py
index fb3fc047..74b96e13 100644
--- a/pybamm/solvers/scikits_ode_solver.py
+++ b/pybamm/solvers/scikits_ode_solver.py
@@ -7,6 +7,7 @@ import pybamm
 
 import numpy as np
 import importlib
+
 scikits_odes_spec = importlib.util.find_spec("scikits")
 if scikits_odes_spec is not None:
     scikits_odes_spec = importlib.util.find_spec("scikits.odes")
@@ -42,7 +43,7 @@ class ScikitsOdeSolver(pybamm.OdeSolver):
     def method(self, value):
         self._method = value
 
-    def integrate(self, derivs, y0, t_eval):
+    def integrate(self, derivs, y0, t_eval, events=None):
         """
         Solve a model defined by dydt with initial conditions y0.
 
@@ -54,17 +55,21 @@ class ScikitsOdeSolver(pybamm.OdeSolver):
             The initial conditions
         t_eval : numeric type
             The times at which to compute the solution
+        events : method, optional
+            A function that takes in t and y and returns conditions for the solver to
+            stop
 
         """
 
         def eqsydot(t, y, return_ydot):
             return_ydot[:] = derivs(t, y)
 
-        extra_options = {
-            'old_api': False,
-            'rtol': self.tol,
-            'atol': self.tol,
-        }
+        def rootfn(t, y, return_root):
+            return_root[:] = [event(t, y) for event in events]
+
+        extra_options = {"old_api": False, "rtol": self.tol, "atol": self.tol}
+        if events:
+            extra_options.update({"rootfn": rootfn, "nr_rootfns": len(events)})
 
         ode_solver = scikits_odes.ode(self.method, eqsydot, **extra_options)
         sol = ode_solver.solve(t_eval, y0)
diff --git a/pybamm/solvers/scipy_solver.py b/pybamm/solvers/scipy_solver.py
index 580dceb2..e47de6d7 100644
--- a/pybamm/solvers/scipy_solver.py
+++ b/pybamm/solvers/scipy_solver.py
@@ -32,7 +32,7 @@ class ScipySolver(pybamm.OdeSolver):
     def method(self, value):
         self._method = value
 
-    def integrate(self, derivs, y0, t_eval):
+    def integrate(self, derivs, y0, t_eval, events=None):
         """
         Solve a model defined by dydt with initial conditions y0.
 
@@ -45,6 +45,9 @@ class ScipySolver(pybamm.OdeSolver):
             The initial conditions
         t_eval : :class:`numpy.array`, size (k,)
             The times at which to compute the solution
+        events : method, optional
+            A function that takes in t and y and returns conditions for the solver to
+            stop
 
         Returns
         -------
@@ -52,6 +55,11 @@ class ScipySolver(pybamm.OdeSolver):
             An object containing the times and values of the solution, as well as
             various diagnostic messages.
         """
+        # make events terminal so that the solver stops when they are reached
+        if events:
+            for event in events:
+                event.terminal = True
+
         sol = it.solve_ivp(
             derivs,
             (t_eval[0], t_eval[-1]),
@@ -60,6 +68,7 @@ class ScipySolver(pybamm.OdeSolver):
             method=self.method,
             rtol=self.tol,
             atol=self.tol,
+            events=events,
         )
-        # TODO: implement concentration cut-off event
+
         return sol.t, sol.y


Test patch:
diff --git a/tests/test_discretisations/test_discretisation.py b/tests/test_discretisations/test_discretisation.py
index a08c90d2..7327f60b 100644
--- a/tests/test_discretisations/test_discretisation.py
+++ b/tests/test_discretisations/test_discretisation.py
@@ -38,7 +38,7 @@ class TestDiscretise(unittest.TestCase):
     def test_discretise_slicing(self):
         # create discretisation
         mesh = get_mesh_for_testing()
-        spatial_methods = {"macroscale": pybamm.FiniteVolume}
+        spatial_methods = {}
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
@@ -71,29 +71,6 @@ class TestDiscretise(unittest.TestCase):
         np.testing.assert_array_equal(y[disc._y_slices[d.id]], d_true)
         np.testing.assert_array_equal(y[disc._y_slices[jn.id]], jn_true)
 
-        # Variables with a concatenation
-        js = pybamm.Variable("js", domain=["separator"])
-        jp = pybamm.Variable("jp", domain=["positive electrode"])
-        j = pybamm.Concatenation(jn, js, jp)
-        variables = [c, d, j]
-        disc.set_variable_slices(variables)
-        self.assertEqual(
-            disc._y_slices,
-            {
-                c.id: slice(0, 100),
-                d.id: slice(100, 200),
-                jn.id: slice(200, 240),
-                js.id: slice(240, 265),
-                jp.id: slice(265, 300),
-            },
-        )
-        d_true = 4 * combined_submesh.nodes
-        jn_true = mesh["negative electrode"].nodes ** 3
-        y = np.concatenate([c_true, d_true, jn_true])
-        np.testing.assert_array_equal(y[disc._y_slices[c.id]], c_true)
-        np.testing.assert_array_equal(y[disc._y_slices[d.id]], d_true)
-        np.testing.assert_array_equal(y[disc._y_slices[jn.id]], jn_true)
-
     def test_process_symbol_base(self):
         # create discretisation
         mesh = get_mesh_for_testing()
@@ -460,36 +437,6 @@ class TestDiscretise(unittest.TestCase):
         with self.assertRaises(pybamm.ModelError):
             disc.process_model(model)
 
-    def test_process_model_concatenation(self):
-        # concatenation of variables as the key
-        cn = pybamm.Variable("c", domain=["negative electrode"])
-        cs = pybamm.Variable("c", domain=["separator"])
-        cp = pybamm.Variable("c", domain=["positive electrode"])
-        c = pybamm.Concatenation(cn, cs, cp)
-        N = pybamm.grad(c)
-        model = pybamm.BaseModel()
-        model.rhs = {c: pybamm.div(N)}
-        model.initial_conditions = {c: pybamm.Scalar(3)}
-
-        model.boundary_conditions = {N: {"left": 0, "right": 0}}
-        model.check_well_posedness()
-
-        # create discretisation
-        disc = get_discretisation_for_testing()
-        mesh = disc.mesh
-
-        combined_submesh = mesh.combine_submeshes(
-            "negative electrode", "separator", "positive electrode"
-        )
-
-        disc.process_model(model)
-        y0 = model.concatenated_initial_conditions
-        np.testing.assert_array_equal(y0, 3 * np.ones_like(combined_submesh.nodes))
-
-        # grad and div are identity operators here
-        np.testing.assert_array_equal(y0, model.concatenated_rhs.evaluate(None, y0))
-        model.check_well_posedness()
-
     def test_broadcast(self):
         whole_cell = ["negative electrode", "separator", "positive electrode"]
 
diff --git a/tests/test_expression_tree/test_concatenations.py b/tests/test_expression_tree/test_concatenations.py
index 6c55b316..e25515f9 100644
--- a/tests/test_expression_tree/test_concatenations.py
+++ b/tests/test_expression_tree/test_concatenations.py
@@ -77,13 +77,6 @@ class TestConcatenations(unittest.TestCase):
             conc.evaluate(None, y), np.concatenate([y, np.array([16]), np.array([3])])
         )
 
-        # with time
-        b = pybamm.t
-        conc = pybamm.NumpyConcatenation(a, b, c)
-        np.testing.assert_array_equal(
-            conc.evaluate(16, y), np.concatenate([y, np.array([16]), np.array([3])])
-        )
-
     def test_numpy_domain_concatenation(self):
         # create discretisation
         disc = get_discretisation_for_testing()
diff --git a/tests/test_models/standard_model_tests.py b/tests/test_models/standard_model_tests.py
index 3917eba8..6c5d4100 100644
--- a/tests/test_models/standard_model_tests.py
+++ b/tests/test_models/standard_model_tests.py
@@ -40,15 +40,16 @@ class StandardModelTest(object):
         # Model should still be well-posed after processing
         self.model.check_well_posedness()
 
-    def test_solving(self, solver=None):
+    def test_solving(self, solver=None, t_eval=None):
         # Overwrite solver if given
         if solver is not None:
             self.solver = solver
-        t_eval = np.linspace(0, 1, 100)
+        if t_eval is None:
+            t_eval = np.linspace(0, 1, 100)
         self.solver.solve(self.model, t_eval)
 
-    def test_all(self, param=None, disc=None, solver=None):
+    def test_all(self, param=None, disc=None, solver=None, t_eval=None):
         self.model.check_well_posedness()
         self.test_processing_parameters(param)
         self.test_processing_disc(disc)
-        self.test_solving(solver)
+        self.test_solving(solver, t_eval)
diff --git a/tests/test_models/test_lead_acid.py b/tests/test_models/test_lead_acid.py
index b58460d9..d90cdde5 100644
--- a/tests/test_models/test_lead_acid.py
+++ b/tests/test_models/test_lead_acid.py
@@ -20,7 +20,7 @@ class TestLeadAcidLOQS(unittest.TestCase):
     def test_solution(self):
         model = pybamm.lead_acid.LOQS()
         modeltest = tests.StandardModelTest(model)
-        modeltest.test_all()
+        modeltest.test_all(t_eval=np.linspace(0, 2))
         T, Y = modeltest.solver.t, modeltest.solver.y
 
         # check output
@@ -34,6 +34,8 @@ class TestLeadAcidLOQS(unittest.TestCase):
             model.variables["V"].evaluate(T, Y)[1:],
             model.variables["V"].evaluate(T, Y)[:-1],
         )
+        # Make sure the concentration is always positive (cut-off event working)
+        np.testing.assert_array_less(0, model.variables["c"].evaluate(T, Y))
 
 
 if __name__ == "__main__":
diff --git a/tests/test_solvers/test_scikits_solvers.py b/tests/test_solvers/test_scikits_solvers.py
index c2bbcbd6..c0ebaee6 100644
--- a/tests/test_solvers/test_scikits_solvers.py
+++ b/tests/test_solvers/test_scikits_solvers.py
@@ -37,6 +37,46 @@ class TestScikitsSolver(unittest.TestCase):
         t_sol, y_sol = solver.integrate(exponential_decay, y0, t_eval)
         np.testing.assert_allclose(y_sol[0], np.exp(-0.1 * t_sol))
 
+    def test_ode_integrate_with_event(self):
+        # Constant
+        solver = pybamm.ScikitsOdeSolver(tol=1e-8)
+
+        def constant_decay(t, y):
+            return -2 * np.ones_like(y)
+
+        def y_equal_0(t, y):
+            return y[0]
+
+        y0 = np.array([1])
+        t_eval = np.linspace(0, 1, 100)
+        t_sol, y_sol = solver.integrate(constant_decay, y0, t_eval, events=[y_equal_0])
+        np.testing.assert_allclose(1 - 2 * t_sol, y_sol[0])
+        self.assertLess(len(t_sol), len(t_eval))
+        np.testing.assert_array_less(0, y_sol[0])
+        np.testing.assert_array_less(t_sol, 0.5)
+
+        # Expnonential growth
+        solver = pybamm.ScikitsOdeSolver(tol=1e-8)
+
+        def exponential_growth(t, y):
+            return y
+
+        def y_eq_9(t, y):
+            return y - 9
+
+        def ysq_eq_7(t, y):
+            return y ** 2 - 7
+
+        y0 = np.array([1])
+        t_eval = np.linspace(0, 3, 100)
+        t_sol, y_sol = solver.integrate(
+            exponential_growth, y0, t_eval, events=[y_eq_9, ysq_eq_7]
+        )
+        self.assertLess(len(t_sol), len(t_eval))
+        np.testing.assert_allclose(np.exp(t_sol), y_sol[0], rtol=1e-4)
+        np.testing.assert_array_less(y_sol, 9)
+        np.testing.assert_array_less(y_sol ** 2, 7)
+
     def test_dae_integrate(self):
         # Constant
         solver = pybamm.ScikitsDaeSolver(tol=1e-8)
@@ -65,6 +105,56 @@ class TestScikitsSolver(unittest.TestCase):
         np.testing.assert_allclose(y_sol[0], np.exp(-0.1 * t_sol))
         np.testing.assert_allclose(y_sol[1], 2 * np.exp(-0.1 * t_sol))
 
+    def test_dae_integrate_with_event(self):
+        # Constant
+        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+
+        def constant_growth_dae(t, y, ydot):
+            return [0.5 * np.ones_like(y[0]) - ydot[0], 2 * y[0] - y[1]]
+
+        def y0_eq_2(t, y):
+            return y[0] - 2
+
+        def y1_eq_5(t, y):
+            return y[1] - 5
+
+        y0 = np.array([0, 0])
+        ydot0 = np.array([0.5, 1.0])
+        t_eval = np.linspace(0, 7, 100)
+        t_sol, y_sol = solver.integrate(
+            constant_growth_dae, y0, ydot0, t_eval, events=[y0_eq_2, y1_eq_5]
+        )
+        self.assertLess(len(t_sol), len(t_eval))
+        np.testing.assert_allclose(0.5 * t_sol, y_sol[0])
+        np.testing.assert_allclose(1.0 * t_sol, y_sol[1])
+        np.testing.assert_array_less(y_sol[0], 2)
+        np.testing.assert_array_less(y_sol[1], 5)
+
+        # Exponential decay
+        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+
+        def exponential_decay_dae(t, y, ydot):
+            return [-0.1 * y[0] - ydot[0], 2 * y[0] - y[1]]
+
+        def y0_eq_0pt9(t, y):
+            return y[0] - 0.9
+
+        def t_eq_0pt5(t, y):
+            return t - 0.5
+
+        y0 = np.array([1, 2])
+        ydot0 = np.array([-0.1, -0.2])
+        t_eval = np.linspace(0, 1, 100)
+        t_sol, y_sol = solver.integrate(
+            exponential_decay_dae, y0, ydot0, t_eval, events=[y0_eq_0pt9, t_eq_0pt5]
+        )
+
+        self.assertLess(len(t_sol), len(t_eval))
+        np.testing.assert_allclose(y_sol[0], np.exp(-0.1 * t_sol))
+        np.testing.assert_allclose(y_sol[1], 2 * np.exp(-0.1 * t_sol))
+        np.testing.assert_array_less(0.9, y_sol[0])
+        np.testing.assert_array_less(t_sol, 0.5)
+
     def test_model_solver_ode(self):
         # Create model
         model = pybamm.BaseModel()
@@ -83,6 +173,25 @@ class TestScikitsSolver(unittest.TestCase):
         np.testing.assert_array_equal(solver.t, t_eval)
         np.testing.assert_allclose(solver.y[0], np.exp(0.1 * solver.t))
 
+    def test_model_solver_ode_events(self):
+        # Create model
+        model = pybamm.BaseModel()
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=whole_cell)
+        model.rhs = {var: 0.1 * var}
+        model.initial_conditions = {var: 1}
+        model.initial_conditions_ydot = {var: 0.1}
+        model.events = [pybamm.Function(np.min, var - 1.5)]
+        disc = StandardModelTest(model).disc
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.ScikitsOdeSolver(tol=1e-9)
+        t_eval = np.linspace(0, 10, 100)
+        solver.solve(model, t_eval)
+        np.testing.assert_allclose(solver.y[0], np.exp(0.1 * solver.t))
+        np.testing.assert_array_less(solver.y[0], 1.5)
+
     def test_model_solver_dae(self):
         # Create model
         model = pybamm.BaseModel()
@@ -104,6 +213,32 @@ class TestScikitsSolver(unittest.TestCase):
         np.testing.assert_allclose(solver.y[0], np.exp(0.1 * solver.t))
         np.testing.assert_allclose(solver.y[-1], 2 * np.exp(0.1 * solver.t))
 
+    def test_model_solver_dae_events(self):
+        # Create model
+        model = pybamm.BaseModel()
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        var1 = pybamm.Variable("var1", domain=whole_cell)
+        var2 = pybamm.Variable("var2", domain=whole_cell)
+        model.rhs = {var1: 0.1 * var1}
+        model.algebraic = {var2: 2 * var1 - var2}
+        model.initial_conditions = {var1: 1, var2: 2}
+        model.initial_conditions_ydot = {var1: 0.1, var2: 0.2}
+        model.events = [
+            pybamm.Function(np.min, var1 - 1.5),
+            pybamm.Function(np.min, var2 - 2.5),
+        ]
+        disc = StandardModelTest(model).disc
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        t_eval = np.linspace(0, 1, 100)
+        solver.solve(model, t_eval)
+        np.testing.assert_array_less(solver.y[0], 1.5)
+        np.testing.assert_array_less(solver.y[-1], 2.5)
+        np.testing.assert_allclose(solver.y[0], np.exp(0.1 * solver.t))
+        np.testing.assert_allclose(solver.y[-1], 2 * np.exp(0.1 * solver.t))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/test_solvers/test_scipy_solver.py b/tests/test_solvers/test_scipy_solver.py
index b7901174..813208dd 100644
--- a/tests/test_solvers/test_scipy_solver.py
+++ b/tests/test_solvers/test_scipy_solver.py
@@ -35,13 +35,51 @@ class TestScipySolver(unittest.TestCase):
         t_sol, y_sol = solver.integrate(exponential_decay, y0, t_eval)
         np.testing.assert_allclose(y_sol[0], np.exp(-0.1 * t_sol))
 
+    def test_integrate_with_event(self):
+        # Constant
+        solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
+
+        def constant_growth(t, y):
+            return 0.5 * np.ones_like(y)
+
+        def y_eq_2(t, y):
+            return y - 2
+
+        y0 = np.array([0])
+        t_eval = np.linspace(0, 10, 100)
+        t_sol, y_sol = solver.integrate(constant_growth, y0, t_eval, events=[y_eq_2])
+        self.assertLess(len(t_sol), len(t_eval))
+        np.testing.assert_allclose(0.5 * t_sol, y_sol[0])
+
+        # Exponential decay
+        solver = pybamm.ScipySolver(tol=1e-8, method="BDF")
+
+        def exponential_growth(t, y):
+            return y
+
+        def y_eq_5(t, y):
+            return np.max(y - 5)
+
+        def t_eq_6(t, y):
+            return t - 6
+
+        y0 = np.array([1, 2])
+        t_eval = np.linspace(0, 7, 100)
+        t_sol, y_sol = solver.integrate(
+            exponential_growth, y0, t_eval, events=[y_eq_5, t_eq_6]
+        )
+        np.testing.assert_allclose(y_sol[0], np.exp(t_sol), rtol=1e-6)
+        np.testing.assert_allclose(y_sol[1], 2 * np.exp(t_sol), rtol=1e-6)
+        np.testing.assert_array_less(t_sol, 6)
+        np.testing.assert_array_less(y_sol, 5)
+
     def test_model_solver(self):
         # Create model
         model = pybamm.BaseModel()
         domain = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=domain)
-        model.rhs = {var: pybamm.Scalar(0.1) * var}
-        model.initial_conditions = {var: pybamm.Scalar(1)}
+        model.rhs = {var: 0.1 * var}
+        model.initial_conditions = {var: 1}
         # No need to set parameters; can use base discretisation (no spatial operators)
 
         # create discretisation
@@ -56,6 +94,29 @@ class TestScipySolver(unittest.TestCase):
         np.testing.assert_array_equal(solver.t, t_eval)
         np.testing.assert_allclose(solver.y[0], np.exp(0.1 * solver.t))
 
+    def test_model_solver_with_event(self):
+        # Create model
+        model = pybamm.BaseModel()
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=domain)
+        model.rhs = {var: -0.1 * var}
+        model.initial_conditions = {var: 1}
+        model.events = [pybamm.Function(np.min, var - 0.5)]
+        # No need to set parameters; can use base discretisation (no spatial operators)
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+        # Solve
+        solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
+        t_eval = np.linspace(0, 10, 100)
+        solver.solve(model, t_eval)
+        self.assertLess(len(solver.t), len(t_eval))
+        np.testing.assert_array_equal(solver.t, t_eval[: len(solver.t)])
+        np.testing.assert_allclose(solver.y[0], np.exp(-0.1 * solver.t))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")