Repository: pybamm-team/PyBaMM (python)

Issue:
Initialize with casadi
**Summary**
Use casadi root-finder to initialize DAEs

 **Motivation**
 Should be faster




Pull request description:
Issue 804 initialize casadi
# Description

Add the (now default) option to use casadi's rootfinder to initialize the DAEs

Fixes #804 

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/master/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [x] New feature (non-breaking change which adds functionality)
- [ ] Optimization (back-end change that speeds up the code)
- [ ] Bug fix (non-breaking change which fixes an issue)


# Key checklist:

- [x] No style issues: `$ flake8`
- [x] All tests pass: `$ python run-tests.py --unit`
- [x] The documentation builds: `$ cd docs` and then `$ make clean; make html`

You can run all three at once, using `$ python run-tests.py --quick`.

## Further checks:

- [x] Code is commented, particularly in hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3e6cf962d..6f2e22249 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Changed rootfinding algorithm to CasADi, scipy.optimize.root still accessible as an option ([#844](https://github.com/pybamm-team/PyBaMM/pull/844))
 -   Added NCA parameter set ([#824](https://github.com/pybamm-team/PyBaMM/pull/824))
 -   Added functionality to `Solution` that automatically gets `t_eval` from the data when simulating drive cycles and performs checks to ensure the output has the required resolution to accurately capture the input current ([#819](https://github.com/pybamm-team/PyBaMM/pull/819))
 -   Added options to export a solution to matlab or csv ([#811](https://github.com/pybamm-team/PyBaMM/pull/811))
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
index fd01abfa2..16d17eef2 100644
--- a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
@@ -177,17 +177,11 @@ class BasicDFN(BaseModel):
         # Boundary conditions must be provided for equations with spatial derivatives
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (
-                -param.C_n * j_n / param.a_n / param.D_n(c_s_surf_n, T),
-                "Neumann",
-            ),
+            "right": (-param.C_n * j_n / param.a_n, "Neumann"),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (
-                -param.C_p * j_p / param.a_p / param.gamma_p / param.D_p(c_s_surf_p, T),
-                "Neumann",
-            ),
+            "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"),
         }
         # c_n_init and c_p_init can in general be functions of x
         # Note the broadcasting, for domains
@@ -201,22 +195,14 @@ class BasicDFN(BaseModel):
         self.initial_conditions[c_s_p] = param.c_p_init(x_p)
         # Events specify points at which a solution should terminate
         self.events += [
-            pybamm.Event(
-                "Minimum negative particle surface concentration",
-                pybamm.min(c_s_surf_n) - 0.01,
-            ),
-            pybamm.Event(
-                "Maximum negative particle surface concentration",
-                (1 - 0.01) - pybamm.max(c_s_surf_n),
-            ),
-            pybamm.Event(
-                "Minimum positive particle surface concentration",
-                pybamm.min(c_s_surf_p) - 0.01,
-            ),
-            pybamm.Event(
-                "Maximum positive particle surface concentration",
-                (1 - 0.01) - pybamm.max(c_s_surf_p),
-            ),
+            pybamm.Event("Minimum negative particle surface concentration",
+                         pybamm.min(c_s_surf_n) - 0.01),
+            pybamm.Event("Maximum negative particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_n)),
+            pybamm.Event("Minimum positive particle surface concentration",
+                         pybamm.min(c_s_surf_p) - 0.01),
+            pybamm.Event("Maximum positive particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_p)),
         ]
         ######################
         # Current in the solid
@@ -271,11 +257,8 @@ class BasicDFN(BaseModel):
             "right": (pybamm.Scalar(0), "Neumann"),
         }
         self.initial_conditions[c_e] = param.c_e_init
-        self.events.append(
-            pybamm.Event(
-                "Zero electrolyte concentration cut-off", pybamm.min(c_e) - 0.002
-            )
-        )
+        self.events.append(pybamm.Event("Zero electrolyte concentration cut-off",
+                                        pybamm.min(c_e) - 0.002))
 
         ######################
         # (Some) variables
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
index 60956d210..15d8e2ab9 100644
--- a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
@@ -80,48 +80,34 @@ class BasicSPM(BaseModel):
         N_s_p = -param.D_p(c_s_p, T) * pybamm.grad(c_s_p)
         self.rhs[c_s_n] = -(1 / param.C_n) * pybamm.div(N_s_n)
         self.rhs[c_s_p] = -(1 / param.C_p) * pybamm.div(N_s_p)
-        # Surf takes the surface value of a variable, i.e. its boundary value on the
-        # right side. This is also accessible via `boundary_value(x, "right")`, with
-        # "left" providing the boundary value of the left side
-        c_s_surf_n = pybamm.surf(c_s_n)
-        c_s_surf_p = pybamm.surf(c_s_p)
         # Boundary conditions must be provided for equations with spatial derivatives
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (
-                -param.C_n * j_n / param.a_n / param.D_n(c_s_surf_n, T),
-                "Neumann",
-            ),
+            "right": (-param.C_n * j_n / param.a_n, "Neumann"),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (
-                -param.C_p * j_p / param.a_p / param.gamma_p / param.D_p(c_s_surf_p, T),
-                "Neumann",
-            ),
+            "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"),
         }
         # c_n_init and c_p_init are functions, but for the SPM we evaluate them at x=0
         # and x=1 since there is no x-dependence in the particles
         self.initial_conditions[c_s_n] = param.c_n_init(0)
         self.initial_conditions[c_s_p] = param.c_p_init(1)
+        # Surf takes the surface value of a variable, i.e. its boundary value on the
+        # right side. This is also accessible via `boundary_value(x, "right")`, with
+        # "left" providing the boundary value of the left side
+        c_s_surf_n = pybamm.surf(c_s_n)
+        c_s_surf_p = pybamm.surf(c_s_p)
         # Events specify points at which a solution should terminate
         self.events += [
-            pybamm.Event(
-                "Minimum negative particle surface concentration",
-                pybamm.min(c_s_surf_n) - 0.01,
-            ),
-            pybamm.Event(
-                "Maximum negative particle surface concentration",
-                (1 - 0.01) - pybamm.max(c_s_surf_n),
-            ),
-            pybamm.Event(
-                "Minimum positive particle surface concentration",
-                pybamm.min(c_s_surf_p) - 0.01,
-            ),
-            pybamm.Event(
-                "Maximum positive particle surface concentration",
-                (1 - 0.01) - pybamm.max(c_s_surf_p),
-            ),
+            pybamm.Event("Minimum negative particle surface concentration",
+                         pybamm.min(c_s_surf_n) - 0.01),
+            pybamm.Event("Maximum negative particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_n)),
+            pybamm.Event("Minimum positive particle surface concentration",
+                         pybamm.min(c_s_surf_p) - 0.01),
+            pybamm.Event("Maximum positive particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_p)),
         ]
 
         # Note that the SPM does not have any algebraic equations, so the `algebraic`
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 1c031b8b8..806015e90 100644
--- a/pybamm/solvers/base_solver.py
+++ b/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
@@ -36,7 +39,7 @@ class BaseSolver(object):
         method=None,
         rtol=1e-6,
         atol=1e-6,
-        root_method="lm",
+        root_method="casadi",
         root_tol=1e-6,
         max_steps=1000,
     ):
@@ -125,10 +128,11 @@ class BaseSolver(object):
                 "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)
-            and model.convert_to_format != "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"
             )
@@ -268,6 +272,16 @@ class BaseSolver(object):
         model.terminate_events_eval = terminate_events_eval
         model.discontinuity_events_eval = discontinuity_events_eval
 
+        # Save CasADi functions for the CasADi solver
+        # Note: when we pass to casadi the ode part of the problem must be in explicit
+        # form so we pre-multiply by the inverse of the mass matrix
+        if self.root_method == "casadi" or isinstance(self, pybamm.CasadiSolver):
+            mass_matrix_inv = casadi.MX(model.mass_matrix_inv.entries)
+            explicit_rhs = mass_matrix_inv @ rhs(t_casadi, y_casadi, u_casadi_stacked)
+            model.casadi_rhs = casadi.Function(
+                "rhs", [t_casadi, y_casadi, u_casadi_stacked], [explicit_rhs]
+            )
+            model.casadi_algebraic = algebraic
         # Calculate consistent initial conditions for the algebraic equations
         if len(model.algebraic) > 0:
             all_states = pybamm.NumpyConcatenation(
@@ -278,26 +292,13 @@ class BaseSolver(object):
             model.residuals_eval = residuals_eval
             model.jacobian_eval = jacobian_eval
             y0_guess = y0.flatten()
-            model.y0 = self.calculate_consistent_state(model, 0, y0_guess)
+            model.y0 = self.calculate_consistent_state(model, 0, y0_guess, inputs)
         else:
             # can use DAE solver to solve ODE model
             model.residuals_eval = Residuals(rhs, "residuals", model)
             model.jacobian_eval = jac_rhs
             model.y0 = y0.flatten()
 
-        # Save CasADi functions for the CasADi solver
-        # Note: when we pass to casadi the ode part of the problem must be in explicit
-        # form so we pre-multiply by the inverse of the mass matrix
-        if model.convert_to_format == "casadi" and isinstance(
-            self, pybamm.CasadiSolver
-        ):
-            mass_matrix_inv = casadi.MX(model.mass_matrix_inv.entries)
-            explicit_rhs = mass_matrix_inv @ rhs(t_casadi, y_casadi, u_casadi_stacked)
-            model.casadi_rhs = casadi.Function(
-                "rhs", [t_casadi, y_casadi, u_casadi_stacked], [explicit_rhs]
-            )
-            model.casadi_algebraic = algebraic
-
         pybamm.logger.info("Finish solver set-up")
 
     def set_inputs(self, model, ext_and_inputs):
@@ -319,7 +320,7 @@ class BaseSolver(object):
         if model.jacobian_eval:
             model.jacobian_eval.set_inputs(ext_and_inputs)
 
-    def calculate_consistent_state(self, model, time=0, y0_guess=None):
+    def calculate_consistent_state(self, model, time=0, y0_guess=None, inputs=None):
         """
         Calculate consistent state for the algebraic equations through
         root-finding
@@ -328,6 +329,12 @@ class BaseSolver(object):
         ----------
         model : :class:`pybamm.BaseModel`
             The model for which to calculate initial conditions.
+        time : float
+            The time at which to calculate the states
+        y0_guess : :class:`np.array`
+            Guess for the rootfinding
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         Returns
         -------
@@ -336,65 +343,101 @@ class BaseSolver(object):
             of the algebraic equations)
         """
         pybamm.logger.info("Start calculating consistent states")
-        rhs = model.rhs_eval
-        algebraic = model.algebraic_eval
-        jac = model.jac_algebraic_eval
         if y0_guess is None:
             y0_guess = model.concatenated_initial_conditions.flatten()
 
         # Split y0_guess into differential and algebraic
-        len_rhs = rhs(time, y0_guess).shape[0]
+        len_rhs = model.rhs_eval(time, y0_guess).shape[0]
         y0_diff, y0_alg_guess = np.split(y0_guess, [len_rhs])
+        inputs = inputs or {}
 
-        def root_fun(y0_alg):
-            "Evaluates algebraic using y0_diff (fixed) and y0_alg (changed by algo)"
-            y0 = np.concatenate([y0_diff, y0_alg])
-            out = algebraic(time, y0)
-            pybamm.logger.debug(
-                "Evaluating algebraic equations at t={}, L2-norm is {}".format(
-                    time * model.timescale, np.linalg.norm(out)
-                )
+        # Solve using casadi or scipy
+        if self.root_method == "casadi":
+            # Set up
+            u_stacked = casadi.vertcat(*[x for x in inputs.values()])
+            u = casadi.MX.sym("u", u_stacked.shape[0])
+            y_alg = casadi.MX.sym("y_alg", y0_alg_guess.shape[0])
+            y = casadi.vertcat(y0_diff, y_alg)
+            alg_root = model.casadi_algebraic(time, y, u)
+            # Solve
+            # 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),
+                {"abstol": self.root_tol},
             )
-            return out
-
-        if jac:
-            if issparse(jac(0, y0_guess)):
+            try:
+                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
+
+            def root_fun(y0_alg):
+                "Evaluates algebraic using y0_diff (fixed) and y0_alg (changed by algo)"
+                y0 = np.concatenate([y0_diff, y0_alg])
+                out = algebraic(time, y0)
+                pybamm.logger.debug(
+                    "Evaluating algebraic equations at t={}, L2-norm is {}".format(
+                        time * model.timescale, np.linalg.norm(out)
+                    )
+                )
+                return out
 
-                def jac_fn(y0_alg):
-                    """
-                    Evaluates jacobian using y0_diff (fixed) and y0_alg (varying)
-                    """
-                    y0 = np.concatenate([y0_diff, y0_alg])
-                    return jac(0, y0)[:, len_rhs:].toarray()
+            if jac:
+                if issparse(jac(0, y0_guess)):
 
-            else:
+                    def jac_fn(y0_alg):
+                        """
+                        Evaluates jacobian using y0_diff (fixed) and y0_alg (varying)
+                        """
+                        y0 = np.concatenate([y0_diff, y0_alg])
+                        return jac(0, y0)[:, len_rhs:].toarray()
 
-                def jac_fn(y0_alg):
-                    """
-                    Evaluates jacobian using y0_diff (fixed) and y0_alg (varying)
-                    """
-                    y0 = np.concatenate([y0_diff, y0_alg])
-                    return jac(0, y0)[:, len_rhs:]
+                else:
 
-        else:
-            jac_fn = None
-        # Find the values of y0_alg that are roots of the algebraic equations
-        sol = optimize.root(
-            root_fun,
-            y0_alg_guess,
-            jac=jac_fn,
-            method=self.root_method,
-            tol=self.root_tol,
-        )
-        # Return full set of consistent initial conditions (y0_diff unchanged)
-        y0_consistent = np.concatenate([y0_diff, sol.x])
+                    def jac_fn(y0_alg):
+                        """
+                        Evaluates jacobian using y0_diff (fixed) and y0_alg (varying)
+                        """
+                        y0 = np.concatenate([y0_diff, y0_alg])
+                        return jac(0, y0)[:, len_rhs:]
 
-        if sol.success and np.all(sol.fun < self.root_tol * len(sol.x)):
+            else:
+                jac_fn = None
+            # Find the values of y0_alg that are roots of the algebraic equations
+            sol = optimize.root(
+                root_fun,
+                y0_alg_guess,
+                jac=jac_fn,
+                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, y0_alg])
             pybamm.logger.info("Finish calculating consistent initial conditions")
             return y0_consistent
-        elif not sol.success:
+        elif not success:
             raise pybamm.SolverError(
-                "Could not find consistent initial conditions: {}".format(sol.message)
+                "Could not find consistent initial conditions: {}".format(message)
             )
         else:
             raise pybamm.SolverError(
@@ -402,7 +445,7 @@ class BaseSolver(object):
                 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)
+                    np.max(fun), self.root_tol * len(y0_alg)
                 )
             )
 
@@ -514,7 +557,7 @@ class BaseSolver(object):
                 )
             else:
                 t_eval_dimensionless = np.insert(
-                    t_eval_dimensionless, dindex, [dtime - eps, dtime + eps],
+                    t_eval_dimensionless, dindex, [dtime - eps, dtime + eps]
                 )
         end_indices.append(len(t_eval_dimensionless))
 
@@ -548,7 +591,10 @@ class BaseSolver(object):
                 last_state = solution.y[:, -1]
                 if len(model.algebraic) > 0:
                     model.y0 = self.calculate_consistent_state(
-                        model, t_eval_dimensionless[end_index], last_state
+                        model,
+                        t_eval_dimensionless[end_index],
+                        last_state,
+                        ext_and_inputs,
                     )
                 else:
                     model.y0 = last_state
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index ea809fda7..34fd74043 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -47,7 +47,7 @@ class CasadiSolver(pybamm.BaseSolver):
         mode="safe",
         rtol=1e-6,
         atol=1e-6,
-        root_method="lm",
+        root_method="casadi",
         root_tol=1e-6,
         max_step_decrease_count=5,
         **extra_options,
diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py
index 795f3361b..be4a1932c 100644
--- a/pybamm/solvers/idaklu_solver.py
+++ b/pybamm/solvers/idaklu_solver.py
@@ -36,7 +36,7 @@ class IDAKLUSolver(pybamm.BaseSolver):
     """
 
     def __init__(
-        self, rtol=1e-6, atol=1e-6, root_method="lm", root_tol=1e-6, max_steps=1000
+        self, rtol=1e-6, atol=1e-6, root_method="casadi", root_tol=1e-6, max_steps=1000
     ):
 
         if idaklu_spec is None:
diff --git a/pybamm/solvers/scikits_dae_solver.py b/pybamm/solvers/scikits_dae_solver.py
index 5f9999cbb..74620cafe 100644
--- a/pybamm/solvers/scikits_dae_solver.py
+++ b/pybamm/solvers/scikits_dae_solver.py
@@ -40,7 +40,7 @@ class ScikitsDaeSolver(pybamm.BaseSolver):
         method="ida",
         rtol=1e-6,
         atol=1e-6,
-        root_method="lm",
+        root_method="casadi",
         root_tol=1e-6,
         max_steps=1000,
     ):
diff --git a/setup.py b/setup.py
index db141d4f4..fa13bd203 100644
--- a/setup.py
+++ b/setup.py
@@ -492,7 +492,7 @@ setup(
     # List of dependencies
     install_requires=[
         "numpy>=1.16",
-        "scipy>=1.0",
+        "scipy>=1.3",
         "pandas>=0.24",
         "anytree>=2.4.3",
         "autograd>=1.2",


Test patch:
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index ce44197a3..7ab105982 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/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,9 +50,16 @@ class TestBaseSolver(unittest.TestCase):
     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([])
@@ -59,18 +67,30 @@ class TestBaseSolver(unittest.TestCase):
             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", root_tol=1e-12)
+        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 @@ class TestBaseSolver(unittest.TestCase):
         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 @@ class TestBaseSolver(unittest.TestCase):
 
     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,12 +150,19 @@ class TestBaseSolver(unittest.TestCase):
             "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: .../casadi",
+        ):
+            solver.calculate_consistent_state(Model())
 
     def test_time_too_short(self):
         solver = pybamm.BaseSolver()
diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py
index 9683b5cbc..703673ae8 100644
--- a/tests/unit/test_solvers/test_idaklu_solver.py
+++ b/tests/unit/test_solvers/test_idaklu_solver.py
@@ -28,7 +28,7 @@ class TestIDAKLUSolver(unittest.TestCase):
             disc = pybamm.Discretisation()
             disc.process_model(model)
 
-            solver = pybamm.IDAKLUSolver()
+            solver = pybamm.IDAKLUSolver(root_method="lm")
 
             t_eval = np.linspace(0, 3, 100)
             solution = solver.solve(model, t_eval)
@@ -58,7 +58,7 @@ class TestIDAKLUSolver(unittest.TestCase):
         mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
         disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
         disc.process_model(model)
-        solver = pybamm.IDAKLUSolver()
+        solver = pybamm.IDAKLUSolver(root_method="lm")
 
         variable_tols = {"Electrolyte concentration": 1e-3}
         solver.set_atol_by_variable(variable_tols, model)
@@ -76,7 +76,7 @@ class TestIDAKLUSolver(unittest.TestCase):
         disc = pybamm.Discretisation()
         disc.process_model(model)
 
-        solver = pybamm.IDAKLUSolver()
+        solver = pybamm.IDAKLUSolver(root_method="lm")
 
         t_eval = np.linspace(0, 3, 100)
         with self.assertRaisesRegex(pybamm.SolverError, "KLU requires the Jacobian"):
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index e3b39791a..905507c69 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -194,7 +194,7 @@ class TestScikitsSolvers(unittest.TestCase):
         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 @@ class TestScikitsSolvers(unittest.TestCase):
         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 @@ class TestScikitsSolvers(unittest.TestCase):
         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 @@ class TestScikitsSolvers(unittest.TestCase):
         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 @@ class TestScikitsSolvers(unittest.TestCase):
 
         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 @@ class TestScikitsSolvers(unittest.TestCase):
         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 @@ class TestScikitsSolvers(unittest.TestCase):
         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 @@ class TestScikitsSolvers(unittest.TestCase):
             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)