#644 set atol and rtol

pybamm/models/full_battery_models/base_battery_model.py

@@ -128,9 +128,9 @@ def default_geometry(self):
     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,

pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py

@@ -60,6 +60,11 @@ def default_geometry(self):
         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

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 @@ def method(self, value):
         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 @@ def solve(self, model, t_eval):
         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),

pybamm/solvers/dae_solver.py

@@ -12,21 +12,29 @@ 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
         (default is 1000).
     """
 
     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

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.

pybamm/solvers/scikits_dae_solver.py

@@ -22,25 +22,32 @@ 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
         (default is 1000).
     """
 
     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 @@ def rootfn(t, y, ydot, return_root):
 
         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 @@ def jacfn(t, y, ydot, residuals, cj, J):
                 np.transpose(sol.values.y),
                 sol.roots.t,
                 np.transpose(sol.roots.y),
-                termination
+                termination,
             )
         else:
             raise pybamm.SolverError(sol.message)

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 @@ def jac_times_setupfn(t, y, fy, userdata):
 
         extra_options = {
             "old_api": False,
-            "rtol": self.tol,
-            "atol": self.tol,
+            "rtol": self.rtol,
+            "atol": self.atol,
             "linsolver": self.linsolver,
         }
 
@@ -134,7 +135,7 @@ def jac_times_setupfn(t, y, fy, userdata):
                 np.transpose(sol.values.y),
                 sol.roots.t,
                 np.transpose(sol.roots.y),
-                termination
+                termination,
             )
         else:
             raise pybamm.SolverError(sol.message)

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 @@ def integrate(
             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 @@ def integrate(
                 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)

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()