#804 add root_tol as an option

Author: valentinsulzer

pybamm/solvers/base_solver.py

@@ -131,8 +131,7 @@ def set_up(self, model, inputs=None):
         if self.ode_solver is True:
             self.root_method = None
         if (
-            isinstance(self, pybamm.CasadiSolver)
-            or self.root_method == "casadi"
+            isinstance(self, pybamm.CasadiSolver) or self.root_method == "casadi"
         ) and model.convert_to_format != "casadi":
             pybamm.logger.warning(
                 f"Converting {model.name} to CasADi for solving with CasADi solver"
@@ -367,15 +366,20 @@ def calculate_consistent_state(self, model, time=0, y0_guess=None, inputs=None):
                 "roots",
                 "newton",
                 dict(x=y_alg, p=u, g=alg_root),
-                {"error_on_fail": False},
-            )
-            y0_alg = roots(y0_alg_guess, u_stacked).full().flatten()
-            success = True
-            message = None
-            # Check final output
-            fun = model.casadi_algebraic(
-                time, casadi.vertcat(y0_diff, y0_alg), u_stacked
+                {"abstol": self.root_tol},
             )
+            try:
+                y0_alg = roots(y0_alg_guess, u_stacked).full().flatten()
+                success = True
+                message = None
+                # Check final output
+                fun = model.casadi_algebraic(
+                    time, casadi.vertcat(y0_diff, y0_alg), u_stacked
+                )
+            except RuntimeError as err:
+                success = False
+                message = err.args[0]
+                fun = None
         else:
             algebraic = model.algebraic_eval
             jac = model.jac_algebraic_eval
@@ -553,7 +557,7 @@ def solve(self, model, t_eval, external_variables=None, inputs=None):
                 )
             else:
                 t_eval_dimensionless = np.insert(
-                    t_eval_dimensionless, dindex, [dtime - eps, dtime + eps],
+                    t_eval_dimensionless, dindex, [dtime - eps, dtime + eps]
                 )
         end_indices.append(len(t_eval_dimensionless))
 

tests/unit/test_solvers/test_base_solver.py

@@ -72,7 +72,7 @@ def algebraic_eval(self, t, y):
         init_cond = solver.calculate_consistent_state(model)
         np.testing.assert_array_equal(init_cond, -2)
         # with casadi
-        solver_with_casadi = pybamm.BaseSolver(root_method="casadi")
+        solver_with_casadi = pybamm.BaseSolver(root_method="casadi", root_tol=1e-12)
         model = ScalarModel()
         init_cond = solver_with_casadi.calculate_consistent_state(model)
         np.testing.assert_array_equal(init_cond, -2)
@@ -160,7 +160,7 @@ def algebraic_eval(self, t, y):
         solver = pybamm.BaseSolver(root_method="casadi")
         with self.assertRaisesRegex(
             pybamm.SolverError,
-            "Could not find consistent initial conditions: solver terminated",
+            "Could not find consistent initial conditions: .../casadi",
         ):
             solver.calculate_consistent_state(Model())