Merge pull request #840 from pybamm-team/issue-837-car-current

Issue 837 car current

pybamm/simulation.py

@@ -393,7 +393,10 @@ def solve(
             # to correspond to a single discharge
             elif t_eval is None:
                 C_rate = self._parameter_values["C-rate"]
-                t_end = 3600 / C_rate
+                try:
+                    t_end = 3600 / C_rate
+                except TypeError:
+                    t_end = 3600
                 t_eval = np.linspace(0, t_end, 100)
 
             self.t_eval = t_eval

pybamm/solvers/base_solver.py

@@ -219,21 +219,31 @@ def report(string):
 
         # Check for heaviside functions in rhs and algebraic and add discontinuity
         # events if these exist.
-        # Note: only checks for the case of t < X, t <= X, X < t, or X <= t
-        for symbol in itertools.chain(
-            model.concatenated_rhs.pre_order(), model.concatenated_algebraic.pre_order()
-        ):
-            if isinstance(symbol, pybamm.Heaviside):
-                if symbol.right.id == pybamm.t.id:
-                    expr = symbol.left
-                elif symbol.left.id == pybamm.t.id:
-                    expr = symbol.right
-
-                model.events.append(
-                    pybamm.Event(
-                        str(symbol), expr.new_copy(), pybamm.EventType.DISCONTINUITY
+        # Note: only checks for the case of t < X, t <= X, X < t, or X <= t, but also
+        # accounts for the fact that t might be dimensional
+        # Only do this for DAE models as ODE models can deal with discontinuities fine
+        if len(model.algebraic) > 0:
+            for symbol in itertools.chain(
+                model.concatenated_rhs.pre_order(),
+                model.concatenated_algebraic.pre_order(),
+            ):
+                if isinstance(symbol, pybamm.Heaviside):
+                    # Dimensionless
+                    if symbol.right.id == pybamm.t.id:
+                        expr = symbol.left
+                    elif symbol.left.id == pybamm.t.id:
+                        expr = symbol.right
+                    # Dimensional
+                    elif symbol.right.id == (pybamm.t * model.timescale).id:
+                        expr = symbol.left.new_copy() / symbol.right.right.new_copy()
+                    elif symbol.left.id == (pybamm.t * model.timescale).id:
+                        expr = symbol.right.new_copy() / symbol.left.right.new_copy()
+
+                    model.events.append(
+                        pybamm.Event(
+                            str(symbol), expr.new_copy(), pybamm.EventType.DISCONTINUITY
+                        )
                     )
-                )
 
         # Process rhs, algebraic and event expressions
         rhs, rhs_eval, jac_rhs = process(model.concatenated_rhs, "RHS")
@@ -328,24 +338,24 @@ def calculate_consistent_state(self, model, time=0, y0_guess=None):
             Initial conditions that are consistent with the algebraic equations (roots
             of the algebraic equations)
         """
-        pybamm.logger.info("Start calculating consistent initial conditions")
+        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(0, y0_guess).shape[0]
+        len_rhs = rhs(time, y0_guess).shape[0]
         y0_diff, y0_alg_guess = np.split(y0_guess, [len_rhs])
 
         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=0, L2-norm is {}".format(
-                    np.linalg.norm(out)
+                "Evaluating algebraic equations at t={}, L2-norm is {}".format(
+                    time * model.timescale, np.linalg.norm(out)
                 )
             )
             return out
@@ -472,43 +482,46 @@ def solve(self, model, t_eval, external_variables=None, inputs=None):
 
         # make sure they are increasing in time
         discontinuities = sorted(discontinuities)
-        if len(discontinuities) > 0:
-            pybamm.logger.info(
-                "Discontinuity events found at t = {}".format(discontinuities)
-            )
-        else:
-            pybamm.logger.info("No discontinuity events found")
 
         # remove any identical discontinuities
         discontinuities = [
             v
             for i, v in enumerate(discontinuities)
-            if i == len(discontinuities) - 1
-            or discontinuities[i] < discontinuities[i + 1]
+            if (
+                i == len(discontinuities) - 1
+                or discontinuities[i] < discontinuities[i + 1]
+            )
+            and v > 0
         ]
 
+        if len(discontinuities) > 0:
+            pybamm.logger.info(
+                "Discontinuity events found at t = {}".format(discontinuities)
+            )
+        else:
+            pybamm.logger.info("No discontinuity events found")
+
         # insert time points around discontinuities in t_eval
         # keep track of sub sections to integrate by storing start and end indices
         start_indices = [0]
         end_indices = []
+        eps = sys.float_info.epsilon
         for dtime in discontinuities:
             dindex = np.searchsorted(t_eval_dimensionless, dtime, side="left")
             end_indices.append(dindex + 1)
             start_indices.append(dindex + 1)
-            if t_eval_dimensionless[dindex] == dtime:
-                t_eval_dimensionless[dindex] += sys.float_info.epsilon
+            if dtime - eps < t_eval_dimensionless[dindex] < dtime + eps:
+                t_eval_dimensionless[dindex] += eps
                 t_eval_dimensionless = np.insert(
-                    t_eval_dimensionless, dindex, dtime - sys.float_info.epsilon
+                    t_eval_dimensionless, dindex, dtime - eps
                 )
             else:
                 t_eval_dimensionless = np.insert(
-                    t_eval_dimensionless,
-                    dindex,
-                    [dtime - sys.float_info.epsilon, dtime + sys.float_info.epsilon],
+                    t_eval_dimensionless, dindex, [dtime - eps, dtime + eps],
                 )
         end_indices.append(len(t_eval_dimensionless))
 
-        # integrate separatly over each time segment and accumulate into the solution
+        # integrate separately over each time segment and accumulate into the solution
         # object, restarting the solver at each discontinuity (and recalculating a
         # consistent state afterwards if a dae)
         old_y0 = model.y0
@@ -521,31 +534,20 @@ def solve(self, model, t_eval, external_variables=None, inputs=None):
                 )
             )
             timer.reset()
+            new_solution = self._integrate(
+                model, t_eval_dimensionless[start_index:end_index], ext_and_inputs
+            )
+            new_solution.solve_time = timer.time()
             if solution is None:
-                solution = self._integrate(
-                    model, t_eval_dimensionless[start_index:end_index], ext_and_inputs
-                )
-                solution.solve_time = timer.time()
+                solution = new_solution
             else:
-                new_solution = self._integrate(
-                    model, t_eval_dimensionless[start_index:end_index], ext_and_inputs
-                )
-                new_solution.solve_time = timer.time()
                 solution.append(new_solution, start_index=0)
 
             if solution.termination != "final time":
                 break
 
             if end_index != len(t_eval_dimensionless):
                 # setup for next integration subsection
-                y0_guess = solution.y[:, -1]
-                if model.algebraic:
-                    model.y0 = self.calculate_consistent_state(
-                        model, t_eval_dimensionless[end_index], y0_guess
-                    )
-                else:
-                    model.y0 = y0_guess
-
                 last_state = solution.y[:, -1]
                 if len(model.algebraic) > 0:
                     model.y0 = self.calculate_consistent_state(

pybamm/solvers/casadi_solver.py

@@ -100,13 +100,13 @@ def _integrate(self, model, t_eval, inputs=None):
             return solution
         elif self.mode == "safe":
             # Step-and-check
+            t = t_eval[0]
             init_event_signs = np.sign(
                 np.concatenate(
-                    [event(0, model.y0) for event in model.terminate_events_eval]
+                    [event(t, model.y0) for event in model.terminate_events_eval]
                 )
             )
             pybamm.logger.info("Start solving {} with {}".format(model.name, self.name))
-            t = t_eval[0]
             y0 = model.y0
             # Initialize solution
             solution = pybamm.Solution(np.array([t]), y0[:, np.newaxis])
@@ -144,7 +144,7 @@ def _integrate(self, model, t_eval, inputs=None):
                 new_event_signs = np.sign(
                     np.concatenate(
                         [
-                            event(0, current_step_sol.y[:, -1])
+                            event(t, current_step_sol.y[:, -1])
                             for event in model.terminate_events_eval
                         ]
                     )
@@ -186,15 +186,15 @@ def get_integrator(self, model, t_eval, inputs):
             # set up and solve
             t = casadi.MX.sym("t")
             u = casadi.MX.sym("u", u_stacked.shape[0])
-            y_diff = casadi.MX.sym("y_diff", rhs(0, y0, u).shape[0])
+            y_diff = casadi.MX.sym("y_diff", rhs(t_eval[0], y0, u).shape[0])
             problem = {"t": t, "x": y_diff, "p": u}
-            if algebraic(0, y0, u).is_empty():
+            if algebraic(t_eval[0], y0, u).is_empty():
                 method = "cvodes"
                 problem.update({"ode": rhs(t, y_diff, u)})
             else:
                 options["calc_ic"] = True
                 method = "idas"
-                y_alg = casadi.MX.sym("y_alg", algebraic(0, y0, u).shape[0])
+                y_alg = casadi.MX.sym("y_alg", algebraic(t_eval[0], y0, u).shape[0])
                 y_full = casadi.vertcat(y_diff, y_alg)
                 problem.update(
                     {
@@ -215,7 +215,7 @@ def get_integrator(self, model, t_eval, inputs):
         )
 
     def _run_integrator(self, integrator, model, y0, inputs, t_eval):
-        rhs_size = model.rhs_eval(0, y0).shape[0]
+        rhs_size = model.rhs_eval(t_eval[0], y0).shape[0]
         y0_diff, y0_alg = np.split(y0, [rhs_size])
         try:
             # Try solving

tests/unit/test_simulation.py

@@ -385,6 +385,29 @@ def test_drive_cycle_data(self):
         with self.assertWarns(pybamm.SolverWarning):
             sim.solve(t_eval=np.linspace(0, time_data[-1], 800))
 
+    def test_discontinuous_current(self):
+        def car_current(t):
+            current = (
+                1 * (t >= 0) * (t <= 1000)
+                - 0.5 * (1000 < t) * (t <= 2000)
+                + 0.5 * (2000 < t)
+            )
+            return current
+
+        model = pybamm.lithium_ion.DFN()
+        param = model.default_parameter_values
+        param["Current function [A]"] = car_current
+
+        sim = pybamm.Simulation(
+            model, parameter_values=param, solver=pybamm.CasadiSolver(mode="fast")
+        )
+        sim.solve()
+        current = sim.solution["Current [A]"]
+        tau = model.timescale.evaluate()
+        self.assertEqual(current(0), 1)
+        self.assertEqual(current(1500 / tau), -0.5)
+        self.assertEqual(current(3000 / tau), 0.5)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

tests/unit/test_solvers/test_base_solver.py

@@ -50,6 +50,7 @@ def test_find_consistent_initial_conditions(self):
         class ScalarModel:
             concatenated_initial_conditions = np.array([[2]])
             jac_algebraic_eval = None
+            timescale = 1
 
             def rhs_eval(self, t, y):
                 return np.array([])
@@ -68,6 +69,7 @@ def algebraic_eval(self, t, y):
         class VectorModel:
             concatenated_initial_conditions = np.zeros_like(vec)
             jac_algebraic_eval = None
+            timescale = 1
 
             def rhs_eval(self, t, y):
                 return y[0:1]
@@ -101,6 +103,7 @@ def test_fail_consistent_initial_conditions(self):
         class Model:
             concatenated_initial_conditions = np.array([2])
             jac_algebraic_eval = None
+            timescale = 1
 
             def rhs_eval(self, t, y):
                 return np.array([])
@@ -123,6 +126,16 @@ def algebraic_eval(self, t, y):
         ):
             solver.calculate_consistent_state(Model())
 
+    def test_time_too_short(self):
+        solver = pybamm.BaseSolver()
+        model = pybamm.BaseModel()
+        v = pybamm.StateVector(slice(0, 1))
+        model.rhs = {v: v}
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "It looks like t_eval might be dimensionless"
+        ):
+            solver.solve(model, np.linspace(0, 0.1))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")