Repository: pybamm-team/PyBaMM

Issue:
Solving algebraic models
Fix a few small bugs that make it harder to solve algebraic models:
- casadi solver complains because mass_matrix_inv is empty
- solution complains because len(t) = 1

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9600e4e69..d4a0c03eb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,7 +1,9 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM)
+
 ## Features 
 
-- Add ambient temperature as a function of time ([#872](https://github.com/pybamm-team/PyBaMM/pull/872))
+-   Add ambient temperature as a function of time ([#872](https://github.com/pybamm-team/PyBaMM/pull/872))
+-   Added `CasadiAlgebraicSolver` for solving algebraic systems with casadi ([#868](https://github.com/pybamm-team/PyBaMM/pull/868))
 
 ## Bug fixes
 
diff --git a/docs/source/solvers/algebraic_solvers.rst b/docs/source/solvers/algebraic_solvers.rst
index cb87d0e0f..df241f927 100644
--- a/docs/source/solvers/algebraic_solvers.rst
+++ b/docs/source/solvers/algebraic_solvers.rst
@@ -3,3 +3,6 @@ Algebraic Solvers
 
 .. autoclass:: pybamm.AlgebraicSolver
   :members:
+
+.. autoclass:: pybamm.CasadiAlgebraicSolver
+  :members:
diff --git a/examples/scripts/DFN_ambient_temperature.py b/examples/scripts/DFN_ambient_temperature.py
index a2dfea865..f50e066fc 100644
--- a/examples/scripts/DFN_ambient_temperature.py
+++ b/examples/scripts/DFN_ambient_temperature.py
@@ -1,5 +1,5 @@
 #
-# Example showing how to load and solve the DFN
+# Example showing how to solve the DFN with a varying ambient temperature
 #
 
 import pybamm
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index da9022e60..65419a5e0 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -245,6 +245,7 @@ from .solvers.solution import Solution, _BaseSolution
 from .solvers.base_solver import BaseSolver
 from .solvers.algebraic_solver import AlgebraicSolver
 from .solvers.casadi_solver import CasadiSolver
+from .solvers.casadi_algebraic_solver import CasadiAlgebraicSolver
 from .solvers.scikits_dae_solver import ScikitsDaeSolver
 from .solvers.scikits_ode_solver import ScikitsOdeSolver, have_scikits_odes
 from .solvers.scipy_solver import ScipySolver
diff --git a/pybamm/solvers/algebraic_solver.py b/pybamm/solvers/algebraic_solver.py
index db366fc1f..35ed42f08 100644
--- a/pybamm/solvers/algebraic_solver.py
+++ b/pybamm/solvers/algebraic_solver.py
@@ -6,7 +6,7 @@ import numpy as np
 from scipy import optimize
 
 
-class AlgebraicSolver(object):
+class AlgebraicSolver(pybamm.BaseSolver):
     """Solve a discretised model which contains only (time independent) algebraic
     equations using a root finding algorithm.
     Note: this solver could be extended for quasi-static models, or models in
@@ -17,22 +17,16 @@ class AlgebraicSolver(object):
     ----------
     method : str, optional
         The method to use to solve the system (default is "lm")
-    tolerance : float, optional
+    tol : float, optional
         The tolerance for the solver (default is 1e-6).
     """
 
     def __init__(self, method="lm", tol=1e-6):
-        self.method = method
+        super().__init__(method=method)
         self.tol = tol
         self.name = "Algebraic solver ({})".format(method)
-
-    @property
-    def method(self):
-        return self._method
-
-    @method.setter
-    def method(self, value):
-        self._method = value
+        self.algebraic_solver = True
+        pybamm.citations.register("virtanen2020scipy")
 
     @property
     def tol(self):
@@ -42,178 +36,75 @@ class AlgebraicSolver(object):
     def tol(self, value):
         self._tol = value
 
-    def solve(self, model):
-        """Calculate the solution of the model.
-
-        Parameters
-        ----------
-        model : :class:`pybamm.BaseModel`
-            The model whose solution to calculate. Must only contain algebraic
-            equations.
-
-        """
-        pybamm.logger.info("Start solving {} with {}".format(model.name, self.name))
-
-        # Set up
-        timer = pybamm.Timer()
-        start_time = timer.time()
-        concatenated_algebraic, jac = self.set_up(model)
-        set_up_time = timer.time() - start_time
-
-        # Create function to evaluate algebraic
-        def algebraic(y):
-            return concatenated_algebraic.evaluate(0, y, known_evals={})[0][:, 0]
-
-        # Create function to evaluate jacobian
-        if jac is not None:
-
-            def jacobian(y):
-                # Note: we only use this solver for time independent algebraic
-                # systems, so jac is arbitrarily evaluated at t=0. Also, needs
-                # to be converted from sparse to dense, so in very large
-                # algebraic models it may be best to switch use_jacobian to False
-                # by default.
-                return jac.evaluate(0, y, known_evals={})[0].toarray()
-
-        else:
-            jacobian = None
-
-        # Use "initial conditions" set in model as initial guess
-        y0_guess = model.concatenated_initial_conditions.evaluate(t=0)
-
-        # Solve
-        solve_start_time = timer.time()
-        pybamm.logger.info("Calling root finding algorithm")
-        solution = self.root(algebraic, y0_guess, jacobian=jacobian)
-        solution.model = model
-
-        # Assign times
-        solution.solve_time = timer.time() - solve_start_time
-        solution.set_up_time = set_up_time
-
-        pybamm.logger.info("Finish solving {}".format(model.name))
-        pybamm.logger.info(
-            "Set-up time: {}, Solve time: {}, Total time: {}".format(
-                timer.format(solution.set_up_time),
-                timer.format(solution.solve_time),
-                timer.format(solution.total_time),
-            )
-        )
-        return solution
-
-    def root(self, algebraic, y0_guess, jacobian=None):
+    def _integrate(self, model, t_eval, inputs=None):
         """
         Calculate the solution of the algebraic equations through root-finding
 
         Parameters
         ----------
-        algebraic : method
-            Function that takes in y and returns the value of the algebraic
-            equations
-        y0_guess : array-like
-            Array of the user's guess for the solution, used to initialise
-            the root finding algorithm
-        jacobian : method, optional
-            A function that takes in t and y and returns the Jacobian. If
-            None, the solver will approximate the Jacobian if required.
+        model : :class:`pybamm.BaseModel`
+            The model whose solution to calculate.
+        t_eval : :class:`numpy.array`, size (k,)
+            The times at which to compute the solution
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
         """
+        algebraic = model.algebraic_eval
+        y0 = model.y0
 
-        def root_fun(y0):
-            "Evaluates algebraic using y0"
-            out = algebraic(y0)
-            pybamm.logger.debug(
-                "Evaluating algebraic equations, L2-norm is {}".format(
-                    np.linalg.norm(out)
+        y = np.empty((len(y0), len(t_eval)))
+
+        for idx, t in enumerate(t_eval):
+
+            def root_fun(y):
+                "Evaluates algebraic using y"
+                out = algebraic(t, y)
+                pybamm.logger.debug(
+                    "Evaluating algebraic equations at t={}, L2-norm is {}".format(
+                        t, np.linalg.norm(out)
+                    )
                 )
-            )
-            return out
-
-        if jacobian:
-            sol = optimize.root(
-                root_fun, y0_guess, method=self.method, tol=self.tol, jac=jacobian
-            )
-        else:
-            sol = optimize.root(root_fun, y0_guess, method=self.method, tol=self.tol)
-
-        if sol.success and np.all(sol.fun < self.tol * len(sol.x)):
-            termination = "success"
-            # Return solution object (no events, so pass None to t_event, y_event)
-            return pybamm.Solution([0], sol.x[:, np.newaxis], termination=termination)
-        elif not sol.success:
-            raise pybamm.SolverError(
-                "Could not find acceptable solution: {}".format(sol.message)
-            )
-        else:
-            raise pybamm.SolverError(
-                """
-                Could not find acceptable solution: solver terminated
-                successfully, but maximum solution error ({}) above tolerance ({})
-                """.format(
-                    np.max(sol.fun), self.tol * len(sol.x)
+                return out
+
+            if model.jacobian_eval is not None:
+
+                def jac(y):
+                    return model.jacobian_eval(t, y)
+
+            else:
+                jac = None
+
+            # Evaluate algebraic with new t and previous y0, if it's already close
+            # enough then keep it
+            if np.all(abs(algebraic(t, y0)) < self.tol):
+                pybamm.logger.debug("Keeping same solution at t={}".format(t))
+                y[:, idx] = y0
+            # Otherwise calculate new y0
+            else:
+                sol = optimize.root(
+                    root_fun, y0, method=self.method, tol=self.tol, jac=jac,
                 )
-            )
 
-    def set_up(self, model):
-        """Unpack model, perform checks, simplify and calculate jacobian.
+                if sol.success and np.all(abs(sol.fun) < self.tol):
+                    # update initial guess for the next iteration
+                    y0 = sol.x
+                    # update solution array
+                    y[:, idx] = y0
+                elif not sol.success:
+                    raise pybamm.SolverError(
+                        "Could not find acceptable solution: {}".format(sol.message)
+                    )
+                else:
+                    raise pybamm.SolverError(
+                        """
+                        Could not find acceptable solution: solver terminated
+                        successfully, but maximum solution error ({})
+                        above tolerance ({})
+                        """.format(
+                            np.max(sol.fun), self.tol
+                        )
+                    )
+
+        # Return solution object (no events, so pass None to t_event, y_event)
+        return pybamm.Solution(t_eval, y, termination="success")
 
-        Parameters
-        ----------
-        model : :class:`pybamm.BaseModel`
-            The model whose solution to calculate. Must have attributes rhs and
-            initial_conditions
-
-        Returns
-        -------
-        concatenated_algebraic : :class:`pybamm.Concatenation`
-            Algebraic equations, which should evaluate to zero
-        jac : :class:`pybamm.SparseStack`
-            Jacobian matrix for the differential and algebraic equations
-
-        Raises
-        ------
-        :class:`pybamm.SolverError`
-            If the model contains any time derivatives, i.e. rhs equations (in
-            which case an ODE or DAE solver should be used instead)
-        """
-        if len(model.rhs) > 0:
-            raise pybamm.SolverError(
-                """Cannot use algebraic solver to solve model with time derivatives"""
-            )
-
-        # create simplified algebraic expressions
-        concatenated_algebraic = model.concatenated_algebraic
-
-        if model.use_simplify:
-            # set up simplification object, for re-use of dict
-            simp = pybamm.Simplification()
-            pybamm.logger.info("Simplifying algebraic")
-            concatenated_algebraic = simp.simplify(concatenated_algebraic)
-
-        if model.use_jacobian:
-            # Create Jacobian from concatenated algebraic
-            y = pybamm.StateVector(
-                slice(0, np.size(model.concatenated_initial_conditions))
-            )
-            # set up Jacobian object, for re-use of dict
-            jacobian = pybamm.Jacobian()
-            pybamm.logger.info("Calculating jacobian")
-            jac = jacobian.jac(concatenated_algebraic, y)
-            model.jacobian = jac
-            model.jacobian_algebraic = jac
-
-            if model.use_simplify:
-                pybamm.logger.info("Simplifying jacobian")
-                jac = simp.simplify(jac)
-
-            if model.convert_to_format == "python":
-                pybamm.logger.info("Converting jacobian to python")
-                jac = pybamm.EvaluatorPython(jac)
-
-        else:
-            jac = None
-
-        if model.convert_to_format == "python":
-            pybamm.logger.info("Converting algebraic to python")
-            concatenated_algebraic = pybamm.EvaluatorPython(concatenated_algebraic)
-
-        return concatenated_algebraic, jac
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 648edf646..38405a6aa 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -55,6 +55,7 @@ class BaseSolver(object):
         # Defaults, can be overwritten by specific solver
         self.name = "Base solver"
         self.ode_solver = False
+        self.algebraic_solver = False
 
     @property
     def method(self):
@@ -116,22 +117,29 @@ class BaseSolver(object):
             Any input parameters to pass to the model when solving
 
         """
-        inputs = inputs or {}
-        y0 = model.concatenated_initial_conditions.evaluate(0, None, inputs)
-
-        # Set model timescale
-        model.timescale_eval = model.timescale.evaluate(u=inputs)
 
         # Check model.algebraic for ode solvers
         if self.ode_solver is True and len(model.algebraic) > 0:
             raise pybamm.SolverError(
                 "Cannot use ODE solver '{}' to solve DAE model".format(self.name)
             )
+        # Check model.rhs for algebraic solvers
+        if self.algebraic_solver is True and len(model.rhs) > 0:
+            raise pybamm.SolverError(
+                """Cannot use algebraic solver to solve model with time derivatives"""
+            )
+
+        inputs = inputs or {}
+        y0 = model.concatenated_initial_conditions.evaluate(0, None, inputs)
+
+        # Set model timescale
+        model.timescale_eval = model.timescale.evaluate(u=inputs)
 
         if self.ode_solver is True:
             self.root_method = None
         if (
-            isinstance(self, pybamm.CasadiSolver) or self.root_method == "casadi"
+            isinstance(self, (pybamm.CasadiSolver, pybamm.CasadiAlgebraicSolver))
+            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"
@@ -276,28 +284,43 @@ class BaseSolver(object):
         # 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]
-            )
+            # can use DAE solver to solve model with algebraic equations only
+            if len(model.rhs) > 0:
+                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(
-                model.concatenated_rhs, model.concatenated_algebraic
-            )
-            # Process again, uses caching so should be quick
-            residuals, residuals_eval, jacobian_eval = process(all_states, "residuals")
-            model.residuals_eval = residuals_eval
-            model.jacobian_eval = jacobian_eval
-            y0_guess = y0.flatten()
-            model.y0 = self.calculate_consistent_state(model, 0, y0_guess, inputs)
-        else:
+        if self.algebraic_solver is True:
+            # we don't calculate consistent initial conditions
+            # for an algebraic solver as this will be the job of the algebraic solver
+            model.residuals_eval = Residuals(algebraic, "residuals", model)
+            model.jacobian_eval = jac_algebraic
+            model.y0 = y0.flatten()
+        elif len(model.algebraic) == 0:
             # can use DAE solver to solve ODE model
+            # - no initial condition initialization needed
             model.residuals_eval = Residuals(rhs, "residuals", model)
             model.jacobian_eval = jac_rhs
             model.y0 = y0.flatten()
+        # Calculate consistent initial conditions for the algebraic equations
+        else:
+            if len(model.rhs) > 0:
+                all_states = pybamm.NumpyConcatenation(
+                    model.concatenated_rhs, model.concatenated_algebraic
+                )
+                # Process again, uses caching so should be quick
+                residuals_eval, jacobian_eval = process(all_states, "residuals")[1:]
+                model.residuals_eval = residuals_eval
+                model.jacobian_eval = jacobian_eval
+            else:
+                model.residuals_eval = Residuals(algebraic, "residuals", model)
+                model.jacobian_eval = jac_algebraic
+            y0_guess = y0.flatten()
+            model.y0 = self.calculate_consistent_state(model, 0, y0_guess, inputs)
 
         pybamm.logger.info("Finish solver set-up")
 
@@ -360,8 +383,6 @@ class BaseSolver(object):
             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",
@@ -376,10 +397,13 @@ class BaseSolver(object):
                 fun = model.casadi_algebraic(
                     time, casadi.vertcat(y0_diff, y0_alg), u_stacked
                 )
+                abs_fun = casadi.fabs(fun)
+                max_fun = casadi.mmax(fun)
             except RuntimeError as err:
                 success = False
                 message = err.args[0]
-                fun = None
+                abs_fun = None
+                max_fun = None
         else:
             algebraic = model.algebraic_eval
             jac = model.jac_algebraic_eval
@@ -430,9 +454,11 @@ class BaseSolver(object):
             y0_alg = sol.x
             success = sol.success
             fun = sol.fun
+            abs_fun = np.abs(fun)
+            max_fun = np.max(fun)
             message = sol.message
 
-        if success and np.all(fun < self.root_tol * len(y0_alg)):
+        if success and np.all(abs_fun < self.root_tol):
             # 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")
@@ -447,11 +473,11 @@ class BaseSolver(object):
                 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)
+                    max_fun, self.root_tol
                 )
             )
 
-    def solve(self, model, t_eval, external_variables=None, inputs=None):
+    def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         """
         Execute the solver setup and calculate the solution of the model at
         specified times.
@@ -481,6 +507,14 @@ class BaseSolver(object):
         if len(model.rhs) == 0 and len(model.algebraic) == 0:
             raise pybamm.ModelError("Cannot solve empty model")
 
+        # t_eval can only be None if the solver is an algebraic solver. In that case
+        # set it to 0
+        if t_eval is None:
+            if self.algebraic_solver is True:
+                t_eval = np.array([0])
+            else:
+                raise ValueError("t_eval cannot be None")
+
         # Make sure t_eval is monotonic
         if (np.diff(t_eval) < 0).any():
             raise pybamm.SolverError("t_eval must increase monotonically")
@@ -495,14 +529,6 @@ class BaseSolver(object):
         inputs = inputs or {}
         ext_and_inputs = {**external_variables, **inputs}
 
-        # Raise warning if t_eval looks like it was supposed to be dimensionless
-        # already
-        if t_eval[-1] < 0.5:
-            raise pybamm.SolverError(
-                """It looks like t_eval might be dimensionless.
-                t_eval should now be provided in seconds"""
-            )
-
         # Set up (if not done already)
         if model not in self.models_set_up:
             self.set_up(model, ext_and_inputs)
diff --git a/pybamm/solvers/casadi_algebraic_solver.py b/pybamm/solvers/casadi_algebraic_solver.py
new file mode 100644
index 000000000..c9e4c44af
--- /dev/null
+++ b/pybamm/solvers/casadi_algebraic_solver.py
@@ -0,0 +1,114 @@
+#
+# Casadi algebraic solver class
+#
+import casadi
+import pybamm
+import numpy as np
+
+
+class CasadiAlgebraicSolver(pybamm.BaseSolver):
+    """Solve a discretised model which contains only (time independent) algebraic
+    equations using CasADi's root finding algorithm.
+    Note: this solver could be extended for quasi-static models, or models in
+    which the time derivative is manually discretised and results in a (possibly
+    nonlinear) algebaric system at each time level.
+
+    Parameters
+    ----------
+    tol : float, optional
+        The tolerance for the solver (default is 1e-6).
+    """
+
+    def __init__(self, method="lm", tol=1e-6, **extra_options):
+        super().__init__()
+        self.tol = tol
+        self.name = "CasADi algebraic solver"
+        self.algebraic_solver = True
+        self.extra_options = extra_options
+        pybamm.citations.register("Andersson2019")
+
+    @property
+    def tol(self):
+        return self._tol
+
+    @tol.setter
+    def tol(self, value):
+        self._tol = value
+
+    def _integrate(self, model, t_eval, inputs=None):
+        """
+        Calculate the solution of the algebraic equations through root-finding
+
+        Parameters
+        ----------
+        model : :class:`pybamm.BaseModel`
+            The model whose solution to calculate.
+        t_eval : :class:`numpy.array`, size (k,)
+            The times at which to compute the solution
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
+        """
+        y0 = model.y0
+
+        y = np.empty((len(y0), len(t_eval)))
+
+        # Set up
+        u_stacked = casadi.vertcat(*[x for x in inputs.values()])
+        t_sym = casadi.MX.sym("t")
+        y_sym = casadi.MX.sym("y_alg", y0.shape[0])
+        u_sym = casadi.MX.sym("u", u_stacked.shape[0])
+
+        t_u_sym = casadi.vertcat(t_sym, u_sym)
+        alg = model.casadi_algebraic(t_sym, y_sym, u_sym)
+
+        # Set up rootfinder
+        roots = casadi.rootfinder(
+            "roots",
+            "newton",
+            dict(x=y_sym, p=t_u_sym, g=alg),
+            {**self.extra_options, "abstol": self.tol},
+        )
+        for idx, t in enumerate(t_eval):
+            # Evaluate algebraic with new t and previous y0, if it's already close
+            # enough then keep it
+            if np.all(abs(model.algebraic_eval(t, y0)) < self.tol):
+                pybamm.logger.debug("Keeping same solution at t={}".format(t))
+                y[:, idx] = y0
+            # Otherwise calculate new y0
+            else:
+                t_u_stacked = casadi.vertcat(t, u_stacked)
+                # Solve
+                try:
+                    y_sol = roots(y0, t_u_stacked).full().flatten()
+                    success = True
+                    message = None
+                    # Check final output
+                    fun = model.casadi_algebraic(t, y_sol, u_stacked)
+                except RuntimeError as err:
+                    success = False
+                    message = err.args[0]
+                    fun = None
+
+                if success and np.all(casadi.fabs(fun) < self.tol):
+                    # update initial guess for the next iteration
+                    y0 = y_sol
+                    # update solution array
+                    y[:, idx] = y_sol
+                elif not success:
+                    raise pybamm.SolverError(
+                        "Could not find acceptable solution: {}".format(message)
+                    )
+                else:
+                    raise pybamm.SolverError(
+                        """
+                        Could not find acceptable solution: solver terminated
+                        successfully, but maximum solution error ({})
+                        above tolerance ({})
+                        """.format(
+                            casadi.mmax(fun), self.tol
+                        )
+                    )
+
+        # Return solution object (no events, so pass None to t_event, y_event)
+        return pybamm.Solution(t_eval, y, termination="success")
+
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index 1029a8560..4075732a6 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -88,7 +88,13 @@ class CasadiSolver(pybamm.BaseSolver):
             Any external variables or input parameters to pass to the model when solving
         """
         inputs = inputs or {}
-
+        if len(model.rhs) == 0:
+            # casadi solver won't allow solving algebraic model so we have to raise an
+            # error here
+            raise pybamm.SolverError(
+                "Cannot use CasadiSolver to solve algebraic model,"
+                "use CasadiAlgebraicSolver instead"
+            )
         if self.mode == "fast":
             integrator = self.get_integrator(model, t_eval, inputs)
             solution = self._run_integrator(integrator, model, model.y0, inputs, t_eval)


Test patch:
diff --git a/tests/unit/test_citations.py b/tests/unit/test_citations.py
index 0590120f2..a930f7ac3 100644
--- a/tests/unit/test_citations.py
+++ b/tests/unit/test_citations.py
@@ -110,6 +110,11 @@ class TestCitations(unittest.TestCase):
         pybamm.ScipySolver()
         self.assertIn("virtanen2020scipy", citations._papers_to_cite)
 
+        citations._reset()
+        self.assertNotIn("virtanen2020scipy", citations._papers_to_cite)
+        pybamm.AlgebraicSolver()
+        self.assertIn("virtanen2020scipy", citations._papers_to_cite)
+
         if pybamm.have_scikits_odes():
             citations._reset()
             self.assertNotIn("scikits-odes", citations._papers_to_cite)
diff --git a/tests/unit/test_solvers/test_algebraic_solver.py b/tests/unit/test_solvers/test_algebraic_solver.py
index 772426bb6..5999f4d9d 100644
--- a/tests/unit/test_solvers/test_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_algebraic_solver.py
@@ -35,53 +35,61 @@ class TestAlgebraicSolver(unittest.TestCase):
 
     def test_simple_root_find(self):
         # Simple system: a single algebraic equation
-        def algebraic(y):
-            return y + 2
+        class Model:
+            y0 = np.array([2])
+            jacobian_eval = None
+
+            def algebraic_eval(self, t, y):
+                return y + 2
 
         solver = pybamm.AlgebraicSolver()
-        y0 = np.array([2])
-        solution = solver.root(algebraic, y0)
+        model = Model()
+        solution = solver._integrate(model, np.array([0]))
         np.testing.assert_array_equal(solution.y, -2)
 
     def test_root_find_fail(self):
-        def algebraic(y):
-            # algebraic equation has no real root
-            return y ** 2 + 1
+        class Model:
+            y0 = np.array([2])
+            jacobian_eval = None
 
-        solver = pybamm.AlgebraicSolver(method="hybr")
-        y0 = np.array([2])
+            def algebraic_eval(self, t, y):
+                # algebraic equation has no real root
+                return y ** 2 + 1
 
+        model = Model()
+
+        solver = pybamm.AlgebraicSolver(method="hybr")
         with self.assertRaisesRegex(
             pybamm.SolverError,
             "Could not find acceptable solution: The iteration is not making",
         ):
-            solver.root(algebraic, y0)
+            solver._integrate(model, np.array([0]))
+
         solver = pybamm.AlgebraicSolver()
         with self.assertRaisesRegex(
             pybamm.SolverError, "Could not find acceptable solution: solver terminated"
         ):
-            solver.root(algebraic, y0)
+            solver._integrate(model, np.array([0]))
 
     def test_with_jacobian(self):
         A = np.array([[4, 3], [1, -1]])
         b = np.array([0, 7])
 
-        def algebraic(y):
-            return A @ y - b
+        class Model:
+            y0 = np.zeros(2)
+
+            def algebraic_eval(self, t, y):
+                return A @ y - b
 
-        def jac(y):
-            return A
+            def jacobian_eval(self, t, y):
+                return A
 
-        y0 = np.zeros(2)
+        model = Model()
         sol = np.array([3, -4])[:, np.newaxis]
 
         solver = pybamm.AlgebraicSolver()
-
-        solution_no_jac = solver.root(algebraic, y0)
-        solution_with_jac = solver.root(algebraic, y0, jacobian=jac)
-
-        np.testing.assert_array_almost_equal(solution_no_jac.y, sol)
-        np.testing.assert_array_almost_equal(solution_with_jac.y, sol)
+        solution = solver._integrate(model, np.array([0]))
+        np.testing.assert_array_almost_equal(solution.y, sol)
 
     def test_model_solver(self):
         # Create model
@@ -117,6 +125,50 @@ class TestAlgebraicSolver(unittest.TestCase):
             model.variables["var2"].evaluate(t=None, y=solution_no_jac.y), sol[100:]
         )
 
+    def test_model_solver_with_time(self):
+        # Create model
+        model = pybamm.BaseModel()
+        var1 = pybamm.Variable("var1")
+        var2 = pybamm.Variable("var2")
+        model.algebraic = {var1: var1 - 3 * pybamm.t, var2: 2 * var1 - var2}
+        model.initial_conditions = {var1: pybamm.Scalar(1), var2: pybamm.Scalar(4)}
+        model.variables = {"var1": var1, "var2": var2}
+
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        # Solve
+        t_eval = np.linspace(0, 1)
+        solver = pybamm.AlgebraicSolver()
+        solution = solver.solve(model, t_eval)
+
+        sol = np.vstack((3 * t_eval, 6 * t_eval))
+        np.testing.assert_array_equal(solution.y, sol)
+        np.testing.assert_array_equal(
+            model.variables["var1"].evaluate(t=t_eval, y=solution.y).flatten(),
+            sol[0, :],
+        )
+        np.testing.assert_array_equal(
+            model.variables["var2"].evaluate(t=t_eval, y=solution.y).flatten(),
+            sol[1, :],
+        )
+
+    def test_solve_with_input(self):
+        # Simple system: a single algebraic equation
+        var = pybamm.Variable("var")
+        model = pybamm.BaseModel()
+        model.algebraic = {var: var + pybamm.InputParameter("value")}
+        model.initial_conditions = {var: 2}
+
+        # create discretisation
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.AlgebraicSolver()
+        solution = solver.solve(model, np.linspace(0, 1, 10), inputs={"value": 7})
+        np.testing.assert_array_equal(solution.y, -7)
+
 
 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 7ab105982..7100d18b5 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -28,6 +28,18 @@ class TestBaseSolver(unittest.TestCase):
         with self.assertRaisesRegex(pybamm.ModelError, "Cannot solve empty model"):
             solver.solve(model, None)
 
+    def test_t_eval_none(self):
+        model = pybamm.BaseModel()
+        v = pybamm.Variable("v")
+        model.rhs = {v: 1}
+        model.initial_conditions = {v: 1}
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        solver = pybamm.BaseSolver()
+        with self.assertRaisesRegex(ValueError, "t_eval cannot be None"):
+            solver.solve(model, None)
+
     def test_nonmonotonic_teval(self):
         solver = pybamm.BaseSolver(rtol=1e-2, atol=1e-4)
         model = pybamm.BaseModel()
@@ -164,16 +176,6 @@ 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_casadi_algebraic_solver.py b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
new file mode 100644
index 000000000..0c9e0d332
--- /dev/null
+++ b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
@@ -0,0 +1,111 @@
+#
+# Tests for the Casadi Algebraic Solver class
+#
+import casadi
+import pybamm
+import unittest
+import numpy as np
+
+
+class TestCasadiAlgebraicSolver(unittest.TestCase):
+    def test_algebraic_solver_init(self):
+        solver = pybamm.CasadiAlgebraicSolver(tol=1e-4)
+        self.assertEqual(solver.tol, 1e-4)
+
+        solver.tol = 1e-5
+        self.assertEqual(solver.tol, 1e-5)
+
+    def test_simple_root_find(self):
+        # Simple system: a single algebraic equation
+        var = pybamm.Variable("var")
+        model = pybamm.BaseModel()
+        model.algebraic = {var: var + 2}
+        model.initial_conditions = {var: 2}
+
+        # create discretisation
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.CasadiAlgebraicSolver()
+        solution = solver.solve(model, np.linspace(0, 1, 10))
+        np.testing.assert_array_equal(solution.y, -2)
+
+    def test_root_find_fail(self):
+        class Model:
+            y0 = np.array([2])
+            t = casadi.MX.sym("t")
+            y = casadi.MX.sym("y")
+            u = casadi.MX.sym("u")
+            casadi_algebraic = casadi.Function("alg", [t, y, u], [y ** 2 + 1])
+
+            def algebraic_eval(self, t, y):
+                # algebraic equation has no real root
+                return y ** 2 + 1
+
+        model = Model()
+
+        solver = pybamm.CasadiAlgebraicSolver()
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "Could not find acceptable solution: .../casadi",
+        ):
+            solver._integrate(model, np.array([0]), {})
+        solver = pybamm.CasadiAlgebraicSolver(error_on_fail=False)
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "Could not find acceptable solution: solver terminated",
+        ):
+            solver._integrate(model, np.array([0]), {})
+
+    def test_model_solver_with_time(self):
+        # Create model
+        model = pybamm.BaseModel()
+        var1 = pybamm.Variable("var1")
+        var2 = pybamm.Variable("var2")
+        model.algebraic = {var1: var1 - 3 * pybamm.t, var2: 2 * var1 - var2}
+        model.initial_conditions = {var1: pybamm.Scalar(1), var2: pybamm.Scalar(4)}
+        model.variables = {"var1": var1, "var2": var2}
+
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        # Solve
+        t_eval = np.linspace(0, 1)
+        solver = pybamm.CasadiAlgebraicSolver()
+        solution = solver.solve(model, t_eval)
+
+        sol = np.vstack((3 * t_eval, 6 * t_eval))
+        np.testing.assert_array_almost_equal(solution.y, sol)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=t_eval, y=solution.y).flatten(),
+            sol[0, :],
+        )
+        np.testing.assert_array_almost_equal(
+            model.variables["var2"].evaluate(t=t_eval, y=solution.y).flatten(),
+            sol[1, :],
+        )
+
+    def test_solve_with_input(self):
+        # Simple system: a single algebraic equation
+        var = pybamm.Variable("var")
+        model = pybamm.BaseModel()
+        model.algebraic = {var: var + pybamm.InputParameter("value")}
+        model.initial_conditions = {var: 2}
+
+        # create discretisation
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.CasadiAlgebraicSolver()
+        solution = solver.solve(model, np.linspace(0, 1, 10), inputs={"value": 7})
+        np.testing.assert_array_equal(solution.y, -7)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 731976846..6148a33c6 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -309,6 +309,22 @@ class TestCasadiSolver(unittest.TestCase):
         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_dae_solver_algebraic_model(self):
+        model = pybamm.BaseModel()
+        var = pybamm.Variable("var")
+        model.algebraic = {var: var + 1}
+        model.initial_conditions = {var: 0}
+
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        solver = pybamm.CasadiSolver()
+        t_eval = np.linspace(0, 1)
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "Cannot use CasadiSolver to solve algebraic model"
+        ):
+            solver.solve(model, t_eval)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py
index 703673ae8..b440d29a3 100644
--- a/tests/unit/test_solvers/test_idaklu_solver.py
+++ b/tests/unit/test_solvers/test_idaklu_solver.py
@@ -82,6 +82,20 @@ class TestIDAKLUSolver(unittest.TestCase):
         with self.assertRaisesRegex(pybamm.SolverError, "KLU requires the Jacobian"):
             solver.solve(model, t_eval)
 
+    def test_dae_solver_algebraic_model(self):
+        model = pybamm.BaseModel()
+        var = pybamm.Variable("var")
+        model.algebraic = {var: var + 1}
+        model.initial_conditions = {var: 0}
+
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        solver = pybamm.IDAKLUSolver()
+        t_eval = np.linspace(0, 1)
+        solution = solver.solve(model, t_eval)
+        np.testing.assert_array_equal(solution.y, -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 905507c69..82b452033 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -696,6 +696,20 @@ class TestScikitsSolvers(unittest.TestCase):
         with self.assertRaisesRegex(pybamm.SolverError, "Cannot use ODE solver"):
             solver.set_up(model)
 
+    def test_dae_solver_algebraic_model(self):
+        model = pybamm.BaseModel()
+        var = pybamm.Variable("var")
+        model.algebraic = {var: var + 1}
+        model.initial_conditions = {var: 0}
+
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        solver = pybamm.ScikitsDaeSolver()
+        t_eval = np.linspace(0, 1)
+        solution = solver.solve(model, t_eval)
+        np.testing.assert_array_equal(solution.y, -1)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")