Merge pull request #1300 from pybamm-team/return-event-state

Return event state

CHANGELOG.md

@@ -4,6 +4,7 @@
 
 -   Added new functionality for `Interpolant` ([#1312](https://github.com/pybamm-team/PyBaMM/pull/1312))
 -   Added option to express experiments (and extract solutions) in terms of cycles of operating condition ([#1309](https://github.com/pybamm-team/PyBaMM/pull/1309))
+-   The event time and state are now returned as part of `Solution.t` and `Solution.y` so that the event is accurately captured in the returned solution ([#1300](https://github.com/pybamm-team/PyBaMM/pull/1300))
 -   Reformatted the `BasicDFNHalfCell` to be consistent with the other models ([#1282](https://github.com/pybamm-team/PyBaMM/pull/1282))
 -   Added option to make the total interfacial current density a state ([#1280](https://github.com/pybamm-team/PyBaMM/pull/1280))
 -   Added functionality to initialize a model using the solution from another model ([#1278](https://github.com/pybamm-team/PyBaMM/pull/1278))
@@ -60,7 +61,7 @@ This release introduces a new aging model for particle swelling and cracking, a
 
 ## Breaking changes
 
--    The parameters "Positive/Negative particle distribution in x" and "Positive/Negative surface area to volume ratio distribution in x" have been deprecated. Instead, users can provide "Positive/Negative particle radius [m]" and "Positive/Negative surface area to volume ratio [m-1]" directly as functions of through-cell position (x [m]) ([#1237](https://github.com/pybamm-team/PyBaMM/pull/1237)) 
+-    The parameters "Positive/Negative particle distribution in x" and "Positive/Negative surface area to volume ratio distribution in x" have been deprecated. Instead, users can provide "Positive/Negative particle radius [m]" and "Positive/Negative surface area to volume ratio [m-1]" directly as functions of through-cell position (x [m]) ([#1237](https://github.com/pybamm-team/PyBaMM/pull/1237))
 
 # [v0.2.4](https://github.com/pybamm-team/PyBaMM/tree/v0.2.4) - 2020-09-07
 

pybamm/simulation.py

@@ -330,13 +330,7 @@ def build(self, check_model=True):
                 self._model_with_set_params, inplace=False, check_model=check_model
             )
 
-    def solve(
-        self,
-        t_eval=None,
-        solver=None,
-        check_model=True,
-        **kwargs,
-    ):
+    def solve(self, t_eval=None, solver=None, check_model=True, **kwargs):
         """
         A method to solve the model. This method will automatically build
         and set the model parameters if not already done so.
@@ -518,12 +512,7 @@ def step(self, dt, solver=None, npts=2, save=True, **kwargs):
             solver = self.solver
 
         self._solution = solver.step(
-            self._solution,
-            self.built_model,
-            dt,
-            npts=npts,
-            save=save,
-            **kwargs,
+            self._solution, self.built_model, dt, npts=npts, save=save, **kwargs
         )
 
         return self.solution

pybamm/solvers/base_solver.py

@@ -781,6 +781,10 @@ def step(
                     "Cannot step empty model, use `pybamm.DummySolver` instead"
                 )
 
+        # Make sure dt is positive
+        if dt <= 0:
+            raise pybamm.SolverError("Step time must be positive")
+
         # Set timer
         timer = pybamm.Timer()
 
@@ -901,6 +905,16 @@ def get_termination_reason(self, solution, events):
             termination_event = min(final_event_values, key=final_event_values.get)
             # Add the event to the solution object
             solution.termination = "event: {}".format(termination_event)
+            # Update t, y and inputs to include event time and state
+            # Note: if the final entry of t is equal to the event time to within
+            # the absolute tolerance we skip this (having duplicate entries
+            # causes an error later in ProcessedVariable)
+            if solution.t_event - solution._t[-1] > self.atol:
+                solution._t = np.concatenate((solution._t, solution.t_event))
+                solution._y = np.concatenate((solution._y, solution.y_event), axis=1)
+                for name, inp in solution.inputs.items():
+                    solution._inputs[name] = np.c_[inp, inp[:, -1]]
+
             return "the termination event '{}' occurred".format(termination_event)
 
     def _set_up_ext_and_inputs(self, model, external_variables, inputs):

pybamm/solvers/casadi_solver.py

@@ -282,8 +282,8 @@ def event_fun(t):
                         # append solution from the current step to solution
                         solution.append(current_step_sol)
                     solution.termination = "event"
-                    solution.t_event = t_event
-                    solution.y_event = y_event
+                    solution.t_event = np.array([t_event])
+                    solution.y_event = y_event[:, np.newaxis]
 
                     break
                 else:

pybamm/solvers/idaklu_solver.py

@@ -35,7 +35,7 @@ class IDAKLUSolver(pybamm.BaseSolver):
         conditions. Otherwise, the solver uses 'scipy.optimize.root' with method
         specified by 'root_method' (e.g. "lm", "hybr", ...)
     root_tol : float, optional
-        The tolerance for the initial-condition solver (default is 1e-8).
+        The tolerance for the initial-condition solver (default is 1e-6).
     """
 
     def __init__(

pybamm/solvers/processed_variable.py

@@ -353,15 +353,11 @@ def initialise_2D(self):
             self.second_dimension = "x"
             self.r_sol = first_dim_pts
             self.x_sol = second_dim_pts
-        elif (
-            self.domain[0]
-            in [
-                "negative electrode",
-                "separator",
-                "positive electrode",
-            ]
-            and self.auxiliary_domains["secondary"] == ["current collector"]
-        ):
+        elif self.domain[0] in [
+            "negative electrode",
+            "separator",
+            "positive electrode",
+        ] and self.auxiliary_domains["secondary"] == ["current collector"]:
             self.first_dimension = "x"
             self.second_dimension = "z"
             self.x_sol = first_dim_pts

tests/integration/test_models/standard_output_tests.py

@@ -603,19 +603,23 @@ def test_interfacial_current_average(self):
         multiplied by the interfacial current density is equal to the true
         value."""
 
+        # The final time corresponds to an event (voltage cut-off). At this time
+        # the average property is satisfied but to a lesser degree of accuracy
+        t = self.t[:-1]
+
         np.testing.assert_array_almost_equal(
             np.mean(
-                self.a_n(self.t, self.x_n)
-                * (self.j_n(self.t, self.x_n) + self.j_n_sei(self.t, self.x_n)),
+                self.a_n(t, self.x_n)
+                * (self.j_n(t, self.x_n) + self.j_n_sei(t, self.x_n)),
                 axis=0,
             ),
             self.i_cell / self.l_n,
             decimal=4,
         )
         np.testing.assert_array_almost_equal(
             np.mean(
-                self.a_p(self.t, self.x_p)
-                * (self.j_p(self.t, self.x_p) + self.j_p_sei(self.t, self.x_p)),
+                self.a_p(t, self.x_p)
+                * (self.j_p(t, self.x_p) + self.j_p_sei(t, self.x_p)),
                 axis=0,
             ),
             -self.i_cell / self.l_p,

tests/unit/test_solvers/test_base_solver.py

@@ -71,6 +71,11 @@ def test_nonmonotonic_teval(self):
         ):
             solver.solve(model, np.array([1, 2, 3, 2]))
 
+        # Check stepping with negative step size
+        dt = -1
+        with self.assertRaisesRegex(pybamm.SolverError, "Step time must be positive"):
+            solver.step(None, model, dt)
+
     def test_solution_time_length_fail(self):
         model = pybamm.BaseModel()
         v = pybamm.Scalar(1)

tests/unit/test_solvers/test_scikits_solvers.py

@@ -149,8 +149,8 @@ def test_model_solver_ode_events_python(self):
         t_eval = np.linspace(0, 10, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
-        np.testing.assert_array_less(solution.y[0], 1.5)
-        np.testing.assert_array_less(solution.y[0], 1.25)
+        np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6)
+        np.testing.assert_array_less(solution.y[0], 1.25 + 1e-6)
 
     def test_model_solver_ode_jacobian_python(self):
         model = pybamm.BaseModel()
@@ -249,8 +249,8 @@ def test_model_solver_dae_events_python(self):
         solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
         t_eval = np.linspace(0, 5, 100)
         solution = solver.solve(model, t_eval)
-        np.testing.assert_array_less(solution.y[0], 1.5)
-        np.testing.assert_array_less(solution.y[-1], 2.5)
+        np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6)
+        np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-6)
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
         np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
@@ -321,8 +321,8 @@ def nonsmooth_mult(t):
 
         # check solution
         for solution in [solution1, solution2]:
-            np.testing.assert_array_less(solution.y[0], 1.5)
-            np.testing.assert_array_less(solution.y[-1], 2.5)
+            np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6)
+            np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-6)
             var1_soln = np.exp(0.2 * solution.t)
             y0 = np.exp(0.2 * discontinuity)
             var1_soln[solution.t > discontinuity] = y0 * np.exp(
@@ -388,8 +388,8 @@ def test_model_solver_dae_multiple_nonsmooth_python(self):
 
         # check solution
         for solution in [solution1, solution2]:
-            np.testing.assert_array_less(solution.y[0], 0.55)
-            np.testing.assert_array_less(solution.y[-1], 1.2)
+            np.testing.assert_array_less(solution.y[0], 0.55 + 1e-6)
+            np.testing.assert_array_less(solution.y[-1], 1.2 + 1e-6)
             var1_soln = (solution.t % a) ** 2 / 2 + a ** 2 / 2 * (solution.t // a)
             var2_soln = 2 * var1_soln
             np.testing.assert_allclose(solution.y[0], var1_soln, rtol=1e-06)
@@ -569,8 +569,8 @@ def test_model_solver_ode_events_casadi(self):
         t_eval = np.linspace(0, 10, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
-        np.testing.assert_array_less(solution.y[0], 1.5)
-        np.testing.assert_array_less(solution.y[0], 1.25)
+        np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6)
+        np.testing.assert_array_less(solution.y[0], 1.25 + 1e-6)
 
     def test_model_solver_dae_events_casadi(self):
         # Create model
@@ -595,8 +595,8 @@ def test_model_solver_dae_events_casadi(self):
             solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
             t_eval = np.linspace(0, 5, 100)
             solution = solver.solve(model_disc, t_eval)
-            np.testing.assert_array_less(solution.y[0], 1.5)
-            np.testing.assert_array_less(solution.y[-1], 2.5)
+            np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6)
+            np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-6)
             np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
             np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
@@ -625,8 +625,8 @@ def test_model_solver_dae_inputs_events(self):
                 solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
             t_eval = np.linspace(0, 5, 100)
             solution = solver.solve(model, t_eval, inputs={"rate 1": 0.1, "rate 2": 2})
-            np.testing.assert_array_less(solution.y[0], 1.5)
-            np.testing.assert_array_less(solution.y[-1], 2.5)
+            np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6)
+            np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-6)
             np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
             np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
@@ -730,8 +730,8 @@ def test_model_step_events(self):
         while time < end_time:
             step_solution = step_solver.step(step_solution, model, dt=dt, npts=10)
             time += dt
-        np.testing.assert_array_less(step_solution.y[0], 1.5)
-        np.testing.assert_array_less(step_solution.y[-1], 2.5001)
+        np.testing.assert_array_less(step_solution.y[0], 1.5 + 1e-6)
+        np.testing.assert_array_less(step_solution.y[-1], 2.5 + 1e-6)
         np.testing.assert_array_almost_equal(
             step_solution.y[0], np.exp(0.1 * step_solution.t), decimal=5
         )
@@ -771,8 +771,8 @@ def test_model_step_nonsmooth_events(self):
         while time < end_time:
             step_solution = step_solver.step(step_solution, model, dt=dt, npts=10)
             time += dt
-        np.testing.assert_array_less(step_solution.y[0], 0.55)
-        np.testing.assert_array_less(step_solution.y[-1], 1.2)
+        np.testing.assert_array_less(step_solution.y[0], 0.55 + 1e-6)
+        np.testing.assert_array_less(step_solution.y[-1], 1.2 + 1e-6)
         var1_soln = (step_solution.t % a) ** 2 / 2 + a ** 2 / 2 * (step_solution.t // a)
         var2_soln = 2 * var1_soln
         np.testing.assert_array_almost_equal(step_solution.y[0], var1_soln, decimal=5)
@@ -854,8 +854,8 @@ def nonsmooth_rate(t):
 
             # check solution
             for solution in [solution1, solution2]:
-                np.testing.assert_array_less(solution.y[0], 1.5)
-                np.testing.assert_array_less(solution.y[-1], 2.5)
+                np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6)
+                np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-6)
                 var1_soln = np.exp(0.2 * solution.t)
                 y0 = np.exp(0.2 * discontinuity)
                 var1_soln[solution.t > discontinuity] = y0 * np.exp(

tests/unit/test_solvers/test_scipy_solver.py

@@ -99,7 +99,7 @@ def test_model_solver_with_event_python(self):
         t_eval = np.linspace(0, 10, 100)
         solution = solver.solve(model, t_eval)
         self.assertLess(len(solution.t), len(t_eval))
-        np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+        np.testing.assert_array_equal(solution.t[:-1], t_eval[: len(solution.t) - 1])
         np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
 
     def test_model_solver_ode_with_jacobian_python(self):
@@ -214,7 +214,7 @@ def test_model_solver_with_inputs(self):
         t_eval = np.linspace(0, 10, 100)
         solution = solver.solve(model, t_eval, inputs={"rate": 0.1})
         self.assertLess(len(solution.t), len(t_eval))
-        np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+        np.testing.assert_array_equal(solution.t[:-1], t_eval[: len(solution.t) - 1])
         np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
 
     def test_model_solver_with_external(self):
@@ -272,7 +272,9 @@ def test_model_solver_with_event_with_casadi(self):
             t_eval = np.linspace(0, 10, 100)
             solution = solver.solve(model_disc, t_eval)
             self.assertLess(len(solution.t), len(t_eval))
-            np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+            np.testing.assert_array_equal(
+                solution.t[:-1], t_eval[: len(solution.t) - 1]
+            )
             np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
 
     def test_model_solver_with_inputs_with_casadi(self):
@@ -297,7 +299,7 @@ def test_model_solver_with_inputs_with_casadi(self):
         t_eval = np.linspace(0, 10, 100)
         solution = solver.solve(model, t_eval, inputs={"rate": 0.1})
         self.assertLess(len(solution.t), len(t_eval))
-        np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+        np.testing.assert_array_equal(solution.t[:-1], t_eval[: len(solution.t) - 1])
         np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
 
     def test_model_solver_inputs_in_initial_conditions(self):