Merge branch 'develop' into composite-porosity

.github/workflows/scorecard.yml

@@ -68,6 +68,6 @@ jobs:
       # Upload the results to GitHub's code scanning dashboard (optional).
       # Commenting out will disable upload of results to your repo's Code Scanning dashboard
       - name: "Upload to code-scanning"
-        uses: github/codeql-action/upload-sarif@f779452ac5af1c261dce0346a8f964149f49322b # v3.26.13
+        uses: github/codeql-action/upload-sarif@662472033e021d55d94146f66f6058822b0b39fd # v3.27.0
         with:
           sarif_file: results.sarif

.pre-commit-config.yaml

@@ -4,7 +4,7 @@ ci:
 
 repos:
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: "v0.7.0"
+    rev: "v0.7.1"
     hooks:
       - id: ruff
         args: [--fix, --show-fixes]
@@ -13,7 +13,7 @@ repos:
         types_or: [python, pyi, jupyter]
 
   - repo: https://github.com/adamchainz/blacken-docs
-    rev: "1.19.0"
+    rev: "1.19.1"
     hooks:
        - id: blacken-docs
          additional_dependencies: [black==23.*]

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- Adds support to `pybamm.Experiment` for the `output_variables` option in the `IDAKLUSolver`. ([#4534](https://github.com/pybamm-team/PyBaMM/pull/4534))
 - Adds an option "voltage as a state" that can be "false" (default) or "true". If "true" adds an explicit algebraic equation for the voltage. ([#4507](https://github.com/pybamm-team/PyBaMM/pull/4507))
 - Improved `QuickPlot` accuracy for simulations with Hermite interpolation. ([#4483](https://github.com/pybamm-team/PyBaMM/pull/4483))
 - Added Hermite interpolation to the (`IDAKLUSolver`) that improves the accuracy and performance of post-processing variables. ([#4464](https://github.com/pybamm-team/PyBaMM/pull/4464))
@@ -21,7 +22,7 @@
 - Removed the `start_step_offset` setting and disabled minimum `dt` warnings for drive cycles with the (`IDAKLUSolver`). ([#4416](https://github.com/pybamm-team/PyBaMM/pull/4416))
 
 ## Bug Fixes
-
+- Fixed bug in post-processing solutions with infeasible experiments using the (`IDAKLUSolver`). ([#4541](https://github.com/pybamm-team/PyBaMM/pull/4541))
 - Disabled IREE on MacOS due to compatibility issues and added the CasADI
   path to the environment to resolve issues on MacOS and Linux. Windows
   users may still experience issues with interpolation. ([#4528](https://github.com/pybamm-team/PyBaMM/pull/4528))

src/pybamm/solvers/base_solver.py

@@ -1452,6 +1452,7 @@ def get_termination_reason(solution, events):
                 solution.t_event,
                 solution.y_event,
                 solution.termination,
+                variables_returned=solution.variables_returned,
             )
             event_sol.solve_time = 0
             event_sol.integration_time = 0

src/pybamm/solvers/c_solvers/idaklu/observe.cpp

@@ -100,6 +100,11 @@ class TimeSeriesInterpolator {
         ) {
         for (size_t i = 0; i < ts_data_np.size(); i++) {
             const auto& t_data = ts_data_np[i].unchecked<1>();
+            // Continue if there is no data
+            if (t_data.size() == 0) {
+                continue;
+            }
+
             const realtype t_data_final = t_data(t_data.size() - 1);
             realtype t_interp_next = t_interp(i_interp);
             // Continue if the next interpolation point is beyond the final data point
@@ -227,6 +232,10 @@ class TimeSeriesProcessor {
         int i_entries = 0;
         for (size_t i = 0; i < ts.size(); i++) {
             const auto& t = ts[i].unchecked<1>();
+            // Continue if there is no data
+            if (t.size() == 0) {
+                continue;
+            }
             const auto& y = ys[i].unchecked<2>();
             const auto input = inputs[i].data();
             const auto func = *funcs[i];

src/pybamm/solvers/idaklu_solver.py

@@ -863,6 +863,7 @@ def _post_process_solution(self, sol, model, integration_time, inputs_dict):
             termination,
             all_sensitivities=yS_out,
             all_yps=yp,
+            variables_returned=bool(save_outputs_only),
         )
 
         newsol.integration_time = integration_time

src/pybamm/solvers/processed_variable.py

@@ -133,19 +133,22 @@ def _setup_cpp_inputs(self, t, full_range):
         ys = self.all_ys
         yps = self.all_yps
         inputs = self.all_inputs_casadi
-        # Find the indices of the time points to observe
-        if full_range:
-            idxs = range(len(ts))
-        else:
-            idxs = _find_ts_indices(ts, t)
 
-        if isinstance(idxs, list):
-            # Extract the time points and inputs
-            ts = [ts[idx] for idx in idxs]
-            ys = [ys[idx] for idx in idxs]
-            if self.hermite_interpolation:
-                yps = [yps[idx] for idx in idxs]
-            inputs = [self.all_inputs_casadi[idx] for idx in idxs]
+        # Remove all empty ts
+        idxs = np.where([ti.size > 0 for ti in ts])[0]
+
+        # Find the indices of the time points to observe
+        if not full_range:
+            ts_nonempty = [ts[idx] for idx in idxs]
+            idxs_subset = _find_ts_indices(ts_nonempty, t)
+            idxs = idxs[idxs_subset]
+
+        # Extract the time points and inputs
+        ts = [ts[idx] for idx in idxs]
+        ys = [ys[idx] for idx in idxs]
+        if self.hermite_interpolation:
+            yps = [yps[idx] for idx in idxs]
+        inputs = [self.all_inputs_casadi[idx] for idx in idxs]
 
         is_f_contiguous = _is_f_contiguous(ys)
 
@@ -977,8 +980,4 @@ def _find_ts_indices(ts, t):
     if (t[-1] > ts[-1][-1]) and (len(indices) == 0 or indices[-1] != len(ts) - 1):
         indices.append(len(ts) - 1)
 
-    if len(indices) == len(ts):
-        # All indices are included
-        return range(len(ts))
-
     return indices

src/pybamm/solvers/processed_variable_computed.py

@@ -1,6 +1,7 @@
 #
-# Processed Variable class
+# Processed Variable Computed class
 #
+from __future__ import annotations
 import casadi
 import numpy as np
 import pybamm
@@ -450,3 +451,27 @@ def sensitivities(self):
         if len(self.all_inputs[0]) == 0:
             return {}
         return self._sensitivities
+
+    def _update(
+        self, other: pybamm.ProcessedVariableComputed, new_sol: pybamm.Solution
+    ) -> pybamm.ProcessedVariableComputed:
+        """
+        Returns a new ProcessedVariableComputed object that is the result of appending
+        the data from other to this object. Used exclusively in running experiments, to
+        append the data from one cycle to the next.
+
+        Parameters
+        ----------
+        other : :class:`pybamm.ProcessedVariableComputed`
+            The other ProcessedVariableComputed object to append to this one
+        new_sol : :class:`pybamm.Solution`
+            The new solution object to be used to create the processed variables
+        """
+
+        bv = self.base_variables + other.base_variables
+        bvc = self.base_variables_casadi + other.base_variables_casadi
+        bvd = self.base_variables_data + other.base_variables_data
+
+        new_var = self.__class__(bv, bvc, bvd, new_sol)
+
+        return new_var

src/pybamm/solvers/solution.py

@@ -62,6 +62,10 @@ class Solution:
         True if sensitivities included as the solution of the explicit forwards
         equations.  False if no sensitivities included/wanted. Dict if sensitivities are
         provided as a dict of {parameter: [sensitivities]} pairs.
+    variables_returned: bool
+        Bool to indicate if `all_ys` contains the full state vector, or is empty because
+        only requested variables have been returned. True if `output_variables` is used
+        with a solver, otherwise False.
 
     """
 
@@ -76,6 +80,7 @@ def __init__(
         termination="final time",
         all_sensitivities=False,
         all_yps=None,
+        variables_returned=False,
         check_solution=True,
     ):
         if not isinstance(all_ts, list):
@@ -93,6 +98,8 @@ def __init__(
             all_yps = [all_yps]
         self._all_yps = all_yps
 
+        self.variables_returned = variables_returned
+
         # Set up inputs
         if not isinstance(all_inputs, list):
             all_inputs_copy = dict(all_inputs)
@@ -460,9 +467,15 @@ def first_state(self):
         else:
             all_yps = self.all_yps[0][:, :1]
 
+        if not self.variables_returned:
+            all_ys = self.all_ys[0][:, :1]
+        else:
+            # Get first state from initial conditions as all_ys is empty
+            all_ys = self.all_models[0].y0full.reshape(-1, 1)
+
         new_sol = Solution(
             self.all_ts[0][:1],
-            self.all_ys[0][:, :1],
+            all_ys,
             self.all_models[:1],
             self.all_inputs[:1],
             None,
@@ -500,9 +513,15 @@ def last_state(self):
         else:
             all_yps = self.all_yps[-1][:, -1:]
 
+        if not self.variables_returned:
+            all_ys = self.all_ys[-1][:, -1:]
+        else:
+            # Get last state from y_event as all_ys is empty
+            all_ys = self.y_event.reshape(len(self.y_event), 1)
+
         new_sol = Solution(
             self.all_ts[-1][-1:],
-            self.all_ys[-1][:, -1:],
+            all_ys,
             self.all_models[-1:],
             self.all_inputs[-1:],
             self.t_event,
@@ -583,7 +602,7 @@ def _update_variable(self, variable):
         for i, (model, ys, inputs, var_pybamm) in enumerate(
             zip(self.all_models, self.all_ys, self.all_inputs, vars_pybamm)
         ):
-            if ys.size == 0 and var_pybamm.has_symbol_of_classes(
+            if self.variables_returned and var_pybamm.has_symbol_of_classes(
                 pybamm.expression_tree.state_vector.StateVector
             ):
                 raise KeyError(
@@ -678,7 +697,7 @@ def __getitem__(self, key):
 
         Returns
         -------
-        :class:`pybamm.ProcessedVariable`
+        :class:`pybamm.ProcessedVariable` or :class:`pybamm.ProcessedVariableComputed`
             A variable that can be evaluated at any time or spatial point. The
             underlying data for this variable is available in its attribute ".data"
         """
@@ -946,6 +965,7 @@ def __add__(self, other):
             other.termination,
             all_sensitivities=all_sensitivities,
             all_yps=all_yps,
+            variables_returned=other.variables_returned,
         )
 
         new_sol.closest_event_idx = other.closest_event_idx
@@ -962,6 +982,19 @@ def __add__(self, other):
         # Set sub_solutions
         new_sol._sub_solutions = self.sub_solutions + other.sub_solutions
 
+        # update variables which were derived at the solver stage
+        if other._variables and all(
+            isinstance(v, pybamm.ProcessedVariableComputed)
+            for v in other._variables.values()
+        ):
+            if not self._variables:
+                new_sol._variables = other._variables.copy()
+            else:
+                new_sol._variables = {
+                    v: self._variables[v]._update(other._variables[v], new_sol)
+                    for v in self._variables.keys()
+                }
+
         return new_sol
 
     def __radd__(self, other):
@@ -979,6 +1012,7 @@ def copy(self):
             self.termination,
             self._all_sensitivities,
             self.all_yps,
+            self.variables_returned,
         )
         new_sol._all_inputs_casadi = self.all_inputs_casadi
         new_sol._sub_solutions = self.sub_solutions
@@ -988,6 +1022,13 @@ def copy(self):
         new_sol.integration_time = self.integration_time
         new_sol.set_up_time = self.set_up_time
 
+        # copy over variables which were derived at the solver stage
+        if self._variables and all(
+            isinstance(v, pybamm.ProcessedVariableComputed)
+            for v in self._variables.values()
+        ):
+            new_sol._variables = self._variables.copy()
+
         return new_sol
 
     def plot_voltage_components(
@@ -1090,6 +1131,7 @@ def make_cycle_solution(
         sum_sols.termination,
         sum_sols._all_sensitivities,
         sum_sols.all_yps,
+        sum_sols.variables_returned,
     )
     cycle_solution._all_inputs_casadi = sum_sols.all_inputs_casadi
     cycle_solution._sub_solutions = sum_sols.sub_solutions

tests/integration/test_solvers/test_idaklu.py

@@ -152,3 +152,55 @@ def test_interpolation(self):
             # test that y[1:3] = to true solution
             true_solution = b_value * sol.t
             np.testing.assert_array_almost_equal(sol.y[1:3], true_solution)
+
+    def test_with_experiments(self):
+        summary_vars = []
+        sols = []
+        for out_vars in [True, False]:
+            model = pybamm.lithium_ion.SPM()
+
+            if out_vars:
+                output_variables = [
+                    "Discharge capacity [A.h]",  # 0D variables
+                    "Time [s]",
+                    "Current [A]",
+                    "Voltage [V]",
+                    "Pressure [Pa]",  # 1D variable
+                    "Positive particle effective diffusivity [m2.s-1]",  # 2D variable
+                ]
+            else:
+                output_variables = None
+
+            solver = pybamm.IDAKLUSolver(output_variables=output_variables)
+
+            experiment = pybamm.Experiment(
+                [
+                    (
+                        "Charge at 1C until 4.2 V",
+                        "Hold at 4.2 V until C/50",
+                        "Rest for 1 hour",
+                    )
+                ]
+            )
+
+            sim = pybamm.Simulation(
+                model,
+                experiment=experiment,
+                solver=solver,
+            )
+
+            sol = sim.solve()
+            sols.append(sol)
+            summary_vars.append(sol.summary_variables)
+
+        # check computed variables are propegated sucessfully
+        np.testing.assert_array_equal(
+            sols[0]["Pressure [Pa]"].data, sols[1]["Pressure [Pa]"].data
+        )
+        np.testing.assert_array_almost_equal(
+            sols[0]["Voltage [V]"].data, sols[1]["Voltage [V]"].data
+        )
+
+        # check summary variables are the same if output variables are specified
+        for var in summary_vars[0].keys():
+            assert summary_vars[0][var] == summary_vars[1][var]

tests/unit/test_solvers/test_idaklu_solver.py

@@ -939,6 +939,9 @@ def construct_model():
             with pytest.raises(KeyError):
                 sol[varname].data
 
+        # Check Solution is marked
+        assert sol.variables_returned is True
+
         # Mock a 1D current collector and initialise (none in the model)
         sol["x_s [m]"].domain = ["current collector"]
         sol["x_s [m]"].entries