#941 filter out discon events after the end of the solver time

CHANGELOG.md

@@ -31,6 +31,8 @@ New expression tree node types, models, parameter sets and solvers, as well as g
 
 ## Bug fixes
 
+-   Filter out discontinuities that occur after solve times 
+    ([#941](https://github.com/pybamm-team/PyBaMM/pull/945))
 -   Fixed tight layout for QuickPlot in jupyter notebooks ([#930](https://github.com/pybamm-team/PyBaMM/pull/930))
 -   Fixed bug raised if function returns a scalar ([#919](https://github.com/pybamm-team/PyBaMM/pull/919))
 -   Fixed event handling in `ScipySolver` ([#905](https://github.com/pybamm-team/PyBaMM/pull/905))

pybamm/solvers/base_solver.py

@@ -563,6 +563,13 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
             and v > 0
         ]
 
+        # remove any discontinuities after end of t_eval
+        discontinuities = [
+            v
+            for v in discontinuities
+            if v < t_eval_dimensionless[-1]
+        ]
+
         if len(discontinuities) > 0:
             pybamm.logger.info(
                 "Discontinuity events found at t = {}".format(discontinuities)

tests/unit/test_solvers/test_scikits_solvers.py

@@ -329,6 +329,41 @@ def nonsmooth_mult(t):
             np.testing.assert_allclose(solution.y[0], var1_soln, rtol=1e-06)
             np.testing.assert_allclose(solution.y[-1], var2_soln, rtol=1e-06)
 
+    def test_model_solver_dae_no_nonsmooth_python(self):
+        model = pybamm.BaseModel()
+        model.convert_to_format = "python"
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        var1 = pybamm.Variable("var1", domain=whole_cell)
+        var2 = pybamm.Variable("var2", domain=whole_cell)
+        discontinuity = 5.6
+
+        def nonsmooth_rate(t):
+            return 0.1 * int(t < discontinuity) + 0.1
+
+        def nonsmooth_mult(t):
+            return int(t < discontinuity) + 1.0
+
+        # put in an extra heaviside with no time dependence, this should be ignored by
+        # the solver i.e. no extra discontinuities added
+        model.rhs = {var1: 0.1 * var1}
+        model.algebraic = {var2: 2 * var1 - var2}
+        model.initial_conditions = {var1: 1, var2: 2}
+        disc = get_discretisation_for_testing()
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-9, atol=1e-9, root_method="lm")
+
+        # create two time series, one without a time point on the discontinuity,
+        # and one with
+        t_eval = np.linspace(0, 5, 10)
+        solution = solver.solve(model, t_eval)
+
+        # test solution, discontinuity should not be triggered
+        np.testing.assert_array_equal(solution.t, t_eval)
+        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))
+
     def test_model_solver_dae_with_jacobian_python(self):
         model = pybamm.BaseModel()
         model.convert_to_format = "python"