Repository: pybamm-team/PyBaMM (python)

Issue:
set relative and absolute tolerances separately for solver
At the moment we use the same tol for both. Would be good to have additional control over these. Also, the defaults are pretty tight so we might relax them a bit to get a good balance of speed and accuracy.  

Pull request description:
#644 set atol and rtol
# Description

Allows the user to set `rtol` and `atol` separately. For the defaults I've selected 1e-3 for rtol and 1e-6 for atol. This is the same as other solvers (e.g. scipy.integrate or MATLABs ODE15s) and provides a pretty significant speedup vs 1e-8 for both. I've also reduced the default number of points in x, as again this speed things up quite a bit. Both of these result in a RMS voltage error of around 1e-4 compared with the previous defaults, but the default DFN now runs in around 3 seconds vs 23 seconds (on my computer at least).

I think for most users, this kind of error is ok out of the box (the plots are indistinguishable), particularly if they are new users and just want to quickly run some examples. Users now have more control over both solver tolerances, so I think this is an ok compromise.

EDIT: I've decreased rtol to 1e-6 too. Seems like the thermal model and some other models didn't like rtol as low as 1e-3. Probably safer to have both set as 1e-6 so that the solvers don't immediately  fall over when people add new physics.

Fixes #644 

## Type of change

Please delete options that are not relevant.

- [ ] New feature (non-breaking change which adds functionality)


# Key checklist:

- [ ] No style issues: `$ flake8`
- [ ] All tests pass: `$ python run-tests.py --unit`
- [ ] 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: 

- [ ] Code is commented, particularly in hard-to-understand areas
- [ ] Tests added that prove fix is effective or that feature works
- [ ] Any dependent changes have been merged and published in downstream modules


