Issue 684 casadi events

examples/scripts/compare_lithium_ion.py

@@ -45,7 +45,7 @@
 
 # solve model
 solutions = [None] * len(models)
-t_eval = np.linspace(0, 0.17, 100)
+t_eval = np.linspace(0, 0.3, 100)
 for i, model in enumerate(models):
     solutions[i] = model.default_solver.solve(model, t_eval)
 

pybamm/solvers/algebraic_solver.py

@@ -24,6 +24,7 @@ class AlgebraicSolver(object):
     def __init__(self, method="lm", tol=1e-6):
         self.method = method
         self.tol = tol
+        self.name = "Algebraic solver ({})".format(method)
 
     @property
     def method(self):
@@ -51,7 +52,7 @@ def solve(self, model):
             equations.
 
         """
-        pybamm.logger.info("Start solving {}".format(model.name))
+        pybamm.logger.info("Start solving {} with {}".format(model.name, self.name))
 
         # Set up
         timer = pybamm.Timer()
@@ -87,7 +88,6 @@ def jacobian(y):
 
         # Assign times
         solution.solve_time = timer.time() - solve_start_time
-        solution.total_time = timer.time() - start_time
         solution.set_up_time = set_up_time
 
         pybamm.logger.info("Finish solving {}".format(model.name))

pybamm/solvers/base_solver.py

@@ -20,6 +20,7 @@ def __init__(self, method=None, rtol=1e-6, atol=1e-6):
         self._method = method
         self._rtol = rtol
         self._atol = atol
+        self.name = "Base solver"
 
     @property
     def method(self):
@@ -64,7 +65,7 @@ def solve(self, model, t_eval):
             If an empty model is passed (`model.rhs = {}` and `model.algebraic={}`)
 
         """
-        pybamm.logger.info("Start solving {}".format(model.name))
+        pybamm.logger.info("Start solving {} with {}".format(model.name, self.name))
 
         # Make sure model isn't empty
         if len(model.rhs) == 0 and len(model.algebraic) == 0:
@@ -84,7 +85,6 @@ def solve(self, model, t_eval):
 
         # Assign times
         solution.solve_time = solve_time
-        solution.total_time = timer.time() - start_time
         solution.set_up_time = set_up_time
 
         pybamm.logger.info("Finish solving {} ({})".format(model.name, termination))
@@ -97,7 +97,7 @@ def solve(self, model, t_eval):
         )
         return solution
 
-    def step(self, model, dt, npts=2):
+    def step(self, model, dt, npts=2, log=True):
         """
         Step the solution of the model forward by a given time increment. The
         first time this method is called it executes the necessary setup by
@@ -129,44 +129,48 @@ def step(self, model, dt, npts=2):
 
         # Run set up on first step
         if not hasattr(self, "y0"):
-            start_time = timer.time()
+            pybamm.logger.info(
+                "Start stepping {} with {}".format(model.name, self.name)
+            )
+
             if model.convert_to_format == "casadi" or isinstance(
                 self, pybamm.CasadiSolver
             ):
                 self.set_up_casadi(model)
             else:
+                pybamm.logger.debug(
+                    "Start stepping {} with {}".format(model.name, self.name)
+                )
                 self.set_up(model)
             self.t = 0.0
-            set_up_time = timer.time() - start_time
+            set_up_time = timer.time()
         else:
             set_up_time = None
 
         # Step
-        pybamm.logger.info("Start stepping {}".format(model.name))
         t_eval = np.linspace(self.t, self.t + dt, npts)
         solution, solve_time, termination = self.compute_solution(model, t_eval)
 
         # Assign times
         solution.solve_time = solve_time
         if set_up_time:
-            solution.total_time = timer.time() - start_time
             solution.set_up_time = set_up_time
 
         # Set self.t and self.y0 to their values at the final step
         self.t = solution.t[-1]
         self.y0 = solution.y[:, -1]
 
-        pybamm.logger.info("Finish stepping {} ({})".format(model.name, termination))
+        pybamm.logger.debug("Finish stepping {} ({})".format(model.name, termination))
         if set_up_time:
-            pybamm.logger.info(
+            pybamm.logger.debug(
                 "Set-up time: {}, Step time: {}, Total time: {}".format(
                     timer.format(solution.set_up_time),
                     timer.format(solution.solve_time),
                     timer.format(solution.total_time),
                 )
             )
         else:
-            pybamm.logger.info(
+            pybamm.logger.debug(
                 "Step time: {}".format(timer.format(solution.solve_time))
             )
         return solution

pybamm/solvers/casadi_solver.py

@@ -9,14 +9,29 @@
 class CasadiSolver(pybamm.DaeSolver):
     """Solve a discretised model, using CasADi.
 
