diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS deleted file mode 100644 index 163b26428a..0000000000 --- a/.github/CODEOWNERS +++ /dev/null @@ -1,2 +0,0 @@ -# Automatically request reviews from maintainers -* @pybamm-team/maintainers diff --git a/.github/workflows/benchmark_on_push.yml b/.github/workflows/benchmark_on_push.yml index 18a9e7931d..f70556c782 100644 --- a/.github/workflows/benchmark_on_push.yml +++ b/.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: | diff --git a/.github/workflows/periodic_benchmarks.yml b/.github/workflows/periodic_benchmarks.yml index 86cd1991ab..0a7a7d9c47 100644 --- a/.github/workflows/periodic_benchmarks.yml +++ b/.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: | diff --git a/CHANGELOG.md b/CHANGELOG.md index 29af0fbe72..d03ef3f88d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,13 +2,19 @@ ## 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 +- Fixed a bug where `IDAKLUSolver` errors were not raised correctly. ([#5291](https://github.com/pybamm-team/PyBaMM/pull/5291)) +- 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 - Allow setting initial conditions from `y_slices` of a `Solution` object. ([#5257](https://github.com/pybamm-team/PyBaMM/pull/5257)) - Added docstring to `FuzzyDict.copy` explaining its return value and behavior. ([#5242](https://github.com/pybamm-team/PyBaMM/pull/5242)) +- Porosity and active material fractions are now `FunctionParameters` of y and z, as well as x ([#5214](https://github.com/pybamm-team/PyBaMM/pull/5214)) ## Bug fixes @@ -19,6 +25,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)) diff --git a/docs/source/examples/notebooks/creating_models/4-comparing-full-and-reduced-order-models.ipynb b/docs/source/examples/notebooks/creating_models/4-comparing-full-and-reduced-order-models.ipynb index 071fd54f95..19780371c0 100644 --- a/docs/source/examples/notebooks/creating_models/4-comparing-full-and-reduced-order-models.ipynb +++ b/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", diff --git a/docs/source/examples/notebooks/models/graded-electrodes.ipynb b/docs/source/examples/notebooks/models/graded-electrodes.ipynb index 272248d7e8..7954c0fbd6 100644 --- a/docs/source/examples/notebooks/models/graded-electrodes.ipynb +++ b/docs/source/examples/notebooks/models/graded-electrodes.ipynb @@ -22,13 +22,6 @@ "text": [ "Note: you may need to restart the kernel to use updated packages.\n" ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.\n" - ] } ], "source": [ @@ -59,7 +52,7 @@ "source": [ "We will vary the porosity in both electrodes and we will try three different scenarios: constant porosity, one where lower porosity occurs near the separator and one where lower porosity occurs near the current collector. All other parameters are kept constant. The varying porosity is defined to be linear centered around the default value and with a variation of $\\pm$ 10%.\n", "\n", - "We define the varying porosities and store them in a list so we can loop over when solving the model." + "We define the varying porosities and store them in a list so we can loop over when solving the model. Note that the porosities must be specified as functions of x, y, and z, even if y and z are not present in the model and not referred to at any other point, as in this example." ] }, { @@ -77,14 +70,14 @@ "\n", "eps_ns = [\n", " eps_n_0,\n", - " lambda x: eps_n_0 * (1.1 - 0.2 * (x / L_n)),\n", - " lambda x: eps_n_0 * (0.9 + 0.2 * (x / L_n)),\n", + " lambda x, y, z: eps_n_0 * (1.1 - 0.2 * (x / L_n)),\n", + " lambda x, y, z: eps_n_0 * (0.9 + 0.2 * (x / L_n)),\n", "]\n", "\n", "eps_ps = [\n", " eps_p_0,\n", - " lambda x: eps_p_0 * (0.9 - 0.2 / L_p * (L_n + L_s) + 0.2 * (x / L_p)),\n", - " lambda x: eps_p_0 * (1.1 + 0.2 / L_p * (L_n + L_s) - 0.2 * (x / L_p)),\n", + " lambda x, y, z: eps_p_0 * (0.9 - 0.2 / L_p * (L_n + L_s) + 0.2 * (x / L_p)),\n", + " lambda x, y, z: eps_p_0 * (1.1 + 0.2 / L_p * (L_n + L_s) - 0.2 * (x / L_p)),\n", "]" ] }, @@ -132,12 +125,12 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "99f847ca09da40cba550dd02dd8281a2", + "model_id": "9198b6c974714d07b0a68f1cbc42c988", "version_major": 2, "version_minor": 0 }, "text/plain": [ - "interactive(children=(FloatSlider(value=0.0, description='t', max=673.9136958613059, step=6.7391369586130585),…" + "interactive(children=(FloatSlider(value=0.0, description='t', max=673.9107055793413, step=6.739107055793413), …" ] }, "metadata": {}, @@ -146,7 +139,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 5, @@ -176,18 +169,18 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "829b68c6b3e04e0ebe5537daabec2278", + "model_id": "666bc13bf8cb4ee9b1a110cc828f025b", "version_major": 2, "version_minor": 0 }, "text/plain": [ - "interactive(children=(FloatSlider(value=0.0, description='t', max=673.9136958613059, step=6.7391369586130585),…" + "interactive(children=(FloatSlider(value=0.0, description='t', max=673.9107055793413, step=6.739107055793413), …" ] }, "metadata": {}, @@ -196,10 +189,10 @@ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 8, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } @@ -218,7 +211,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -232,10 +225,7 @@ "[5] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n", "[6] Alan C. Hindmarsh. The PVODE and IDA algorithms. Technical Report, Lawrence Livermore National Lab., CA (US), 2000. doi:10.2172/802599.\n", "[7] Alan C. Hindmarsh, Peter N. Brown, Keith E. Grant, Steven L. Lee, Radu Serban, Dan E. Shumaker, and Carol S. Woodward. SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Transactions on Mathematical Software (TOMS), 31(3):363–396, 2005. doi:10.1145/1089014.1089020.\n", - "[8] Peyman Mohtat, Suhak Lee, Jason B Siegel, and Anna G Stefanopoulou. Towards better estimability of electrode-specific state of health: decoding the cell expansion. Journal of Power Sources, 427:101–111, 2019.\n", - "[9] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n", - "[10] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.\n", - "[11] Andrew Weng, Jason B Siegel, and Anna Stefanopoulou. Differential voltage analysis for battery manufacturing process control. arXiv preprint arXiv:2303.07088, 2023.\n", + "[8] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n", "\n" ] } @@ -254,7 +244,7 @@ ], "metadata": { "kernelspec": { - "display_name": "venv", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -268,9 +258,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.6" + "version": "3.10.12" } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } diff --git a/examples/scripts/3d_examples/dried_out_pouch.py b/examples/scripts/3d_examples/dried_out_pouch.py new file mode 100644 index 0000000000..f795e1676b --- /dev/null +++ b/examples/scripts/3d_examples/dried_out_pouch.py @@ -0,0 +1,48 @@ +import numpy as np + +import pybamm + +model = pybamm.lithium_ion.DFN( + {"current collector": "potential pair", "dimensionality": 2} +) +param = pybamm.ParameterValues("Ecker2015") +Ly = param["Electrode width [m]"] +Lz = param["Electrode height [m]"] + + +def _sigmoid(arg): + return (1 + np.tanh(arg)) / 2 + + +def _top_hat(arg, a, b, k=500): + return _sigmoid(k * (arg - a)) * _sigmoid(k * (b - arg)) + + +# Simulate drying out of the negative electrode edges by reducing porosity +def eps_n(x, y, z): + return 0.329 * ( + _top_hat(arg=y, a=Ly * 0.02, b=Ly * 0.98) + * _top_hat(arg=z, a=Lz * 0.02, b=Lz * 0.98) + ) + + +param.update({"Negative electrode porosity": eps_n}) +var_pts = {"x_n": 8, "x_s": 8, "x_p": 8, "r_n": 8, "r_p": 8, "y": 24, "z": 24} +exp = pybamm.Experiment( + [ + "Discharge at 1C until 2.7 V", + "Charge at 1C until 4.2 V", + "Hold at 4.2 V until C/20", + ] +) +sim = pybamm.Simulation(model, var_pts=var_pts, parameter_values=param, experiment=exp) +sol = sim.solve() +output_variables = [ + "X-averaged negative electrode porosity", + "X-averaged negative particle surface stoichiometry", + "X-averaged negative electrode surface potential difference [V]", + "Current collector current density [A.m-2]", + "Current [A]", + "Voltage [V]", +] +plot = sol.plot(output_variables, variable_limits="tight", shading="auto") diff --git a/pyproject.toml b/pyproject.toml index 5f214ef642..364b31a52e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ classifiers = [ "Topic :: Scientific/Engineering", ] dependencies = [ - "pybammsolvers>=0.3.0,<0.4.0", + "pybammsolvers>=0.3.3,<0.4.0", "black", "numpy", "scipy>=1.11.4", diff --git a/src/pybamm/expression_tree/operations/serialise.py b/src/pybamm/expression_tree/operations/serialise.py index 6b733e0fad..3cc8f3e98f 100644 --- a/src/pybamm/expression_tree/operations/serialise.py +++ b/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"]) diff --git a/src/pybamm/meshes/meshes.py b/src/pybamm/meshes/meshes.py index f717b0a347..ddc295d57b 100644 --- a/src/pybamm/meshes/meshes.py +++ b/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 # 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. diff --git a/src/pybamm/models/base_model.py b/src/pybamm/models/base_model.py index 9fd7483385..67096912cf 100644 --- a/src/pybamm/models/base_model.py +++ b/src/pybamm/models/base_model.py @@ -906,7 +906,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 +923,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 +1088,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 diff --git a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py b/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py index a9322991b2..b9abb4261f 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py +++ b/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 diff --git a/src/pybamm/models/submodels/interface/kinetics/butler_volmer.py b/src/pybamm/models/submodels/interface/kinetics/butler_volmer.py index c3e9ac0fc0..84cee586a3 100644 --- a/src/pybamm/models/submodels/interface/kinetics/butler_volmer.py +++ b/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 ---------- diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/base_hysteresis_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/base_hysteresis_ocp.py index 934ed5e6cf..cdfc47ac5e 100644 --- a/src/pybamm/models/submodels/interface/open_circuit_potential/base_hysteresis_ocp.py +++ b/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") diff --git a/src/pybamm/parameters/lithium_ion_parameters.py b/src/pybamm/parameters/lithium_ion_parameters.py index f64e7a71c5..85f573c7a0 100644 --- a/src/pybamm/parameters/lithium_ion_parameters.py +++ b/src/pybamm/parameters/lithium_ion_parameters.py @@ -233,8 +233,15 @@ def _set_parameters(self): if domain == "separator": x = pybamm.standard_spatial_vars.x_s + y = pybamm.PrimaryBroadcast(pybamm.standard_spatial_vars.y, "separator") + z = pybamm.PrimaryBroadcast(pybamm.standard_spatial_vars.z, "separator") self.epsilon_init = pybamm.FunctionParameter( - "Separator porosity", {"Through-cell distance (x) [m]": x} + "Separator porosity", + { + "Through-cell distance (x) [m]": x, + "Horizontal distance (y) [m]": y, + "Vertical distance (z) [m]": z, + }, ) self.epsilon_inactive = 1 - self.epsilon_init return @@ -248,6 +255,12 @@ def _set_parameters(self): auxiliary_domains={"secondary": "current collector"}, coord_sys="cartesian", ) + y = pybamm.PrimaryBroadcast( + pybamm.standard_spatial_vars.y, f"{domain} electrode" + ) + z = pybamm.PrimaryBroadcast( + pybamm.standard_spatial_vars.z, f"{domain} electrode" + ) # Macroscale geometry self.L_cc = self.geo.L_cc @@ -267,7 +280,12 @@ def _set_parameters(self): ) if main.options.electrode_types[domain] == "porous": self.epsilon_init = pybamm.FunctionParameter( - f"{Domain} electrode porosity", {"Through-cell distance (x) [m]": x} + f"{Domain} electrode porosity", + { + "Through-cell distance (x) [m]": x, + "Horizontal distance (y) [m]": y, + "Vertical distance (z) [m]": z, + }, ) epsilon_s_tot = sum(phase.epsilon_s for phase in self.phase_params.values()) self.epsilon_inactive = 1 - self.epsilon_init - epsilon_s_tot @@ -414,6 +432,12 @@ def _set_parameters(self): auxiliary_domains={"secondary": "current collector"}, coord_sys="cartesian", ) + y = pybamm.PrimaryBroadcast( + pybamm.standard_spatial_vars.y, f"{domain} electrode" + ) + z = pybamm.PrimaryBroadcast( + pybamm.standard_spatial_vars.z, f"{domain} electrode" + ) r = pybamm.SpatialVariable( f"r_{domain[0]}", domain=[f"{domain} {self.phase_name}particle"], @@ -438,7 +462,11 @@ def _set_parameters(self): # Particle properties self.epsilon_s = pybamm.FunctionParameter( f"{pref}{Domain} electrode active material volume fraction", - {"Through-cell distance (x) [m]": x}, + { + "Through-cell distance (x) [m]": x, + "Horizontal distance (y) [m]": y, + "Vertical distance (z) [m]": z, + }, ) self.epsilon_s_av = pybamm.xyz_average(self.epsilon_s) self.c_max = pybamm.Parameter( diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py index 2bec326a93..845be2b0bf 100644 --- a/src/pybamm/parameters/parameter_values.py +++ b/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 ) } diff --git a/src/pybamm/simulation.py b/src/pybamm/simulation.py index 8f3fcaaad1..93d01c5af4 100644 --- a/src/pybamm/simulation.py +++ b/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, @@ -778,7 +822,7 @@ def solve( feasible = False # If none of the cycles worked, raise an error if cycle_num == 1 and step_num == 1: - raise error + raise error from error # Otherwise, just stop this cycle break diff --git a/src/pybamm/solvers/base_solver.py b/src/pybamm/solvers/base_solver.py index 92c7381f43..a1db0c2953 100644 --- a/src/pybamm/solvers/base_solver.py +++ b/src/pybamm/solvers/base_solver.py @@ -1381,7 +1381,9 @@ def step( else: _, concatenated_initial_conditions = model.set_initial_conditions_from( - old_solution, return_type="ics" + old_solution, + inputs=model_inputs, + return_type="ics", ) model.y0 = concatenated_initial_conditions.evaluate(0, inputs=model_inputs) if using_sensitivities: diff --git a/src/pybamm/solvers/idaklu_solver.py b/src/pybamm/solvers/idaklu_solver.py index 4a0d9edb07..147bc0314e 100644 --- a/src/pybamm/solvers/idaklu_solver.py +++ b/src/pybamm/solvers/idaklu_solver.py @@ -76,6 +76,8 @@ class IDAKLUSolver(pybamm.BaseSolver): "increment_factor": 1.0, # Enable or disable linear solution scaling "linear_solution_scaling": True, + # Silence Sundials errors during solve + "silence_sundials_errors": False, ## Main solver # Maximum order of the linear multistep method "max_order_bdf": 5, @@ -176,6 +178,7 @@ def __init__( "epsilon_linear_tolerance": 0.05, "increment_factor": 1.0, "linear_solution_scaling": True, + "silence_sundials_errors": False, "max_order_bdf": 5, "max_num_steps": 100000, "dt_init": 0.0, @@ -681,13 +684,17 @@ def _integrate( atol = self._check_atol_type(atol, y0full.size) timer = pybamm.Timer() - solns = self._setup["solver"].solve( - t_eval, - t_interp, - y0full, - ydot0full, - inputs, - ) + try: + solns = self._setup["solver"].solve( + t_eval, + t_interp, + y0full, + ydot0full, + inputs, + ) + except ValueError as e: + # Return from None to replace the C++ runtime error + raise pybamm.SolverError(str(e)) from None integration_time = timer.time() return [ @@ -734,14 +741,15 @@ def _post_process_solution(self, sol, model, integration_time, inputs_dict): termination = "final time" elif sol.flag < 0: termination = "failure" + msg = idaklu.sundials_error_message(sol.flag) match self._on_failure: case "warn": warnings.warn( - f"FAILURE {self._solver_flag(sol.flag)}, returning a partial solution.", + msg + ", returning a partial solution.", stacklevel=2, ) case "raise": - raise pybamm.SolverError(f"FAILURE {self._solver_flag(sol.flag)}") + raise pybamm.SolverError(msg) if sol.yp.size > 0: yp = sol.yp.reshape((number_of_timesteps, number_of_states)).T @@ -1002,40 +1010,3 @@ def jaxify( t_interp=t_interp, ) return obj - - @staticmethod - def _solver_flag(flag): - flags = { - 99: "IDA_WARNING: IDASolve succeeded but an unusual situation occurred.", - 2: "IDA_ROOT_RETURN: IDASolve succeeded and found one or more roots.", - 1: "IDA_TSTOP_RETURN: IDASolve succeeded by reaching the specified stopping point.", - 0: "IDA_SUCCESS: Successful function return.", - -1: "IDA_TOO_MUCH_WORK: The solver took mxstep internal steps but could not reach tout.", - -2: "IDA_TOO_MUCH_ACC: The solver could not satisfy the accuracy demanded by the user for some internal step.", - -3: "IDA_ERR_FAIL: Error test failures occurred too many times during one internal time step or minimum step size was reached.", - -4: "IDA_CONV_FAIL: Convergence test failures occurred too many times during one internal time step or minimum step size was reached.", - -5: "IDA_LINIT_FAIL: The linear solver's initialization function failed.", - -6: "IDA_LSETUP_FAIL: The linear solver's setup function failed in an unrecoverable manner.", - -7: "IDA_LSOLVE_FAIL: The linear solver's solve function failed in an unrecoverable manner.", - -8: "IDA_RES_FAIL: The user-provided residual function failed in an unrecoverable manner.", - -9: "IDA_REP_RES_FAIL: The user-provided residual function repeatedly returned a recoverable error flag, but the solver was unable to recover.", - -10: "IDA_RTFUNC_FAIL: The rootfinding function failed in an unrecoverable manner.", - -11: "IDA_CONSTR_FAIL: The inequality constraints were violated and the solver was unable to recover.", - -12: "IDA_FIRST_RES_FAIL: The user-provided residual function failed recoverably on the first call.", - -13: "IDA_LINESEARCH_FAIL: The line search failed.", - -14: "IDA_NO_RECOVERY: The residual function, linear solver setup function, or linear solver solve function had a recoverable failure, but IDACalcIC could not recover.", - -15: "IDA_NLS_INIT_FAIL: The nonlinear solver's init routine failed.", - -16: "IDA_NLS_SETUP_FAIL: The nonlinear solver's setup routine failed.", - -20: "IDA_MEM_NULL: The ida mem argument was NULL.", - -21: "IDA_MEM_FAIL: A memory allocation failed.", - -22: "IDA_ILL_INPUT: One of the function inputs is illegal.", - -23: "IDA_NO_MALLOC: The ida memory was not allocated by a call to IDAInit.", - -24: "IDA_BAD_EWT: Zero value of some error weight component.", - -25: "IDA_BAD_K: The k-th derivative is not available.", - -26: "IDA_BAD_T: The time t is outside the last step taken.", - -27: "IDA_BAD_DKY: The vector argument where derivative should be stored is NULL.", - } - - flag_unknown = "Unknown IDA flag." - - return flags.get(flag, flag_unknown) diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py index 2b62b83db3..3c5099dde7 100644 --- a/tests/unit/test_experiments/test_simulation_with_experiment.py +++ b/tests/unit/test_experiments/test_simulation_with_experiment.py @@ -1,3 +1,4 @@ +import logging import os from datetime import datetime @@ -1022,3 +1023,125 @@ def neg_stoich_cutoff(variables): neg_stoich = sol["Negative electrode stoichiometry"].data assert neg_stoich[-1] == pytest.approx(0.5, abs=0.0001) + + def test_simulation_changing_capacity_crate_steps(self): + """Test that C-rate steps are correctly updated when capacity changes""" + model = pybamm.lithium_ion.SPM() + experiment = pybamm.Experiment( + [ + ( + "Discharge at C/5 for 20 minutes", + "Discharge at C/2 for 20 minutes", + "Discharge at 1C for 20 minutes", + ) + ] + ) + param = pybamm.ParameterValues("Chen2020") + sim = pybamm.Simulation(model, experiment=experiment, parameter_values=param) + + # First solve + sol1 = sim.solve(calc_esoh=False) + original_capacity = param["Nominal cell capacity [A.h]"] + + # Check that C-rates correspond to expected currents + I_C5_1 = np.abs(sol1.cycles[0].steps[0]["Current [A]"].data).mean() + I_C2_1 = np.abs(sol1.cycles[0].steps[1]["Current [A]"].data).mean() + I_1C_1 = np.abs(sol1.cycles[0].steps[2]["Current [A]"].data).mean() + + np.testing.assert_allclose(I_C5_1, original_capacity / 5, rtol=1e-2) + np.testing.assert_allclose(I_C2_1, original_capacity / 2, rtol=1e-2) + np.testing.assert_allclose(I_1C_1, original_capacity, rtol=1e-2) + + # Update capacity + new_capacity = 0.9 * original_capacity + sim._parameter_values.update({"Nominal cell capacity [A.h]": new_capacity}) + + # Second solve with updated capacity + sol2 = sim.solve(calc_esoh=False) + + # Check that C-rates now correspond to updated currents + I_C5_2 = np.abs(sol2.cycles[0].steps[0]["Current [A]"].data).mean() + I_C2_2 = np.abs(sol2.cycles[0].steps[1]["Current [A]"].data).mean() + I_1C_2 = np.abs(sol2.cycles[0].steps[2]["Current [A]"].data).mean() + + np.testing.assert_allclose(I_C5_2, new_capacity / 5, rtol=1e-2) + np.testing.assert_allclose(I_C2_2, new_capacity / 2, rtol=1e-2) + np.testing.assert_allclose(I_1C_2, new_capacity, rtol=1e-2) + + # Verify all currents scaled proportionally + np.testing.assert_allclose(I_C5_2 / I_C5_1, 0.9, rtol=1e-2) + np.testing.assert_allclose(I_C2_2 / I_C2_1, 0.9, rtol=1e-2) + np.testing.assert_allclose(I_1C_2 / I_1C_1, 0.9, rtol=1e-2) + + def test_simulation_multiple_cycles_with_capacity_change(self): + """Test capacity changes across multiple experiment cycles""" + model = pybamm.lithium_ion.SPM() + experiment = pybamm.Experiment( + [("Discharge at 1C for 5 minutes", "Charge at 1C for 5 minutes")] * 2 + ) + param = pybamm.ParameterValues("Chen2020") + sim = pybamm.Simulation(model, experiment=experiment, parameter_values=param) + + # First solve + sol1 = sim.solve(calc_esoh=False) + original_capacity = param["Nominal cell capacity [A.h]"] + + # Get discharge currents for both cycles + I_discharge_cycle1 = np.abs(sol1.cycles[0].steps[0]["Current [A]"].data).mean() + I_discharge_cycle2 = np.abs(sol1.cycles[1].steps[0]["Current [A]"].data).mean() + + # Both cycles should use the same capacity initially + np.testing.assert_allclose(I_discharge_cycle1, original_capacity, rtol=1e-2) + np.testing.assert_allclose(I_discharge_cycle2, original_capacity, rtol=1e-2) + + # Update capacity between cycles + new_capacity = 0.85 * original_capacity + sim._parameter_values.update({"Nominal cell capacity [A.h]": new_capacity}) + + # Solve again + sol2 = sim.solve(calc_esoh=False) + + # All cycles in the new solution should use updated capacity + I_discharge_cycle1_new = np.abs( + sol2.cycles[0].steps[0]["Current [A]"].data + ).mean() + I_discharge_cycle2_new = np.abs( + sol2.cycles[1].steps[0]["Current [A]"].data + ).mean() + + np.testing.assert_allclose(I_discharge_cycle1_new, new_capacity, rtol=1e-2) + np.testing.assert_allclose(I_discharge_cycle2_new, new_capacity, rtol=1e-2) + + def test_simulation_logging_with_capacity_change(self, caplog): + """Test that capacity changes are logged appropriately""" + model = pybamm.lithium_ion.SPM() + experiment = pybamm.Experiment([("Discharge at 1C for 10 minutes",)]) + param = pybamm.ParameterValues("Chen2020") + sim = pybamm.Simulation(model, experiment=experiment, parameter_values=param) + + # First solve + sim.solve(calc_esoh=False) + original_capacity = param["Nominal cell capacity [A.h]"] + + # Update capacity + new_capacity = 0.75 * original_capacity + sim._parameter_values.update({"Nominal cell capacity [A.h]": new_capacity}) + + # Set logging level to capture INFO messages + original_log_level = pybamm.logger.level + pybamm.set_logging_level("INFO") + + try: + # Second solve should log capacity change + with caplog.at_level(logging.INFO, logger="pybamm.logger"): + sim.solve(calc_esoh=False) + + # Check that a log message about capacity change was recorded + log_messages = [record.message for record in caplog.records] + capacity_change_logged = any( + "Nominal capacity changed" in msg for msg in log_messages + ) + assert capacity_change_logged + finally: + # Restore original logging level + pybamm.logger.setLevel(original_log_level) diff --git a/tests/unit/test_meshes/test_meshes.py b/tests/unit/test_meshes/test_meshes.py index 8c26a0900f..15ef8317ba 100644 --- a/tests/unit/test_meshes/test_meshes.py +++ b/tests/unit/test_meshes/test_meshes.py @@ -584,6 +584,37 @@ def test_to_json(self): assert mesh_json == expected_json + def test_compute_var_pts_from_thicknesses_cell_size(self): + from pybamm.meshes.meshes import compute_var_pts_from_thicknesses + + electrode_thicknesses = { + "negative electrode": 100e-6, + "separator": 25e-6, + "positive electrode": 100e-6, + } + + cell_size = 5e-6 # 5 micrometres per cell + var_pts = compute_var_pts_from_thicknesses(electrode_thicknesses, cell_size) + + assert isinstance(var_pts, dict) + assert all(isinstance(v, dict) for v in var_pts.values()) + assert var_pts["negative electrode"]["x_n"] == 20 + assert var_pts["separator"]["x_s"] == 5 + assert var_pts["positive electrode"]["x_p"] == 20 + + def test_compute_var_pts_from_thicknesses_invalid_thickness_type(self): + from pybamm.meshes.meshes import compute_var_pts_from_thicknesses + + with pytest.raises(TypeError): + compute_var_pts_from_thicknesses(["not", "a", "dict"], 1e-6) + + def test_compute_var_pts_from_thicknesses_invalid_grid_size(self): + from pybamm.meshes.meshes import compute_var_pts_from_thicknesses + + electrode_thicknesses = {"negative electrode": 100e-6} + with pytest.raises(ValueError): + compute_var_pts_from_thicknesses(electrode_thicknesses, -1e-6) + class TestMeshGenerator: def test_init_name(self): diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py index 769705cf51..50b4c7065d 100644 --- a/tests/unit/test_parameters/test_parameter_values.py +++ b/tests/unit/test_parameters/test_parameter_values.py @@ -1298,6 +1298,31 @@ def test_to_json_with_filename(self): finally: os.remove(temp_path) + def test_roundtrip_with_keyword_args(self): + def func_no_kwargs(x): + return 2 * x + + def func_with_kwargs(x, y=1): + return 2 * x + + x = pybamm.Scalar(2) + func_param = pybamm.FunctionParameter("func", {"x": x}) + + parameter_values = pybamm.ParameterValues({"func": func_no_kwargs}) + assert parameter_values.evaluate(func_param) == 4.0 + + serialized = parameter_values.to_json() + parameter_values_loaded = pybamm.ParameterValues.from_json(serialized) + assert parameter_values_loaded.evaluate(func_param) == 4.0 + + parameter_values = pybamm.ParameterValues({"func": func_with_kwargs}) + assert parameter_values.evaluate(func_param) == 4.0 + + serialized = parameter_values.to_json() + parameter_values_loaded = pybamm.ParameterValues.from_json(serialized) + + assert parameter_values_loaded.evaluate(func_param) == 4.0 + def test_convert_symbols_in_dict_with_interpolator(self): """Test convert_symbols_in_dict with interpolator (covers lines 1154-1170).""" import numpy as np diff --git a/tests/unit/test_serialisation/test_serialisation.py b/tests/unit/test_serialisation/test_serialisation.py index d6ddac10e4..4dafcc5641 100644 --- a/tests/unit/test_serialisation/test_serialisation.py +++ b/tests/unit/test_serialisation/test_serialisation.py @@ -617,6 +617,14 @@ def test_serialise_time(self): t2 = convert_symbol_from_json(j) assert isinstance(t2, pybamm.Time) + def test_serialise_input_parameter(self): + """Test InputParameter serialization and deserialization.""" + ip = pybamm.InputParameter("test_param") + j = convert_symbol_to_json(ip) + ip_restored = convert_symbol_from_json(j) + assert isinstance(ip_restored, pybamm.InputParameter) + assert ip_restored.name == "test_param" + def test_convert_symbol_to_json_with_number_and_list(self): for val in (0, 3.14, -7, True): out = convert_symbol_to_json(val) diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py index a86e18b21d..5b55e5e3d8 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -583,7 +583,7 @@ def test_failures(self): solver = pybamm.IDAKLUSolver() t_eval = [0, 3] - with pytest.raises(ValueError): + with pytest.raises(pybamm.SolverError): solver.solve(model, t_eval) def test_dae_solver_algebraic_model(self): @@ -779,7 +779,7 @@ def test_solver_options(self): options = {option: options_fail[option]} solver = pybamm.IDAKLUSolver(options=options) - with pytest.raises(ValueError): + with pytest.raises(pybamm.SolverError): solver.solve(model, t_eval) def test_with_output_variables(self): @@ -1487,7 +1487,7 @@ def test_on_failure_option(self): model, t_eval=t_eval, t_interp=t_interp, inputs=input_parameters ) assert len(w) > 0 - assert "FAILURE" in str(w[0].message) + assert "_FAIL" in str(w[0].message) def test_no_progress_early_termination(self): # SPM at rest