Patch:
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 1db1276c5..c301255a7 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -128,9 +128,9 @@ class BaseBatteryModel(pybamm.BaseModel):
     def default_var_pts(self):
         var = pybamm.standard_spatial_vars
         return {
-            var.x_n: 40,
-            var.x_s: 25,
-            var.x_p: 35,
+            var.x_n: 20,
+            var.x_s: 20,
+            var.x_p: 20,
             var.r_n: 10,
             var.r_p: 10,
             var.y: 10,
diff --git a/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py b/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py
index 9f223ae9a..cf753ddb4 100644
--- a/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py
+++ b/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py
@@ -60,6 +60,11 @@ class BaseModel(pybamm.BaseBatteryModel):
         elif self.options["dimensionality"] == 2:
             return pybamm.Geometry("2+1D macro")
 
+    @property
+    def default_var_pts(self):
+        var = pybamm.standard_spatial_vars
+        return {var.x_n: 30, var.x_s: 30, var.x_p: 30, var.y: 10, var.z: 10}
+
     def set_standard_output_variables(self):
         super().set_standard_output_variables()
         # Current
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 63dd82ad6..ad1ad2b4a 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -10,13 +10,16 @@ class BaseSolver(object):
 
     Parameters
     ----------
-    tolerance : float, optional
-        The tolerance for the solver (default is 1e-8).
+    rtol : float, optional
+        The relative tolerance for the solver (default is 1e-6).
+    atol : float, optional
+        The absolute tolerance for the solver (default is 1e-6).
     """
 
-    def __init__(self, method=None, tol=1e-8):
+    def __init__(self, method=None, rtol=1e-6, atol=1e-6):
         self._method = method
-        self._tol = tol
+        self._rtol = rtol
+        self._atol = atol
 
     @property
     def method(self):
@@ -27,12 +30,20 @@ class BaseSolver(object):
         self._method = value
 
     @property
-    def tol(self):
-        return self._tol
+    def rtol(self):
+        return self._rtol
 
-    @tol.setter
-    def tol(self, value):
-        self._tol = value
+    @rtol.setter
+    def rtol(self, value):
+        self._rtol = value
+
+    @property
+    def atol(self):
+        return self._atol
+
+    @atol.setter
+    def atol(self, value):
+        self._atol = value
 
     def solve(self, model, t_eval):
         """
@@ -73,7 +84,7 @@ class BaseSolver(object):
         solution.total_time = timer.time() - start_time
         solution.set_up_time = set_up_time
 
-        pybamm.logger.warning("Finish solving {} ({})".format(model.name, termination))
+        pybamm.logger.info("Finish solving {} ({})".format(model.name, termination))
         pybamm.logger.info(
             "Set-up time: {}, Solve time: {}, Total time: {}".format(
                 timer.format(solution.set_up_time),
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 035a31f05..003403254 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -12,11 +12,13 @@ class DaeSolver(pybamm.BaseSolver):
 
     Parameters
     ----------
-    tolerance : float, optional
-        The tolerance for the solver (default is 1e-8).
+    rtol : float, optional
+        The relative tolerance for the solver (default is 1e-6).
+    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")
-    tolerance : float, optional
+    root_tol : float, optional
         The tolerance for the initial-condition solver (default is 1e-8).
     max_steps: int, optional
         The maximum number of steps the solver will take before terminating
@@ -24,9 +26,15 @@ class DaeSolver(pybamm.BaseSolver):
     """
 
     def __init__(
-        self, method=None, tol=1e-8, root_method="lm", root_tol=1e-6, max_steps=1000
+        self,
+        method=None,
+        rtol=1e-6,
+        atol=1e-6,
+        root_method="lm",
+        root_tol=1e-6,
+        max_steps=1000,
     ):
-        super().__init__(method, tol)
+        super().__init__(method, rtol, atol)
         self.root_method = root_method
         self.root_tol = root_tol
         self.max_steps = max_steps
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 7d16dcb99..db279b257 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -10,12 +10,14 @@ class OdeSolver(pybamm.BaseSolver):
 
     Parameters
     ----------
-    tolerance : float, optional
-        The tolerance for the solver (default is 1e-8).
+    rtol : float, optional
+        The relative tolerance for the solver (default is 1e-6).
+    atol : float, optional
+        The absolute tolerance for the solver (default is 1e-6).
     """
 
-    def __init__(self, method=None, tol=1e-8):
-        super().__init__(method, tol)
+    def __init__(self, method=None, rtol=1e-6, atol=1e-6):
+        super().__init__(method, rtol, atol)
 
     def compute_solution(self, model, t_eval):
         """Calculate the solution of the model at specified times.
diff --git a/pybamm/solvers/scikits_dae_solver.py b/pybamm/solvers/scikits_dae_solver.py
index 0ed22a103..ee99d4602 100644
--- a/pybamm/solvers/scikits_dae_solver.py
+++ b/pybamm/solvers/scikits_dae_solver.py
@@ -22,12 +22,13 @@ class ScikitsDaeSolver(pybamm.DaeSolver):
     ----------
     method : str, optional
         The method to use in solve_ivp (default is "BDF")
-    tolerance : float, optional
-        The tolerance for the solver (default is 1e-8). Set as the both reltol and
-        abstol in solve_ivp.
+    rtol : float, optional
+        The relative tolerance for the solver (default is 1e-6).
+    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")
-    tolerance : float, optional
+    root_tol : float, optional
         The tolerance for the initial-condition solver (default is 1e-8).
     max_steps: int, optional
         The maximum number of steps the solver will take before terminating
@@ -35,12 +36,18 @@ class ScikitsDaeSolver(pybamm.DaeSolver):
     """
 
     def __init__(
-        self, method="ida", tol=1e-8, root_method="lm", root_tol=1e-6, max_steps=1000
+        self,
+        method="ida",
+        rtol=1e-6,
+        atol=1e-6,
+        root_method="lm",
+        root_tol=1e-6,
+        max_steps=1000,
     ):
         if scikits_odes_spec is None:
             raise ImportError("scikits.odes is not installed")
 
-        super().__init__(method, tol, root_method, root_tol, max_steps)
+        super().__init__(method, rtol, atol, root_method, root_tol, max_steps)
 
     def integrate(
         self, residuals, y0, t_eval, events=None, mass_matrix=None, jacobian=None
@@ -76,8 +83,8 @@ class ScikitsDaeSolver(pybamm.DaeSolver):
 
         extra_options = {
             "old_api": False,
-            "rtol": self.tol,
-            "atol": self.tol,
+            "rtol": self.rtol,
+            "atol": self.atol,
             "max_steps": self.max_steps,
         }
 
@@ -120,7 +127,7 @@ class ScikitsDaeSolver(pybamm.DaeSolver):
                 np.transpose(sol.values.y),
                 sol.roots.t,
                 np.transpose(sol.roots.y),
-                termination
+                termination,
             )
         else:
             raise pybamm.SolverError(sol.message)
diff --git a/pybamm/solvers/scikits_ode_solver.py b/pybamm/solvers/scikits_ode_solver.py
index 17cfff62a..542b15d19 100644
--- a/pybamm/solvers/scikits_ode_solver.py
+++ b/pybamm/solvers/scikits_ode_solver.py
@@ -26,18 +26,19 @@ class ScikitsOdeSolver(pybamm.OdeSolver):
     ----------
     method : str, optional
         The method to use in solve_ivp (default is "BDF")
-    tolerance : float, optional
-        The tolerance for the solver (default is 1e-8). Set as the both reltol and
-        abstol in solve_ivp.
+    rtol : float, optional
+        The relative tolerance for the solver (default is 1e-6).
+    atol : float, optional
+        The absolute tolerance for the solver (default is 1e-6).
     linsolver : str, optional
             Can be 'dense' (= default), 'lapackdense', 'spgmr', 'spbcgs', 'sptfqmr'
     """
 
-    def __init__(self, method="cvode", tol=1e-8, linsolver="dense"):
+    def __init__(self, method="cvode", rtol=1e-6, atol=1e-6, linsolver="dense"):
         if scikits_odes_spec is None:
             raise ImportError("scikits.odes is not installed")
 
-        super().__init__(method, tol)
+        super().__init__(method, rtol, atol)
         self.linsolver = linsolver
 
     def integrate(
@@ -98,8 +99,8 @@ class ScikitsOdeSolver(pybamm.OdeSolver):
 
         extra_options = {
             "old_api": False,
-            "rtol": self.tol,
-            "atol": self.tol,
+            "rtol": self.rtol,
+            "atol": self.atol,
             "linsolver": self.linsolver,
         }
 
@@ -134,7 +135,7 @@ class ScikitsOdeSolver(pybamm.OdeSolver):
                 np.transpose(sol.values.y),
                 sol.roots.t,
                 np.transpose(sol.roots.y),
-                termination
+                termination,
             )
         else:
             raise pybamm.SolverError(sol.message)
diff --git a/pybamm/solvers/scipy_solver.py b/pybamm/solvers/scipy_solver.py
index 60e3db9ac..fd7afbf49 100644
--- a/pybamm/solvers/scipy_solver.py
+++ b/pybamm/solvers/scipy_solver.py
@@ -14,13 +14,14 @@ class ScipySolver(pybamm.OdeSolver):
     ----------
     method : str, optional
         The method to use in solve_ivp (default is "BDF")
-    tolerance : float, optional
-        The tolerance for the solver (default is 1e-8). Set as the both reltol and
-        abstol in solve_ivp.
+    rtol : float, optional
+        The relative tolerance for the solver (default is 1e-6).
+    atol : float, optional
+        The absolute tolerance for the solver (default is 1e-6).
     """
 
-    def __init__(self, method="BDF", tol=1e-8):
-        super().__init__(method, tol)
+    def __init__(self, method="BDF", rtol=1e-6, atol=1e-6):
+        super().__init__(method, rtol, atol)
 
     def integrate(
         self, derivs, y0, t_eval, events=None, mass_matrix=None, jacobian=None
@@ -52,7 +53,7 @@ class ScipySolver(pybamm.OdeSolver):
             various diagnostic messages.
 
         """
-        extra_options = {"rtol": self.tol, "atol": self.tol}
+        extra_options = {"rtol": self.rtol, "atol": self.atol}
 
         # check for user-supplied Jacobian
         implicit_methods = ["Radau", "BDF", "LSODA"]
@@ -90,12 +91,6 @@ class ScipySolver(pybamm.OdeSolver):
                 termination = "final time"
                 t_event = None
                 y_event = np.array(None)
-            return pybamm.Solution(
-                sol.t,
-                sol.y,
-                t_event,
-                y_event,
-                termination
-            )
+            return pybamm.Solution(sol.t, sol.y, t_event, y_event, termination)
         else:
             raise pybamm.SolverError(sol.message)
diff --git a/results/change_settings/change_solver_tolerances.py b/results/change_settings/change_solver_tolerances.py
new file mode 100644
index 000000000..81261e796
--- /dev/null
+++ b/results/change_settings/change_solver_tolerances.py
@@ -0,0 +1,52 @@
+#
+# Compare solution of li-ion battery models when changing solver tolerances
+#
+import numpy as np
+import pybamm
+
+pybamm.set_logging_level("INFO")
+
+# load model
+model = pybamm.lithium_ion.DFN()
+
+
+# process and discretise
+param = model.default_parameter_values
+param.process_model(model)
+geometry = model.default_geometry
+param.process_geometry(geometry)
+mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
+disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+disc.process_model(model)
+
+# tolerances (rtol, atol)
+tols = [[1e-8, 1e-8], [1e-6, 1e-6], [1e-3, 1e-6], [1e-3, 1e-3]]
+
+# solve model
+solutions = [None] * len(tols)
+voltages = [None] * len(tols)
+voltage_rmse = [None] * len(tols)
+labels = [None] * len(tols)
+t_eval = np.linspace(0, 0.17, 100)
+for i, tol in enumerate(tols):
+    solver = pybamm.ScikitsDaeSolver(rtol=tol[0], atol=tol[1])
+    solutions[i] = solver.solve(model, t_eval)
+    voltages[i] = pybamm.ProcessedVariable(
+        model.variables["Terminal voltage [V]"],
+        solutions[i].t,
+        solutions[i].y,
+        mesh=mesh,
+    )(solutions[i].t)
+    voltage_rmse[i] = pybamm.rmse(voltages[0], voltages[i])
+    labels[i] = "rtol = {}, atol = {}".format(tol[0], tol[1])
+
+# print RMSE voltage errors vs tighest tolerance
+for i, tol in enumerate(tols):
+    print(
+        "rtol = {}, atol = {}, solve time = {} s, Voltage RMSE = {}".format(
+            tol[0], tol[1], solutions[i].solve_time, voltage_rmse[i]
+        )
+    )
+# plot
+plot = pybamm.QuickPlot([model] * len(solutions), mesh, solutions, labels=labels)
+plot.dynamic_plot()
diff --git a/results/change_settings/compare_var_pts.py b/results/change_settings/compare_var_pts.py
new file mode 100644
index 000000000..2153a96f6
--- /dev/null
+++ b/results/change_settings/compare_var_pts.py
@@ -0,0 +1,68 @@
+#
+# Compare solution of li-ion battery models when varying the number of grid points
+#
+import numpy as np
+import pybamm
+import matplotlib.pyplot as plt
+
+pybamm.set_logging_level("INFO")
+
+# choose number of points per domain (all domains will have same npts)
+Npts = [30, 20, 10, 5]
+
+# create models
+models = [None] * len(Npts)
+for i, npts in enumerate(Npts):
+    models[i] = pybamm.lithium_ion.DFN(name="Npts = {}".format(npts))
+
+# load parameter values and process models and geometry
+param = models[0].default_parameter_values
+for model in models:
+    param.process_model(model)
+
+# set mesh
+meshes = [None] * len(models)
+
+# create geometry and discretise models
+var = pybamm.standard_spatial_vars
+for i, model in enumerate(models):
+    geometry = model.default_geometry
+    param.process_geometry(geometry)
+    var_pts = {
+        var.x_n: Npts[i],
+        var.x_s: Npts[i],
+        var.x_p: Npts[i],
+        var.r_n: Npts[i],
+        var.r_p: Npts[i],
+    }
+    meshes[i] = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts)
+    disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods)
+    disc.process_model(model)
+
+# solve model and plot voltage
+solutions = [None] * len(models)
+voltages = [None] * len(models)
+voltage_rmse = [None] * len(models)
+t_eval = np.linspace(0, 0.17, 100)
+for i, model in enumerate(models):
+    solutions[i] = model.default_solver.solve(model, t_eval)
+    voltages[i] = pybamm.ProcessedVariable(
+        model.variables["Terminal voltage [V]"],
+        solutions[i].t,
+        solutions[i].y,
+        mesh=meshes[i],
+    )(solutions[i].t)
+    voltage_rmse[i] = pybamm.rmse(voltages[0], voltages[i])
+    plt.plot(solutions[i].t, voltages[i], label=model.name)
+
+for i, npts in enumerate(Npts):
+    print(
+        "npts = {}, solve time = {} s, Voltage RMSE = {}".format(
+            npts, solutions[i].solve_time, voltage_rmse[i]
+        )
+    )
+
+plt.xlabel(r"$t$")
+plt.ylabel("Voltage [V]")
+plt.legend()
+plt.show()


Test patch:
diff --git a/tests/integration/test_models/standard_model_tests.py b/tests/integration/test_models/standard_model_tests.py
index 09fcd23a2..2b5afa124 100644
--- a/tests/integration/test_models/standard_model_tests.py
+++ b/tests/integration/test_models/standard_model_tests.py
@@ -71,6 +71,10 @@ class StandardModelTest(object):
         # Overwrite solver if given
         if solver is not None:
             self.solver = solver
+        # Use tighter default tolerances for testing
+        self.solver.rtol = 1e-8
+        self.solver.atol = 1e-8
+
         if t_eval is None:
             t_eval = np.linspace(0, 1, 100)
 
diff --git a/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_asymptotics_convergence.py b/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_asymptotics_convergence.py
index 4d2405fac..f7158abd8 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_asymptotics_convergence.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_asymptotics_convergence.py
@@ -47,13 +47,19 @@ class TestAsymptoticConvergence(unittest.TestCase):
             param.update_model(leading_order_model, loqs_disc)
             param.update_model(composite_model, comp_disc)
             param.update_model(full_model, full_disc)
-            # Solve, make sure times are the same
+            # Solve, make sure times are the same and use tight tolerances
             t_eval = np.linspace(0, 0.6)
             solver_loqs = leading_order_model.default_solver
+            solver_loqs.rtol = 1e-8
+            solver_loqs.atol = 1e-8
             solution_loqs = solver_loqs.solve(leading_order_model, t_eval)
             solver_comp = composite_model.default_solver
+            solver_comp.rtol = 1e-8
+            solver_comp.atol = 1e-8
             solution_comp = solver_comp.solve(composite_model, t_eval)
             solver_full = full_model.default_solver
+            solver_full.rtol = 1e-8
+            solver_full.atol = 1e-8
             solution_full = solver_full.solve(full_model, t_eval)
 
             # Post-process variables
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index 000b01c6b..9d0075e3e 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -8,11 +8,14 @@ import unittest
 
 class TestBaseSolver(unittest.TestCase):
     def test_base_solver_init(self):
-        solver = pybamm.BaseSolver(tol=1e-4)
-        self.assertEqual(solver.tol, 1e-4)
-
-        solver.tol = 1e-5
-        self.assertEqual(solver.tol, 1e-5)
+        solver = pybamm.BaseSolver(rtol=1e-2, atol=1e-4)
+        self.assertEqual(solver.rtol, 1e-2)
+        self.assertEqual(solver.atol, 1e-4)
+
+        solver.rtol = 1e-5
+        self.assertEqual(solver.rtol, 1e-5)
+        solver.rtol = 1e-7
+        self.assertEqual(solver.rtol, 1e-7)
 
         with self.assertRaises(NotImplementedError):
             solver.compute_solution(None, None)
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 597ecf12e..36b3cb21f 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -13,7 +13,7 @@ from tests import get_mesh_for_testing, get_discretisation_for_testing
 class TestScikitsSolvers(unittest.TestCase):
     def test_ode_integrate(self):
         # Constant
-        solver = pybamm.ScikitsOdeSolver(tol=1e-8)
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-8, atol=1e-8)
 
         def constant_growth(t, y):
             return 0.5 * np.ones_like(y)
@@ -25,7 +25,7 @@ class TestScikitsSolvers(unittest.TestCase):
         np.testing.assert_allclose(0.5 * solution.t, solution.y[0])
 
         # Exponential decay
-        solver = pybamm.ScikitsOdeSolver(tol=1e-8)
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-8, atol=1e-8)
 
         def exponential_decay(t, y):
             return -0.1 * y
@@ -55,7 +55,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
     def test_ode_integrate_with_event(self):
         # Constant
-        solver = pybamm.ScikitsOdeSolver(tol=1e-8)
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-8, atol=1e-8)
 
         def constant_decay(t, y):
             return -2 * np.ones_like(y)
@@ -75,7 +75,7 @@ class TestScikitsSolvers(unittest.TestCase):
         self.assertEqual(solution.termination, "event")
 
         # Expnonential growth
-        solver = pybamm.ScikitsOdeSolver(tol=1e-8)
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-8, atol=1e-8)
 
         def exponential_growth(t, y):
             return y
@@ -101,7 +101,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
     def test_ode_integrate_with_jacobian(self):
         # Linear
-        solver = pybamm.ScikitsOdeSolver(tol=1e-8)
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-8, atol=1e-8)
 
         def linear_ode(t, y):
             return np.array([0.5, 2 - y[0]])
@@ -134,7 +134,7 @@ class TestScikitsSolvers(unittest.TestCase):
         )
         np.testing.assert_allclose(0.5 * solution.t, solution.y[0])
 
-        solver = pybamm.ScikitsOdeSolver(tol=1e-8, linsolver="spgmr")
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-8, atol=1e-8, linsolver="spgmr")
 
         solution = solver.integrate(linear_ode, y0, t_eval, jacobian=jacobian)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -152,7 +152,7 @@ class TestScikitsSolvers(unittest.TestCase):
         np.testing.assert_allclose(0.5 * solution.t, solution.y[0])
 
         # Nonlinear exponential grwoth
-        solver = pybamm.ScikitsOdeSolver(tol=1e-8)
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-8, atol=1e-8)
 
         def exponential_growth(t, y):
             return np.array([y[0], (1.0 - y[0]) * y[1]])
@@ -182,7 +182,7 @@ class TestScikitsSolvers(unittest.TestCase):
             np.exp(1 + solution.t - np.exp(solution.t)), solution.y[1], rtol=1e-4
         )
 
-        solver = pybamm.ScikitsOdeSolver(tol=1e-8, linsolver="spgmr")
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-8, atol=1e-8, linsolver="spgmr")
 
         solution = solver.integrate(exponential_growth, y0, t_eval, jacobian=jacobian)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -202,7 +202,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
     def test_dae_integrate(self):
         # Constant
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
 
         def constant_growth_dae(t, y, ydot):
             return [0.5 * np.ones_like(y[0]) - ydot[0], 2 * y[0] - y[1]]
@@ -215,7 +215,7 @@ class TestScikitsSolvers(unittest.TestCase):
         np.testing.assert_allclose(1.0 * solution.t, solution.y[1])
 
         # Exponential decay
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
 
         def exponential_decay_dae(t, y, ydot):
             return [-0.1 * y[0] - ydot[0], 2 * y[0] - y[1]]
@@ -228,7 +228,7 @@ class TestScikitsSolvers(unittest.TestCase):
         self.assertEqual(solution.termination, "final time")
 
     def test_dae_integrate_failure(self):
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
 
         def constant_growth_dae(t, y, ydot):
             return [0.5 * np.ones_like(y[0]) - ydot[0], 2 * y[0] - y[1]]
@@ -240,7 +240,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
     def test_dae_integrate_bad_ics(self):
         # Constant
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
 
         def constant_growth_dae(t, y, ydot):
             return [0.5 * np.ones_like(y[0]) - ydot[0], 2 * y[0] - y[1]]
@@ -265,7 +265,7 @@ class TestScikitsSolvers(unittest.TestCase):
         np.testing.assert_allclose(1.0 * solution.t, solution.y[1])
 
         # Exponential decay
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
 
         def exponential_decay_dae(t, y, ydot):
             return [-0.1 * y[0] - ydot[0], 2 * y[0] - y[1]]
@@ -278,7 +278,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
     def test_dae_integrate_with_event(self):
         # Constant
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
 
         def constant_growth_dae(t, y, ydot):
             return [0.5 * np.ones_like(y[0]) - ydot[0], 2 * y[0] - y[1]]
@@ -305,7 +305,7 @@ class TestScikitsSolvers(unittest.TestCase):
         self.assertEqual(solution.termination, "event")
 
         # Exponential decay
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
 
         def exponential_decay_dae(t, y, ydot):
             return np.array([-0.1 * y[0] - ydot[0], 2 * y[0] - y[1]])
@@ -334,7 +334,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
     def test_dae_integrate_with_jacobian(self):
         # Constant
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
 
         def constant_growth_dae(t, y, ydot):
             return np.array([0.5 * np.ones_like(y[0]) - ydot[0], 2.0 * y[0] - y[1]])
@@ -354,7 +354,7 @@ class TestScikitsSolvers(unittest.TestCase):
         np.testing.assert_allclose(1.0 * solution.t, solution.y[1])
 
         # Nonlinear (tests when Jacobian a function of y)
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
 
         def nonlinear_dae(t, y, ydot):
             return np.array([0.5 * np.ones_like(y[0]) - ydot[0], 2 * y[0] ** 2 - y[1]])
@@ -375,7 +375,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
     def test_dae_integrate_with_non_unity_mass(self):
         # Constant
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
 
         def constant_growth_dae(t, y, ydot):
             return np.array([0.5 * np.ones_like(y[0]) - 4 * ydot[0], 2.0 * y[0] - y[1]])
@@ -405,7 +405,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsOdeSolver(tol=1e-9)
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-9, atol=1e-9)
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -431,7 +431,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsOdeSolver(tol=1e-9)
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-9, atol=1e-9)
         t_eval = np.linspace(0, 10, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
@@ -474,7 +474,7 @@ class TestScikitsSolvers(unittest.TestCase):
         model.jacobian = jacobian
 
         # Solve
-        solver = pybamm.ScikitsOdeSolver(tol=1e-9)
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-9, atol=1e-9)
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -503,7 +503,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -528,7 +528,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -552,7 +552,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
         t_eval = np.linspace(0, 5, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_less(solution.y[0], 1.5)
@@ -591,7 +591,7 @@ class TestScikitsSolvers(unittest.TestCase):
         model.jacobian = jacobian
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -608,7 +608,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -624,7 +624,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc = get_discretisation_for_testing()
         disc.process_model(model)
 
-        solver = pybamm.ScikitsOdeSolver(tol=1e-9)
+        solver = pybamm.ScikitsOdeSolver(rtol=1e-9, atol=1e-9)
 
         # Step once
         dt = 0.1
@@ -656,7 +656,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc = get_discretisation_for_testing()
         disc.process_model(model)
 
-        solver = pybamm.ScikitsDaeSolver(tol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
 
         # Step once
         dt = 0.1
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index cadc2d12b..cef29c65e 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -11,7 +11,7 @@ import warnings
 class TestScipySolver(unittest.TestCase):
     def test_integrate(self):
         # Constant
-        solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
 
         def constant_growth(t, y):
             return 0.5 * np.ones_like(y)
@@ -23,7 +23,7 @@ class TestScipySolver(unittest.TestCase):
         np.testing.assert_allclose(0.5 * solution.t, solution.y[0])
 
         # Exponential decay
-        solver = pybamm.ScipySolver(tol=1e-8, method="BDF")
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="BDF")
 
         def exponential_decay(t, y):
             return -0.1 * y
@@ -43,7 +43,7 @@ class TestScipySolver(unittest.TestCase):
 
         y0 = np.array([1])
         t_eval = np.linspace(0, 3, 100)
-        solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
         # Expect solver to fail when y goes negative
         with self.assertRaises(pybamm.SolverError):
             solver.integrate(sqrt_decay, y0, t_eval)
@@ -53,7 +53,7 @@ class TestScipySolver(unittest.TestCase):
 
     def test_integrate_with_event(self):
         # Constant
-        solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
 
         def constant_growth(t, y):
             return 0.5 * np.ones_like(y)
@@ -71,7 +71,7 @@ class TestScipySolver(unittest.TestCase):
         self.assertEqual(solution.termination, "event")
 
         # Exponential decay
-        solver = pybamm.ScipySolver(tol=1e-8, method="BDF")
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="BDF")
 
         def exponential_growth(t, y):
             return y
@@ -97,7 +97,7 @@ class TestScipySolver(unittest.TestCase):
 
     def test_ode_integrate_with_jacobian(self):
         # Linear
-        solver = pybamm.ScipySolver(tol=1e-8, method="BDF")
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="BDF")
 
         def linear_ode(t, y):
             return np.array([0.5 * np.ones_like(y[0]), 2.0 - y[0]])
@@ -115,7 +115,7 @@ class TestScipySolver(unittest.TestCase):
         )
 
         # Nonlinear exponential grwoth
-        solver = pybamm.ScipySolver(tol=1e-8, method="BDF")
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="BDF")
 
         def exponential_growth(t, y):
             return np.array([y[0], (1.0 - y[0]) * y[1]])
@@ -147,7 +147,7 @@ class TestScipySolver(unittest.TestCase):
         disc = pybamm.Discretisation(mesh, spatial_methods)
         disc.process_model(model)
         # Solve
-        solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -174,7 +174,7 @@ class TestScipySolver(unittest.TestCase):
         disc = pybamm.Discretisation(mesh, spatial_methods)
         disc.process_model(model)
         # Solve
-        solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
         t_eval = np.linspace(0, 10, 100)
         solution = solver.solve(model, t_eval)
         self.assertLess(len(solution.t), len(t_eval))
@@ -219,7 +219,7 @@ class TestScipySolver(unittest.TestCase):
         model.jacobian = jacobian
 
         # Solve
-        solver = pybamm.ScipySolver(tol=1e-9)
+        solver = pybamm.ScipySolver(rtol=1e-9, atol=1e-9)
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -249,7 +249,7 @@ class TestScipySolver(unittest.TestCase):
         disc = pybamm.Discretisation(mesh, spatial_methods)
         disc.process_model(model)
 
-        solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
 
         # Step once
         dt = 0.1