+    **Extends**: :class:`pybamm.DaeSolver`
+
     Parameters
     ----------
+    method : str, optional
+        The method to use for solving the system ('cvodes', for ODEs, or 'idas', for
+        DAEs). Default is 'idas'.
     rtol : float, optional
         The relative tolerance for the solver (default is 1e-6).
     atol : float, optional
         The absolute tolerance for the solver (default is 1e-6).
+    root_method : str, optional
+        The method to use for finding consistend initial conditions. Default is 'lm'.
+    root_tol : float, optional
+        The tolerance for root-finding. Default is 1e-6.
+    max_step_decrease_counts : float, optional
+        The maximum number of times step size can be decreased before an error is
+        raised. Default is 10.
+    extra_options : keyword arguments, optional
+        Any extra keyword-arguments; these are passed directly to the CasADi integrator.
+        Please consult `CasADi documentation <https://tinyurl.com/y5rk76os>`_ for
+        details.
 
-    **Extends**: :class:`pybamm.DaeSolver`
     """
 
     def __init__(
@@ -26,11 +41,138 @@ def __init__(
         atol=1e-6,
         root_method="lm",
         root_tol=1e-6,
-        max_steps=1000,
+        max_step_decrease_count=10,
         **extra_options,
     ):
-        super().__init__(method, rtol, atol, root_method, root_tol, max_steps)
+        super().__init__(method, rtol, atol, root_method, root_tol)
+        self.max_step_decrease_count = max_step_decrease_count
         self.extra_options = extra_options
+        self.name = "CasADi solver ({})".format(method)
+
+    def solve(self, model, t_eval, mode="safe"):
+        """
+        Execute the solver setup and calculate the solution of the model at
+        specified times.
+
+        Parameters
+        ----------
+        model : :class:`pybamm.BaseModel`
+            The model whose solution to calculate. Must have attributes rhs and
+            initial_conditions
+        t_eval : numeric type
+            The times at which to compute the solution
+        mode : str
+            How to solve the model (default is "safe"):
+
+            - "fast": perform direct integration, without accounting for events. \
+            Recommended when simulating a drive cycle or other simulation where \
+            no events should be triggered.
+            - "safe": perform step-and-check integration, checking whether events have \
+            been triggered. Recommended for simulations of a full charge or discharge.
+
+        Raises
+        ------
+        :class:`pybamm.ValueError`
+            If an invalid mode is passed.
+        :class:`pybamm.ModelError`
+            If an empty model is passed (`model.rhs = {}` and `model.algebraic={}`)
+
+        """
+        if mode == "fast":
+            # Solve model normally by calling the solve method from parent class
+            return super().solve(model, t_eval)
+        elif model.events == {}:
+            pybamm.logger.info("No events found, running fast mode")
+            # Solve model normally by calling the solve method from parent class
+            return super().solve(model, t_eval)
+        elif mode == "safe":
+            # Step-and-check
+            timer = pybamm.Timer()
+            self.set_up_casadi(model)
+            set_up_time = timer.time()
+            init_event_signs = np.sign(
+                np.concatenate([event(0, self.y0) for event in self.event_funs])
+            )
+            solution = None
+            pybamm.logger.info(
+                "Start solving {} with {} in 'safe' mode".format(model.name, self.name)
+            )
+            self.t = 0.0
+            for dt in np.diff(t_eval):
+                # Step
+                solved = False
+                count = 0
+                while not solved:
+                    # Try to solve with the current step, if it fails then halve the
+                    # step size and try again. This will make solution.t slightly
+                    # different to t_eval, but shouldn't matter too much as it should
+                    # only happen near events.
+                    try:
+                        current_step_sol = self.step(model, dt)
+                        solved = True
+                    except pybamm.SolverError:
+                        dt /= 2
+                    count += 1
+                    if count >= self.max_step_decrease_count:
+                        if solution is None:
+                            t = 0
+                        else:
+                            t = solution.t
+                        raise pybamm.SolverError(
+                            """
+                            Maximum number of decreased steps occurred at t={}. Try
+                            solving the model up to this time only
+                            """.format(
+                                t
+                            )
+                        )
+                # Check most recent y
+                new_event_signs = np.sign(
+                    np.concatenate(
+                        [
+                            event(0, current_step_sol.y[:, -1])
+                            for event in self.event_funs
+                        ]
+                    )
+                )
+                # Exit loop if the sign of an event changes
+                if (new_event_signs != init_event_signs).any():
+                    solution.termination = "event"
+                    solution.t_event = solution.t[-1]
+                    solution.y_event = solution.y[:, -1]
+                    break
+                else:
+                    if not solution:
+                        # create solution object on first step
+                        solution = current_step_sol
+                        solution.set_up_time = set_up_time
+                    else:
+                        # append solution from the current step to solution
+                        solution.append(current_step_sol)
+            if not hasattr(solution, "termination"):
+                solution.termination = "final time"
+
+            # Calculate more exact termination reason
+            solution.termination = self.get_termination_reason(solution, self.events)
+            pybamm.logger.info(
+                "Finish solving {} ({})".format(model.name, solution.termination)
+            )
+            pybamm.logger.info(
+                "Set-up time: {}, Solve time: {}, Total time: {}".format(
+                    timer.format(solution.set_up_time),
+                    timer.format(solution.solve_time),
+                    timer.format(solution.total_time),
+                )
+            )
+            return solution
+        else:
+            raise ValueError(
+                """
+                invalid mode '{}'. Must be either 'safe', for solving with events,
+                or 'fast', for solving quickly without events""".format(
+                    mode
+                )
+            )
 
     def compute_solution(self, model, t_eval):
         """Calculate the solution of the model at specified times. In this class, we
