diff --git a/CHANGELOG.md b/CHANGELOG.md index 4c1ba6a8a3..29a5f3a55d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- Changed rootfinding algorithm to CasADi, scipy.optimize.root still accessible as an option ([#844](https://github.com/pybamm-team/PyBaMM/pull/844)) - Added capacitance effects to lithium-ion models ([#842](https://github.com/pybamm-team/PyBaMM/pull/842)) - Added NCA parameter set ([#824](https://github.com/pybamm-team/PyBaMM/pull/824)) - Added functionality to `Solution` that automatically gets `t_eval` from the data when simulating drive cycles and performs checks to ensure the output has the required resolution to accurately capture the input current ([#819](https://github.com/pybamm-team/PyBaMM/pull/819)) diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py index 16d17eef28..fd01abfa2c 100644 --- a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py +++ b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py @@ -177,11 +177,17 @@ def __init__(self, name="Doyle-Fuller-Newman model"): # Boundary conditions must be provided for equations with spatial derivatives self.boundary_conditions[c_s_n] = { "left": (pybamm.Scalar(0), "Neumann"), - "right": (-param.C_n * j_n / param.a_n, "Neumann"), + "right": ( + -param.C_n * j_n / param.a_n / param.D_n(c_s_surf_n, T), + "Neumann", + ), } self.boundary_conditions[c_s_p] = { "left": (pybamm.Scalar(0), "Neumann"), - "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"), + "right": ( + -param.C_p * j_p / param.a_p / param.gamma_p / param.D_p(c_s_surf_p, T), + "Neumann", + ), } # c_n_init and c_p_init can in general be functions of x # Note the broadcasting, for domains @@ -195,14 +201,22 @@ def __init__(self, name="Doyle-Fuller-Newman model"): self.initial_conditions[c_s_p] = param.c_p_init(x_p) # Events specify points at which a solution should terminate self.events += [ - pybamm.Event("Minimum negative particle surface concentration", - pybamm.min(c_s_surf_n) - 0.01), - pybamm.Event("Maximum negative particle surface concentration", - (1 - 0.01) - pybamm.max(c_s_surf_n)), - pybamm.Event("Minimum positive particle surface concentration", - pybamm.min(c_s_surf_p) - 0.01), - pybamm.Event("Maximum positive particle surface concentration", - (1 - 0.01) - pybamm.max(c_s_surf_p)), + pybamm.Event( + "Minimum negative particle surface concentration", + pybamm.min(c_s_surf_n) - 0.01, + ), + pybamm.Event( + "Maximum negative particle surface concentration", + (1 - 0.01) - pybamm.max(c_s_surf_n), + ), + pybamm.Event( + "Minimum positive particle surface concentration", + pybamm.min(c_s_surf_p) - 0.01, + ), + pybamm.Event( + "Maximum positive particle surface concentration", + (1 - 0.01) - pybamm.max(c_s_surf_p), + ), ] ###################### # Current in the solid @@ -257,8 +271,11 @@ def __init__(self, name="Doyle-Fuller-Newman model"): "right": (pybamm.Scalar(0), "Neumann"), } self.initial_conditions[c_e] = param.c_e_init - self.events.append(pybamm.Event("Zero electrolyte concentration cut-off", - pybamm.min(c_e) - 0.002)) + self.events.append( + pybamm.Event( + "Zero electrolyte concentration cut-off", pybamm.min(c_e) - 0.002 + ) + ) ###################### # (Some) variables diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py index 15d8e2ab92..60956d210d 100644 --- a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py +++ b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py @@ -80,34 +80,48 @@ def __init__(self, name="Single Particle Model"): N_s_p = -param.D_p(c_s_p, T) * pybamm.grad(c_s_p) self.rhs[c_s_n] = -(1 / param.C_n) * pybamm.div(N_s_n) self.rhs[c_s_p] = -(1 / param.C_p) * pybamm.div(N_s_p) + # Surf takes the surface value of a variable, i.e. its boundary value on the + # right side. This is also accessible via `boundary_value(x, "right")`, with + # "left" providing the boundary value of the left side + c_s_surf_n = pybamm.surf(c_s_n) + c_s_surf_p = pybamm.surf(c_s_p) # Boundary conditions must be provided for equations with spatial derivatives self.boundary_conditions[c_s_n] = { "left": (pybamm.Scalar(0), "Neumann"), - "right": (-param.C_n * j_n / param.a_n, "Neumann"), + "right": ( + -param.C_n * j_n / param.a_n / param.D_n(c_s_surf_n, T), + "Neumann", + ), } self.boundary_conditions[c_s_p] = { "left": (pybamm.Scalar(0), "Neumann"), - "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"), + "right": ( + -param.C_p * j_p / param.a_p / param.gamma_p / param.D_p(c_s_surf_p, T), + "Neumann", + ), } # c_n_init and c_p_init are functions, but for the SPM we evaluate them at x=0 # and x=1 since there is no x-dependence in the particles self.initial_conditions[c_s_n] = param.c_n_init(0) self.initial_conditions[c_s_p] = param.c_p_init(1) - # Surf takes the surface value of a variable, i.e. its boundary value on the - # right side. This is also accessible via `boundary_value(x, "right")`, with - # "left" providing the boundary value of the left side - c_s_surf_n = pybamm.surf(c_s_n) - c_s_surf_p = pybamm.surf(c_s_p) # Events specify points at which a solution should terminate self.events += [ - pybamm.Event("Minimum negative particle surface concentration", - pybamm.min(c_s_surf_n) - 0.01), - pybamm.Event("Maximum negative particle surface concentration", - (1 - 0.01) - pybamm.max(c_s_surf_n)), - pybamm.Event("Minimum positive particle surface concentration", - pybamm.min(c_s_surf_p) - 0.01), - pybamm.Event("Maximum positive particle surface concentration", - (1 - 0.01) - pybamm.max(c_s_surf_p)), + pybamm.Event( + "Minimum negative particle surface concentration", + pybamm.min(c_s_surf_n) - 0.01, + ), + pybamm.Event( + "Maximum negative particle surface concentration", + (1 - 0.01) - pybamm.max(c_s_surf_n), + ), + pybamm.Event( + "Minimum positive particle surface concentration", + pybamm.min(c_s_surf_p) - 0.01, + ), + pybamm.Event( + "Maximum positive particle surface concentration", + (1 - 0.01) - pybamm.max(c_s_surf_p), + ), ] # Note that the SPM does not have any algebraic equations, so the `algebraic` diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index 1c031b8b8f..806015e900 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -23,7 +23,10 @@ class BaseSolver(object): 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") + The method to use to find initial conditions (default is "casadi"). If "casadi", + the solver uses casadi's Newton rootfinding algorithm to find initial + conditions. Otherwise, the solver uses 'scipy.optimize.root' with method + specified by 'root_method' (e.g. "lm", "hybr", ...) root_tol : float, optional The tolerance for the initial-condition solver (default is 1e-6). max_steps: int, optional @@ -36,7 +39,7 @@ def __init__( method=None, rtol=1e-6, atol=1e-6, - root_method="lm", + root_method="casadi", root_tol=1e-6, max_steps=1000, ): @@ -125,10 +128,11 @@ def set_up(self, model, inputs=None): "Cannot use ODE solver '{}' to solve DAE model".format(self.name) ) + if self.ode_solver is True: + self.root_method = None if ( - isinstance(self, pybamm.CasadiSolver) - and model.convert_to_format != "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" ) @@ -268,6 +272,16 @@ def report(string): model.terminate_events_eval = terminate_events_eval model.discontinuity_events_eval = discontinuity_events_eval + # Save CasADi functions for the CasADi solver + # 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] + ) + model.casadi_algebraic = algebraic # Calculate consistent initial conditions for the algebraic equations if len(model.algebraic) > 0: all_states = pybamm.NumpyConcatenation( @@ -278,26 +292,13 @@ def report(string): model.residuals_eval = residuals_eval model.jacobian_eval = jacobian_eval y0_guess = y0.flatten() - model.y0 = self.calculate_consistent_state(model, 0, y0_guess) + model.y0 = self.calculate_consistent_state(model, 0, y0_guess, inputs) else: # can use DAE solver to solve ODE model model.residuals_eval = Residuals(rhs, "residuals", model) model.jacobian_eval = jac_rhs model.y0 = y0.flatten() - # Save CasADi functions for the CasADi solver - # 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 model.convert_to_format == "casadi" and 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] - ) - model.casadi_algebraic = algebraic - pybamm.logger.info("Finish solver set-up") def set_inputs(self, model, ext_and_inputs): @@ -319,7 +320,7 @@ def set_inputs(self, model, ext_and_inputs): if model.jacobian_eval: model.jacobian_eval.set_inputs(ext_and_inputs) - def calculate_consistent_state(self, model, time=0, y0_guess=None): + def calculate_consistent_state(self, model, time=0, y0_guess=None, inputs=None): """ Calculate consistent state for the algebraic equations through root-finding @@ -328,6 +329,12 @@ def calculate_consistent_state(self, model, time=0, y0_guess=None): ---------- model : :class:`pybamm.BaseModel` The model for which to calculate initial conditions. + time : float + The time at which to calculate the states + y0_guess : :class:`np.array` + Guess for the rootfinding + inputs : dict, optional + Any input parameters to pass to the model when solving Returns ------- @@ -336,65 +343,101 @@ def calculate_consistent_state(self, model, time=0, y0_guess=None): of the algebraic equations) """ pybamm.logger.info("Start calculating consistent states") - rhs = model.rhs_eval - algebraic = model.algebraic_eval - jac = model.jac_algebraic_eval if y0_guess is None: y0_guess = model.concatenated_initial_conditions.flatten() # Split y0_guess into differential and algebraic - len_rhs = rhs(time, y0_guess).shape[0] + len_rhs = model.rhs_eval(time, y0_guess).shape[0] y0_diff, y0_alg_guess = np.split(y0_guess, [len_rhs]) + inputs = inputs or {} - def root_fun(y0_alg): - "Evaluates algebraic using y0_diff (fixed) and y0_alg (changed by algo)" - y0 = np.concatenate([y0_diff, y0_alg]) - out = algebraic(time, y0) - pybamm.logger.debug( - "Evaluating algebraic equations at t={}, L2-norm is {}".format( - time * model.timescale, np.linalg.norm(out) - ) + # Solve using casadi or scipy + if self.root_method == "casadi": + # Set up + u_stacked = casadi.vertcat(*[x for x in inputs.values()]) + u = casadi.MX.sym("u", u_stacked.shape[0]) + y_alg = casadi.MX.sym("y_alg", y0_alg_guess.shape[0]) + 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", + dict(x=y_alg, p=u, g=alg_root), + {"abstol": self.root_tol}, ) - return out - - if jac: - if issparse(jac(0, y0_guess)): + 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 + + def root_fun(y0_alg): + "Evaluates algebraic using y0_diff (fixed) and y0_alg (changed by algo)" + y0 = np.concatenate([y0_diff, y0_alg]) + out = algebraic(time, y0) + pybamm.logger.debug( + "Evaluating algebraic equations at t={}, L2-norm is {}".format( + time * model.timescale, np.linalg.norm(out) + ) + ) + return out - def jac_fn(y0_alg): - """ - Evaluates jacobian using y0_diff (fixed) and y0_alg (varying) - """ - y0 = np.concatenate([y0_diff, y0_alg]) - return jac(0, y0)[:, len_rhs:].toarray() + if jac: + if issparse(jac(0, y0_guess)): - else: + def jac_fn(y0_alg): + """ + Evaluates jacobian using y0_diff (fixed) and y0_alg (varying) + """ + y0 = np.concatenate([y0_diff, y0_alg]) + return jac(0, y0)[:, len_rhs:].toarray() - def jac_fn(y0_alg): - """ - Evaluates jacobian using y0_diff (fixed) and y0_alg (varying) - """ - y0 = np.concatenate([y0_diff, y0_alg]) - return jac(0, y0)[:, len_rhs:] + else: - else: - jac_fn = None - # Find the values of y0_alg that are roots of the algebraic equations - sol = optimize.root( - root_fun, - y0_alg_guess, - jac=jac_fn, - method=self.root_method, - tol=self.root_tol, - ) - # Return full set of consistent initial conditions (y0_diff unchanged) - y0_consistent = np.concatenate([y0_diff, sol.x]) + def jac_fn(y0_alg): + """ + Evaluates jacobian using y0_diff (fixed) and y0_alg (varying) + """ + y0 = np.concatenate([y0_diff, y0_alg]) + return jac(0, y0)[:, len_rhs:] - if sol.success and np.all(sol.fun < self.root_tol * len(sol.x)): + else: + jac_fn = None + # Find the values of y0_alg that are roots of the algebraic equations + sol = optimize.root( + root_fun, + y0_alg_guess, + jac=jac_fn, + method=self.root_method, + tol=self.root_tol, + ) + # Set outputs + y0_alg = sol.x + success = sol.success + fun = sol.fun + message = sol.message + + if success and np.all(fun < self.root_tol * len(y0_alg)): + # 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") return y0_consistent - elif not sol.success: + elif not success: raise pybamm.SolverError( - "Could not find consistent initial conditions: {}".format(sol.message) + "Could not find consistent initial conditions: {}".format(message) ) else: raise pybamm.SolverError( @@ -402,7 +445,7 @@ def jac_fn(y0_alg): Could not find consistent initial conditions: solver terminated successfully, but maximum solution error ({}) above tolerance ({}) """.format( - np.max(sol.fun), self.root_tol * len(sol.x) + np.max(fun), self.root_tol * len(y0_alg) ) ) @@ -514,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)) @@ -548,7 +591,10 @@ def solve(self, model, t_eval, external_variables=None, inputs=None): last_state = solution.y[:, -1] if len(model.algebraic) > 0: model.y0 = self.calculate_consistent_state( - model, t_eval_dimensionless[end_index], last_state + model, + t_eval_dimensionless[end_index], + last_state, + ext_and_inputs, ) else: model.y0 = last_state diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index ea809fda7c..34fd740438 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -47,7 +47,7 @@ def __init__( mode="safe", rtol=1e-6, atol=1e-6, - root_method="lm", + root_method="casadi", root_tol=1e-6, max_step_decrease_count=5, **extra_options, diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py index 795f3361b8..be4a1932cc 100644 --- a/pybamm/solvers/idaklu_solver.py +++ b/pybamm/solvers/idaklu_solver.py @@ -36,7 +36,7 @@ class IDAKLUSolver(pybamm.BaseSolver): """ def __init__( - self, rtol=1e-6, atol=1e-6, root_method="lm", root_tol=1e-6, max_steps=1000 + self, rtol=1e-6, atol=1e-6, root_method="casadi", root_tol=1e-6, max_steps=1000 ): if idaklu_spec is None: diff --git a/pybamm/solvers/scikits_dae_solver.py b/pybamm/solvers/scikits_dae_solver.py index 5f9999cbb7..74620cafed 100644 --- a/pybamm/solvers/scikits_dae_solver.py +++ b/pybamm/solvers/scikits_dae_solver.py @@ -40,7 +40,7 @@ def __init__( method="ida", rtol=1e-6, atol=1e-6, - root_method="lm", + root_method="casadi", root_tol=1e-6, max_steps=1000, ): diff --git a/setup.py b/setup.py index db141d4f4d..fa13bd203a 100644 --- a/setup.py +++ b/setup.py @@ -492,7 +492,7 @@ def load_version(): # List of dependencies install_requires=[ "numpy>=1.16", - "scipy>=1.0", + "scipy>=1.3", "pandas>=0.24", "anytree>=2.4.3", "autograd>=1.2", diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py index ce44197a34..7ab105982f 100644 --- a/tests/unit/test_solvers/test_base_solver.py +++ b/tests/unit/test_solvers/test_base_solver.py @@ -1,6 +1,7 @@ # # Tests for the Base Solver class # +import casadi import pybamm import numpy as np from scipy.sparse import csr_matrix @@ -49,9 +50,16 @@ def test_ode_solver_fail_with_dae(self): def test_find_consistent_initial_conditions(self): # Simple system: a single algebraic equation class ScalarModel: - concatenated_initial_conditions = np.array([[2]]) - jac_algebraic_eval = None - timescale = 1 + def __init__(self): + self.concatenated_initial_conditions = np.array([[2]]) + self.jac_algebraic_eval = None + self.timescale = 1 + t = casadi.MX.sym("t") + y = casadi.MX.sym("y") + u = casadi.MX.sym("u") + self.casadi_algebraic = casadi.Function( + "alg", [t, y, u], [self.algebraic_eval(t, y)] + ) def rhs_eval(self, t, y): return np.array([]) @@ -59,18 +67,30 @@ def rhs_eval(self, t, y): def algebraic_eval(self, t, y): return y + 2 - solver = pybamm.BaseSolver() + solver = pybamm.BaseSolver(root_method="lm") model = ScalarModel() 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", root_tol=1e-12) + model = ScalarModel() + init_cond = solver_with_casadi.calculate_consistent_state(model) + np.testing.assert_array_equal(init_cond, -2) # More complicated system vec = np.array([0.0, 1.0, 1.5, 2.0]) class VectorModel: - concatenated_initial_conditions = np.zeros_like(vec) - jac_algebraic_eval = None - timescale = 1 + def __init__(self): + self.concatenated_initial_conditions = np.zeros_like(vec) + self.jac_algebraic_eval = None + self.timescale = 1 + t = casadi.MX.sym("t") + y = casadi.MX.sym("y", vec.size) + u = casadi.MX.sym("u") + self.casadi_algebraic = casadi.Function( + "alg", [t, y, u], [self.algebraic_eval(t, y)] + ) def rhs_eval(self, t, y): return y[0:1] @@ -81,6 +101,9 @@ def algebraic_eval(self, t, y): model = VectorModel() init_cond = solver.calculate_consistent_state(model) np.testing.assert_array_almost_equal(init_cond, vec) + # with casadi + init_cond = solver_with_casadi.calculate_consistent_state(model) + np.testing.assert_array_almost_equal(init_cond, vec) # With jacobian def jac_dense(t, y): @@ -102,9 +125,16 @@ def jac_sparse(t, y): def test_fail_consistent_initial_conditions(self): class Model: - concatenated_initial_conditions = np.array([2]) - jac_algebraic_eval = None - timescale = 1 + def __init__(self): + self.concatenated_initial_conditions = np.array([2]) + self.jac_algebraic_eval = None + self.timescale = 1 + t = casadi.MX.sym("t") + y = casadi.MX.sym("y") + u = casadi.MX.sym("u") + self.casadi_algebraic = casadi.Function( + "alg", [t, y, u], [self.algebraic_eval(t, y)] + ) def rhs_eval(self, t, y): return np.array([]) @@ -120,12 +150,19 @@ def algebraic_eval(self, t, y): "Could not find consistent initial conditions: The iteration is not making", ): solver.calculate_consistent_state(Model()) - solver = pybamm.BaseSolver() + solver = pybamm.BaseSolver(root_method="lm") with self.assertRaisesRegex( pybamm.SolverError, "Could not find consistent initial conditions: solver terminated", ): solver.calculate_consistent_state(Model()) + # with casadi + solver = pybamm.BaseSolver(root_method="casadi") + with self.assertRaisesRegex( + pybamm.SolverError, + "Could not find consistent initial conditions: .../casadi", + ): + solver.calculate_consistent_state(Model()) def test_time_too_short(self): solver = pybamm.BaseSolver() diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py index 9683b5cbc6..703673ae8e 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -28,7 +28,7 @@ def test_ida_roberts_klu(self): disc = pybamm.Discretisation() disc.process_model(model) - solver = pybamm.IDAKLUSolver() + solver = pybamm.IDAKLUSolver(root_method="lm") t_eval = np.linspace(0, 3, 100) solution = solver.solve(model, t_eval) @@ -58,7 +58,7 @@ def test_set_atol(self): mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts) disc = pybamm.Discretisation(mesh, model.default_spatial_methods) disc.process_model(model) - solver = pybamm.IDAKLUSolver() + solver = pybamm.IDAKLUSolver(root_method="lm") variable_tols = {"Electrolyte concentration": 1e-3} solver.set_atol_by_variable(variable_tols, model) @@ -76,7 +76,7 @@ def test_failures(self): disc = pybamm.Discretisation() disc.process_model(model) - solver = pybamm.IDAKLUSolver() + solver = pybamm.IDAKLUSolver(root_method="lm") t_eval = np.linspace(0, 3, 100) with self.assertRaisesRegex(pybamm.SolverError, "KLU requires the Jacobian"): diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py index e3b39791a9..905507c692 100644 --- a/tests/unit/test_solvers/test_scikits_solvers.py +++ b/tests/unit/test_solvers/test_scikits_solvers.py @@ -194,7 +194,7 @@ def test_model_solver_dae_python(self): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") t_eval = np.linspace(0, 1, 100) solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) @@ -214,7 +214,7 @@ def test_model_solver_dae_bad_ics_python(self): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") t_eval = np.linspace(0, 1, 100) solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) @@ -238,7 +238,7 @@ def test_model_solver_dae_events_python(self): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") t_eval = np.linspace(0, 5, 100) solution = solver.solve(model, t_eval) np.testing.assert_array_less(solution.y[0], 1.5) @@ -283,7 +283,7 @@ def nonsmooth_mult(t): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") # create two time series, one without a time point on the discontinuity, # and one with @@ -355,7 +355,7 @@ def jacobian(t, y): model.jacobian = jacobian # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") t_eval = np.linspace(0, 1, 100) solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) @@ -372,7 +372,7 @@ def test_solve_ode_model_with_dae_solver_python(self): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") t_eval = np.linspace(0, 1, 100) solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) @@ -421,7 +421,7 @@ def test_model_step_dae_python(self): disc = get_discretisation_for_testing() disc.process_model(model) - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") # Step once dt = 1 @@ -514,7 +514,10 @@ def test_model_solver_dae_inputs_events(self): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + if form == "python": + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") + else: + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) t_eval = np.linspace(0, 5, 100) solution = solver.solve(model, t_eval, inputs={"rate 1": 0.1, "rate 2": 2}) np.testing.assert_array_less(solution.y[0], 1.5)