#1219 start script on speeding up solvers

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added `Solution.integration_time`, which is the time taken just by the integration subroutine, without extra setups ([#1223](https://github.com/pybamm-team/PyBaMM/pull/1223))
 -   Added parameter set for an A123 LFP cell ([#1209](https://github.com/pybamm-team/PyBaMM/pull/1209))
 -   Added variables related to equivalent circuit models ([#1204](https://github.com/pybamm-team/PyBaMM/pull/1204))
 -   Added an example script to check conservation of lithium ([#1186](https://github.com/pybamm-team/PyBaMM/pull/1186))
@@ -17,7 +18,7 @@
 
 ## Bug fixes
 
--   Raise error if saving to matlab with variable names that matlab can't read, and give option of providing alternative variable names ([#1206](https://github.com/pybamm-team/PyBaMM/pull/1206))
+-   Raise error if saving to MATLAB with variable names that MATLAB can't read, and give option of providing alternative variable names ([#1206](https://github.com/pybamm-team/PyBaMM/pull/1206))
 -   Raise error if the boundary condition at the origin in a spherical domain is other than no-flux ([#1175](https://github.com/pybamm-team/PyBaMM/pull/1175))
 -   Fix boundary conditions at r = 0 for Creating Models notebooks ([#1173](https://github.com/pybamm-team/PyBaMM/pull/1173))
 

examples/notebooks/change-input-current.ipynb

@@ -327,7 +327,20 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.9"
+   "version": "3.8.5"
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": true
   }
  },
  "nbformat": 4,

examples/scripts/cycling_ageing_yang.py

@@ -1,51 +1,29 @@
 import pybamm as pb
 
 pb.set_logging_level("INFO")
-options = {"sei": "ec reaction limited", "sei porosity change": True}
+options = {
+    "sei": "ec reaction limited",
+    "sei porosity change": True,
+    "thermal": "x-lumped",
+}
 param = pb.ParameterValues(chemistry=pb.parameter_sets.Ramadass2004)
+param.update(
+    {
+        "Separator density [kg.m-3]": 397,
+        "Separator specific heat capacity [J.kg-1.K-1]": 700,
+        "Separator thermal conductivity [W.m-1.K-1]": 0.16,
+    },
+    check_already_exists=False,
+)
 model = pb.lithium_ion.DFN(options)
 experiment = pb.Experiment(
     [
-        "Charge at 1 C until 4.2 V",
-        "Hold at 4.2 V until C/10",
-        "Rest for 5 minutes",
-        "Discharge at 2 C until 2.8 V",
+        "Charge at 0.3 C until 4.2 V",
         "Rest for 5 minutes",
-    ]
-    * 2
-    + [
-        "Charge at 1 C until 4.2 V",
-        "Hold at 4.2 V until C/20",
-        "Rest for 30 minutes",
-        "Discharge at C/3 until 2.8 V",
-        "Charge at 1 C until 4.2 V",
-        "Hold at 4.2 V until C/20",
-        "Rest for 30 minutes",
         "Discharge at 1 C until 2.8 V",
-        "Charge at 1 C until 4.2 V",
-        "Hold at 4.2 V until C/20",
-        "Rest for 30 minutes",
-        "Discharge at 2 C until 2.8 V",
-        "Charge at 1 C until 4.2 V",
-        "Hold at 4.2 V until C/20",
-        "Rest for 30 minutes",
-        "Discharge at 3 C until 2.8 V",
+        "Rest for 5 minutes",
     ]
+    * 5
 )
 sim = pb.Simulation(model, experiment=experiment, parameter_values=param)
-sim.solve(solver=pb.CasadiSolver(mode="safe", dt_max=120))
-sim.plot(
-    [
-        "Current [A]",
-        "Total current density [A.m-2]",
-        "Terminal voltage [V]",
-        "Discharge capacity [A.h]",
-        "Electrolyte potential [V]",
-        "Electrolyte concentration [mol.m-3]",
-        "Total negative electrode sei thickness",
-        "Negative electrode porosity",
-        "X-averaged negative electrode porosity",
-        "Negative electrode sei interfacial current density [A.m-2]",
-        "X-averaged total negative electrode sei thickness [m]",
-    ]
-)
+sim.solve(solver=pb.CasadiSolver(mode="safe", dt_max=120))

pybamm/expression_tree/binary_operators.py

@@ -138,7 +138,23 @@ def format(self, left, right):
 
     def __str__(self):
         """ See :meth:`pybamm.Symbol.__str__()`. """
-        return "{!s} {} {!s}".format(self.left, self.name, self.right)
+        # Possibly add brackets for clarity
+        if isinstance(self.left, pybamm.BinaryOperator) and not (
+            (self.left.name == self.name)
+            or (self.left.name == "*" and self.name == "/")
+        ):
+            left_str = "({!s})".format(self.left)
+        else:
+            left_str = "{!s}".format(self.left)
+        if isinstance(self.right, pybamm.BinaryOperator) and not (
+            (self.name == "*" and self.right.name == "*")
+            or (self.name == "+" and self.right.name == "+")
+            or (self.name == "*" and self.right.name == "/")
+        ):
+            right_str = "({!s})".format(self.right)
+        else:
+            right_str = "{!s}".format(self.right)
+        return "{} {} {}".format(left_str, self.name, right_str)
 
     def get_children_domains(self, ldomain, rdomain):
         "Combine domains from children in appropriate way"

pybamm/models/submodels/interface/sei/ec_reaction_limited.py

@@ -93,7 +93,8 @@ def set_algebraic(self, variables):
                 + self.domain.lower()
                 + " electrode interfacial current density"
             ]
-        except KeyError:
+        except KeyError as e:
+            print(e)
             j = variables[
                 "X-averaged "
                 + self.domain.lower()

pybamm/solvers/algebraic_solver.py

@@ -82,6 +82,8 @@ def _integrate(self, model, t_eval, inputs=None):
 
         y_alg = np.empty((len(y0_alg), len(t_eval)))
 
+        timer = pybamm.Timer()
+        integration_time = 0
         for idx, t in enumerate(t_eval):
 
             def root_fun(y_alg):
@@ -135,6 +137,7 @@ def jac_fn(y_alg):
                         method = self.method[5:]
                     if jac_fn is None:
                         jac_fn = "2-point"
+                    timer.reset()
                     sol = optimize.least_squares(
                         root_fun,
                         y0_alg,
@@ -144,6 +147,7 @@ def jac_fn(y_alg):
                         bounds=model.bounds,
                         **self.extra_options,
                     )
+                    integration_time += timer.time()
                 # Methods which use minimize are specified as either "minimize", which
                 # uses the default method, or with "minimize__methodname"
                 elif self.method.startswith("minimize"):
@@ -170,6 +174,7 @@ def jac_norm(y):
                             (lb, ub) for lb, ub in zip(model.bounds[0], model.bounds[1])
                         ]
                         extra_options["bounds"] = bounds
+                    timer.reset()
                     sol = optimize.minimize(
                         root_norm,
                         y0_alg,
@@ -178,7 +183,9 @@ def jac_norm(y):
                         jac=jac_norm,
                         **extra_options,
                     )
+                    integration_time += timer.time()
                 else:
+                    timer.reset()
                     sol = optimize.root(
                         root_fun,
                         y0_alg,
@@ -187,6 +194,7 @@ def jac_norm(y):
                         jac=jac_fn,
                         options=self.extra_options,
                     )
+                    integration_time += timer.time()
 
                 if sol.success and np.all(abs(sol.fun) < self.tol):
                     # update initial guess for the next iteration
@@ -210,4 +218,6 @@ def jac_norm(y):
         y_diff = np.r_[[y0_diff] * len(t_eval)].T
         y_sol = np.r_[y_diff, y_alg]
         # Return solution object (no events, so pass None to t_event, y_event)
-        return pybamm.Solution(t_eval, y_sol, termination="success")
+        sol = pybamm.Solution(t_eval, y_sol, termination="success")
+        sol.integration_time = integration_time
+        return sol

pybamm/solvers/base_solver.py

@@ -578,6 +578,7 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
                     "initial conditions"
                 ] = model.concatenated_initial_conditions
         set_up_time = timer.time()
+        timer.reset()
 
         # (Re-)calculate consistent initial conditions
         self._set_initial_conditions(model, ext_and_inputs, update_rhs=True)
@@ -647,7 +648,6 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
                     t_eval_dimensionless[end_index - 1] * model.timescale_eval,
                 )
             )
-            timer.reset()
             new_solution = self._integrate(
                 model, t_eval_dimensionless[start_index:end_index], ext_and_inputs
             )
@@ -671,13 +671,13 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
                         model, t_eval_dimensionless[end_index], ext_and_inputs
                     )
 
-        # restore old y0
-        model.y0 = old_y0
-
         # Assign times
         solution.set_up_time = set_up_time
         solution.solve_time = timer.time()
 
+        # restore old y0
+        model.y0 = old_y0
+
         # Add model and inputs to solution
         solution.model = model
         solution.inputs = ext_and_inputs

pybamm/solvers/casadi_algebraic_solver.py

@@ -118,6 +118,8 @@ def _integrate(self, model, t_eval, inputs=None):
                 "constraints": list(constraints[len_rhs:]),
             },
         )
+        timer = pybamm.Timer()
+        integration_time = 0
         for idx, t in enumerate(t_eval):
             # Evaluate algebraic with new t and previous y0, if it's already close
             # enough then keep it
@@ -137,7 +139,9 @@ def _integrate(self, model, t_eval, inputs=None):
                 t_eval_inputs_sym = casadi.vertcat(t, symbolic_inputs)
                 # Solve
                 try:
+                    timer.reset()
                     y_alg_sol = roots(y0_alg, t_eval_inputs_sym)
+                    integration_time += timer.time()
                     success = True
                     message = None
                     # Check final output
@@ -179,4 +183,6 @@ def _integrate(self, model, t_eval, inputs=None):
         y_diff = casadi.horzcat(*[y0_diff] * len(t_eval))
         y_sol = casadi.vertcat(y_diff, y_alg)
         # Return solution object (no events, so pass None to t_event, y_event)
-        return pybamm.Solution(t_eval, y_sol, termination="success")
+        sol = pybamm.Solution(t_eval, y_sol, termination="success")
+        sol.integration_time = integration_time
+        return sol

pybamm/solvers/casadi_solver.py

@@ -396,11 +396,15 @@ def _run_integrator(self, model, y0, inputs, t_eval):
             # Try solving
             if use_grid is True:
                 # Call the integrator once, with the grid
+                timer = pybamm.Timer()
                 sol = integrator(
                     x0=y0_diff, z0=y0_alg, p=inputs, **self.extra_options_call
                 )
+                integration_time = timer.time()
                 y_sol = np.concatenate([sol["xf"].full(), sol["zf"].full()])
-                return pybamm.Solution(t_eval, y_sol)
+                sol = pybamm.Solution(t_eval, y_sol)
+                sol.integration_time = integration_time
+                return sol
             else:
                 # Repeated calls to the integrator
                 x = y0_diff
@@ -411,19 +415,24 @@ def _run_integrator(self, model, y0, inputs, t_eval):
                     t_min = t_eval[i]
                     t_max = t_eval[i + 1]
                     inputs_with_tlims = casadi.vertcat(inputs, t_min, t_max)
+                    timer = pybamm.Timer()
                     sol = integrator(
                         x0=x, z0=z, p=inputs_with_tlims, **self.extra_options_call
                     )
+                    integration_time = timer.time()
                     x = sol["xf"]
                     z = sol["zf"]
                     y_diff = casadi.horzcat(y_diff, x)
                     if not z.is_empty():
                         y_alg = casadi.horzcat(y_alg, z)
                 if z.is_empty():
-                    return pybamm.Solution(t_eval, y_diff)
+                    sol = pybamm.Solution(t_eval, y_diff)
                 else:
                     y_sol = casadi.vertcat(y_diff, y_alg)
-                    return pybamm.Solution(t_eval, y_sol)
+                    sol = pybamm.Solution(t_eval, y_sol)
+
+                sol.integration_time = integration_time
+                return sol
         except RuntimeError as e:
             # If it doesn't work raise error
             raise pybamm.SolverError(e.args[0])

pybamm/solvers/dummy_solver.py

@@ -33,4 +33,6 @@ def _integrate(self, model, t_eval, inputs=None):
 
         """
         y_sol = np.zeros((1, t_eval.size))
-        return pybamm.Solution(t_eval, y_sol, termination="final time")
+        sol = pybamm.Solution(t_eval, y_sol, termination="final time")
+        sol.integration_time = 0
+        return sol

pybamm/solvers/idaklu_solver.py

@@ -234,6 +234,7 @@ def rootfn(t, y):
         ids = np.concatenate((rhs_ids, alg_ids))
 
         # solve
+        timer = pybamm.Timer()
         sol = idaklu.solve(
             t_eval,
             y0,
@@ -251,6 +252,7 @@ def rootfn(t, y):
             atol,
             rtol,
         )
+        integration_time = timer.time()
 
         t = sol.t
         number_of_timesteps = t.size
@@ -266,12 +268,14 @@ def rootfn(t, y):
             elif sol.flag == 2:
                 termination = "event"
 
-            return pybamm.Solution(
+            sol = pybamm.Solution(
                 sol.t,
                 np.transpose(y_out),
                 t[-1],
                 np.transpose(y_out[-1])[:, np.newaxis],
                 termination,
             )
+            sol.integration_time = integration_time
+            return sol
         else:
             raise pybamm.SolverError(sol.message)

pybamm/solvers/jax_solver.py

@@ -188,15 +188,19 @@ def _integrate(self, model, t_eval, inputs=None):
             various diagnostic messages.
 
         """
+        timer = pybamm.Timer()
         if model not in self._cached_solves:
             self._cached_solves[model] = self.create_solve(model, t_eval)
 
         y = self._cached_solves[model](inputs)
+        integration_time = timer.time()
 
         # note - the actual solve is not done until this line!
         y = onp.array(y)
 
         termination = "final time"
         t_event = None
         y_event = onp.array(None)
-        return pybamm.Solution(t_eval, y, t_event, y_event, termination)
+        sol = pybamm.Solution(t_eval, y, t_event, y_event, termination)
+        sol.integration_time = integration_time
+        return sol

pybamm/solvers/scikits_dae_solver.py

@@ -130,7 +130,9 @@ def jacfn(t, y, ydot, residuals, cj, J):
 
         # set up and solve
         dae_solver = scikits_odes.dae(self.method, eqsres, **extra_options)
+        timer = pybamm.Timer()
         sol = dae_solver.solve(t_eval, y0, ydot0)
+        integration_time = timer.time()
 
         # return solution, we need to tranpose y to match scipy's interface
         if sol.flag in [0, 2]:
@@ -144,12 +146,14 @@ def jacfn(t, y, ydot, residuals, cj, J):
                 t_root = None
             else:
                 t_root = sol.roots.t
-            return pybamm.Solution(
+            sol = pybamm.Solution(
                 sol.values.t,
                 np.transpose(sol.values.y),
                 t_root,
                 np.transpose(sol.roots.y),
                 termination,
             )
+            sol.integration_time = integration_time
+            return sol
         else:
             raise pybamm.SolverError(sol.message)