Merge branch 'develop' into output_vars_extrapolation

.github/CODEOWNERS

@@ -10,7 +10,7 @@ src/pybamm/meshes/ @martinjrobins @rtimms @valentinsulzer @rtimms
 src/pybamm/models/ @brosaplanella @DrSOKane @rtimms @valentinsulzer @TomTranter @rtimms
 src/pybamm/parameters/ @brosaplanella @DrSOKane @rtimms @valentinsulzer @TomTranter @rtimms  @kratman
 src/pybamm/plotting/ @martinjrobins @rtimms @Saransh-cpp @valentinsulzer @rtimms  @kratman @agriyakhetarpal
-src/pybamm/solvers/ @martinjrobins @rtimms @valentinsulzer @TomTranter @rtimms
+src/pybamm/solvers/ @martinjrobins @rtimms @valentinsulzer @TomTranter @rtimms @MarcBerliner
 src/pybamm/spatial_methods/ @martinjrobins @rtimms @valentinsulzer @rtimms
 src/pybamm/* @pybamm-team/maintainers # the files directly under /pybamm/, will not recurse
 

CHANGELOG.md

@@ -1,20 +1,24 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
 ## Features
+
 - 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))
+- Added a lithium ion equivalent circuit model with split open circuit voltages for each electrode (`SplitOCVR`). ([#4330](https://github.com/pybamm-team/PyBaMM/pull/4330))
 
 ## Optimizations
+
 - 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))
 
 ## Features
 
 - Added phase-dependent particle options to LAM #4369
 
 ## Breaking changes
+
 - The parameters "... electrode OCP entropic change [V.K-1]" and "... electrode volume change" are now expected to be functions of stoichiometry only instead of functions of both stoichiometry and maximum concentration ([#4427](https://github.com/pybamm-team/PyBaMM/pull/4427))
 - Renamed `set_events` function to `add_events_from` to better reflect its purpose. ([#4421](https://github.com/pybamm-team/PyBaMM/pull/4421))
 
-
 # [v24.9.0](https://github.com/pybamm-team/PyBaMM/tree/v24.9.0) - 2024-09-03
 
 ## Features

CMakeLists.txt

@@ -19,7 +19,7 @@ endif()
 
 project(idaklu)
 
-set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 set(CMAKE_EXPORT_COMPILE_COMMANDS 1)
@@ -82,6 +82,8 @@ pybind11_add_module(idaklu
   src/pybamm/solvers/c_solvers/idaklu/idaklu_solver.hpp
   src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolver.cpp
   src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolver.hpp
+  src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.cpp
+  src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.hpp
   src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl
   src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp
   src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP_solvers.cpp
@@ -94,6 +96,8 @@ pybind11_add_module(idaklu
   src/pybamm/solvers/c_solvers/idaklu/common.cpp
   src/pybamm/solvers/c_solvers/idaklu/Solution.cpp
   src/pybamm/solvers/c_solvers/idaklu/Solution.hpp
+  src/pybamm/solvers/c_solvers/idaklu/SolutionData.cpp
+  src/pybamm/solvers/c_solvers/idaklu/SolutionData.hpp
   src/pybamm/solvers/c_solvers/idaklu/Options.hpp
   src/pybamm/solvers/c_solvers/idaklu/Options.cpp
   # IDAKLU expressions / function evaluation [abstract]
@@ -138,6 +142,23 @@ set_target_properties(
   INSTALL_RPATH_USE_LINK_PATH TRUE
 )
 
+# openmp
+if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
+  execute_process(
+      COMMAND "brew" "--prefix"
+      OUTPUT_VARIABLE HOMEBREW_PREFIX
+      OUTPUT_STRIP_TRAILING_WHITESPACE)
+  if (OpenMP_ROOT)
+    set(OpenMP_ROOT "${OpenMP_ROOT}:${HOMEBREW_PREFIX}/opt/libomp")
+  else()
+    set(OpenMP_ROOT "${HOMEBREW_PREFIX}/opt/libomp")
+  endif()
+endif()
+find_package(OpenMP)
+if(OpenMP_CXX_FOUND)
+    target_link_libraries(idaklu PRIVATE OpenMP::OpenMP_CXX)
+endif()
+
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR})
 # Sundials
 find_package(SUNDIALS REQUIRED)

docs/source/api/models/lithium_ion/ecm_split_ocv.rst

@@ -0,0 +1,7 @@
+Equivalent Circuit Model with Split OCV (SplitOCVR)
+=====================================================
+
+.. autoclass:: pybamm.lithium_ion.SplitOCVR
+    :members:
+
+.. footbibliography::

docs/source/api/models/lithium_ion/index.rst

@@ -12,3 +12,4 @@ Lithium-ion Models
   msmr
   yang2017
   electrode_soh
+  ecm_split_ocv

src/pybamm/models/full_battery_models/lithium_ion/__init__.py

@@ -24,8 +24,9 @@
 from .Yang2017 import Yang2017
 from .mpm import MPM
 from .msmr import MSMR
+from .basic_splitOCVR import SplitOCVR
 
 __all__ = ['Yang2017', 'base_lithium_ion_model', 'basic_dfn',
            'basic_dfn_composite', 'basic_dfn_half_cell', 'basic_spm', 'dfn',
            'electrode_soh', 'electrode_soh_half_cell', 'mpm', 'msmr',
-           'newman_tobias', 'spm', 'spme']
+           'newman_tobias', 'spm', 'spme', 'basic_splitOCVR']

src/pybamm/models/full_battery_models/lithium_ion/basic_splitOCVR.py

@@ -0,0 +1,100 @@
+#
+# Equivalent Circuit Model with split OCV
+#
+import pybamm
+
+
+class SplitOCVR(pybamm.BaseModel):
+    """Basic Equivalent Circuit Model that uses two OCV functions
+    for each electrode. This model is easily parameterizable with minimal parameters.
+    This class differs from the :class: pybamm.equivalent_circuit.Thevenin() due
+    to dual OCV functions to make up the voltage from each electrode.
+
+    Parameters
+    ----------
+    name: str, optional
+        The name of the model.
+    """
+
+    def __init__(self, name="ECM with split OCV"):
+        super().__init__(name)
+
+        ######################
+        # Variables
+        ######################
+        # All variables are only time-dependent
+        # No domain definition needed
+
+        theta_n = pybamm.Variable("Negative particle stoichiometry")
+        theta_p = pybamm.Variable("Positive particle stoichiometry")
+        Q = pybamm.Variable("Discharge capacity [A.h]")
+        V = pybamm.Variable("Voltage [V]")
+
+        # model is isothermal
+        I = pybamm.FunctionParameter("Current function [A]", {"Time [s]": pybamm.t})
+
+        # Capacity equation
+        self.rhs[Q] = I / 3600
+        self.initial_conditions[Q] = pybamm.Scalar(0)
+
+        # Capacity in each electrode
+        Q_n = pybamm.Parameter("Negative electrode capacity [A.h]")
+        Q_p = pybamm.Parameter("Positive electrode capacity [A.h]")
+
+        # State of charge electrode equations
+        theta_n_0 = pybamm.Parameter("Negative electrode initial stoichiometry")
+        theta_p_0 = pybamm.Parameter("Positive electrode initial stoichiometry")
+        self.rhs[theta_n] = -I / Q_n / 3600
+        self.rhs[theta_p] = I / Q_p / 3600
+        self.initial_conditions[theta_n] = theta_n_0
+        self.initial_conditions[theta_p] = theta_p_0
+
+        # Resistance for IR expression
+        R = pybamm.Parameter("Ohmic resistance [Ohm]")
+
+        # Open-circuit potential for each electrode
+        Un = pybamm.FunctionParameter(
+            "Negative electrode OCP [V]", {"Negative particle stoichiometry": theta_n}
+        )
+        Up = pybamm.FunctionParameter(
+            "Positive electrode OCP [V]", {"Positive particle stoichiometry": theta_p}
+        )
+
+        # Voltage expression
+        V = Up - Un - I * R
+
+        # Parameters for Voltage cutoff
+        voltage_high_cut = pybamm.Parameter("Upper voltage cut-off [V]")
+        voltage_low_cut = pybamm.Parameter("Lower voltage cut-off [V]")
+
+        self.variables = {
+            "Negative particle stoichiometry": theta_n,
+            "Positive particle stoichiometry": theta_p,
+            "Current [A]": I,
+            "Discharge capacity [A.h]": Q,
+            "Voltage [V]": V,
+            "Times [s]": pybamm.t,
+            "Positive electrode OCP [V]": Up,
+            "Negative electrode OCP [V]": Un,
+            "Current function [A]": I,
+        }
+
+        # Events specify points at which a solution should terminate
+        self.events += [
+            pybamm.Event("Minimum voltage [V]", V - voltage_low_cut),
+            pybamm.Event("Maximum voltage [V]", voltage_high_cut - V),
+            pybamm.Event("Maximum Negative Electrode stoichiometry", 0.999 - theta_n),
+            pybamm.Event("Maximum Positive Electrode stoichiometry", 0.999 - theta_p),
+            pybamm.Event("Minimum Negative Electrode stoichiometry", theta_n - 0.0001),
+            pybamm.Event("Minimum Positive Electrode stoichiometry", theta_p - 0.0001),
+        ]
+
+    @property
+    def default_quick_plot_variables(self):
+        return [
+            "Voltage [V]",
+            ["Negative particle stoichiometry", "Positive particle stoichiometry"],
+            "Negative electrode OCP [V]",
+            "Positive electrode OCP [V]",
+            "Current [A]",
+        ]

src/pybamm/solvers/base_solver.py

@@ -86,6 +86,10 @@ def supports_interp(self):
     def root_method(self):
         return self._root_method
 
+    @property
+    def supports_parallel_solve(self):
+        return False
+
     @root_method.setter
     def root_method(self, method):
         if method == "casadi":
@@ -896,36 +900,37 @@ def solve(
             pybamm.logger.verbose(
                 f"Calling solver for {t_eval[start_index]} < t < {t_eval[end_index - 1]}"
             )
-            ninputs = len(model_inputs_list)
-            if ninputs == 1:
-                new_solution = self._integrate(
-                    model,
-                    t_eval[start_index:end_index],
-                    model_inputs_list[0],
-                    t_interp=t_interp,
-                )
-                new_solutions = [new_solution]
-            elif model.convert_to_format == "jax":
-                # Jax can parallelize over the inputs efficiently
+            if self.supports_parallel_solve:
+                # Jax and IDAKLU solver can accept a list of inputs
                 new_solutions = self._integrate(
                     model,
                     t_eval[start_index:end_index],
                     model_inputs_list,
                     t_interp,
                 )
             else:
-                with mp.get_context(self._mp_context).Pool(processes=nproc) as p:
-                    new_solutions = p.starmap(
-                        self._integrate,
-                        zip(
-                            [model] * ninputs,
-                            [t_eval[start_index:end_index]] * ninputs,
-                            model_inputs_list,
-                            [t_interp] * ninputs,
-                        ),
+                ninputs = len(model_inputs_list)
+                if ninputs == 1:
+                    new_solution = self._integrate(
+                        model,
+                        t_eval[start_index:end_index],
+                        model_inputs_list[0],
+                        t_interp=t_interp,
                     )
-                    p.close()
-                    p.join()
+                    new_solutions = [new_solution]
+                else:
+                    with mp.get_context(self._mp_context).Pool(processes=nproc) as p:
+                        new_solutions = p.starmap(
+                            self._integrate,
+                            zip(
+                                [model] * ninputs,
+                                [t_eval[start_index:end_index]] * ninputs,
+                                model_inputs_list,
+                                [t_interp] * ninputs,
+                            ),
+                        )
+                        p.close()
+                        p.join()
             # Setting the solve time for each segment.
             # pybamm.Solution.__add__ assumes attribute solve_time.
             solve_time = timer.time()
@@ -995,7 +1000,7 @@ def solve(
             )
 
         # Return solution(s)
-        if ninputs == 1:
+        if len(solutions) == 1:
             return solutions[0]
         else:
             return solutions
@@ -1350,7 +1355,13 @@ def step(
         # Step
         pybamm.logger.verbose(f"Stepping for {t_start_shifted:.0f} < t < {t_end:.0f}")
         timer.reset()
-        solution = self._integrate(model, t_eval, model_inputs, t_interp)
+
+        # API for _integrate is different for JaxSolver and IDAKLUSolver
+        if self.supports_parallel_solve:
+            solutions = self._integrate(model, t_eval, [model_inputs], t_interp)
+            solution = solutions[0]
+        else:
+            solution = self._integrate(model, t_eval, model_inputs, t_interp)
         solution.solve_time = timer.time()
 
         # Check if extrapolation occurred

src/pybamm/solvers/c_solvers/idaklu.cpp

@@ -9,6 +9,7 @@
 #include <pybind11/stl_bind.h>
 
 #include "idaklu/idaklu_solver.hpp"
+#include "idaklu/IDAKLUSolverGroup.hpp"
 #include "idaklu/IdakluJax.hpp"
 #include "idaklu/common.hpp"
 #include "idaklu/Expressions/Casadi/CasadiFunctions.hpp"
@@ -26,15 +27,17 @@ casadi::Function generate_casadi_function(const std::string &data)
 namespace py = pybind11;
 
 PYBIND11_MAKE_OPAQUE(std::vector<np_array>);
+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<Solution>>(m, "VectorSolution");
 
-  py::class_<IDAKLUSolver>(m, "IDAKLUSolver")
-  .def("solve", &IDAKLUSolver::solve,
+  py::class_<IDAKLUSolverGroup>(m, "IDAKLUSolverGroup")
+  .def("solve", &IDAKLUSolverGroup::solve,
     "perform a solve",
     py::arg("t_eval"),
     py::arg("t_interp"),
@@ -43,8 +46,8 @@ PYBIND11_MODULE(idaklu, m)
     py::arg("inputs"),
     py::return_value_policy::take_ownership);
 
-  m.def("create_casadi_solver", &create_idaklu_solver<CasadiFunctions>,
-    "Create a casadi idaklu solver object",
+  m.def("create_casadi_solver_group", &create_idaklu_solver_group<CasadiFunctions>,
+    "Create a group of casadi idaklu solver objects",
     py::arg("number_of_states"),
     py::arg("number_of_parameters"),
     py::arg("rhs_alg"),
@@ -70,8 +73,8 @@ PYBIND11_MODULE(idaklu, m)
     py::return_value_policy::take_ownership);
 
 #ifdef IREE_ENABLE
-  m.def("create_iree_solver", &create_idaklu_solver<IREEFunctions>,
-    "Create a iree idaklu solver object",
+  m.def("create_iree_solver_group", &create_idaklu_solver_group<IREEFunctions>,
+    "Create a group of iree idaklu solver objects",
     py::arg("number_of_states"),
     py::arg("number_of_parameters"),
     py::arg("rhs_alg"),

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

@@ -2,7 +2,8 @@
 #define PYBAMM_IDAKLU_CASADI_SOLVER_HPP
 
 #include "common.hpp"
-#include "Solution.hpp"
+#include "SolutionData.hpp"
+
 
 /**
  * Abstract base class for solutions that can use different solvers and vector
@@ -24,14 +25,17 @@ class IDAKLUSolver
   ~IDAKLUSolver() = default;
 
   /**
-   * @brief Abstract solver method that returns a Solution class
+   * @brief Abstract solver method that executes the solver
    */
-  virtual Solution solve(
-    np_array t_eval_np,
-    np_array t_interp_np,
-    np_array y0_np,
-    np_array yp0_np,
-    np_array_dense inputs) = 0;
+  virtual SolutionData solve(
+    const std::vector<realtype> &t_eval,
+    const std::vector<realtype> &t_interp,
+    const realtype *y0,
+    const realtype *yp0,
+    const realtype *inputs,
+    bool save_adaptive_steps,
+    bool save_interp_steps
+  ) = 0;
 
   /**
    * Abstract method to initialize the solver, once vectors and solver classes