diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a6c42188b..a3659b4742 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,14 +2,14 @@ ## Bug fixes -- For simulations with events that cause the simulation to stop early, the sensitivities could be evaluated incorrectly to zero ([#2331](https://github.com/pybamm-team/PyBaMM/pull/2337)) +- For simulations with events that cause the simulation to stop early, the sensitivities could be evaluated incorrectly to zero ([#2337](https://github.com/pybamm-team/PyBaMM/pull/2337)) ## Optimizations +- Added small perturbation to initial conditions for casadi solver. This seems to help the solver converge better in some cases ([#2356](https://github.com/pybamm-team/PyBaMM/pull/2356)) +- Added `ExplicitTimeIntegral` functionality to move variables which do not appear anywhere on the rhs to a new location, and to integrate those variables explicitly when `get` is called by the solution object. ([#2348](https://github.com/pybamm-team/PyBaMM/pull/2348)) - Added more rules for simplifying expressions ([#2211](https://github.com/pybamm-team/PyBaMM/pull/2211)) - - Sped up calculations of Electrode SOH variables for summary variables ([#2210](https://github.com/pybamm-team/PyBaMM/pull/2210)) -- Added `ExplicitTimeIntegral` functionality to move variables which do not appear anywhere on the rhs to a new location, and to integrate those variables explicitly when `get` is called by the solution object. ([#2348](https://github.com/pybamm-team/PyBaMM/pull/2348)) ## Breaking change diff --git a/README.md b/README.md index eaea0cf618..8027d4fa04 100644 --- a/README.md +++ b/README.md @@ -13,7 +13,9 @@ [![black code style](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/ambv/black) + [![All Contributors](https://img.shields.io/badge/all_contributors-46-orange.svg)](#-contributors) + diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 87b80d4fde..7b4ef000bf 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -66,6 +66,10 @@ class CasadiSolver(pybamm.BaseSolver): return_solution_if_failed_early : bool, optional Whether to return a Solution object if the solver fails to reach the end of the simulation, but managed to take some successful steps. Default is False. + perturb_algebraic_initial_conditions : bool, optional + Whether to perturb algebraic initial conditions to avoid a singularity. This + can sometimes slow down the solver, but is kept True as default for "safe" mode + as it seems to be more robust (False by default for other modes). """ def __init__( @@ -81,6 +85,7 @@ def __init__( extra_options_setup=None, extra_options_call=None, return_solution_if_failed_early=False, + perturb_algebraic_initial_conditions=None, ): super().__init__( "problem dependent", @@ -106,6 +111,17 @@ def __init__( self.extrap_tol = extrap_tol self.return_solution_if_failed_early = return_solution_if_failed_early + # Decide whether to perturb algebraic initial conditions, True by default for + # "safe" mode, False by default for other modes + if perturb_algebraic_initial_conditions is None: + if mode == "safe": + self.perturb_algebraic_initial_conditions = True + else: + self.perturb_algebraic_initial_conditions = False + else: + self.perturb_algebraic_initial_conditions = ( + perturb_algebraic_initial_conditions + ) self.name = "CasADi solver with '{}' mode".format(mode) # Initialize @@ -658,14 +674,22 @@ def _run_integrator( integrator = self.integrators[model]["no grid"] len_rhs = model.concatenated_rhs.size + len_alg = model.concatenated_algebraic.size # Check y0 to see if it includes sensitivities if explicit_sensitivities: num_parameters = model.len_rhs_sens // model.len_rhs len_rhs = len_rhs * (num_parameters + 1) + len_alg = len_alg * (num_parameters + 1) y0_diff = y0[:len_rhs] y0_alg = y0[len_rhs:] + if self.perturb_algebraic_initial_conditions and len_alg > 0: + # Add a tiny perturbation to the algebraic initial conditions + # For some reason this helps with convergence + # The actual value of the initial conditions for the algebraic variables + # doesn't matter + y0_alg = y0_alg * (1 + 1e-6 * casadi.DM(np.random.rand(len_alg))) pybamm.logger.spam("Finished preliminary setup for integrator run") # Solve @@ -680,9 +704,10 @@ def _run_integrator( casadi_sol = integrator( x0=y0_diff, z0=y0_alg, p=inputs_with_tmin, **self.extra_options_call ) - except RuntimeError as e: + except RuntimeError as error: # If it doesn't work raise error - raise pybamm.SolverError(e.args[0]) + pybamm.logger.debug(f"Casadi integrator failed with error {error}") + raise pybamm.SolverError(error.args[0]) pybamm.logger.debug("Finished casadi integrator") integration_time = timer.time() y_sol = casadi.vertcat(casadi_sol["xf"], casadi_sol["zf"]) @@ -711,9 +736,10 @@ def _run_integrator( casadi_sol = integrator( x0=x, z0=z, p=inputs_with_tlims, **self.extra_options_call ) - except RuntimeError as e: + except RuntimeError as error: # If it doesn't work raise error - raise pybamm.SolverError(e.args[0]) + pybamm.logger.debug(f"Casadi integrator failed with error {error}") + raise pybamm.SolverError(error.args[0]) integration_time = timer.time() x = casadi_sol["xf"] z = casadi_sol["zf"] diff --git a/tests/integration/test_models/standard_model_tests.py b/tests/integration/test_models/standard_model_tests.py index b36a9ca16a..9d8040e516 100644 --- a/tests/integration/test_models/standard_model_tests.py +++ b/tests/integration/test_models/standard_model_tests.py @@ -72,11 +72,6 @@ def test_solving( self.solver.rtol = 1e-8 self.solver.atol = 1e-8 - # Somehow removing an equation makes the solver fail at - # the low tolerances - if isinstance(self.model, pybamm.lithium_ion.NewmanTobias): - self.solver.rtol = 1e-7 - Crate = abs( self.parameter_values["Current function [A]"] / self.parameter_values["Nominal cell capacity [A.h]"] diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py index 8c24bac17b..a473842fee 100644 --- a/tests/unit/test_experiments/test_simulation_with_experiment.py +++ b/tests/unit/test_experiments/test_simulation_with_experiment.py @@ -476,11 +476,11 @@ def test_run_experiment_skip_steps(self): model, parameter_values=parameter_values, experiment=experiment2 ) sol2 = sim2.solve() - np.testing.assert_array_equal( + np.testing.assert_array_almost_equal( sol["Terminal voltage [V]"].data, sol2["Terminal voltage [V]"].data ) for idx1, idx2 in [(1, 0), (2, 1), (4, 2)]: - np.testing.assert_array_equal( + np.testing.assert_array_almost_equal( sol.cycles[0].steps[idx1]["Terminal voltage [V]"].data, sol2.cycles[0].steps[idx2]["Terminal voltage [V]"].data, ) diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py index 30923c40c2..6b9d073ab6 100644 --- a/tests/unit/test_solvers/test_casadi_solver.py +++ b/tests/unit/test_solvers/test_casadi_solver.py @@ -24,7 +24,12 @@ def test_model_solver(self): disc = pybamm.Discretisation() model_disc = disc.process_model(model, inplace=False) # Solve - solver = pybamm.CasadiSolver(mode="fast", rtol=1e-8, atol=1e-8) + solver = pybamm.CasadiSolver( + mode="fast", + rtol=1e-8, + atol=1e-8, + perturb_algebraic_initial_conditions=False, # added for coverage + ) t_eval = np.linspace(0, 1, 100) solution = solver.solve(model_disc, t_eval) np.testing.assert_array_equal(solution.t, t_eval)