Merge pull request #1089 from pybamm-team/issue-1082-casadi

Issue 1082 casadi

CHANGELOG.md

@@ -2,7 +2,7 @@
 
 ## Features
 
--   Reformatted Getting Starte notebooks ([#1083](https://github.com/pybamm-team/PyBaMM/pull/1083))
+-   Reformatted Getting Started notebooks ([#1083](https://github.com/pybamm-team/PyBaMM/pull/1083))
 -   Reformatted Landesfeind electrolytes ([#1064](https://github.com/pybamm-team/PyBaMM/pull/1064))
 -   Adapted examples to be run in Google Colab ([#1061](https://github.com/pybamm-team/PyBaMM/pull/1061))
 -   Added some new solvers for algebraic models ([#1059](https://github.com/pybamm-team/PyBaMM/pull/1059))
@@ -12,6 +12,8 @@
 
 ## Optimizations
 
+-   Reformatted CasADi "safe" mode to deal with events better ([#1089](https://github.com/pybamm-team/PyBaMM/pull/1089))
+
 ## Bug fixes
 
 -   2D processed variables can now be evaluated at the domain boundaries ([#1088](https://github.com/pybamm-team/PyBaMM/pull/1088))

examples/scripts/compare_lithium_ion.py

@@ -3,7 +3,7 @@
 #
 import pybamm
 
-# pybamm.set_logging_level("INFO")
+pybamm.set_logging_level("INFO")
 
 # load models
 models = [

pybamm/solvers/casadi_solver.py

@@ -151,7 +151,8 @@ def _integrate(self, model, t_eval, inputs=None):
                 # Non-dimensionalise provided dt_max
                 dt_max = self.dt_max / model.timescale_eval
             else:
-                dt_max = 0.01
+                # t_f is the dimensionless final time (scaled with the timescale)
+                dt_max = 0.1 * min(1, t_f)
             dt_eval_max = np.max(np.diff(t_eval)) * 1.01
             dt_max = np.max([dt_max, dt_eval_max])
             while t < t_f:
@@ -186,11 +187,11 @@ def _integrate(self, model, t_eval, inputs=None):
                     count += 1
                     if count >= self.max_step_decrease_count:
                         raise pybamm.SolverError(
-                            """
-                            Maximum number of decreased steps occurred at t={}. Try
-                            solving the model up to this time only or reducing dt_max.
-                            """.format(
-                                t
+                            "Maximum number of decreased steps occurred at t={}. Try "
+                            "solving the model up to this time only or reducing dt_max "
+                            "(currently, dt_max={})."
+                            "".format(
+                                t * model.timescale_eval, dt_max * model.timescale_eval
                             )
                         )
                 # Check most recent y to see if any events have been crossed
@@ -246,18 +247,27 @@ def event_fun(t):
                     t_event = np.nanmin(t_events)
                     y_event = y_sol(t_event)
 
-                    # return truncated solution
-                    t_truncated = current_step_sol.t[current_step_sol.t < t_event]
-                    y_trunctaed = current_step_sol.y[:, 0 : len(t_truncated)]
-                    truncated_step_sol = pybamm.Solution(t_truncated, y_trunctaed)
+                    # solve again until the event time
+                    # See comments above on creating t_window
+                    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])
+
+                    integrator = self.get_integrator(model, t_window, inputs)
+                    current_step_sol = self._run_integrator(
+                        integrator, model, y0, inputs, t_window
+                    )
+
                     # assign temporary solve time
-                    truncated_step_sol.solve_time = np.nan
+                    current_step_sol.solve_time = np.nan
                     # append solution from the current step to solution
-                    solution.append(truncated_step_sol)
-
+                    solution.append(current_step_sol)
                     solution.termination = "event"
                     solution.t_event = t_event
                     solution.y_event = y_event
+
                     break
                 else:
                     # assign temporary solve time

tests/integration/test_models/standard_output_comparison.py

@@ -25,7 +25,7 @@ def get_output_times(self):
 
         # Assign common time
         solution0 = self.solutions[0]
-        max_index = np.where(solution0.t == max_t)[0][0]
+        max_index = np.where(solution0.t >= max_t)[0][0]
         t_common = solution0.t[:max_index]
 
         # Check times

tests/unit/test_solvers/test_casadi_solver.py

@@ -119,7 +119,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[0], 1.5)
-        np.testing.assert_array_less(solution.y[-1], 2.5)
+        np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-10)
         np.testing.assert_array_almost_equal(
             solution.y[0], np.exp(0.1 * solution.t), decimal=5
         )
@@ -133,7 +133,9 @@ 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[0], 1.5)
-        np.testing.assert_array_less(solution.y[-1], 2.5)
+        np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-10)
+        # test the last entry is exactly 2.5
+        np.testing.assert_array_almost_equal(solution.y[-1, -1], 2.5, decimal=2)
         np.testing.assert_array_almost_equal(
             solution.y[0], np.exp(0.1 * solution.t), decimal=5
         )
@@ -168,7 +170,8 @@ def test_model_solver_events(self):
         disc.process_model(model)
         solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8)
         solution = solver.solve(model, t_eval)
-        np.testing.assert_array_less(solution.y[0], 1.02)
+        np.testing.assert_array_less(solution.y[0], 1.02 + 1e-10)
+        np.testing.assert_array_almost_equal(solution.y[0, -1], 1.02, decimal=2)
 
     def test_model_step(self):
         # Create model
@@ -299,8 +302,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_allclose(solution.y[0], np.exp(-0.1 * solution.t), rtol=1e-06)
+        np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t), rtol=1e-04)
 
     def test_model_solver_dae_inputs_in_initial_conditions(self):
         # Create model