Merge pull request #2773 from pybamm-team/issue-2760-step

Issue 2760 step

Author: valentinsulzer

CHANGELOG.md

@@ -18,6 +18,7 @@
 
 ## Breaking changes
 
+- When using `solver.step()`, the first time point in the step is shifted by `pybamm.settings.step_start_offset` (default 1 ns) to avoid having duplicate times in the solution steps from the end of one step and the start of the next. ([#2773](https://github.com/pybamm-team/PyBaMM/pull/2773))
 - Renamed "Measured open circuit voltage [V]" to "Surface open-circuit voltage [V]". This variable was calculated from surface particle concentrations, and hence "hid" the overpotential from particle gradients. The new variable "Open-circuit voltage [V]" is calculated from bulk particle concentrations instead. ([#2740](https://github.com/pybamm-team/PyBaMM/pull/2740))
 - Renamed all references to "open circuit" to be "open-circuit" instead. ([#2740](https://github.com/pybamm-team/PyBaMM/pull/2740))
 - Renamed parameter "1 + dlnf/dlnc" to "Thermodynamic factor". ([#2727](https://github.com/pybamm-team/PyBaMM/pull/2727))

pybamm/settings.py

@@ -12,6 +12,7 @@ class Settings(object):
     _abs_smoothing = "exact"
     max_words_in_line = 4
     max_y_value = 1e5
+    step_start_offset = 1e-9
     tolerances = {
         "D_e__c_e": 10,  # dimensional
         "kappa_e__c_e": 10,  # dimensional

pybamm/solvers/base_solver.py

@@ -1085,18 +1085,34 @@ 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")
+        # Make sure dt is greater than the offset
+        step_start_offset = pybamm.settings.step_start_offset
+        if dt <= step_start_offset:
+            raise pybamm.SolverError(
+                f"Step time must be at least {pybamm.TimerTime(step_start_offset)}"
+            )
+
+        t_start = old_solution.t[-1]
+        t_end = t_start + dt
+        # Calculate t_eval
+        t_eval = np.linspace(t_start, t_end, npts)
+
+        if t_start == 0:
+            t_start_shifted = t_start
+        else:
+            # offset t_start by t_start_offset (default 1 ns)
+            # to avoid repeated times in the solution
+            # from having the same time at the end of the previous step and
+            # the start of the next step
+            t_start_shifted = t_start + step_start_offset
+            t_eval[0] = t_start_shifted
 
         # Set timer
         timer = pybamm.Timer()
 
         # Set up inputs
         model_inputs = self._set_up_model_inputs(model, inputs)
 
-        t = old_solution.t[-1]
-
         first_step_this_model = False
         if model not in self._model_set_up:
             first_step_this_model = True
@@ -1139,16 +1155,17 @@ def step(
         set_up_time = timer.time()
 
         # (Re-)calculate consistent initial conditions
-        self._set_initial_conditions(model, t, model_inputs, update_rhs=False)
-
-        # Calculate t_eval
-        t_eval = np.linspace(t, t + dt, npts)
+        self._set_initial_conditions(
+            model, t_start_shifted, model_inputs, update_rhs=False
+        )
 
         # Check initial conditions don't violate events
         self._check_events_with_initial_conditions(t_eval, model, model_inputs)
 
         # Step
-        pybamm.logger.verbose("Stepping for {:.0f} < t < {:.0f}".format(t, (t + dt)))
+        pybamm.logger.verbose(
+            "Stepping for {:.0f} < t < {:.0f}".format(t_start_shifted, t_end)
+        )
         timer.reset()
         solution = self._integrate(model, t_eval, model_inputs)
         solution.solve_time = timer.time()

tests/integration/test_experiments.py

@@ -14,20 +14,21 @@ def test_discharge_rest_charge(self):
                 "Rest for 1 hour",
                 "Charge at C/2 for 1 hour",
             ],
-            period="0.25 hours",
+            period="0.5 hours",
         )
         model = pybamm.lithium_ion.SPM()
         sim = pybamm.Simulation(
             model, experiment=experiment, solver=pybamm.CasadiSolver()
         )
         sim.solve()
         np.testing.assert_array_almost_equal(
-            sim._solution["Time [h]"].entries, np.linspace(0, 3, 13)
+            sim._solution["Time [h]"].entries,
+            np.array([0, 0.5, 1, 1 + 1e-9, 1.5, 2, 2 + 1e-9, 2.5, 3]),
         )
         cap = model.default_parameter_values["Nominal cell capacity [A.h]"]
         np.testing.assert_array_almost_equal(
             sim._solution["Current [A]"].entries,
-            [cap / 2] * 5 + [0] * 4 + [-cap / 2] * 4,
+            [cap / 2] * 3 + [0] * 3 + [-cap / 2] * 3,
         )
 
     def test_rest_discharge_rest(self):
@@ -57,13 +58,10 @@ def test_gitt(self):
             model, experiment=experiment, solver=pybamm.CasadiSolver()
         )
         sim.solve()
-        np.testing.assert_array_almost_equal(
-            sim._solution["Time [h]"].entries, np.arange(0, 20.01, 0.1)
-        )
         cap = model.default_parameter_values["Nominal cell capacity [A.h]"]
         np.testing.assert_array_almost_equal(
             sim._solution["Current [A]"].entries,
-            [cap / 20] * 11 + [0] * 10 + ([cap / 20] * 10 + [0] * 10) * 9,
+            [cap / 20] * 11 + [0] * 11 + ([cap / 20] * 11 + [0] * 11) * 9,
         )
 
     def test_infeasible(self):

tests/unit/test_experiments/test_simulation_with_experiment.py

@@ -176,22 +176,18 @@ def test_run_experiment_cccv_ode(self):
             solution = sim.solve(solver=pybamm.CasadiSolver("fast with events"))
             solutions.append(solution)
 
+        t = solutions[1]["Time [s]"].data
         np.testing.assert_array_almost_equal(
-            solutions[0]["Voltage [V]"].data,
-            solutions[1]["Voltage [V]"].data,
+            solutions[0]["Voltage [V]"](t=t),
+            solutions[1]["Voltage [V]"](t=t),
             decimal=1,
         )
         np.testing.assert_array_almost_equal(
-            solutions[0]["Current [A]"].data,
-            solutions[1]["Current [A]"].data,
+            solutions[0]["Current [A]"](t=t),
+            solutions[1]["Current [A]"](t=t),
             decimal=0,
         )
 
-        np.testing.assert_array_equal(
-            solutions[0]["Ambient temperature [C]"].data,
-            solutions[1]["Ambient temperature [C]"].data,
-        )
-
         self.assertEqual(solutions[1].termination, "final time")
 
     @unittest.skipIf(not pybamm.have_idaklu(), "idaklu solver is not installed")
@@ -516,12 +512,13 @@ def test_run_experiment_skip_steps(self):
         )
         sol2 = sim2.solve()
         np.testing.assert_array_almost_equal(
-            sol["Voltage [V]"].data, sol2["Voltage [V]"].data
+            sol["Voltage [V]"].data, sol2["Voltage [V]"].data, decimal=5
         )
         for idx1, idx2 in [(1, 0), (2, 1), (4, 2)]:
             np.testing.assert_array_almost_equal(
                 sol.cycles[0].steps[idx1]["Voltage [V]"].data,
                 sol2.cycles[0].steps[idx2]["Voltage [V]"].data,
+                decimal=5,
             )
 
     def test_all_empty_solution_errors(self):

tests/unit/test_simulation.py

@@ -148,34 +148,27 @@ def test_step(self):
         sim = pybamm.Simulation(model)
 
         sim.step(dt)  # 1 step stores first two points
-        self.assertEqual(sim.solution.t.size, 2)
         self.assertEqual(sim.solution.y.full()[0, :].size, 2)
-        self.assertEqual(sim.solution.t[0], 0)
-        self.assertEqual(sim.solution.t[1], dt)
+        np.testing.assert_array_almost_equal(sim.solution.t, np.array([0, dt]))
         saved_sol = sim.solution
 
         sim.step(dt)  # automatically append the next step
-        self.assertEqual(sim.solution.t.size, 3)
-        self.assertEqual(sim.solution.y.full()[0, :].size, 3)
-        self.assertEqual(sim.solution.t[0], 0)
-        self.assertEqual(sim.solution.t[1], dt)
-        self.assertEqual(sim.solution.t[2], 2 * dt)
+        self.assertEqual(sim.solution.y.full()[0, :].size, 4)
+        np.testing.assert_array_almost_equal(
+            sim.solution.t, np.array([0, dt, dt + 1e-9, 2 * dt])
+        )
 
         sim.step(dt, save=False)  # now only store the two end step points
-        self.assertEqual(sim.solution.t.size, 2)
         self.assertEqual(sim.solution.y.full()[0, :].size, 2)
-        self.assertEqual(sim.solution.t[0], 2 * dt)
-        self.assertEqual(sim.solution.t[1], 3 * dt)
-
+        np.testing.assert_array_almost_equal(
+            sim.solution.t, np.array([2 * dt + 1e-9, 3 * dt])
+        )
         # Start from saved solution
-        sim.step(
-            dt, starting_solution=saved_sol
-        )  # now only store the two end step points
-        self.assertEqual(sim.solution.t.size, 3)
-        self.assertEqual(sim.solution.y.full()[0, :].size, 3)
-        self.assertEqual(sim.solution.t[0], 0)
-        self.assertEqual(sim.solution.t[1], dt)
-        self.assertEqual(sim.solution.t[2], 2 * dt)
+        sim.step(dt, starting_solution=saved_sol)
+        self.assertEqual(sim.solution.y.full()[0, :].size, 4)
+        np.testing.assert_array_almost_equal(
+            sim.solution.t, np.array([0, dt, dt + 1e-9, 2 * dt])
+        )
 
     def test_solve_with_initial_soc(self):
         model = pybamm.lithium_ion.SPM()
@@ -238,11 +231,10 @@ def test_step_with_inputs(self):
         sim.step(
             dt, inputs={"Current function [A]": 2}
         )  # automatically append the next step
-        self.assertEqual(sim.solution.t.size, 3)
-        self.assertEqual(sim.solution.y.full()[0, :].size, 3)
-        self.assertEqual(sim.solution.t[0], 0)
-        self.assertEqual(sim.solution.t[1], dt)
-        self.assertEqual(sim.solution.t[2], 2 * dt)
+        self.assertEqual(sim.solution.y.full()[0, :].size, 4)
+        np.testing.assert_array_almost_equal(
+            sim.solution.t, np.array([0, dt, dt + 1e-9, 2 * dt])
+        )
         np.testing.assert_array_equal(
             sim.solution.all_inputs[1]["Current function [A]"], 2
         )

tests/unit/test_solvers/test_base_solver.py

@@ -67,9 +67,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"):
+        # Check stepping with step size too small
+        dt = 1e-9
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "Step time must be at least 1.000 ns"
+        ):
             solver.step(None, model, dt)
 
     def test_solution_time_length_fail(self):

tests/unit/test_solvers/test_casadi_solver.py

@@ -292,7 +292,7 @@ def test_model_step(self):
         # Step again (return 5 points)
         step_sol_2 = solver.step(step_sol, model, dt, npts=5)
         np.testing.assert_array_equal(
-            step_sol_2.t, np.concatenate([np.array([0]), np.linspace(dt, 2 * dt, 5)])
+            step_sol_2.t, np.array([0, 1, 1 + 1e-9, 1.25, 1.5, 1.75, 2])
         )
         np.testing.assert_array_almost_equal(
             step_sol_2.y.full()[0], np.exp(0.1 * step_sol_2.t)
@@ -322,16 +322,21 @@ def test_model_step_with_input(self):
 
         # Step again with different inputs
         step_sol_2 = solver.step(step_sol, model, dt, npts=5, inputs={"a": -1})
-        np.testing.assert_array_equal(step_sol_2.t, np.linspace(0, 2 * dt, 9))
+        np.testing.assert_array_almost_equal(
+            step_sol_2.t,
+            np.array([0, 0.025, 0.05, 0.075, 0.1, 0.1 + 1e-9, 0.125, 0.15, 0.175, 0.2]),
+        )
         np.testing.assert_array_equal(
-            step_sol_2["a"].entries, np.array([0.1, 0.1, 0.1, 0.1, 0.1, -1, -1, -1, -1])
+            step_sol_2["a"].entries,
+            np.array([0.1, 0.1, 0.1, 0.1, 0.1, -1, -1, -1, -1, -1]),
         )
         np.testing.assert_allclose(
             step_sol_2.y.full()[0],
             np.concatenate(
                 [
                     np.exp(0.1 * step_sol_2.t[:5]),
-                    np.exp(0.1 * step_sol_2.t[4]) * np.exp(-(step_sol_2.t[5:] - dt)),
+                    np.exp(0.1 * step_sol_2.t[4])
+                    * np.exp(-(step_sol_2.t[5:] - step_sol_2.t[5])),
                 ]
             ),
         )

tests/unit/test_solvers/test_dummy_solver.py

@@ -38,9 +38,11 @@ def test_dummy_solver_step(self):
         for dt in np.diff(t_eval):
             sol = solver.step(sol, model, dt)
 
-        np.testing.assert_array_equal(sol.t, t_eval)
-        np.testing.assert_array_equal(sol.y, np.zeros((1, t_eval.size)))
-        np.testing.assert_array_equal(sol["v"].data, np.ones(t_eval.size))
+        # if t_eval is [0,1,2,3]
+        # sol.t should be [0,1,1+eps,2,2+eps,3] i.e. size 2n-2
+        np.testing.assert_array_equal(len(sol.t), t_eval.size * 2 - 2)
+        np.testing.assert_array_equal(sol.y, np.zeros((1, sol.t.size)))
+        np.testing.assert_array_equal(sol["v"].data, np.ones(sol.t.size))
 
 
 if __name__ == "__main__":

tests/unit/test_solvers/test_scikits_solvers.py

@@ -532,14 +532,14 @@ def test_model_step_ode_python(self):
         # Step again (return 5 points)
         step_sol_2 = solver.step(step_sol, model, dt, npts=5)
         np.testing.assert_array_equal(
-            step_sol_2.t, np.concatenate([np.array([0]), np.linspace(dt, 2 * dt, 5)])
+            step_sol_2.t, np.array([0, 1, 1 + 1e-9, 1.25, 1.5, 1.75, 2])
         )
         np.testing.assert_allclose(step_sol_2.y[0], np.exp(-0.1 * step_sol_2.t))
 
         # Check steps give same solution as solve
         t_eval = step_sol.t
         solution = solver.solve(model, t_eval)
-        np.testing.assert_allclose(solution.y[0], step_sol.y[0])
+        np.testing.assert_allclose(solution.y[0], step_sol.y[0], atol=1e-6, rtol=1e-6)
 
     def test_model_step_dae_python(self):
         model = pybamm.BaseModel()
@@ -560,22 +560,22 @@ def test_model_step_dae_python(self):
         dt = 1
         step_sol = solver.step(None, model, dt)
         np.testing.assert_array_equal(step_sol.t, [0, dt])
-        np.testing.assert_allclose(step_sol.y[0], np.exp(0.1 * step_sol.t))
-        np.testing.assert_allclose(step_sol.y[-1], 2 * np.exp(0.1 * step_sol.t))
+        np.testing.assert_allclose(step_sol.y[0, :], np.exp(0.1 * step_sol.t))
+        np.testing.assert_allclose(step_sol.y[-1, :], 2 * np.exp(0.1 * step_sol.t))
 
         # Step again (return 5 points)
         step_sol_2 = solver.step(step_sol, model, dt, npts=5)
         np.testing.assert_array_equal(
-            step_sol_2.t, np.concatenate([np.array([0]), np.linspace(dt, 2 * dt, 5)])
+            step_sol_2.t, np.array([0, 1, 1 + 1e-9, 1.25, 1.5, 1.75, 2])
         )
-        np.testing.assert_allclose(step_sol_2.y[0], np.exp(0.1 * step_sol_2.t))
-        np.testing.assert_allclose(step_sol_2.y[-1], 2 * np.exp(0.1 * step_sol_2.t))
+        np.testing.assert_allclose(step_sol_2.y[0, :], np.exp(0.1 * step_sol_2.t))
+        np.testing.assert_allclose(step_sol_2.y[-1, :], 2 * np.exp(0.1 * step_sol_2.t))
 
         # Check steps give same solution as solve
         t_eval = step_sol.t
         solution = solver.solve(model, t_eval)
-        np.testing.assert_allclose(solution.y[0], step_sol.y[0, :])
-        np.testing.assert_allclose(solution.y[-1], step_sol.y[-1, :])
+        np.testing.assert_allclose(solution.y[0, :], step_sol.y[0, :])
+        np.testing.assert_allclose(solution.y[-1, :], step_sol.y[-1, :])
 
     def test_model_solver_ode_events_casadi(self):
         # Create model
@@ -598,7 +598,7 @@ def test_model_solver_ode_events_casadi(self):
         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], 1.5)
-        np.testing.assert_array_less(solution.y[0:, -1], 1.25 + 1e-6)
+        np.testing.assert_array_less(solution.y[0:, -1], 1.25 + 1e-9)
         np.testing.assert_equal(solution.t_event[0], solution.t[-1])
         np.testing.assert_array_equal(solution.y_event[:, 0], solution.y[:, -1])
 
@@ -794,8 +794,8 @@ def test_model_step_nonsmooth_events(self):
         )
         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)
-        np.testing.assert_array_almost_equal(step_solution.y[-1], var2_soln, decimal=5)
+        np.testing.assert_array_almost_equal(step_solution.y[0], var1_soln, decimal=4)
+        np.testing.assert_array_almost_equal(step_solution.y[-1], var2_soln, decimal=4)
 
     def test_model_solver_dae_nonsmooth(self):
         whole_cell = ["negative electrode", "separator", "positive electrode"]