Merge pull request pybamm-team#1331 from pybamm-team/issue-1329-add-s…

…olutions

pybamm-team#1329 reformat Solution to avoid concatenating y

CHANGELOG.md

@@ -19,7 +19,8 @@
 
 ## Optimizations
 
--   If solver method `solve()` is passed a list of inputs as the `inputs` keyword argument, the resolution of the model for each input set is spread acrosss several Python processes, usually running in parallel on different processors. The default number of processes is the number of processors available. `solve()` takes a new keyword argument `nproc` which can be used to set this number a manually.
+-   The `Solution` class now only creates the concatenated `y` when the user asks for it. This is an optimization step as the concatenation can be slow, especially with larger experiments ([#1331](https://github.com/pybamm-team/PyBaMM/pull/1331))
+-   If solver method `solve()` is passed a list of inputs as the `inputs` keyword argument, the resolution of the model for each input set is spread across several Python processes, usually running in parallel on different processors. The default number of processes is the number of processors available. `solve()` takes a new keyword argument `nproc` which can be used to set this number a manually.
 -   Variables are now post-processed using CasADi ([#1316](https://github.com/pybamm-team/PyBaMM/pull/1316))
 -   Operations such as `1*x` and `0+x` now directly return `x` ([#1252](https://github.com/pybamm-team/PyBaMM/pull/1252))
 

docs/source/solvers/solution.rst

@@ -1,8 +1,5 @@
 Solutions
 =========
 
-.. autoclass:: pybamm._BaseSolution
-  :members:
-
 .. autoclass:: pybamm.Solution
   :members:

examples/notebooks/Getting Started/Tutorial 9 - Changing the mesh.ipynb

@@ -22,6 +22,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
+      "\u001b[33mWARNING: You are using pip version 20.2.4; however, version 20.3.3 is available.\n",
+      "You should consider upgrading via the '/Users/vsulzer/Documents/Energy_storage/PyBaMM/.tox/dev/bin/python -m pip install --upgrade pip' command.\u001b[0m\n",
       "Note: you may need to restart the kernel to use updated packages.\n"
      ]
     }
@@ -64,13 +66,13 @@
     {
      "data": {
       "text/plain": [
-       "{SpatialVariable(0x62eb5d16885394c9, x_n, children=[], domain=['negative electrode'], auxiliary_domains={'secondary': \"['current collector']\"}): 20,\n",
-       " SpatialVariable(-0x14459c29151bda45, x_s, children=[], domain=['separator'], auxiliary_domains={'secondary': \"['current collector']\"}): 20,\n",
-       " SpatialVariable(-0x4a38fdaaecee173f, x_p, children=[], domain=['positive electrode'], auxiliary_domains={'secondary': \"['current collector']\"}): 20,\n",
-       " SpatialVariable(-0xf387938c410c357, r_n, children=[], domain=['negative particle'], auxiliary_domains={'secondary': \"['negative electrode']\", 'tertiary': \"['current collector']\"}): 30,\n",
-       " SpatialVariable(0x2f93da3e03d4dcbe, r_p, children=[], domain=['positive particle'], auxiliary_domains={'secondary': \"['positive electrode']\", 'tertiary': \"['current collector']\"}): 30,\n",
-       " SpatialVariable(0x14b0200ccfc1f5cc, y, children=[], domain=['current collector'], auxiliary_domains={}): 10,\n",
-       " SpatialVariable(-0x4938b510c47f7708, z, children=[], domain=['current collector'], auxiliary_domains={}): 10}"
+       "{SpatialVariable(0x3452d5d06e5d6beb, x_n, children=[], domain=['negative electrode'], auxiliary_domains={'secondary': \"['current collector']\"}): 20,\n",
+       " SpatialVariable(-0x561e00d7cc14367b, x_s, children=[], domain=['separator'], auxiliary_domains={'secondary': \"['current collector']\"}): 20,\n",
+       " SpatialVariable(-0x71763e2ccac3ade7, x_p, children=[], domain=['positive electrode'], auxiliary_domains={'secondary': \"['current collector']\"}): 20,\n",
+       " SpatialVariable(-0x19c83c2dfa4ac578, r_n, children=[], domain=['negative particle'], auxiliary_domains={'secondary': \"['negative electrode']\", 'tertiary': \"['current collector']\"}): 30,\n",
+       " SpatialVariable(0x5c2de4f3839b67be, r_p, children=[], domain=['positive particle'], auxiliary_domains={'secondary': \"['positive electrode']\", 'tertiary': \"['current collector']\"}): 30,\n",
+       " SpatialVariable(-0x5746d0678b95cd9a, y, children=[], domain=['current collector'], auxiliary_domains={}): 10,\n",
+       " SpatialVariable(-0x7e1ee8ded1c165e8, z, children=[], domain=['current collector'], auxiliary_domains={}): 10}"
       ]
      },
      "execution_count": 3,
@@ -123,7 +125,7 @@
     {
      "data": {
       "text/plain": [
-       "<pybamm.solvers.solution.Solution at 0x7f8beed56e80>"
+       "<pybamm.solvers.solution.Solution at 0x1453f6880>"
       ]
      },
      "execution_count": 5,
@@ -151,7 +153,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "0613bdd52da440c9925d1dcb26558e71",
+       "model_id": "fbd0fd9401444435940b9933716cb07d",
        "version_major": 2,
        "version_minor": 0
       },
@@ -238,7 +240,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "99cad0a3137b4ee1ba353ae02b02b43f",
+       "model_id": "1ef54c1f3c1941f988d9ce0459149078",
        "version_major": 2,
        "version_minor": 0
       },
@@ -252,7 +254,7 @@
     {
      "data": {
       "text/plain": [
-       "<pybamm.plotting.quick_plot.QuickPlot at 0x7f8beb334e80>"
+       "<pybamm.plotting.quick_plot.QuickPlot at 0x146458400>"
       ]
      },
      "execution_count": 9,
@@ -288,7 +290,20 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.8"
+   "version": "3.8.6"
+  },
+  "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/notebooks/models/pouch-cell-model.ipynb

@@ -353,11 +353,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "comsol_solution = pybamm.Solution(solutions[\"1+1D DFN\"].t, solutions[\"1+1D DFN\"].y)\n",
     "comsol_model.timescale = simulations[\"1+1D DFN\"].model.timescale\n",
     "comsol_model.length_scales = simulations[\"1+1D DFN\"].model.length_scales\n",
-    "comsol_solution.model = comsol_model\n",
-    "comsol_solution.inputs = {}"
+    "comsol_solution = pybamm.Solution(solutions[\"1+1D DFN\"].t, solutions[\"1+1D DFN\"].y, comsol_model, {})"
    ]
   },
   {

examples/scripts/compare_comsol/compare_comsol_DFN.py

@@ -123,11 +123,12 @@ def get_interp_fun(variable_name, domain):
 }
 
 # Make new solution with same t and y
-comsol_solution = pybamm.Solution(pybamm_solution.t, pybamm_solution.y)
 # Update solution scales to match the pybamm model
 comsol_model.timescale_eval = pybamm_model.timescale_eval
 comsol_model.length_scales_eval = pybamm_model.length_scales_eval
-comsol_solution.model = comsol_model
+comsol_solution = pybamm.Solution(
+    pybamm_solution.t, pybamm_solution.y, comsol_model, {}
+)
 
 # plot
 output_variables = [

examples/scripts/experimental_protocols/cccv.py

@@ -8,7 +8,7 @@
 experiment = pybamm.Experiment(
     [
         (
-            "Discharge at C/10 for 10 hours or until 3.3 V",
+            "Discharge at C/5 for 10 hours or until 3.3 V",
             "Rest for 1 hour",
             "Charge at 1 A until 4.1 V",
             "Hold at 4.1 V until 50 mA",
@@ -25,16 +25,16 @@
 fig, ax = plt.subplots()
 for i in range(3):
     # Extract sub solutions
-    sol = sim.solution.cycles[i][0]
+    sol = sim.solution.cycles[i]
     # Extract variables
     t = sol["Time [h]"].entries
     V = sol["Terminal voltage [V]"].entries
     # Plot
     ax.plot(t - t[0], V, label="Discharge {}".format(i + 1))
     ax.set_xlabel("Time [h]")
     ax.set_ylabel("Voltage [V]")
-    ax.set_xlim([0, 13])
-ax.legend()
+    ax.set_xlim([0, 10])
+ax.legend(loc="lower left")
 
 # Save time, voltage, current, discharge capacity, temperature, and electrolyte
 # concentration to csv and matlab formats

pybamm/__init__.py

@@ -214,7 +214,7 @@ def version(formatted=False):
 #
 # Solver classes
 #
-from .solvers.solution import Solution, _BaseSolution
+from .solvers.solution import Solution
 from .solvers.processed_variable import ProcessedVariable
 from .solvers.processed_symbolic_variable import ProcessedSymbolicVariable
 from .solvers.base_solver import BaseSolver

pybamm/parameters/parameter_values.py

@@ -565,14 +565,14 @@ def _process_symbol(self, symbol):
                 self.parameter_events.append(
                     pybamm.Event(
                         "Interpolant {} lower bound".format(name),
-                        new_children[0] - min(data[:, 0]),
+                        pybamm.min(new_children[0] - min(data[:, 0])),
                         pybamm.EventType.INTERPOLANT_EXTRAPOLATION,
                     )
                 )
                 self.parameter_events.append(
                     pybamm.Event(
                         "Interpolant {} upper bound".format(name),
-                        max(data[:, 0]) - new_children[0],
+                        pybamm.min(max(data[:, 0]) - new_children[0]),
                         pybamm.EventType.INTERPOLANT_EXTRAPOLATION,
                     )
                 )

pybamm/plotting/quick_plot.py

@@ -428,6 +428,12 @@ def reset_axis(self):
                 self.axis_limits[key] = [x_min, x_max, var_min, var_max]
             else:
                 self.variable_limits[key] = (var_min, var_max)
+            if (
+                var_min is not None
+                and var_max is not None
+                and (np.isnan(var_min) or np.isnan(var_max))
+            ):  # pragma: no cover
+                raise ValueError(f"Axis limits cannot be NaN for variables '{key}'")
 
     def plot(self, t):
         """Produces a quick plot with the internal states at time t.
@@ -676,7 +682,7 @@ def slider_update(self, t):
                         var = variable(
                             time_in_seconds,
                             **self.spatial_variable_dict[key],
-                            warn=False
+                            warn=False,
                         )
                         plot[i][j].set_ydata(var)
                         var_min = min(var_min, np.nanmin(var))

pybamm/simulation.py

@@ -428,7 +428,6 @@ def solve(self, t_eval=None, solver=None, check_model=True, **kwargs):
                         )
 
             self._solution = solver.solve(self.built_model, t_eval, **kwargs)
-            self.t_eval = self._solution.t * self._solution.timescale_eval
 
         elif self.operating_mode == "with experiment":
             if t_eval is not None:
@@ -438,10 +437,13 @@ def solve(self, t_eval=None, solver=None, check_model=True, **kwargs):
             # Re-initialize solution, e.g. for solving multiple times with different
             # inputs without having to build the simulation again
             self._solution = None
+            previous_num_subsolutions = 0
             # Step through all experimental conditions
             inputs = kwargs.get("inputs", {})
             pybamm.logger.info("Start running experiment")
             timer = pybamm.Timer()
+
+            steps = []
             for idx, (exp_inputs, dt) in enumerate(
                 zip(self._experiment_inputs, self._experiment_times)
             ):
@@ -451,6 +453,25 @@ def solve(self, t_eval=None, solver=None, check_model=True, **kwargs):
                 # Make sure we take at least 2 timesteps
                 npts = max(int(round(dt / exp_inputs["period"])) + 1, 2)
                 self.step(dt, solver=solver, npts=npts, **kwargs)
+
+                # Extract the new parts of the solution to construct the entire "step"
+                sol = self.solution
+                new_num_subsolutions = len(sol.sub_solutions)
+                diff_num_subsolutions = new_num_subsolutions - previous_num_subsolutions
+                previous_num_subsolutions = new_num_subsolutions
+
+                step_solution = pybamm.Solution(
+                    sol.all_ts[-diff_num_subsolutions:],
+                    sol.all_ys[-diff_num_subsolutions:],
+                    sol.model,
+                    sol.all_inputs[-diff_num_subsolutions:],
+                    sol.t_event,
+                    sol.y_event,
+                    sol.termination,
+                )
+                step_solution.solve_time = 0
+                step_solution.integration_time = 0
+                steps.append(step_solution)
                 # Only allow events specified by experiment
                 if not (
                     self._solution.termination == "final time"
@@ -467,19 +488,18 @@ def solve(self, t_eval=None, solver=None, check_model=True, **kwargs):
                         "or reducing the period.\n\n"
                     )
                     break
-            if hasattr(self.solution, "_sub_solutions"):
-                # Construct solution.cycles (a list of tuples) from sub_solutions
-                self.solution.cycles = []
-                for cycle_num, cycle_length in enumerate(self.experiment.cycle_lengths):
-                    cycle_start_idx = sum(self.experiment.cycle_lengths[0:cycle_num])
-                    self.solution.cycles.append(
-                        tuple(
-                            [
-                                self.solution.sub_solutions[cycle_start_idx + idx]
-                                for idx in range(cycle_length)
-                            ]
-                        )
-                    )
+            # Construct solution.cycles (a list of solutions corresponding to
+            # cycles) from sub_solutions
+            self.solution.cycles = []
+            for cycle_num, cycle_length in enumerate(self.experiment.cycle_lengths):
+                cycle_start_idx = sum(self.experiment.cycle_lengths[0:cycle_num])
+                cycle_solution = steps[cycle_start_idx]
+                for idx in range(cycle_length - 1):
+                    cycle_solution = cycle_solution + steps[cycle_start_idx + idx + 1]
+                cycle_solution.steps = steps[
+                    cycle_start_idx : cycle_start_idx + cycle_length
+                ]
+                self.solution.cycles.append(cycle_solution)
             pybamm.logger.info(
                 "Finish experiment simulation, took {}".format(timer.time())
             )

pybamm/solvers/algebraic_solver.py

@@ -46,7 +46,7 @@ def tol(self):
     def tol(self, value):
         self._tol = value
 
-    def _integrate(self, model, t_eval, inputs=None):
+    def _integrate(self, model, t_eval, inputs_dict=None):
         """
         Calculate the solution of the algebraic equations through root-finding
 
@@ -56,12 +56,14 @@ def _integrate(self, model, t_eval, inputs=None):
             The model whose solution to calculate.
         t_eval : :class:`numpy.array`, size (k,)
             The times at which to compute the solution
-        inputs : dict, optional
+        inputs_dict : dict, optional
             Any input parameters to pass to the model when solving
         """
-        inputs = inputs or {}
+        inputs_dict = inputs_dict or {}
         if model.convert_to_format == "casadi":
-            inputs = casadi.vertcat(*[x for x in inputs.values()])
+            inputs = casadi.vertcat(*[x for x in inputs_dict.values()])
+        else:
+            inputs = inputs_dict
 
         y0 = model.y0
         if isinstance(y0, casadi.DM):
@@ -218,6 +220,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)
-        sol = pybamm.Solution(t_eval, y_sol, termination="success")
+        sol = pybamm.Solution(t_eval, y_sol, model, inputs_dict, termination="success")
         sol.integration_time = integration_time
         return sol