Return event state

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+-   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))
 -   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))
 -   Added submodels for active material ([#1262](https://github.com/pybamm-team/PyBaMM/pull/1262))
@@ -14,8 +15,10 @@
 ## Bug fixes
 
 -   Fixed `Simulation` and `model.new_copy` to fix a bug where changes to the model were overwritten ([#1278](https://github.com/pybamm-team/PyBaMM/pull/1278))
+
 ## Breaking changes
 
+-   An extra sub-solution is created when an event is triggered during an experiment  ([#1300](https://github.com/pybamm-team/PyBaMM/pull/1300))
 -   Boolean model options ('sei porosity change', 'convection') must now be given in string format ('true' or 'false' instead of True or False) ([#1280](https://github.com/pybamm-team/PyBaMM/pull/1280))
 -   Operations such as `1*x` and `0+x` now directly return `x`. This can be bypassed by explicitly creating the binary operators, e.g. `pybamm.Multiplication(1, x)` ([#1252](https://github.com/pybamm-team/PyBaMM/pull/1252))
 
@@ -55,7 +58,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
 

examples/scripts/experimental_protocols/cccv.py

@@ -20,10 +20,13 @@
 sim.solve()
 
 # Plot voltages from the discharge segments only
+# Note: an additional sub-solution is created each time an event is triggered
+# so we need to specify the index for each discharge
+discharge_idx = [0, 7, 14]
 fig, ax = plt.subplots()
-for i in range(3):
+for i, idx in enumerate(discharge_idx):
     # Extract sub solutions
-    sol = sim.solution.sub_solutions[i * 5]
+    sol = sim.solution.sub_solutions[idx]
     # Extract variables
     t = sol["Time [h]"].entries
     V = sol["Terminal voltage [V]"].entries

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.
@@ -454,6 +448,20 @@ def solve(
                 pybamm.logger.info(self.experiment.operating_conditions_strings[idx])
                 inputs.update(exp_inputs)
                 kwargs["inputs"] = inputs
+                # If the last time entry was due to an event take a single step
+                # to get the solution back on to the correct reporting period.
+                # e.g. if the period is 30s and an event is triggered at 32s we
+                # return the times [0, 30, 32, 60] instead of [0, 30, 32, 62]
+                if hasattr(self._solution, "t_event"):
+                    if self._solution.t_event is not None and (
+                        self._solution.t_event[0] == self.solution.t[-1]
+                    ):
+                        dt_event = (
+                            exp_inputs["period"]
+                            - (self._solution.t[-1] - self._solution.t[-2])
+                            * self._solution.timescale_eval
+                        )
+                        self.step(dt_event, solver=solver, **kwargs)
                 # Make sure we take at least 2 timesteps
                 npts = max(int(round(dt / exp_inputs["period"])) + 1, 2)
                 self.step(dt, solver=solver, npts=npts, **kwargs)
@@ -505,12 +513,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

@@ -901,6 +901,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

@@ -351,15 +351,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_experiments/test_simulation_with_experiment.py

@@ -115,6 +115,26 @@ def test_inputs(self):
         sim.solve(inputs={"Dsn": 2})
         np.testing.assert_array_equal(sim.solution.inputs["Dsn"], 2)
 
+    def test_experiment_period_with_event(self):
+        experiment = pybamm.Experiment(
+            ["Discharge at 1C for 1 hour or until 3.6V", "Rest for 1 hour"],
+            period="5 minutes",
+        )
+        model = pybamm.lithium_ion.SPM()
+        param = pybamm.ParameterValues(chemistry=pybamm.parameter_sets.Marquis2019)
+        param["Negative electrode diffusivity [m2.s-1]"] = (
+            pybamm.InputParameter("Dsn") * 3.9e-14
+        )
+
+        sim = pybamm.Simulation(model, experiment=experiment, parameter_values=param)
+        sim.solve(solver=pybamm.CasadiSolver(), inputs={"Dsn": 1})
+
+        # check all times apart the event time are divisible by 300 (5m=300s)
+        times = sim.solution["Time [s]"].entries.tolist()
+        times = np.array([time % 300 for time in times])
+        times[times < 1e-9] = 0
+        self.assertEqual(1, np.count_nonzero(times))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

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):