#804 fix tests

Author: valentinsulzer

pybamm/solvers/base_solver.py

@@ -23,7 +23,10 @@ class BaseSolver(object):
     atol : float, optional
         The absolute tolerance for the solver (default is 1e-6).
     root_method : str, optional
-        The method to use to find initial conditions (default is "lm")
+        The method to use to find initial conditions (default is "casadi"). If "casadi",
+        the solver uses casadi's Newton rootfinding algorithm to find initial
+        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-6).
     max_steps: int, optional
@@ -125,8 +128,11 @@ def set_up(self, model, inputs=None):
                 "Cannot use ODE solver '{}' to solve DAE model".format(self.name)
             )
 
+        if self.ode_solver is True:
+            self.root_method = None
         if (
-            isinstance(self, pybamm.CasadiSolver) or self.root_method == "casadi"
+            isinstance(self, pybamm.CasadiSolver)
+            or self.root_method == "casadi"
         ) and model.convert_to_format != "casadi":
             pybamm.logger.warning(
                 f"Converting {model.name} to CasADi for solving with CasADi solver"
@@ -349,17 +355,32 @@ def calculate_consistent_state(self, model, time=0, y0_guess=None, inputs=None):
         # Solve using casadi or scipy
         if self.root_method == "casadi":
             # Set up
-            print("yp")
             u_stacked = casadi.vertcat(*[x for x in inputs.values()])
             u = casadi.MX.sym("u", u_stacked.shape[0])
-            alg = model.casadi_algebraic
             y_alg = casadi.MX.sym("y_alg", y0_alg_guess.shape[0])
             y = casadi.vertcat(y0_diff, y_alg)
-            alg_root = alg(time, y, u)
+            alg_root = model.casadi_algebraic(time, y, u)
             # Solve
-            roots = casadi.rootfinder("roots", "newton", dict(x=y_alg, p=u, g=alg_root))
-            y0_alg = roots(y0_alg_guess, u_stacked).full().flatten()
-            return np.concatenate([y0_diff, y0_alg])
+            try:
+                # set error_on_fail to False and just check the final output is small
+                # enough
+                roots = casadi.rootfinder(
+                    "roots",
+                    "newton",
+                    dict(x=y_alg, p=u, g=alg_root),
+                    {"error_on_fail": False},
+                )
+                y0_alg = roots(y0_alg_guess, u_stacked).full().flatten()
+                success = True
+                message = None
+                # Check final output
+                fun = model.casadi_algebraic(
+                    time, casadi.vertcat(y0_diff, y0_alg), u_stacked
+                )
+            except RuntimeError as err:
+                success = False
+                message = err.args[0]
+                fun = None
         else:
             algebraic = model.algebraic_eval
             jac = model.jac_algebraic_eval
@@ -404,27 +425,30 @@ def jac_fn(y0_alg):
                 method=self.root_method,
                 tol=self.root_tol,
             )
+            # Set outputs
+            y0_alg = sol.x
+            success = sol.success
+            fun = sol.fun
+            message = sol.message
+
+        if success and np.all(fun < self.root_tol * len(y0_alg)):
             # Return full set of consistent initial conditions (y0_diff unchanged)
