Repository: pybamm-team/PyBaMM

Issue:
Problem with simulating step changes in current
Car current simulation doesn't work. Seems to not pick up the right discontinuities because time is dimensional

Discussion:

Patch:
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 75d7f1427..993d4d37e 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -393,7 +393,10 @@ class Simulation:
             # to correspond to a single discharge
             elif t_eval is None:
                 C_rate = self._parameter_values["C-rate"]
-                t_end = 3600 / C_rate
+                try:
+                    t_end = 3600 / C_rate
+                except TypeError:
+                    t_end = 3600
                 t_eval = np.linspace(0, t_end, 100)
 
             self.t_eval = t_eval
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 50099b7e4..8dd2ce273 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -219,21 +219,31 @@ class BaseSolver(object):
 
         # Check for heaviside functions in rhs and algebraic and add discontinuity
         # events if these exist.
-        # Note: only checks for the case of t < X, t <= X, X < t, or X <= t
-        for symbol in itertools.chain(
-            model.concatenated_rhs.pre_order(), model.concatenated_algebraic.pre_order()
-        ):
-            if isinstance(symbol, pybamm.Heaviside):
-                if symbol.right.id == pybamm.t.id:
-                    expr = symbol.left
-                elif symbol.left.id == pybamm.t.id:
-                    expr = symbol.right
-
-                model.events.append(
-                    pybamm.Event(
-                        str(symbol), expr.new_copy(), pybamm.EventType.DISCONTINUITY
+        # Note: only checks for the case of t < X, t <= X, X < t, or X <= t, but also
+        # accounts for the fact that t might be dimensional
+        # Only do this for DAE models as ODE models can deal with discontinuities fine
+        if len(model.algebraic) > 0:
+            for symbol in itertools.chain(
+                model.concatenated_rhs.pre_order(),
+                model.concatenated_algebraic.pre_order(),
+            ):
+                if isinstance(symbol, pybamm.Heaviside):
+                    # Dimensionless
+                    if symbol.right.id == pybamm.t.id:
+                        expr = symbol.left
+                    elif symbol.left.id == pybamm.t.id:
+                        expr = symbol.right
+                    # Dimensional
+                    elif symbol.right.id == (pybamm.t * model.timescale).id:
+                        expr = symbol.left.new_copy() / symbol.right.right.new_copy()
+                    elif symbol.left.id == (pybamm.t * model.timescale).id:
+                        expr = symbol.right.new_copy() / symbol.left.right.new_copy()
+
+                    model.events.append(
+                        pybamm.Event(
+                            str(symbol), expr.new_copy(), pybamm.EventType.DISCONTINUITY
+                        )
                     )
-                )
 
         # Process rhs, algebraic and event expressions
         rhs, rhs_eval, jac_rhs = process(model.concatenated_rhs, "RHS")
@@ -328,7 +338,7 @@ class BaseSolver(object):
             Initial conditions that are consistent with the algebraic equations (roots
             of the algebraic equations)
         """
-        pybamm.logger.info("Start calculating consistent initial conditions")
+        pybamm.logger.info("Start calculating consistent states")
         rhs = model.rhs_eval
         algebraic = model.algebraic_eval
         jac = model.jac_algebraic_eval
@@ -336,7 +346,7 @@ class BaseSolver(object):
             y0_guess = model.concatenated_initial_conditions.flatten()
 
         # Split y0_guess into differential and algebraic
-        len_rhs = rhs(0, y0_guess).shape[0]
+        len_rhs = rhs(time, y0_guess).shape[0]
         y0_diff, y0_alg_guess = np.split(y0_guess, [len_rhs])
 
         def root_fun(y0_alg):
@@ -344,8 +354,8 @@ class BaseSolver(object):
             y0 = np.concatenate([y0_diff, y0_alg])
             out = algebraic(time, y0)
             pybamm.logger.debug(
-                "Evaluating algebraic equations at t=0, L2-norm is {}".format(
-                    np.linalg.norm(out)
+                "Evaluating algebraic equations at t={}, L2-norm is {}".format(
+                    time * model.timescale, np.linalg.norm(out)
                 )
             )
             return out
@@ -472,43 +482,46 @@ class BaseSolver(object):
 
         # make sure they are increasing in time
         discontinuities = sorted(discontinuities)
-        if len(discontinuities) > 0:
-            pybamm.logger.info(
-                "Discontinuity events found at t = {}".format(discontinuities)
-            )
-        else:
-            pybamm.logger.info("No discontinuity events found")
 
         # remove any identical discontinuities
         discontinuities = [
             v
             for i, v in enumerate(discontinuities)
-            if i == len(discontinuities) - 1
-            or discontinuities[i] < discontinuities[i + 1]
+            if (
+                i == len(discontinuities) - 1
+                or discontinuities[i] < discontinuities[i + 1]
+            )
+            and v > 0
         ]
 
+        if len(discontinuities) > 0:
+            pybamm.logger.info(
+                "Discontinuity events found at t = {}".format(discontinuities)
+            )
+        else:
+            pybamm.logger.info("No discontinuity events found")
+
         # insert time points around discontinuities in t_eval
         # keep track of sub sections to integrate by storing start and end indices
         start_indices = [0]
         end_indices = []
+        eps = sys.float_info.epsilon
         for dtime in discontinuities:
             dindex = np.searchsorted(t_eval_dimensionless, dtime, side="left")
             end_indices.append(dindex + 1)
             start_indices.append(dindex + 1)
-            if t_eval_dimensionless[dindex] == dtime:
-                t_eval_dimensionless[dindex] += sys.float_info.epsilon
+            if dtime - eps < t_eval_dimensionless[dindex] < dtime + eps:
+                t_eval_dimensionless[dindex] += eps
                 t_eval_dimensionless = np.insert(
-                    t_eval_dimensionless, dindex, dtime - sys.float_info.epsilon
+                    t_eval_dimensionless, dindex, dtime - eps
                 )
             else:
                 t_eval_dimensionless = np.insert(
-                    t_eval_dimensionless,
-                    dindex,
-                    [dtime - sys.float_info.epsilon, dtime + sys.float_info.epsilon],
+                    t_eval_dimensionless, dindex, [dtime - eps, dtime + eps],
                 )
         end_indices.append(len(t_eval_dimensionless))
 
-        # integrate separatly over each time segment and accumulate into the solution
+        # integrate separately over each time segment and accumulate into the solution
         # object, restarting the solver at each discontinuity (and recalculating a
         # consistent state afterwards if a dae)
         old_y0 = model.y0
@@ -521,16 +534,13 @@ class BaseSolver(object):
                 )
             )
             timer.reset()
+            new_solution = self._integrate(
+                model, t_eval_dimensionless[start_index:end_index], ext_and_inputs
+            )
+            new_solution.solve_time = timer.time()
             if solution is None:
-                solution = self._integrate(
-                    model, t_eval_dimensionless[start_index:end_index], ext_and_inputs
-                )
-                solution.solve_time = timer.time()
+                solution = new_solution
             else:
-                new_solution = self._integrate(
-                    model, t_eval_dimensionless[start_index:end_index], ext_and_inputs
-                )
-                new_solution.solve_time = timer.time()
                 solution.append(new_solution, start_index=0)
 
             if solution.termination != "final time":
@@ -538,14 +548,6 @@ class BaseSolver(object):
 
             if end_index != len(t_eval_dimensionless):
                 # setup for next integration subsection
-                y0_guess = solution.y[:, -1]
-                if model.algebraic:
-                    model.y0 = self.calculate_consistent_state(
-                        model, t_eval_dimensionless[end_index], y0_guess
-                    )
-                else:
-                    model.y0 = y0_guess
-
                 last_state = solution.y[:, -1]
                 if len(model.algebraic) > 0:
                     model.y0 = self.calculate_consistent_state(
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index 141364ac2..ea809fda7 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -100,13 +100,13 @@ class CasadiSolver(pybamm.BaseSolver):
             return solution
         elif self.mode == "safe":
             # Step-and-check
+            t = t_eval[0]
             init_event_signs = np.sign(
                 np.concatenate(
-                    [event(0, model.y0) for event in model.terminate_events_eval]
+                    [event(t, model.y0) for event in model.terminate_events_eval]
                 )
             )
             pybamm.logger.info("Start solving {} with {}".format(model.name, self.name))
-            t = t_eval[0]
             y0 = model.y0
             # Initialize solution
             solution = pybamm.Solution(np.array([t]), y0[:, np.newaxis])
@@ -144,7 +144,7 @@ class CasadiSolver(pybamm.BaseSolver):
                 new_event_signs = np.sign(
                     np.concatenate(
                         [
-                            event(0, current_step_sol.y[:, -1])
+                            event(t, current_step_sol.y[:, -1])
                             for event in model.terminate_events_eval
                         ]
                     )
@@ -186,15 +186,15 @@ class CasadiSolver(pybamm.BaseSolver):
             # set up and solve
             t = casadi.MX.sym("t")
             u = casadi.MX.sym("u", u_stacked.shape[0])
-            y_diff = casadi.MX.sym("y_diff", rhs(0, y0, u).shape[0])
+            y_diff = casadi.MX.sym("y_diff", rhs(t_eval[0], y0, u).shape[0])
             problem = {"t": t, "x": y_diff, "p": u}
-            if algebraic(0, y0, u).is_empty():
+            if algebraic(t_eval[0], y0, u).is_empty():
                 method = "cvodes"
                 problem.update({"ode": rhs(t, y_diff, u)})
             else:
                 options["calc_ic"] = True
                 method = "idas"
-                y_alg = casadi.MX.sym("y_alg", algebraic(0, y0, u).shape[0])
+                y_alg = casadi.MX.sym("y_alg", algebraic(t_eval[0], y0, u).shape[0])
                 y_full = casadi.vertcat(y_diff, y_alg)
                 problem.update(
                     {
@@ -215,7 +215,7 @@ class CasadiSolver(pybamm.BaseSolver):
         )
 
     def _run_integrator(self, integrator, model, y0, inputs, t_eval):
-        rhs_size = model.rhs_eval(0, y0).shape[0]
+        rhs_size = model.rhs_eval(t_eval[0], y0).shape[0]
         y0_diff, y0_alg = np.split(y0, [rhs_size])
         try:
             # Try solving


Test patch:
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index af9c6a179..0d87c632e 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -385,6 +385,29 @@ class TestSimulation(unittest.TestCase):
         with self.assertWarns(pybamm.SolverWarning):
             sim.solve(t_eval=np.linspace(0, time_data[-1], 800))
 
+    def test_discontinuous_current(self):
+        def car_current(t):
+            current = (
+                1 * (t >= 0) * (t <= 1000)
+                - 0.5 * (1000 < t) * (t <= 2000)
+                + 0.5 * (2000 < t)
+            )
+            return current
+
+        model = pybamm.lithium_ion.DFN()
+        param = model.default_parameter_values
+        param["Current function [A]"] = car_current
+
+        sim = pybamm.Simulation(
+            model, parameter_values=param, solver=pybamm.CasadiSolver(mode="fast")
+        )
+        sim.solve()
+        current = sim.solution["Current [A]"]
+        tau = model.timescale.evaluate()
+        self.assertEqual(current(0), 1)
+        self.assertEqual(current(1500 / tau), -0.5)
+        self.assertEqual(current(3000 / tau), 0.5)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index 3921f2031..62120bcde 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -50,6 +50,7 @@ class TestBaseSolver(unittest.TestCase):
         class ScalarModel:
             concatenated_initial_conditions = np.array([[2]])
             jac_algebraic_eval = None
+            timescale = 1
 
             def rhs_eval(self, t, y):
                 return np.array([])
@@ -68,6 +69,7 @@ class TestBaseSolver(unittest.TestCase):
         class VectorModel:
             concatenated_initial_conditions = np.zeros_like(vec)
             jac_algebraic_eval = None
+            timescale = 1
 
             def rhs_eval(self, t, y):
                 return y[0:1]
@@ -101,6 +103,7 @@ class TestBaseSolver(unittest.TestCase):
         class Model:
             concatenated_initial_conditions = np.array([2])
             jac_algebraic_eval = None
+            timescale = 1
 
             def rhs_eval(self, t, y):
                 return np.array([])
@@ -123,6 +126,16 @@ class TestBaseSolver(unittest.TestCase):
         ):
             solver.calculate_consistent_state(Model())
 
+    def test_time_too_short(self):
+        solver = pybamm.BaseSolver()
+        model = pybamm.BaseModel()
+        v = pybamm.StateVector(slice(0, 1))
+        model.rhs = {v: v}
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "It looks like t_eval might be dimensionless"
+        ):
+            solver.solve(model, np.linspace(0, 0.1))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 7e6974c57..da2e6a3a6 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -259,6 +259,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
         def nonsmooth_mult(t):
             return int(t < discontinuity) + 1.0
+
         rate = pybamm.Function(nonsmooth_rate, pybamm.t)
         mult = pybamm.Function(nonsmooth_mult, pybamm.t)
         model.rhs = {var1: rate * var1}
@@ -267,14 +268,16 @@ class TestScikitsSolvers(unittest.TestCase):
         model.events = [
             pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
             pybamm.Event("var2 = 2.5", pybamm.min(var2 - 2.5)),
-            pybamm.Event("nonsmooth rate",
-                         pybamm.Scalar(discontinuity),
-                         pybamm.EventType.DISCONTINUITY
-                         ),
-            pybamm.Event("nonsmooth mult",
-                         pybamm.Scalar(discontinuity),
-                         pybamm.EventType.DISCONTINUITY
-                         )
+            pybamm.Event(
+                "nonsmooth rate",
+                pybamm.Scalar(discontinuity),
+                pybamm.EventType.DISCONTINUITY,
+            ),
+            pybamm.Event(
+                "nonsmooth mult",
+                pybamm.Scalar(discontinuity),
+                pybamm.EventType.DISCONTINUITY,
+            ),
         ]
         disc = get_discretisation_for_testing()
         disc.process_model(model)
@@ -285,9 +288,9 @@ class TestScikitsSolvers(unittest.TestCase):
         # create two time series, one without a time point on the discontinuity,
         # and one with
         t_eval1 = np.linspace(0, 5, 10)
-        t_eval2 = np.insert(t_eval1,
-                            np.searchsorted(t_eval1, discontinuity),
-                            discontinuity)
+        t_eval2 = np.insert(
+            t_eval1, np.searchsorted(t_eval1, discontinuity), discontinuity
+        )
         solution1 = solver.solve(model, t_eval1)
         solution2 = solver.solve(model, t_eval2)
 
@@ -312,13 +315,13 @@ class TestScikitsSolvers(unittest.TestCase):
             np.testing.assert_array_less(solution.y[-1], 2.5)
             var1_soln = np.exp(0.2 * solution.t)
             y0 = np.exp(0.2 * discontinuity)
-            var1_soln[solution.t > discontinuity] = \
-                y0 * np.exp(
+            var1_soln[solution.t > discontinuity] = y0 * np.exp(
                 0.1 * (solution.t[solution.t > discontinuity] - discontinuity)
             )
             var2_soln = 2 * var1_soln
-            var2_soln[solution.t > discontinuity] = \
-                var1_soln[solution.t > discontinuity]
+            var2_soln[solution.t > discontinuity] = var1_soln[
+                solution.t > discontinuity
+            ]
             np.testing.assert_allclose(solution.y[0], var1_soln, rtol=1e-06)
             np.testing.assert_allclose(solution.y[-1], var2_soln, rtol=1e-06)
 
@@ -592,6 +595,95 @@ class TestScikitsSolvers(unittest.TestCase):
             step_solution.y[-1], 2 * np.exp(0.1 * step_solution.t), decimal=5
         )
 
+    def test_model_solver_dae_nonsmooth(self):
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        var1 = pybamm.Variable("var1", domain=whole_cell)
+        var2 = pybamm.Variable("var2")
+        discontinuity = 0.6
+
+        # Create three different models with the same solution, each expressing the
+        # discontinuity in a different way
+
+        # first model explicitly adds a discontinuity event
+        def nonsmooth_rate(t):
+            return 0.1 * (t < discontinuity) + 0.1
+
+        rate = pybamm.Function(nonsmooth_rate, pybamm.t)
+        model1 = pybamm.BaseModel()
+        model1.rhs = {var1: rate * var1}
+        model1.algebraic = {var2: var2}
+        model1.initial_conditions = {var1: 1, var2: 0}
+        model1.events = [
+            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
+            pybamm.Event(
+                "nonsmooth rate",
+                pybamm.Scalar(discontinuity),
+                pybamm.EventType.DISCONTINUITY,
+            ),
+        ]
+
+        # second model implicitly adds a discontinuity event via a heaviside function
+        model2 = pybamm.BaseModel()
+        model2.rhs = {var1: (0.1 * (pybamm.t < discontinuity) + 0.1) * var1}
+        model2.algebraic = {var2: var2}
+        model2.initial_conditions = {var1: 1, var2: 0}
+        model2.events = [
+            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
+        ]
+
+        # third model implicitly adds a discontinuity event via another heaviside
+        # function
+        model3 = pybamm.BaseModel()
+        model3.rhs = {var1: (-0.1 * (discontinuity < pybamm.t) + 0.2) * var1}
+        model3.algebraic = {var2: var2}
+        model3.initial_conditions = {var1: 1, var2: 0}
+        model3.events = [
+            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
+        ]
+
+        for model in [model1, model2, model3]:
+
+            disc = get_discretisation_for_testing()
+            disc.process_model(model)
+
+            # Solve
+            solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+
+            # create two time series, one without a time point on the discontinuity,
+            # and one with
+            t_eval1 = np.linspace(0, 5, 10)
+            t_eval2 = np.insert(
+                t_eval1, np.searchsorted(t_eval1, discontinuity), discontinuity
+            )
+            solution1 = solver.solve(model, t_eval1)
+            solution2 = solver.solve(model, t_eval2)
+
+            # check time vectors
+            for solution in [solution1, solution2]:
+                # time vectors are ordered
+                self.assertTrue(np.all(solution.t[:-1] <= solution.t[1:]))
+
+                # time value before and after discontinuity is an epsilon away
+                dindex = np.searchsorted(solution.t, discontinuity)
+                value_before = solution.t[dindex - 1]
+                value_after = solution.t[dindex]
+                self.assertEqual(value_before + sys.float_info.epsilon, discontinuity)
+                self.assertEqual(value_after - sys.float_info.epsilon, discontinuity)
+
+            # both solution time vectors should have same number of points
+            self.assertEqual(len(solution1.t), len(solution2.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)
+                var1_soln = np.exp(0.2 * solution.t)
+                y0 = np.exp(0.2 * discontinuity)
+                var1_soln[solution.t > discontinuity] = y0 * np.exp(
+                    0.1 * (solution.t[solution.t > discontinuity] - discontinuity)
+                )
+                np.testing.assert_allclose(solution.y[0], var1_soln, rtol=1e-06)
+
     def test_ode_solver_fail_with_dae(self):
         model = pybamm.BaseModel()
         a = pybamm.Scalar(1)
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index dbd2aabc0..fe064b64c 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -6,7 +6,6 @@ import unittest
 import numpy as np
 from tests import get_mesh_for_testing
 import warnings
-from tests import get_discretisation_for_testing
 import sys
 
 
@@ -138,91 +137,6 @@ class TestScipySolver(unittest.TestCase):
             np.ones((N, T.size)) * (T[np.newaxis, :] - np.exp(T[np.newaxis, :])),
         )
 
-    def test_model_solver_ode_nonsmooth(self):
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
-        var1 = pybamm.Variable("var1", domain=whole_cell)
-        discontinuity = 0.6
-
-        # Create three different models with the same solution, each expressing the
-        # discontinuity in a different way
-
-        # first model explicitly adds a discontinuity event
-        def nonsmooth_rate(t):
-            return 0.1 * (t < discontinuity) + 0.1
-
-        rate = pybamm.Function(nonsmooth_rate, pybamm.t)
-        model1 = pybamm.BaseModel()
-        model1.rhs = {var1: rate * var1}
-        model1.initial_conditions = {var1: 1}
-        model1.events = [
-            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
-            pybamm.Event(
-                "nonsmooth rate",
-                pybamm.Scalar(discontinuity),
-                pybamm.EventType.DISCONTINUITY,
-            ),
-        ]
-
-        # second model implicitly adds a discontinuity event via a heaviside function
-        model2 = pybamm.BaseModel()
-        model2.rhs = {var1: (0.1 * (pybamm.t < discontinuity) + 0.1) * var1}
-        model2.initial_conditions = {var1: 1}
-        model2.events = [
-            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
-        ]
-
-        # third model implicitly adds a discontinuity event via another heaviside
-        # function
-        model3 = pybamm.BaseModel()
-        model3.rhs = {var1: (-0.1 * (discontinuity < pybamm.t) + 0.2) * var1}
-        model3.initial_conditions = {var1: 1}
-        model3.events = [
-            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
-        ]
-
-        for model in [model1, model2, model3]:
-
-            disc = get_discretisation_for_testing()
-            disc.process_model(model)
-
-            # Solve
-            solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8)
-
-            # create two time series, one without a time point on the discontinuity,
-            # and one with
-            t_eval1 = np.linspace(0, 5, 10)
-            t_eval2 = np.insert(
-                t_eval1, np.searchsorted(t_eval1, discontinuity), discontinuity
-            )
-            solution1 = solver.solve(model, t_eval1)
-            solution2 = solver.solve(model, t_eval2)
-
-            # check time vectors
-            for solution in [solution1, solution2]:
-                # time vectors are ordered
-                self.assertTrue(np.all(solution.t[:-1] <= solution.t[1:]))
-
-                # time value before and after discontinuity is an epsilon away
-                dindex = np.searchsorted(solution.t, discontinuity)
-                value_before = solution.t[dindex - 1]
-                value_after = solution.t[dindex]
-                self.assertEqual(value_before + sys.float_info.epsilon, discontinuity)
-                self.assertEqual(value_after - sys.float_info.epsilon, discontinuity)
-
-            # both solution time vectors should have same number of points
-            self.assertEqual(len(solution1.t), len(solution2.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)
-                var1_soln = np.exp(0.2 * solution.t)
-                y0 = np.exp(0.2 * discontinuity)
-                var1_soln[solution.t > discontinuity] = y0 * np.exp(
-                    0.1 * (solution.t[solution.t > discontinuity] - discontinuity)
-                )
-                np.testing.assert_allclose(solution.y[0], var1_soln, rtol=1e-06)
-
     def test_model_step_python(self):
         # Create model
         model = pybamm.BaseModel()