merge develop

.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@294a9d92911152fe08befb9ec03e240add280cb3 # v3.26.8
+        uses: github/codeql-action/upload-sarif@e2b3eafc8d227b0241d48be5f425d47c2d750a13 # v3.26.10
         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.6.7"
+    rev: "v0.6.8"
     hooks:
       - id: ruff
         args: [--fix, --show-fixes]

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- 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))
 - Added basic telemetry to record which functions are being run. See [Telemetry section in the User Guide](https://docs.pybamm.org/en/latest/source/user_guide/index.html#telemetry) for more information. ([#4441](https://github.com/pybamm-team/PyBaMM/pull/4441))
 - Added sensitivity calculation support for `pybamm.Simulation` and `pybamm.Experiment` ([#4415](https://github.com/pybamm-team/PyBaMM/pull/4415))
 - Added OpenMP parallelization to IDAKLU solver for lists of input parameters ([#4449](https://github.com/pybamm-team/PyBaMM/pull/4449))

CMakeLists.txt

@@ -105,6 +105,8 @@ pybind11_add_module(idaklu
   src/pybamm/solvers/c_solvers/idaklu/Expressions/Base/Expression.hpp
   src/pybamm/solvers/c_solvers/idaklu/Expressions/Base/ExpressionSet.hpp
   src/pybamm/solvers/c_solvers/idaklu/Expressions/Base/ExpressionTypes.hpp
+  src/pybamm/solvers/c_solvers/idaklu/observe.hpp
+  src/pybamm/solvers/c_solvers/idaklu/observe.cpp
   # IDAKLU expressions - concrete implementations
   ${IDAKLU_EXPR_CASADI_SOURCE_FILES}
   ${IDAKLU_EXPR_IREE_SOURCE_FILES}

docs/source/user_guide/installation/index.rst

@@ -145,8 +145,8 @@ Dependency
 `pre-commit <https://pre-commit.com/index.html>`__                               \-                 dev                For managing and maintaining multi-language pre-commit hooks.
 `ruff <https://beta.ruff.rs/docs/>`__                                            \-                 dev                For code formatting.
 `nox <https://nox.thea.codes/en/stable/>`__                                      \-                 dev                For running testing sessions in multiple environments.
+`pytest-subtests <https://github.com/pytest-dev/pytest-subtests/>`__             \-                 dev                For subtests pytest fixture.
 `pytest-cov <https://pytest-cov.readthedocs.io/en/stable/>`__                    \-                 dev                For calculating test coverage.
-`parameterized <https://github.com/wolever/parameterized>`__                     \-                 dev                For test parameterization.
 `pytest <https://docs.pytest.org/en/stable/>`__                                  6.0.0              dev                For running the test suites.
 `pytest-doctestplus <https://github.com/scientific-python/pytest-doctestplus>`__ \-                 dev                For running doctests.
 `pytest-xdist <https://pytest-xdist.readthedocs.io/en/latest/>`__                \-                 dev                For running tests in parallel across distributed workers.

pyproject.toml

@@ -85,6 +85,7 @@ plot = [
     "matplotlib>=3.6.0",
 ]
 cite = [
+    "setuptools", # Fix for a pybtex issue
     "pybtex>=0.24.0",
 ]
 # Battery Parameter eXchange format
@@ -106,12 +107,11 @@ dev = [
     "pytest-cov",
     # For doctest
     "pytest-doctestplus",
-    # For test parameterization
-    "parameterized>=0.9",
     # pytest and its plugins
     "pytest>=6",
     "pytest-xdist",
     "pytest-mock",
+    "pytest-subtests",
     # For testing Jupyter notebooks
     "nbmake",
     # To access the metadata for python packages
@@ -198,6 +198,8 @@ extend-select = [
     "YTT",           # flake8-2020
     "TID252",        # relative-imports
     "S101",          # to identify use of assert statement
+    "PT027",         # remove unittest style assertion
+    "PT009",         # Use pytest.raises instead of unittest-style
 ]
 ignore = [
     "E741",        # Ambiguous variable name
@@ -231,6 +233,7 @@ minversion = "8"
 required_plugins = [
     "pytest-xdist",
     "pytest-mock",
+    "pytest-subtests",
 ]
 norecursedirs = 'pybind11*'
 addopts = [

setup.py

@@ -327,6 +327,8 @@ def compile_KLU():
         "src/pybamm/solvers/c_solvers/idaklu/Solution.hpp",
         "src/pybamm/solvers/c_solvers/idaklu/Options.hpp",
         "src/pybamm/solvers/c_solvers/idaklu/Options.cpp",
+        "src/pybamm/solvers/c_solvers/idaklu/observe.hpp",
+        "src/pybamm/solvers/c_solvers/idaklu/observe.cpp",
         "src/pybamm/solvers/c_solvers/idaklu.cpp",
     ],
 )

src/pybamm/__init__.py

@@ -157,7 +157,7 @@
 
 # Solver classes
 from .solvers.solution import Solution, EmptySolution, make_cycle_solution
-from .solvers.processed_variable import ProcessedVariable
+from .solvers.processed_variable import ProcessedVariable, process_variable
 from .solvers.processed_variable_computed import ProcessedVariableComputed
 from .solvers.base_solver import BaseSolver
 from .solvers.dummy_solver import DummySolver

src/pybamm/plotting/quick_plot.py

@@ -419,14 +419,14 @@ def reset_axis(self):
                 spatial_vars = self.spatial_variable_dict[key]
                 var_min = np.min(
                     [
-                        ax_min(var(self.ts_seconds[i], **spatial_vars, warn=False))
+                        ax_min(var(self.ts_seconds[i], **spatial_vars))
                         for i, variable_list in enumerate(variable_lists)
                         for var in variable_list
                     ]
                 )
                 var_max = np.max(
                     [
-                        ax_max(var(self.ts_seconds[i], **spatial_vars, warn=False))
+                        ax_max(var(self.ts_seconds[i], **spatial_vars))
                         for i, variable_list in enumerate(variable_lists)
                         for var in variable_list
                     ]
@@ -512,7 +512,7 @@ def plot(self, t, dynamic=False):
                         full_t = self.ts_seconds[i]
                         (self.plots[key][i][j],) = ax.plot(
                             full_t / self.time_scaling_factor,
-                            variable(full_t, warn=False),
+                            variable(full_t),
                             color=self.colors[i],
                             linestyle=linestyle,
                         )
@@ -548,7 +548,7 @@ def plot(self, t, dynamic=False):
                             linestyle = self.linestyles[j]
                         (self.plots[key][i][j],) = ax.plot(
                             self.first_spatial_variable[key],
-                            variable(t_in_seconds, **spatial_vars, warn=False),
+                            variable(t_in_seconds, **spatial_vars),
                             color=self.colors[i],
                             linestyle=linestyle,
                             zorder=10,
@@ -570,13 +570,13 @@ def plot(self, t, dynamic=False):
                     y_name = next(iter(spatial_vars.keys()))[0]
                     x = self.second_spatial_variable[key]
                     y = self.first_spatial_variable[key]
-                    var = variable(t_in_seconds, **spatial_vars, warn=False)
+                    var = variable(t_in_seconds, **spatial_vars)
                 else:
                     x_name = next(iter(spatial_vars.keys()))[0]
                     y_name = list(spatial_vars.keys())[1][0]
                     x = self.first_spatial_variable[key]
                     y = self.second_spatial_variable[key]
-                    var = variable(t_in_seconds, **spatial_vars, warn=False).T
+                    var = variable(t_in_seconds, **spatial_vars).T
                 ax.set_xlabel(f"{x_name} [{self.spatial_unit}]")
                 ax.set_ylabel(f"{y_name} [{self.spatial_unit}]")
                 vmin, vmax = self.variable_limits[key]
@@ -710,7 +710,6 @@ def slider_update(self, t):
                         var = variable(
                             time_in_seconds,
                             **self.spatial_variable_dict[key],
-                            warn=False,
                         )
                         plot[i][j].set_ydata(var)
                         var_min = min(var_min, ax_min(var))
@@ -729,11 +728,11 @@ def slider_update(self, t):
                 if self.x_first_and_y_second[key] is False:
                     x = self.second_spatial_variable[key]
                     y = self.first_spatial_variable[key]
-                    var = variable(time_in_seconds, **spatial_vars, warn=False)
+                    var = variable(time_in_seconds, **spatial_vars)
                 else:
                     x = self.first_spatial_variable[key]
                     y = self.second_spatial_variable[key]
-                    var = variable(time_in_seconds, **spatial_vars, warn=False).T
+                    var = variable(time_in_seconds, **spatial_vars).T
                 # store the plot and the var data (for testing) as cant access
                 # z data from QuadMesh or QuadContourSet object
                 if self.is_y_z[key] is True:

src/pybamm/solvers/c_solvers/idaklu.cpp

@@ -9,6 +9,7 @@
 #include <pybind11/stl_bind.h>
 
 #include "idaklu/idaklu_solver.hpp"
+#include "idaklu/observe.hpp"
 #include "idaklu/IDAKLUSolverGroup.hpp"
 #include "idaklu/IdakluJax.hpp"
 #include "idaklu/common.hpp"
@@ -27,13 +28,15 @@ casadi::Function generate_casadi_function(const std::string &data)
 namespace py = pybind11;
 
 PYBIND11_MAKE_OPAQUE(std::vector<np_array>);
+PYBIND11_MAKE_OPAQUE(std::vector<np_array_realtype>);
 PYBIND11_MAKE_OPAQUE(std::vector<Solution>);
 
 PYBIND11_MODULE(idaklu, m)
 {
   m.doc() = "sundials solvers"; // optional module docstring
 
   py::bind_vector<std::vector<np_array>>(m, "VectorNdArray");
+  py::bind_vector<std::vector<np_array_realtype>>(m, "VectorRealtypeNdArray");
   py::bind_vector<std::vector<Solution>>(m, "VectorSolution");
 
   py::class_<IDAKLUSolverGroup>(m, "IDAKLUSolverGroup")
@@ -72,6 +75,27 @@ PYBIND11_MODULE(idaklu, m)
     py::arg("options"),
     py::return_value_policy::take_ownership);
 
+  m.def("observe", &observe,
+    "Observe variables",
+    py::arg("ts"),
+    py::arg("ys"),
+    py::arg("inputs"),
+    py::arg("funcs"),
+    py::arg("is_f_contiguous"),
+    py::arg("shape"),
+    py::return_value_policy::take_ownership);
+
+  m.def("observe_hermite_interp", &observe_hermite_interp,
+    "Observe and Hermite interpolate variables",
+    py::arg("t_interp"),
+    py::arg("ts"),
+    py::arg("ys"),
+    py::arg("yps"),
+    py::arg("inputs"),
+    py::arg("funcs"),
+    py::arg("shape"),
+    py::return_value_policy::take_ownership);
+
 #ifdef IREE_ENABLE
   m.def("create_iree_solver_group", &create_idaklu_solver_group<IREEFunctions>,
     "Create a group of iree idaklu solver objects",
@@ -167,7 +191,9 @@ PYBIND11_MODULE(idaklu, m)
   py::class_<Solution>(m, "solution")
     .def_readwrite("t", &Solution::t)
     .def_readwrite("y", &Solution::y)
+    .def_readwrite("yp", &Solution::yp)
     .def_readwrite("yS", &Solution::yS)
+    .def_readwrite("ypS", &Solution::ypS)
     .def_readwrite("y_term", &Solution::y_term)
     .def_readwrite("flag", &Solution::flag);
 }

src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp

@@ -54,9 +54,9 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver
   int const number_of_events;  // cppcheck-suppress unusedStructMember
   int number_of_timesteps;
   int precon_type;  // cppcheck-suppress unusedStructMember
-  N_Vector yy, yp, y_cache, avtol;  // y, y', y cache vector, and absolute tolerance
+  N_Vector yy, yyp, y_cache, avtol;  // y, y', y cache vector, and absolute tolerance
   N_Vector *yyS;  // cppcheck-suppress unusedStructMember
-  N_Vector *ypS;  // cppcheck-suppress unusedStructMember
+  N_Vector *yypS;  // cppcheck-suppress unusedStructMember
   N_Vector id;              // rhs_alg_id
   realtype rtol;
   int const jac_times_cjmass_nnz;  // cppcheck-suppress unusedStructMember
@@ -70,11 +70,14 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver
   vector<realtype> res_dvar_dp;
   bool const sensitivity;  // cppcheck-suppress unusedStructMember
   bool const save_outputs_only; // cppcheck-suppress unusedStructMember
+  bool save_hermite;  // cppcheck-suppress unusedStructMember
   bool is_ODE;  // cppcheck-suppress unusedStructMember
   int length_of_return_vector;  // cppcheck-suppress unusedStructMember
   vector<realtype> t;  // cppcheck-suppress unusedStructMember
   vector<vector<realtype>> y;  // cppcheck-suppress unusedStructMember
+  vector<vector<realtype>> yp;  // cppcheck-suppress unusedStructMember
   vector<vector<vector<realtype>>> yS;  // cppcheck-suppress unusedStructMember
+  vector<vector<vector<realtype>>> ypS;  // cppcheck-suppress unusedStructMember
   SetupOptions const setup_opts;
   SolverOptions const solver_opts;
 
@@ -144,6 +147,11 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver
    */
   void InitializeStorage(int const N);
 
+  /**
+   * @brief Initialize the storage for Hermite interpolation
+   */
+  void InitializeHermiteStorage(int const N);
+
   /**
    * @brief Apply user-configurable IDA options
    */
@@ -190,13 +198,20 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver
    */
   void ExtendAdaptiveArrays();
 
+  /**
+   * @brief Extend the Hermite interpolation info by 1
+   */
+  void ExtendHermiteArrays();
+
   /**
    * @brief Set the step values
    */
   void SetStep(
-    realtype &t_val,
+    realtype &tval,
     realtype *y_val,
+    realtype *yp_val,
     vector<realtype *> const &yS_val,
+    vector<realtype *> const &ypS_val,
     int &i_save
   );
 
@@ -211,7 +226,9 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver
     realtype &t_prev,
     realtype const &t_next,
     realtype *y_val,
+    realtype *yp_val,
     vector<realtype *> const &yS_val,
+    vector<realtype *> const &ypS_val,
     int &i_save
   );
 
@@ -255,6 +272,26 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver
     int &i_save
   );
 
+  /**
+   * @brief Save the output function results at the requested time
+   */
+  void SetStepHermite(
+    realtype &t_val,
+    realtype *yp_val,
+    const vector<realtype*> &ypS_val,
+    int &i_save
+  );
+
+  /**
+   * @brief Save the output function sensitivities at the requested time
+   */
+  void SetStepHermiteSensitivities(
+    realtype &t_val,
+    realtype *yp_val,
+    const vector<realtype*> &ypS_val,
+    int &i_save
+  );
+
 };
 
 #include "IDAKLUSolverOpenMP.inl"