From fba12947a04ca0fbe33d0d98cd38ef6a822f9132 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Tue, 4 Feb 2020 16:09:11 +0000
Subject: [PATCH] #759 add nonsmooth test case

---
 .../unit/test_solvers/test_scikits_solvers.py | 39 +++++++++++++++++++
 1 file changed, 39 insertions(+)

diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 8227b3ec4a..53ef1dcba6 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -8,6 +8,9 @@
 import warnings
 from tests import get_mesh_for_testing, get_discretisation_for_testing
 
+# TODO: remove this
+import matplotlib.pylab as plt
+
 
 @unittest.skipIf(not pybamm.have_scikits_odes(), "scikits.odes not installed")
 class TestScikitsSolvers(unittest.TestCase):
@@ -548,6 +551,42 @@ def test_model_solver_dae_events_python(self):
         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_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)
+
+        def nonsmooth_rate(t):
+            return 0.1 * int(t < 2.5) + 0.1
+        rate = pybamm.Function(nonsmooth_rate, pybamm.t)
+        model.rhs = {var1: rate * var1}
+        model.algebraic = {var2: 2 * var1 - var2}
+        model.initial_conditions = {var1: 1, var2: 2}
+        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(2.5),
+                         pybamm.EventType.DISCONTINUITY
+                         )
+        ]
+        disc = get_discretisation_for_testing()
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        t_eval = np.linspace(0, 5, 100)
+        solution = solver.solve(model, t_eval)
+        plt.plot(solution.y[0])
+        plt.plot(solution.y[1])
+        plt.show()
+        #np.testing.assert_array_less(solution.y[0], 1.5)
+        #np.testing.assert_array_less(solution.y[-1], 2.5)
+        #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"