Allow lists of input parameters as initial conditions

.github/workflows/benchmark_on_push.yml

@@ -40,7 +40,7 @@ jobs:
           enable-cache: true
 
       - name: Install python dependencies
-        run: uv pip install asv[virtualenv]
+        run: uv pip install --system asv
 
       - name: Fetch base branch
         run: |

.github/workflows/periodic_benchmarks.yml

@@ -48,7 +48,7 @@ jobs:
           enable-cache: true
 
       - name: Install python dependencies
-        run: uv pip install asv[virtualenv]
+        run: uv pip install --system asv
 
       - name: Run benchmarks
         run: |

CHANGELOG.md

@@ -2,8 +2,12 @@
 
 ## Features
 
+- Adds `silence_sundials_errors` IDAKLU solver option with `default=False` to match historical output. ([#5290](https://github.com/pybamm-team/PyBaMM/pull/5290))
+
 ## Bug fixes
 
+- Fix a bug with serialising `InputParameter`s. ([#5289](https://github.com/pybamm-team/PyBaMM/pull/5289))
+
 # [v25.10.1](https://github.com/pybamm-team/PyBaMM/tree/v25.10.1) - 2025-11-14
 
 ## Features
@@ -19,6 +23,7 @@
 
 ## Features
 
+- Added uniform grid sizing across subdomains in the x-dimension, ensuring consistent grid spacing when geometries have varying lengths. ([#5253](https://github.com/pybamm-team/PyBaMM/pull/5253))
 - Added the `electrode_phases` kwarg to `plot_voltage_components()` which allows choosing between plotting primary or secondary phase overpotentials. ([#5229](https://github.com/pybamm-team/PyBaMM/pull/5229))
 - Added the `num_steps_no_progress` and `t_no_progress` options in the `IDAKLUSolver` to early terminate the simulation if little progress is detected. ([#5201](https://github.com/pybamm-team/PyBaMM/pull/5201))
 - EvaluateAt symbol: add support for children evaluated at edges ([#5190](https://github.com/pybamm-team/PyBaMM/pull/5190))

docs/source/examples/notebooks/creating_models/4-comparing-full-and-reduced-order-models.ipynb

@@ -24,7 +24,7 @@
     "$$\n",
     "\\left.c\\right\\vert_{t=0} = c_0,\n",
     "$$\n",
-    "where $c$$ is the concentration, $r$ the radial coordinate, $t$ time, $R$ the particle radius, $D$ the diffusion coefficient, $j$ the interfacial current density, $F$ Faraday's constant, and $c_0$ the initial concentration. \n",
+    "where $c$ is the concentration, $r$ the radial coordinate, $t$ time, $R$ the particle radius, $D$ the diffusion coefficient, $j$ the interfacial current density, $F$ Faraday's constant, and $c_0$ the initial concentration. \n",
     "\n",
     "As in the previous example we use the following parameters:\n",
     "\n",

src/pybamm/expression_tree/binary_operators.py

@@ -213,6 +213,10 @@ def to_json(self):
 
         return json_dict
 
+    def _t_discon(self, expr, y0, inputs, num_events):
+        """Returns the discontinuity time-points for a function."""
+        raise NotImplementedError(f"`_t_discon` not implemented for {self.name}")
+
 
 class Power(BinaryOperator):
     """
@@ -673,6 +677,18 @@ def _evaluate_for_shape(self):
         # an array of NaNs
         return self._binary_evaluate(left, right) * np.nan
 
+    def _t_discon(self, expr, y0, inputs, num_events):
+        """Returns the discontinuity time-points for the heaviside function."""
+        value = expr.evaluate(0, y0, inputs=inputs)
+        t_discon = [value]
+        t_discon.append(self._t_discon_next(value))
+        return t_discon
+
+    def _t_discon_next(self, value: float):
+        raise NotImplementedError(
+            "_t_discon_next method should be implemented in subclasses of _Heaviside"
+        )  # pragma: no cover
+
 
 class EqualHeaviside(_Heaviside):
     """A heaviside function with equality (return 1 when left = right)"""
@@ -695,6 +711,16 @@ def _binary_evaluate(self, left, right):
         with np.errstate(invalid="ignore"):
             return left <= right
 
+    def _t_discon_next(self, value: float):
+        if self.left == pybamm.t:
+            # t <= x
+            # Stop at t = x and right after t = x
+            return np.nextafter(value, np.inf)
+        else:
+            # t >= x
+            # Stop at t = x and right before t = x
+            return np.nextafter(value, -np.inf)
+
 
 class NotEqualHeaviside(_Heaviside):
     """A heaviside function without equality (return 0 when left = right)"""
@@ -716,6 +742,16 @@ def _binary_evaluate(self, left, right):
         with np.errstate(invalid="ignore"):
             return left < right
 
+    def _t_discon_next(self, value: float):
+        if self.left == pybamm.t:
+            # t < x
+            # Stop at t = x and right before t = x
+            return np.nextafter(value, -np.inf)
+        else:
+            # t > x
+            # Stop at t = x and right after t = x
+            return np.nextafter(value, np.inf)
+
 
 class Modulo(BinaryOperator):
     """Calculates the remainder of an integer division."""
@@ -758,6 +794,17 @@ def _binary_evaluate(self, left, right):
         """See :meth:`pybamm.BinaryOperator._binary_evaluate()`."""
         return left % right
 
+    def _t_discon(self, expr, y0, inputs, num_events):
+        value = expr.evaluate(0, y0, inputs=inputs)
+        t_discon = []
+
+        for i in np.arange(num_events):
+            t = value * (i + 1)
+            # Stop right before t and at t
+            t_discon.append(np.nextafter(t, -np.inf))
+            t_discon.append(t)
+        return t_discon
+
 
 class Minimum(BinaryOperator):
     """Returns the smaller of two objects."""

src/pybamm/expression_tree/operations/serialise.py

@@ -1625,6 +1625,8 @@ def convert_symbol_from_json(json_data):
     elif json_data["type"] == "Parameter":
         # Convert stored parameters back to PyBaMM Parameter objects
         return pybamm.Parameter(json_data["name"])
+    elif json_data["type"] == "InputParameter":
+        return pybamm.InputParameter(json_data["name"])
     elif json_data["type"] == "Scalar":
         # Convert stored numerical values back to PyBaMM Scalar objects
         return pybamm.Scalar(json_data["value"])

src/pybamm/meshes/meshes.py

@@ -9,6 +9,39 @@
 import pybamm
 
 
+def compute_var_pts_from_thicknesses(electrode_thicknesses, grid_size):
+    """
+    Compute a ``var_pts`` dictionary using electrode thicknesses and a target cell size (dx).
+
+    Added as per maintainer feedback in issue #<your-issue-number> to make mesh generation
+    explicit — ``grid_size`` now represents the mesh cell size in metres.
+
+    Parameters
+    ----------
+    electrode_thicknesses : dict
+        Domain thicknesses in metres.
+    grid_size : float
+        Desired uniform mesh cell size (m).
+
+    Returns
+    -------
+    dict
+        Mapping of each domain to its computed grid points.
+    """
+    if not isinstance(electrode_thicknesses, dict):
+        raise TypeError("electrode_thicknesses must be a dictionary")
+
+    if not isinstance(grid_size, (int | float)) or grid_size <= 0:
+        raise ValueError("grid_size must be a positive number")
+
+    var_pts = {}
+    for domain, thickness in electrode_thicknesses.items():
+        npts = max(round(thickness / grid_size), 2)
+        var_pts[domain] = {f"x_{domain[0]}": npts}
+
+    return var_pts
+
+
 class Mesh(dict):
     """
     Mesh contains a list of submeshes on each subdomain.

src/pybamm/models/base_model.py

@@ -468,6 +468,17 @@ def algebraic_root_solver(self):
     def algebraic_root_solver(self, algebraic_root_solver):
         self._algebraic_root_solver = algebraic_root_solver
 
+    @property
+    def y0(self):
+        if not hasattr(self, "y0_list") or self.y0_list is None:
+            return None
+        elif len(self.y0_list) == 1:
+            return self.y0_list[0]
+        else:
+            raise ValueError(
+                "Model contains multiple initial states. Access using y0_list instead."
+            )
+
     def get_parameter_info(self, by_submodel=False):
         """
         Extracts the parameter information and returns it as a dictionary.
@@ -906,7 +917,12 @@ def _build_model(self):
         self.build_model_equations()
 
     def set_initial_conditions_from(
-        self, solution, inplace=True, return_type="model", mesh=None
+        self,
+        solution,
+        inputs=None,
+        inplace=True,
+        return_type="model",
+        mesh=None,
     ):
         """
         Update initial conditions with the final states from a Solution object or from
@@ -918,6 +934,8 @@ def set_initial_conditions_from(
         ----------
         solution : :class:`pybamm.Solution`, or dict
             The solution to use to initialize the model
+        inputs : dict
+            The dictionary of model input parameters.
         inplace : bool, optional
             Whether to modify the model inplace or create a new model (default True)
         return_type : str, optional
@@ -1081,7 +1099,7 @@ def get_variable_state(var):
                 scale, reference = pybamm.Scalar(1), pybamm.Scalar(0)
             initial_conditions[var] = (
                 pybamm.Vector(final_state_eval) - reference
-            ) / scale.evaluate()
+            ) / scale.evaluate(inputs=inputs)
 
         # Also update the concatenated initial conditions if the model is already
         # discretised

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

@@ -561,14 +561,14 @@ def _set_up_solve(self, inputs, direction):
     def _solve_full(self, inputs, ics, direction):
         sim = self._get_electrode_soh_sims_full(direction)
         sim.build()
-        sim.built_model.set_initial_conditions_from(ics)
+        sim.built_model.set_initial_conditions_from(ics, inputs=inputs)
         sol = sim.solve([0], inputs=inputs)
         return sol
 
     def _solve_split(self, inputs, ics, direction):
         x100_sim, x0_sim = self._get_electrode_soh_sims_split(direction)
         x100_sim.build()
-        x100_sim.built_model.set_initial_conditions_from(ics)
+        x100_sim.built_model.set_initial_conditions_from(ics, inputs=inputs)
         x100_sol = x100_sim.solve([0], inputs=inputs)
         if self.options["open-circuit potential"] == "MSMR":
             inputs["Un(x_100)"] = x100_sol["Un(x_100)"].data[0]
@@ -577,7 +577,7 @@ def _solve_split(self, inputs, ics, direction):
             inputs["x_100"] = x100_sol["x_100"].data[0]
             inputs["y_100"] = x100_sol["y_100"].data[0]
         x0_sim.build()
-        x0_sim.built_model.set_initial_conditions_from(ics)
+        x0_sim.built_model.set_initial_conditions_from(ics, inputs=inputs)
         x0_sol = x0_sim.solve([0], inputs=inputs)
 
         return x0_sol

src/pybamm/models/submodels/interface/kinetics/butler_volmer.py

@@ -12,7 +12,7 @@ class SymmetricButlerVolmer(BaseKinetics):
     Submodel which implements the symmetric forward Butler-Volmer equation:
 
     .. math::
-        j = 2 * j_0(c) * \\sinh(ne * F * \\eta_r(c) / RT)
+        j = 2 * j_0(c) * \\sinh(ne * F * \\eta_r(c) / 2RT)
 
     Parameters
     ----------

src/pybamm/models/submodels/interface/open_circuit_potential/base_hysteresis_ocp.py

@@ -84,7 +84,7 @@ def _get_coupled_variables(self, variables):
             U_eq = self.phase_param.U(sto_surf, T)
             U_eq_x_av = self.phase_param.U(sto_surf, T)
             U_lith = self.phase_param.U(sto_surf, T, "lithiation")
-            U_lith_bulk = self.phase_param.U(sto_bulk, T_bulk)
+            U_lith_bulk = self.phase_param.U(sto_bulk, T_bulk, "lithiation")
             U_delith = self.phase_param.U(sto_surf, T, "delithiation")
             U_delith_bulk = self.phase_param.U(sto_bulk, T_bulk, "delithiation")
 

src/pybamm/parameters/parameter_values.py

@@ -866,12 +866,21 @@ def _process_function_parameter(self, symbol):
                 else:
                     new_children.append(self.process_symbol(child))
 
-            # Get the expression and inputs for the function
+            # Get the expression and inputs for the function.
+            # func_args may include arguments that were not explicitly wired up
+            # in this FunctionParameter (e.g., kwargs with default values). After
+            # serialisation/deserialisation, we only recover the children that were
+            # actually connected.
+            #
+            # Using strict=True here therefore raises a ValueError when there are
+            # more args than children. We allow func_args to be longer than
+            # symbol.children and only build the mapping for the args for which we
+            # actually have children.
             expression = function_parameter.child
             inputs = {
                 arg: child
                 for arg, child in zip(
-                    function_parameter.func_args, symbol.children, strict=True
+                    function_parameter.func_args, symbol.children, strict=False
                 )
             }
 

src/pybamm/simulation.py

@@ -128,6 +128,7 @@ def __init__(
         self._model_with_set_params = None
         self._built_model = None
         self._built_initial_soc = None
+        self._built_nominal_capacity = None
         self.steps_to_built_models = None
         self.steps_to_built_solvers = None
         self._mesh = None
@@ -163,6 +164,42 @@ def set_up_and_parameterise_experiment(self, solve_kwargs=None):
         warnings.warn(msg, DeprecationWarning, stacklevel=2)
         self._set_up_and_parameterise_experiment(solve_kwargs=solve_kwargs)
 
+    def _update_experiment_models_for_capacity(self, solve_kwargs=None):
+        """
+        Check if the nominal capacity has changed and update the experiment models
+        if needed. This re-processes the models without rebuilding the mesh and
+        discretisation.
+        """
+        current_capacity = self._parameter_values.get(
+            "Nominal cell capacity [A.h]", None
+        )
+
+        if self._built_nominal_capacity == current_capacity:
+            return
+
+        # Capacity has changed, need to re-process the models
+        pybamm.logger.info(
+            f"Nominal capacity changed from {self._built_nominal_capacity} to "
+            f"{current_capacity}. Re-processing experiment models."
+        )
+
+        # Re-parameterise the experiment with the new capacity
+        self._set_up_and_parameterise_experiment(solve_kwargs)
+
+        # Re-discretise the models
+        self.steps_to_built_models = {}
+        self.steps_to_built_solvers = {}
+        for (
+            step,
+            model_with_set_params,
+        ) in self.experiment_unique_steps_to_model.items():
+            built_model = self._disc.process_model(model_with_set_params, inplace=True)
+            solver = self._solver.copy()
+            self.steps_to_built_solvers[step] = solver
+            self.steps_to_built_models[step] = built_model
+
+        self._built_nominal_capacity = current_capacity
+
     def _set_up_and_parameterise_experiment(self, solve_kwargs=None):
         """
         Create and parameterise the models for each step in the experiment.
@@ -266,6 +303,7 @@ def set_initial_state(self, initial_soc, direction=None, inputs=None):
             # reset
             self._model_with_set_params = None
             self._built_model = None
+            self._built_nominal_capacity = None
             self.steps_to_built_models = None
             self.steps_to_built_solvers = None
 
@@ -338,6 +376,8 @@ def build_for_experiment(
             self.set_initial_state(initial_soc, direction=direction, inputs=inputs)
 
         if self.steps_to_built_models:
+            # Check if we need to update the models due to capacity change
+            self._update_experiment_models_for_capacity(solve_kwargs)
             return
         else:
             self._set_up_and_parameterise_experiment(solve_kwargs)
@@ -366,6 +406,10 @@ def build_for_experiment(
                 self.steps_to_built_solvers[step] = solver
                 self.steps_to_built_models[step] = built_model
 
+            self._built_nominal_capacity = self._parameter_values.get(
+                "Nominal cell capacity [A.h]", None
+            )
+
     def solve(
         self,
         t_eval=None,
@@ -379,7 +423,6 @@ def solve(
         showprogress=False,
         inputs=None,
         t_interp=None,
-        initial_conditions=None,
         **kwargs,
     ):
         """
@@ -542,7 +585,6 @@ def solve(
                 inputs=inputs,
                 t_interp=t_interp,
                 **kwargs,
-                initial_conditions=initial_conditions,
             )
 
         elif self.operating_mode == "with experiment":