Repository: pybamm-team/PyBaMM (python)

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

Pull request:
Issue 830 algebraic
# Description

Reformat algebraic solver to be more in line with the other solvers, and add casadi algebraic solver.
We could call the algebraic solvers to initialize the dae models, would reduce code duplication, but might be a little bit confusing?

Fixes #830

## 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 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")