-            y0_consistent = np.concatenate([y0_diff, sol.x])
-
-            if sol.success and np.all(sol.fun < self.root_tol * len(sol.x)):
-                pybamm.logger.info("Finish calculating consistent initial conditions")
-                return y0_consistent
-            elif not sol.success:
-                raise pybamm.SolverError(
-                    "Could not find consistent initial conditions: {}".format(
-                        sol.message
-                    )
-                )
-            else:
-                raise pybamm.SolverError(
-                    """
-                    Could not find consistent initial conditions: solver terminated
-                    successfully, but maximum solution error ({}) above tolerance ({})
-                    """.format(
-                        np.max(sol.fun), self.root_tol * len(sol.x)
-                    )
+            y0_consistent = np.concatenate([y0_diff, y0_alg])
+            pybamm.logger.info("Finish calculating consistent initial conditions")
+            return y0_consistent
+        elif not success:
+            raise pybamm.SolverError(
+                "Could not find consistent initial conditions: {}".format(message)
+            )
+        else:
+            raise pybamm.SolverError(
+                """
+                Could not find consistent initial conditions: solver terminated
+                successfully, but maximum solution error ({}) above tolerance ({})
+                """.format(
+                    np.max(fun), self.root_tol * len(y0_alg)
                 )
+            )
 
     def solve(self, model, t_eval, external_variables=None, inputs=None):
         """

setup.py

@@ -492,7 +492,7 @@ def load_version():
     # List of dependencies
     install_requires=[
         "numpy>=1.16",
-        "scipy>=1.0",
+        "scipy>=1.3",
         "pandas>=0.24",
         "anytree>=2.4.3",
         "autograd>=1.2",

tests/unit/test_solvers/test_base_solver.py

@@ -1,6 +1,7 @@
 #
 # Tests for the Base Solver class
 #
+import casadi
 import pybamm
 import numpy as np
 from scipy.sparse import csr_matrix
@@ -49,28 +50,47 @@ def test_ode_solver_fail_with_dae(self):
     def test_find_consistent_initial_conditions(self):
         # Simple system: a single algebraic equation
         class ScalarModel:
-            concatenated_initial_conditions = np.array([[2]])
-            jac_algebraic_eval = None
-            timescale = 1
+            def __init__(self):
+                self.concatenated_initial_conditions = np.array([[2]])
+                self.jac_algebraic_eval = None
+                self.timescale = 1
+                t = casadi.MX.sym("t")
+                y = casadi.MX.sym("y")
+                u = casadi.MX.sym("u")
+                self.casadi_algebraic = casadi.Function(
+                    "alg", [t, y, u], [self.algebraic_eval(t, y)]
+                )
 
             def rhs_eval(self, t, y):
                 return np.array([])
 
             def algebraic_eval(self, t, y):
                 return y + 2
 
-        solver = pybamm.BaseSolver()
+        solver = pybamm.BaseSolver(root_method="lm")
         model = ScalarModel()
         init_cond = solver.calculate_consistent_state(model)
         np.testing.assert_array_equal(init_cond, -2)
+        # with casadi
+        solver_with_casadi = pybamm.BaseSolver(root_method="casadi")
+        model = ScalarModel()
+        init_cond = solver_with_casadi.calculate_consistent_state(model)
+        np.testing.assert_array_equal(init_cond, -2)
 
         # More complicated system
         vec = np.array([0.0, 1.0, 1.5, 2.0])
 
         class VectorModel:
-            concatenated_initial_conditions = np.zeros_like(vec)
-            jac_algebraic_eval = None
-            timescale = 1
+            def __init__(self):
+                self.concatenated_initial_conditions = np.zeros_like(vec)
+                self.jac_algebraic_eval = None
+                self.timescale = 1
+                t = casadi.MX.sym("t")
+                y = casadi.MX.sym("y", vec.size)
+                u = casadi.MX.sym("u")
+                self.casadi_algebraic = casadi.Function(
+                    "alg", [t, y, u], [self.algebraic_eval(t, y)]
+                )
 
             def rhs_eval(self, t, y):
                 return y[0:1]
@@ -81,6 +101,9 @@ def algebraic_eval(self, t, y):
         model = VectorModel()
         init_cond = solver.calculate_consistent_state(model)
         np.testing.assert_array_almost_equal(init_cond, vec)
+        # with casadi
+        init_cond = solver_with_casadi.calculate_consistent_state(model)
+        np.testing.assert_array_almost_equal(init_cond, vec)
 
         # With jacobian
         def jac_dense(t, y):
@@ -102,9 +125,16 @@ def jac_sparse(t, y):
 
     def test_fail_consistent_initial_conditions(self):
         class Model:
-            concatenated_initial_conditions = np.array([2])
-            jac_algebraic_eval = None
-            timescale = 1
+            def __init__(self):
+                self.concatenated_initial_conditions = np.array([2])
+                self.jac_algebraic_eval = None
+                self.timescale = 1
+                t = casadi.MX.sym("t")
+                y = casadi.MX.sym("y")
+                u = casadi.MX.sym("u")
+                self.casadi_algebraic = casadi.Function(
+                    "alg", [t, y, u], [self.algebraic_eval(t, y)]
+                )
 
             def rhs_eval(self, t, y):
                 return np.array([])
@@ -120,7 +150,14 @@ def algebraic_eval(self, t, y):
             "Could not find consistent initial conditions: The iteration is not making",
         ):
             solver.calculate_consistent_state(Model())
-        solver = pybamm.BaseSolver()
+        solver = pybamm.BaseSolver(root_method="lm")
+        with self.assertRaisesRegex(
+            pybamm.SolverError,
+            "Could not find consistent initial conditions: solver terminated",
+        ):
+            solver.calculate_consistent_state(Model())
+        # with casadi
+        solver = pybamm.BaseSolver(root_method="casadi")
         with self.assertRaisesRegex(
             pybamm.SolverError,
             "Could not find consistent initial conditions: solver terminated",

tests/unit/test_solvers/test_scikits_solvers.py

@@ -194,7 +194,7 @@ def test_model_solver_dae_python(self):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -214,7 +214,7 @@ def test_model_solver_dae_bad_ics_python(self):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -238,7 +238,7 @@ def test_model_solver_dae_events_python(self):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        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)
@@ -283,7 +283,7 @@ def nonsmooth_mult(t):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
 
         # create two time series, one without a time point on the discontinuity,
         # and one with
@@ -355,7 +355,7 @@ def jacobian(t, y):
 
         model.jacobian = jacobian
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -372,7 +372,7 @@ def test_solve_ode_model_with_dae_solver_python(self):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -421,7 +421,7 @@ def test_model_step_dae_python(self):
         disc = get_discretisation_for_testing()
         disc.process_model(model)
 
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
 
         # Step once
         dt = 1
@@ -514,7 +514,10 @@ def test_model_solver_dae_inputs_events(self):
             disc.process_model(model)
 
             # Solve
-            solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+            if form == "python":
+                solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
+            else:
+                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)