created a function to sort and clean discontinuties, and a corresponding test

Author: isaacbasil

src/pybamm/solvers/base_solver.py

@@ -963,38 +963,45 @@ def solve(
             return solutions
 
     @staticmethod
-    def _get_discontinuity_start_end_indices(model, inputs, t_eval):
-        if not model.discontinuity_events_eval:
-            pybamm.logger.verbose("No discontinuity events found")
-            return [0], [len(t_eval)], t_eval
-
-        # Calculate discontinuities
-        discontinuities = [
-            # Assuming that discontinuities do not depend on
-            # input parameters when len(input_list) > 1, only
-            # `inputs` is passed to `evaluate`.
-            # See https://github.com/pybamm-team/PyBaMM/pull/1261
-            event.expression.evaluate(inputs=inputs)
-            for event in model.discontinuity_events_eval
-        ]
-
+    def _sort_and_clean_discontinuities(discontinuities, t_eval):
         # make sure they are increasing in time
         discontinuities = sorted(discontinuities)
 
-        # remove any identical discontinuities
+        # remove any identical discontinuities, and also unwanted discontinuities
+        # at the beginning of the integration (see https://github.com/pybamm-team/PyBaMM/pull/5075)
         discontinuities = [
             v
             for i, v in enumerate(discontinuities)
             if (
                 i == len(discontinuities) - 1
                 or discontinuities[i] < discontinuities[i + 1]
             )
-            and v > t_eval[0]
+            and v > t_eval[0] 
         ]
 
         # remove any discontinuities after end of t_eval
         discontinuities = [v for v in discontinuities if v < t_eval[-1]]
 
+        return discontinuities
+    
+    def _get_discontinuity_start_end_indices(self, model, inputs, t_eval):
+        if not model.discontinuity_events_eval:
+            pybamm.logger.verbose("No discontinuity events found")
+            return [0], [len(t_eval)], t_eval
+
+        # Calculate discontinuities
+        discontinuities = [
+            # Assuming that discontinuities do not depend on
+            # input parameters when len(input_list) > 1, only
+            # `inputs` is passed to `evaluate`.
+            # See https://github.com/pybamm-team/PyBaMM/pull/1261
+            event.expression.evaluate(inputs=inputs)
+            for event in model.discontinuity_events_eval
+        ]
+
+        # sort and remove unwanted discontinuities
+        discontinuities = self._sort_and_clean_discontinuities(discontinuities, t_eval)
+
         pybamm.logger.verbose(f"Discontinuity events found at t = {discontinuities}")
         if isinstance(inputs, list):
             raise pybamm.SolverError(

tests/unit/test_solvers/test_base_solver.py

@@ -435,18 +435,11 @@ def test_on_extrapolation_settings(self):
             base_solver.on_extrapolation = "invalid"
 
     def test_discontinuity_removed_at_nonzero_initial_time(self):
-        # Test that discontinuity caused by Heaviside(t) is removed when solver called with non-zero initial time
-        model = pybamm.BaseModel()
-        u = pybamm.Variable("u")
-        v = pybamm.Variable("v")
-        model.rhs = {v: -1 * (pybamm.t < 1)}
-        model.algebraic = {u: v - 1 - u}
-        model.initial_conditions = {v: 1, u: 0}
-        disc = pybamm.Discretisation()
-        disc.process_model(model)
-        solver = pybamm.IDAKLUSolver()
-        sol1 = solver.solve(model, t_eval=[0, 1])
+        # Test that the discontinuity caused by Heaviside(t) will be removed when 
+        # the solver is called with t_eval[0] != 0 
+        base_solver = pybamm.BaseSolver()
+        t_eval = [1, 2]
+        uncleaned_discontinuities = [1]
+        cleaned_discontinuities = base_solver._sort_and_clean_discontinuities(uncleaned_discontinuities, t_eval)
 
-        model.set_initial_conditions_from(sol1)
-        sol2 = solver.solve(model, t_eval=[1, 2])
-        assert sol2 is not None
+        assert t_eval[0] not in cleaned_discontinuities