Repository: pybamm-team/PyBaMM

Issue:
Termination events not consistent
**Describe the bug**
I was running simulations on the LG M50 and I noticed a strange behaviour related to the termination events (strange behaviour for this case described below as it might not be that relevant). Trying to isolate what the problem was I ran it with the DFN and SPMe examples (Moura's parameters) and noticed that the termination events depend on t_eval which should not happen. Also, the solution does not provide the time and value of the event (probably because of the interpolation to t_eval).

Probably the reason is that the three events (max cp, min cn and min voltage) happen very close together so a high time resolution is needed to distinguish them.

**To Reproduce**
Steps to reproduce the behaviour:
1. Go to '/examples/scripts/DFN.py' or '/examples/scripts/SPMe.py'
2. Run the model for different number of points in the t_eval linspace.
3. **DFN:** if N = 100 the termination event is 'Maximum positive particle surface concentration' while if N = 1E4 the termination event is 'Minimum voltage'.
4. **SPMe:** if N = 100 the termination event is 'Maximum positive particle surface concentration' while if N = 5000 the termination event is 'Minimum voltage'.

**Additional context**
When running the LG M50 parameters for the SPMe I noticed that the terminal voltage can fall way below the cut-off voltage (set to 2.5V while the simulation went down to 2V). More interestingly, in other cases , when changing the value of the cut-off voltage it stopped at different points (which weren't the cut-off voltage) and the termination event was 'Minimum negative particle surface concentration'.


Discussion:

Patch:
diff --git a/examples/scripts/DFN.py b/examples/scripts/DFN.py
index a2133bd7a..412d3a6fc 100644
--- a/examples/scripts/DFN.py
+++ b/examples/scripts/DFN.py
@@ -1,6 +1,8 @@
 import pybamm
 import numpy as np
 
+pybamm.set_logging_level("INFO")
+
 # load model
 model = pybamm.lithium_ion.DFN()
 
diff --git a/examples/scripts/SPMe.py b/examples/scripts/SPMe.py
index 85acfcf3c..a3d32a290 100644
--- a/examples/scripts/SPMe.py
+++ b/examples/scripts/SPMe.py
@@ -1,6 +1,8 @@
 import pybamm
 import numpy as np
 
+pybamm.set_logging_level("INFO")
+
 # load model
 model = pybamm.lithium_ion.SPMe()
 
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 8accb1a5e..48023269f 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -69,7 +69,7 @@ class BaseSolver(object):
             final_event_values = {}
             for name, event in events.items():
                 final_event_values[name] = abs(
-                    event.evaluate(solution.t[-1], solution.y[:, -1])
+                    event.evaluate(solution.t_event, solution.y_event)
                 )
             termination_event = min(final_event_values, key=final_event_values.get)
             return "the termination event '{}' occurred".format(termination_event)
diff --git a/pybamm/solvers/scikits_dae_solver.py b/pybamm/solvers/scikits_dae_solver.py
index b50b1b217..4d42fcf1e 100644
--- a/pybamm/solvers/scikits_dae_solver.py
+++ b/pybamm/solvers/scikits_dae_solver.py
@@ -116,7 +116,11 @@ class ScikitsDaeSolver(pybamm.DaeSolver):
             elif sol.flag == 2:
                 termination = "event"
             return pybamm.Solution(
-                sol.values.t, np.transpose(sol.values.y), termination
+                sol.values.t,
+                np.transpose(sol.values.y),
+                sol.roots.t,
+                np.transpose(sol.roots.y),
+                termination
             )
         else:
             raise pybamm.SolverError(sol.message)
diff --git a/pybamm/solvers/scikits_ode_solver.py b/pybamm/solvers/scikits_ode_solver.py
index 9ad6cf7ff..17cfff62a 100644
--- a/pybamm/solvers/scikits_ode_solver.py
+++ b/pybamm/solvers/scikits_ode_solver.py
@@ -130,7 +130,11 @@ class ScikitsOdeSolver(pybamm.OdeSolver):
             elif sol.flag == 2:
                 termination = "event"
             return pybamm.Solution(
-                sol.values.t, np.transpose(sol.values.y), termination
+                sol.values.t,
+                np.transpose(sol.values.y),
+                sol.roots.t,
+                np.transpose(sol.roots.y),
+                termination
             )
         else:
             raise pybamm.SolverError(sol.message)
diff --git a/pybamm/solvers/scipy_solver.py b/pybamm/solvers/scipy_solver.py
index 5910c63d0..60e3db9ac 100644
--- a/pybamm/solvers/scipy_solver.py
+++ b/pybamm/solvers/scipy_solver.py
@@ -72,6 +72,7 @@ class ScipySolver(pybamm.OdeSolver):
             y0,
             t_eval=t_eval,
             method=self.method,
+            dense_output=True,
             **extra_options
         )
 
@@ -79,8 +80,22 @@ class ScipySolver(pybamm.OdeSolver):
             # Set the reason for termination
             if sol.message == "A termination event occurred.":
                 termination = "event"
+                t_event = []
+                for time in sol.t_events:
+                    if time:
+                        t_event = np.append(t_event, np.max(time))
+                t_event = np.array([np.max(t_event)])
+                y_event = sol.sol(t_event)
             elif sol.message.startswith("The solver successfully reached the end"):
                 termination = "final time"
-            return pybamm.Solution(sol.t, sol.y, termination)
+                t_event = None
+                y_event = np.array(None)
+            return pybamm.Solution(
+                sol.t,
+                sol.y,
+                t_event,
+                y_event,
+                termination
+            )
         else:
             raise pybamm.SolverError(sol.message)
diff --git a/pybamm/solvers/solution.py b/pybamm/solvers/solution.py
index 5856d1ea3..4b5f7b515 100644
--- a/pybamm/solvers/solution.py
+++ b/pybamm/solvers/solution.py
@@ -15,14 +15,21 @@ class Solution(object):
     y : :class:`numpy.array`, size (m, n)
         A two-dimensional array containing the values of the solution. y[i, :] is the
         vector of solutions at time t[i].
+    t_event : :class:`numpy.array`, size (1,)
+        A zero-dimensional array containing the time at which the event happens.
+    y_event : :class:`numpy.array`, size (m,)
+        A one-dimensional array containing the value of the solution at the time when
+        the event happens.
     termination : str
         String to indicate why the solution terminated
 
     """
 
-    def __init__(self, t, y, termination):
+    def __init__(self, t, y, t_event, y_event, termination):
         self.t = t
         self.y = y
+        self.t_event = t_event
+        self.y_event = y_event
         self.termination = termination
 
     @property
@@ -45,6 +52,26 @@ class Solution(object):
         "Updates the solution values"
         self._y = value
 
+    @property
+    def t_event(self):
+        "Time at which the event happens"
+        return self._t_event
+
+    @t_event.setter
+    def t_event(self, value):
+        "Updates the event time"
+        self._t_event = value
+
+    @property
+    def y_event(self):
+        "Value of the solution at the time of the event"
+        return self._y_event
+
+    @y_event.setter
+    def y_event(self, value):
+        "Updates the solution at the time of the event"
+        self._y_event = value
+
     @property
     def termination(self):
         "Reason for termination"


Test patch:
diff --git a/tests/unit/test_quick_plot.py b/tests/unit/test_quick_plot.py
index dee5a9fa6..2ff5c042a 100644
--- a/tests/unit/test_quick_plot.py
+++ b/tests/unit/test_quick_plot.py
@@ -104,7 +104,7 @@ class TestQuickPlot(unittest.TestCase):
             pybamm.QuickPlot(
                 pybamm.BaseModel(),
                 None,
-                [pybamm.Solution(0, 0, ""), pybamm.Solution(0, 0, "")],
+                [pybamm.Solution(0, 0, 0, 0, ""), pybamm.Solution(0, 0, 0, 0, "")],
             )
 
 
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 9ed921d3c..912af5da1 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -70,6 +70,8 @@ class TestScikitsSolvers(unittest.TestCase):
         self.assertLess(len(solution.t), len(t_eval))
         np.testing.assert_array_less(0, solution.y[0])
         np.testing.assert_array_less(solution.t, 0.5)
+        np.testing.assert_allclose(solution.t_event, 0.5)
+        np.testing.assert_allclose(solution.y_event, 0)
         self.assertEqual(solution.termination, "event")
 
         # Expnonential growth
@@ -93,6 +95,8 @@ class TestScikitsSolvers(unittest.TestCase):
         np.testing.assert_allclose(np.exp(solution.t), solution.y[0], rtol=1e-4)
         np.testing.assert_array_less(solution.y, 9)
         np.testing.assert_array_less(solution.y ** 2, 7)
+        np.testing.assert_allclose(solution.t_event, np.log(7) / 2, rtol=1e-4)
+        np.testing.assert_allclose(solution.y_event ** 2, 7, rtol=1e-4)
         self.assertEqual(solution.termination, "event")
 
     def test_ode_integrate_with_jacobian(self):
@@ -295,6 +299,9 @@ class TestScikitsSolvers(unittest.TestCase):
         np.testing.assert_allclose(1.0 * solution.t, solution.y[1])
         np.testing.assert_array_less(solution.y[0], 2)
         np.testing.assert_array_less(solution.y[1], 5)
+        np.testing.assert_allclose(solution.t_event, 4)
+        np.testing.assert_allclose(solution.y_event[0], 2)
+        np.testing.assert_allclose(solution.y_event[1], 4)
         self.assertEqual(solution.termination, "event")
 
         # Exponential decay
@@ -320,6 +327,9 @@ class TestScikitsSolvers(unittest.TestCase):
         np.testing.assert_allclose(solution.y[1], 2 * np.exp(-0.1 * solution.t))
         np.testing.assert_array_less(0.9, solution.y[0])
         np.testing.assert_array_less(solution.t, 0.5)
+        np.testing.assert_allclose(solution.t_event, 0.5)
+        np.testing.assert_allclose(solution.y_event[0], np.exp(-0.1 * 0.5))
+        np.testing.assert_allclose(solution.y_event[1], 2 * np.exp(-0.1 * 0.5))
         self.assertEqual(solution.termination, "event")
 
     def test_dae_integrate_with_jacobian(self):
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 68ca6f314..6749d61d3 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -66,6 +66,8 @@ class TestScipySolver(unittest.TestCase):
         solution = solver.integrate(constant_growth, y0, t_eval, events=[y_eq_2])
         self.assertLess(len(solution.t), len(t_eval))
         np.testing.assert_allclose(0.5 * solution.t, solution.y[0])
+        np.testing.assert_allclose(solution.t_event, 4)
+        np.testing.assert_allclose(solution.y_event, 2)
         self.assertEqual(solution.termination, "event")
 
         # Exponential decay
@@ -89,6 +91,9 @@ class TestScipySolver(unittest.TestCase):
         np.testing.assert_allclose(solution.y[1], 2 * np.exp(solution.t), rtol=1e-6)
         np.testing.assert_array_less(solution.t, 6)
         np.testing.assert_array_less(solution.y, 5)
+        np.testing.assert_allclose(solution.t_event, np.log(5 / 2), rtol=1e-6)
+        np.testing.assert_allclose(solution.y_event[0], 5 / 2, rtol=1e-6)
+        np.testing.assert_allclose(solution.y_event[1], 5, rtol=1e-6)
 
     def test_ode_integrate_with_jacobian(self):
         # Linear