Merge pull request #1328 from pybamm-team/issue-1318-casadi-event

Re-solve before root finding when handling events in casadi

CHANGELOG.md

@@ -2,6 +2,8 @@
 
 ## Features
 
+
+-   Updated the way events are handled in `CasadiSolver` for more accurate event location ([#1328](https://github.com/pybamm-team/PyBaMM/pull/1328))
 -   Added error message if initial conditions are outside the bounds of a variable ([#1326](https://github.com/pybamm-team/PyBaMM/pull/1326))
 -   Added temperature dependence to density, heat capacity and thermal conductivity ([#1323](https://github.com/pybamm-team/PyBaMM/pull/1323))
 -   Updated solvers' method `solve()` so it can take a list of inputs dictionaries as the `inputs` keyword argument. In this case the model is solved for each input set in the list, and a list of solutions mapping the set of inputs to the solutions is returned. Note that `solve()` can still take a single dictionary as the `inputs` keyword argument. In this case the behaviour is unchanged compared to previous versions.

pybamm/solvers/base_solver.py

@@ -1028,10 +1028,9 @@ def get_termination_reason(self, solution, events):
             # 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:
+            # Note: if the final entry of t is equal to the event time we skip
+            # this (having duplicate entries causes an error later in ProcessedVariable)
+            if solution.t_event != solution._t[-1]:
                 solution._t = np.concatenate((solution._t, solution.t_event))
                 if isinstance(solution.y, casadi.DM):
                     solution._y = casadi.horzcat(solution.y, solution.y_event)

pybamm/solvers/casadi_solver.py

@@ -165,7 +165,7 @@ def _integrate(self, model, t_eval, inputs=None):
                                 == pybamm.EventType.INTERPOLANT_EXTRAPOLATION
                                 and (
                                     event.expression.evaluate(
-                                        t, y0.full(), inputs=inputs,
+                                        t, y0.full(), inputs=inputs
                                     )
                                     < self.extrap_tol
                                 ).any()
@@ -297,9 +297,21 @@ def _integrate(self, model, t_eval, inputs=None):
                     event_ind = np.where(new_event_signs != init_event_signs)[0]
                     active_events = [model.terminate_events_eval[i] for i in event_ind]
 
+                    # solve again with a more dense t_window
+                    if len(t_window) < 10:
+                        t_window_dense = np.linspace(t_window[0], t_window[-1], 10)
+                        if self.mode == "safe":
+                            self.create_integrator(model, inputs, t_window_dense)
+                        current_step_sol = self._run_integrator(
+                            model, y0, inputs, t_window_dense
+                        )
+                    integration_time = current_step_sol.integration_time
+
                     # create interpolant to evaluate y in the current integration
                     # window
-                    y_sol = interp1d(current_step_sol.t, current_step_sol.y)
+                    y_sol = interp1d(
+                        current_step_sol.t, current_step_sol.y, kind="cubic"
+                    )
 
                     # loop over events to compute the time at which they were triggered
                     t_events = [None] * len(active_events)
@@ -332,21 +344,16 @@ def event_fun(t):
                     t_event = np.nanmin(t_events)
                     y_event = y_sol(t_event)
 
-                    # solve again until the event time
-                    # See comments above on creating t_window
+                    # call interpolant at times in t_eval and create solution
                     t_window = np.concatenate(
                         ([t], t_eval[(t_eval > t) & (t_eval < t_event)])
                     )
-                    if len(t_window) == 1:
-                        t_window = np.array([t, t_event])
-
-                    if self.mode == "safe":
-                        self.create_integrator(model, inputs, t_window)
-                    current_step_sol = self._run_integrator(model, y0, inputs, t_window)
-
+                    y_sol_casadi = casadi.DM(y_sol(t_window))
+                    current_step_sol = pybamm.Solution(t_window, y_sol_casadi)
+                    current_step_sol.integration_time = integration_time
                     # assign temporary solve time
                     current_step_sol.solve_time = np.nan
-                    # append solution from the current step to solution
+
                     if solution is None:
                         solution = current_step_sol
                     else:

tests/integration/test_models/standard_output_tests.py

@@ -603,23 +603,19 @@ 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(t, self.x_n)
-                * (self.j_n(t, self.x_n) + self.j_n_sei(t, self.x_n)),
+                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)),
                 axis=0,
             ),
             self.i_cell / self.l_n,
             decimal=4,
         )
         np.testing.assert_array_almost_equal(
             np.mean(
-                self.a_p(t, self.x_p)
-                * (self.j_p(t, self.x_p) + self.j_p_sei(t, self.x_p)),
+                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)),
                 axis=0,
             ),
             -self.i_cell / self.l_p,

tests/unit/test_solvers/test_casadi_solver.py

@@ -152,7 +152,7 @@ def test_model_solver_events(self):
         t_eval = np.linspace(0, 5, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_less(solution.y.full()[0], 1.5)
-        np.testing.assert_array_less(solution.y.full()[-1], 2.5)
+        np.testing.assert_array_less(solution.y.full()[-1], 2.5 + 1e-10)
         np.testing.assert_array_almost_equal(
             solution.y.full()[0], np.exp(0.1 * solution.t), decimal=5
         )