diff --git a/CHANGELOG.md b/CHANGELOG.md index d34578b7b2..ae2ab2e198 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Features - Added option to express experiments (and extract solutions) in terms of cycles of operating condition ([#1309](https://github.com/pybamm-team/PyBaMM/pull/1309)) +- The event time and state are now returned as part of `Solution.t` and `Solution.y` so that the event is accurately captured in the returned solution ([#1300](https://github.com/pybamm-team/PyBaMM/pull/1300)) - Reformatted the `BasicDFNHalfCell` to be consistent with the other models ([#1282](https://github.com/pybamm-team/PyBaMM/pull/1282)) - Added option to make the total interfacial current density a state ([#1280](https://github.com/pybamm-team/PyBaMM/pull/1280)) - Added functionality to initialize a model using the solution from another model ([#1278](https://github.com/pybamm-team/PyBaMM/pull/1278)) @@ -58,7 +59,7 @@ This release introduces a new aging model for particle swelling and cracking, a ## Breaking changes -- The parameters "Positive/Negative particle distribution in x" and "Positive/Negative surface area to volume ratio distribution in x" have been deprecated. Instead, users can provide "Positive/Negative particle radius [m]" and "Positive/Negative surface area to volume ratio [m-1]" directly as functions of through-cell position (x [m]) ([#1237](https://github.com/pybamm-team/PyBaMM/pull/1237)) +- The parameters "Positive/Negative particle distribution in x" and "Positive/Negative surface area to volume ratio distribution in x" have been deprecated. Instead, users can provide "Positive/Negative particle radius [m]" and "Positive/Negative surface area to volume ratio [m-1]" directly as functions of through-cell position (x [m]) ([#1237](https://github.com/pybamm-team/PyBaMM/pull/1237)) # [v0.2.4](https://github.com/pybamm-team/PyBaMM/tree/v0.2.4) - 2020-09-07 diff --git a/pybamm/simulation.py b/pybamm/simulation.py index c99d227754..a1e31f4be7 100644 --- a/pybamm/simulation.py +++ b/pybamm/simulation.py @@ -330,13 +330,7 @@ def build(self, check_model=True): self._model_with_set_params, inplace=False, check_model=check_model ) - def solve( - self, - t_eval=None, - solver=None, - check_model=True, - **kwargs, - ): + def solve(self, t_eval=None, solver=None, check_model=True, **kwargs): """ A method to solve the model. This method will automatically build and set the model parameters if not already done so. @@ -518,12 +512,7 @@ def step(self, dt, solver=None, npts=2, save=True, **kwargs): solver = self.solver self._solution = solver.step( - self._solution, - self.built_model, - dt, - npts=npts, - save=save, - **kwargs, + self._solution, self.built_model, dt, npts=npts, save=save, **kwargs ) return self.solution diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index 60533452f3..29a45843b2 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -781,6 +781,10 @@ def step( "Cannot step empty model, use `pybamm.DummySolver` instead" ) + # Make sure dt is positive + if dt <= 0: + raise pybamm.SolverError("Step time must be positive") + # Set timer timer = pybamm.Timer() @@ -901,6 +905,16 @@ def get_termination_reason(self, solution, events): termination_event = min(final_event_values, key=final_event_values.get) # Add the event to the solution object solution.termination = "event: {}".format(termination_event) + # Update t, y and inputs to include event time and state + # Note: if the final entry of t is equal to the event time to within + # the absolute tolerance we skip this (having duplicate entries + # causes an error later in ProcessedVariable) + if solution.t_event - solution._t[-1] > self.atol: + solution._t = np.concatenate((solution._t, solution.t_event)) + solution._y = np.concatenate((solution._y, solution.y_event), axis=1) + for name, inp in solution.inputs.items(): + solution._inputs[name] = np.c_[inp, inp[:, -1]] + return "the termination event '{}' occurred".format(termination_event) def _set_up_ext_and_inputs(self, model, external_variables, inputs): diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index d05e6283f9..e23f1580b4 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -282,8 +282,8 @@ def event_fun(t): # append solution from the current step to solution solution.append(current_step_sol) solution.termination = "event" - solution.t_event = t_event - solution.y_event = y_event + solution.t_event = np.array([t_event]) + solution.y_event = y_event[:, np.newaxis] break else: diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py index 93a8a3731b..34476d7672 100644 --- a/pybamm/solvers/idaklu_solver.py +++ b/pybamm/solvers/idaklu_solver.py @@ -35,7 +35,7 @@ class IDAKLUSolver(pybamm.BaseSolver): 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-8). + The tolerance for the initial-condition solver (default is 1e-6). """ def __init__( diff --git a/pybamm/solvers/processed_variable.py b/pybamm/solvers/processed_variable.py index ca447f3195..966c5026a6 100644 --- a/pybamm/solvers/processed_variable.py +++ b/pybamm/solvers/processed_variable.py @@ -353,15 +353,11 @@ def initialise_2D(self): self.second_dimension = "x" self.r_sol = first_dim_pts self.x_sol = second_dim_pts - elif ( - self.domain[0] - in [ - "negative electrode", - "separator", - "positive electrode", - ] - and self.auxiliary_domains["secondary"] == ["current collector"] - ): + elif self.domain[0] in [ + "negative electrode", + "separator", + "positive electrode", + ] and self.auxiliary_domains["secondary"] == ["current collector"]: self.first_dimension = "x" self.second_dimension = "z" self.x_sol = first_dim_pts diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py index 6562e9c3ae..97d9acbc5b 100644 --- a/tests/integration/test_models/standard_output_tests.py +++ b/tests/integration/test_models/standard_output_tests.py @@ -603,10 +603,14 @@ def test_interfacial_current_average(self): multiplied by the interfacial current density is equal to the true value.""" + # The final time corresponds to an event (voltage cut-off). At this time + # the average property is satisfied but to a lesser degree of accuracy + t = self.t[:-1] + np.testing.assert_array_almost_equal( np.mean( - self.a_n(self.t, self.x_n) - * (self.j_n(self.t, self.x_n) + self.j_n_sei(self.t, self.x_n)), + self.a_n(t, self.x_n) + * (self.j_n(t, self.x_n) + self.j_n_sei(t, self.x_n)), axis=0, ), self.i_cell / self.l_n, @@ -614,8 +618,8 @@ def test_interfacial_current_average(self): ) np.testing.assert_array_almost_equal( np.mean( - self.a_p(self.t, self.x_p) - * (self.j_p(self.t, self.x_p) + self.j_p_sei(self.t, self.x_p)), + self.a_p(t, self.x_p) + * (self.j_p(t, self.x_p) + self.j_p_sei(t, self.x_p)), axis=0, ), -self.i_cell / self.l_p, diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py index c12f7bf6c3..1fb2c81c13 100644 --- a/tests/unit/test_solvers/test_base_solver.py +++ b/tests/unit/test_solvers/test_base_solver.py @@ -71,6 +71,11 @@ def test_nonmonotonic_teval(self): ): solver.solve(model, np.array([1, 2, 3, 2])) + # Check stepping with negative step size + dt = -1 + with self.assertRaisesRegex(pybamm.SolverError, "Step time must be positive"): + solver.step(None, model, dt) + def test_solution_time_length_fail(self): model = pybamm.BaseModel() v = pybamm.Scalar(1) diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py index 050805029a..e32ce3f69d 100644 --- a/tests/unit/test_solvers/test_scikits_solvers.py +++ b/tests/unit/test_solvers/test_scikits_solvers.py @@ -149,8 +149,8 @@ def test_model_solver_ode_events_python(self): 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)) - np.testing.assert_array_less(solution.y[0], 1.5) - np.testing.assert_array_less(solution.y[0], 1.25) + np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6) + np.testing.assert_array_less(solution.y[0], 1.25 + 1e-6) def test_model_solver_ode_jacobian_python(self): model = pybamm.BaseModel() @@ -249,8 +249,8 @@ def test_model_solver_dae_events_python(self): 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) - np.testing.assert_array_less(solution.y[-1], 2.5) + np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6) + np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-6) np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t)) np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t)) @@ -321,8 +321,8 @@ def nonsmooth_mult(t): # check solution for solution in [solution1, solution2]: - np.testing.assert_array_less(solution.y[0], 1.5) - np.testing.assert_array_less(solution.y[-1], 2.5) + np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6) + np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-6) var1_soln = np.exp(0.2 * solution.t) y0 = np.exp(0.2 * discontinuity) var1_soln[solution.t > discontinuity] = y0 * np.exp( @@ -388,8 +388,8 @@ def test_model_solver_dae_multiple_nonsmooth_python(self): # check solution for solution in [solution1, solution2]: - np.testing.assert_array_less(solution.y[0], 0.55) - np.testing.assert_array_less(solution.y[-1], 1.2) + np.testing.assert_array_less(solution.y[0], 0.55 + 1e-6) + np.testing.assert_array_less(solution.y[-1], 1.2 + 1e-6) var1_soln = (solution.t % a) ** 2 / 2 + a ** 2 / 2 * (solution.t // a) var2_soln = 2 * var1_soln np.testing.assert_allclose(solution.y[0], var1_soln, rtol=1e-06) @@ -569,8 +569,8 @@ def test_model_solver_ode_events_casadi(self): 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)) - np.testing.assert_array_less(solution.y[0], 1.5) - np.testing.assert_array_less(solution.y[0], 1.25) + np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6) + np.testing.assert_array_less(solution.y[0], 1.25 + 1e-6) def test_model_solver_dae_events_casadi(self): # Create model @@ -595,8 +595,8 @@ def test_model_solver_dae_events_casadi(self): solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) t_eval = np.linspace(0, 5, 100) solution = solver.solve(model_disc, t_eval) - np.testing.assert_array_less(solution.y[0], 1.5) - np.testing.assert_array_less(solution.y[-1], 2.5) + np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6) + np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-6) np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t)) np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t)) @@ -625,8 +625,8 @@ def test_model_solver_dae_inputs_events(self): 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) - np.testing.assert_array_less(solution.y[-1], 2.5) + np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6) + np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-6) np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t)) np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t)) @@ -730,8 +730,8 @@ def test_model_step_events(self): while time < end_time: step_solution = step_solver.step(step_solution, model, dt=dt, npts=10) time += dt - np.testing.assert_array_less(step_solution.y[0], 1.5) - np.testing.assert_array_less(step_solution.y[-1], 2.5001) + np.testing.assert_array_less(step_solution.y[0], 1.5 + 1e-6) + np.testing.assert_array_less(step_solution.y[-1], 2.5 + 1e-6) np.testing.assert_array_almost_equal( step_solution.y[0], np.exp(0.1 * step_solution.t), decimal=5 ) @@ -771,8 +771,8 @@ def test_model_step_nonsmooth_events(self): while time < end_time: step_solution = step_solver.step(step_solution, model, dt=dt, npts=10) time += dt - np.testing.assert_array_less(step_solution.y[0], 0.55) - np.testing.assert_array_less(step_solution.y[-1], 1.2) + np.testing.assert_array_less(step_solution.y[0], 0.55 + 1e-6) + np.testing.assert_array_less(step_solution.y[-1], 1.2 + 1e-6) var1_soln = (step_solution.t % a) ** 2 / 2 + a ** 2 / 2 * (step_solution.t // a) var2_soln = 2 * var1_soln np.testing.assert_array_almost_equal(step_solution.y[0], var1_soln, decimal=5) @@ -854,8 +854,8 @@ def nonsmooth_rate(t): # check solution for solution in [solution1, solution2]: - np.testing.assert_array_less(solution.y[0], 1.5) - np.testing.assert_array_less(solution.y[-1], 2.5) + np.testing.assert_array_less(solution.y[0], 1.5 + 1e-6) + np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-6) var1_soln = np.exp(0.2 * solution.t) y0 = np.exp(0.2 * discontinuity) var1_soln[solution.t > discontinuity] = y0 * np.exp( diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py index 7f8d30dcae..7c091b532c 100644 --- a/tests/unit/test_solvers/test_scipy_solver.py +++ b/tests/unit/test_solvers/test_scipy_solver.py @@ -99,7 +99,7 @@ def test_model_solver_with_event_python(self): t_eval = np.linspace(0, 10, 100) solution = solver.solve(model, t_eval) self.assertLess(len(solution.t), len(t_eval)) - np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)]) + np.testing.assert_array_equal(solution.t[:-1], t_eval[: len(solution.t) - 1]) np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t)) def test_model_solver_ode_with_jacobian_python(self): @@ -214,7 +214,7 @@ def test_model_solver_with_inputs(self): t_eval = np.linspace(0, 10, 100) solution = solver.solve(model, t_eval, inputs={"rate": 0.1}) self.assertLess(len(solution.t), len(t_eval)) - np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)]) + np.testing.assert_array_equal(solution.t[:-1], t_eval[: len(solution.t) - 1]) np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t)) def test_model_solver_with_external(self): @@ -272,7 +272,9 @@ def test_model_solver_with_event_with_casadi(self): t_eval = np.linspace(0, 10, 100) solution = solver.solve(model_disc, t_eval) self.assertLess(len(solution.t), len(t_eval)) - np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)]) + np.testing.assert_array_equal( + solution.t[:-1], t_eval[: len(solution.t) - 1] + ) np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t)) def test_model_solver_with_inputs_with_casadi(self): @@ -297,7 +299,7 @@ def test_model_solver_with_inputs_with_casadi(self): t_eval = np.linspace(0, 10, 100) solution = solver.solve(model, t_eval, inputs={"rate": 0.1}) self.assertLess(len(solution.t), len(t_eval)) - np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)]) + np.testing.assert_array_equal(solution.t[:-1], t_eval[: len(solution.t) - 1]) np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t)) def test_model_solver_inputs_in_initial_conditions(self):