@@ -49,14 +191,14 @@ def compute_solution(self, model, t_eval):
         timer = pybamm.Timer()
 
         solve_start_time = timer.time()
-        pybamm.logger.info("Calling DAE solver")
+        pybamm.logger.debug("Calling DAE solver")
         solution = self.integrate_casadi(
             self.casadi_problem, self.y0, t_eval, mass_matrix=model.mass_matrix.entries
         )
         solve_time = timer.time() - solve_start_time
 
-        # Identify the event that caused termination
-        termination = self.get_termination_reason(solution, self.events)
+        # Events not implemented, termination is always 'final time'
+        termination = "final time"
 
         return solution, solve_time, termination
 

pybamm/solvers/dae_solver.py

@@ -39,6 +39,7 @@ def __init__(
         self.root_method = root_method
         self.root_tol = root_tol
         self.max_steps = max_steps
+        self.name = "Base DAE solver"
 
     @property
     def root_method(self):

pybamm/solvers/idaklu_solver.py

@@ -43,6 +43,7 @@ def __init__(
             raise ImportError("KLU is not installed")
 
         super().__init__("ida", rtol, atol, root_method, root_tol, max_steps)
+        self.name = "IDA KLU solver"
 
     def integrate(self, residuals, y0, t_eval, events, mass_matrix, jacobian):
         """

pybamm/solvers/ode_solver.py

@@ -19,6 +19,7 @@ class OdeSolver(pybamm.BaseSolver):
 
     def __init__(self, method=None, rtol=1e-6, atol=1e-6):
         super().__init__(method, rtol, atol)
+        self.name = "Base ODE solver"
 
     def compute_solution(self, model, t_eval):
         """Calculate the solution of the model at specified times.

pybamm/solvers/scikits_dae_solver.py

@@ -48,6 +48,7 @@ def __init__(
             raise ImportError("scikits.odes is not installed")
 
         super().__init__(method, rtol, atol, root_method, root_tol, max_steps)
+        self.name = "Scikits DAE solver ({})".format(method)
 
     def integrate(
         self, residuals, y0, t_eval, events=None, mass_matrix=None, jacobian=None

pybamm/solvers/scikits_ode_solver.py

@@ -40,6 +40,7 @@ def __init__(self, method="cvode", rtol=1e-6, atol=1e-6, linsolver="dense"):
 
         super().__init__(method, rtol, atol)
         self.linsolver = linsolver
+        self.name = "Scikits ODE solver ({})".format(method)
 
     def integrate(
         self, derivs, y0, t_eval, events=None, mass_matrix=None, jacobian=None

pybamm/solvers/scipy_solver.py

@@ -22,6 +22,7 @@ class ScipySolver(pybamm.OdeSolver):
 
     def __init__(self, method="BDF", rtol=1e-6, atol=1e-6):
         super().__init__(method, rtol, atol)
+        self.name = "Scipy solver ({})".format(method)
 
     def integrate(
         self, derivs, y0, t_eval, events=None, mass_matrix=None, jacobian=None

pybamm/solvers/solution.py

@@ -93,3 +93,8 @@ def append(self, solution):
         """
         self.t = np.concatenate((self.t, solution.t[1:]))
         self.y = np.concatenate((self.y, solution.y[:, 1:]), axis=1)
+        self.solve_time += solution.solve_time
+
+    @property
+    def total_time(self):
+        return self.set_up_time + self.solve_time

tests/unit/test_solvers/test_algebraic_solver.py

@@ -107,11 +107,6 @@ def test_model_solver(self):
             model.variables["var2"].evaluate(t=None, y=solution.y), sol[100:]
         )
 
-        # Test time
-        self.assertGreater(
-            solution.total_time, solution.solve_time + solution.set_up_time
-        )
-
         # Test without jacobian
         model.use_jacobian = False
         solution_no_jac = solver.solve(model)