diff --git a/CHANGELOG.md b/CHANGELOG.md index a17810776c..b8d3cd75e1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- `ElectrodeSOHComposite` now solves each stoichiometry by bracketed rootfind, converging from any target. ([#5730](https://github.com/pybamm-team/PyBaMM/pull/5730)) - Added unstructured mesh support (`UnstructuredSubMesh`, generators, and interface coupling) for arbitrary 2D/3D domains. Hexahedra must have planar faces (warped hexes raise a `GeometryError`), and `UserSuppliedUnstructuredMesh` accepts tetrahedral, triangular, and quadrilateral cells only. ([#5687](https://github.com/pybamm-team/PyBaMM/pull/5687)) - Generalised `VectorField` to N components and added `Component`/`Norm` operators for multi-dimensional vector fields. ([#5686](https://github.com/pybamm-team/PyBaMM/pull/5686)) - Removed the left sidebar from the documentation home page for a cleaner landing experience. ([#5699](https://github.com/pybamm-team/PyBaMM/pull/5699)) diff --git a/packages/pybamm/src/pybamm/__init__.py b/packages/pybamm/src/pybamm/__init__.py index 740fca3e26..405783a81f 100644 --- a/packages/pybamm/src/pybamm/__init__.py +++ b/packages/pybamm/src/pybamm/__init__.py @@ -11,6 +11,7 @@ from .util import ( get_parameters_filepath, has_jax, + is_windows, raise_jax_not_found, import_optional_dependency, ) diff --git a/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_composite.py b/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_composite.py index 2eb9d7fa9b..51c5adce3c 100644 --- a/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_composite.py +++ b/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_composite.py @@ -6,6 +6,9 @@ import warnings from typing import Any +import casadi +import numpy as np + import pybamm from .electrode_soh import _ElectrodeSOH, get_esoh_default_solver @@ -15,6 +18,14 @@ get_lithiation_delithiation, ) +# `LithiumIonParameters.U` clips its stoichiometry and adds an asymptote, so every OCP +# diverges outside [0, 1] and both brackets straddle a root. The inner one is wider so +# an inversion can always reach the potential the outer solve asks it for. +_STOICH_LO, _STOICH_HI = -10.0, 10.0 +_BRACKET_LO, _BRACKET_HI = -5.0, 5.0 + +_ABSTOL, _MAX_ITER = 1e-14, 200 + def _get_primary_only_options( options: dict | pybamm.BatteryModelOptions | None, @@ -74,6 +85,19 @@ def _get_primary_only_options( return options_dict +def _initialization_method(initial_value: float | str) -> str: + """``"voltage"`` for a string ending in ``"V"``, ``"SOC"`` for a float.""" + if isinstance(initial_value, str) and initial_value.endswith("V"): + return "voltage" + if isinstance(initial_value, float): + return "SOC" + raise ValueError( + "Invalid initial value. Expected a float between 0 and 1 " + "(for SOC) or a string ending in 'V' (for voltage), got " + f"{initial_value!r} of type {type(initial_value).__name__}" + ) + + def _get_stoich_variables(options): """Create stoichiometry variables for composite electrodes.""" variables = { @@ -300,234 +324,195 @@ def __init__( pybamm.citations.register("Mohtat2019") super().__init__(name) param = pybamm.LithiumIonParameters(options) + neg_composite = check_if_composite(options, "negative") + pos_composite = check_if_composite(options, "positive") Q_Li = pybamm.InputParameter("Q_Li") - is_negative_composite = check_if_composite(options, "negative") - is_positive_composite = check_if_composite(options, "positive") - variables = _get_stoich_variables(options) - x_100_1 = variables["x_100_1"] - y_100_1 = variables["y_100_1"] - x_0_1 = variables["x_0_1"] - y_0_1 = variables["y_0_1"] - V_max = param.ocp_soc_100 - V_min = param.ocp_soc_0 - if is_negative_composite: - x_100_2 = variables["x_100_2"] - x_0_2 = variables["x_0_2"] - self.algebraic[x_100_2] = param.n.sec.U( - x_100_2, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("100", "negative", options, "secondary"), - "negative", - options, - phase="secondary", - ), - ) - param.n.prim.U( - x_100_1, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("100", "negative", options, "primary"), - "negative", - options, - phase="primary", - ), - ) - self.algebraic[x_0_2] = param.n.sec.U( - x_0_2, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("0", "negative", options, "secondary"), - "negative", - options, - phase="secondary", - ), - ) - param.n.prim.U( - x_0_1, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("0", "negative", options, "primary"), - "negative", - options, - phase="primary", - ), - ) - if is_positive_composite: - y_100_2 = variables["y_100_2"] - y_0_2 = variables["y_0_2"] - self.algebraic[y_100_2] = param.p.sec.U( - y_100_2, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("100", "positive", options, "secondary"), - "positive", - options, - phase="secondary", - ), - ) - param.p.prim.U( - y_100_1, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("100", "positive", options, "primary"), - "positive", - options, - phase="primary", - ), - ) - self.algebraic[y_0_2] = param.p.prim.U( - y_0_1, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("0", "positive", options, "primary"), - "positive", - options, - phase="primary", - ), - ) - param.p.sec.U( - y_0_2, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("0", "positive", options, "secondary"), - "positive", - options, - phase="secondary", - ), - ) - self.algebraic[x_100_1] = ( - param.p.prim.U( - y_100_1, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("100", "positive", options, "primary"), - "positive", - options, - phase="primary", - ), + Q_n = [pybamm.InputParameter("Q_n_1")] + Q_p = [pybamm.InputParameter("Q_p_1")] + if neg_composite: + Q_n.append(pybamm.InputParameter("Q_n_2")) + if pos_composite: + Q_p.append(pybamm.InputParameter("Q_p_2")) + + T_ref = param.T_ref + T_init = param.T_init if initialization_method == "voltage" else T_ref + + def ocps(electrode, soc): + """(open-circuit potential, branch, name) per phase, at an SOC level.""" + side = param.n if electrode == "negative" else param.p + phases = ["primary"] + if check_if_composite(options, electrode): + phases.append("secondary") + out = [] + for phase in phases: + U = (side.prim if phase == "primary" else side.sec).U + if soc == "init": + branch = get_lithiation_delithiation( + direction, electrode, options, phase=phase + ) + else: + branch = get_lithiation_delithiation( + get_equilibrium_direction(soc, electrode, options, phase), + electrode, + options, + phase=phase, + ) + out.append((U, branch, f"{electrode} {phase} {soc}")) + return out + + def invert(U, branch, target, T, name): + """Stoichiometry at a given potential, by Brent.""" + sto = pybamm._BrentUnknown(name) + return pybamm._Brent( + U(sto, T, branch) - target, + sto, + (_STOICH_LO, _STOICH_HI), + abstol=_ABSTOL, + max_iter=_MAX_ITER, ) - - param.n.prim.U( - x_100_1, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("100", "negative", options, "primary"), - "negative", - options, - phase="primary", - ), + + def bracket(pairs, offsets, T): + """Potentials for which every inversion in `pairs` brackets.""" + los = [ + U(pybamm.Scalar(_BRACKET_HI), T, b) - off + for (U, b, _), off in zip(pairs, offsets, strict=True) + ] + his = [ + U(pybamm.Scalar(_BRACKET_LO), T, b) - off + for (U, b, _), off in zip(pairs, offsets, strict=True) + ] + lo, hi = los[0], his[0] + for value in los[1:]: + lo = pybamm.maximum(lo, value) + for value in his[1:]: + hi = pybamm.minimum(hi, value) + return lo, hi + + def limits(soc, V_target, closure): + """A limit state: one unknown, the shared negative potential.""" + neg, pos = ocps("negative", soc), ocps("positive", soc) + + # built twice: once against the unknown potential, once against the solved + def states(U_n, suffix): + return ( + [invert(U, b, U_n, T_ref, f"{n} {suffix}") for U, b, n in neg], + [ + invert(U, b, U_n + V_target, T_ref, f"{n} {suffix}") + for U, b, n in pos + ], + ) + + U_n = pybamm._BrentUnknown(f"negative potential {soc}") + lo, hi = bracket(neg + pos, [0] * len(neg) + [V_target] * len(pos), T_ref) + solved = pybamm._Brent( + closure(*states(U_n, "guess")), + U_n, + (lo, hi), + abstol=_ABSTOL, + max_iter=_MAX_ITER, ) - - V_max + return states(solved, "solved") + + def total_lithium(x, y): + return sum(q * s for q, s in zip(Q_n + Q_p, x + y, strict=True)) + + x_100, y_100 = limits( + "100", param.ocp_soc_100, lambda x, y: Q_Li - total_lithium(x, y) ) - self.algebraic[x_0_1] = ( - param.p.prim.U( - y_0_1, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("0", "positive", options, "primary"), - "positive", - options, - phase="primary", - ), - ) - - param.n.prim.U( - x_0_1, - param.T_ref, - get_lithiation_delithiation( - get_equilibrium_direction("0", "negative", options, "primary"), - "negative", - options, - phase="primary", - ), - ) - - V_min + x_0, y_0 = limits( + "0", + param.ocp_soc_0, + lambda x, y: ( + -sum(q * (a - b) for q, a, b in zip(Q_p, y_100, y, strict=True)) + - sum(q * (a - b) for q, a, b in zip(Q_n, x_100, x, strict=True)) + ), ) - self.algebraic[y_0_1] = _get_electrode_capacity_equation( - options, "positive" - ) - _get_electrode_capacity_equation(options, "negative") - self.algebraic[y_100_1] = Q_Li - _get_cyclable_lithium_equation(options) - x_init_1 = variables["x_init_1"] - y_init_1 = variables["y_init_1"] - if initialization_method == "voltage": - V_init = pybamm.InputParameter("V_init") - self.algebraic[x_init_1] = ( - param.p.prim.U( - y_init_1, - param.T_init, - get_lithiation_delithiation( - direction, "positive", options, phase="primary" - ), - ) - - param.n.prim.U( - x_init_1, - param.T_init, - get_lithiation_delithiation( - direction, "negative", options, phase="primary" - ), + neg_init, pos_init = ocps("negative", "init"), ocps("positive", "init") + + def init_states(x_init_1, suffix): + U_n = neg_init[0][0](x_init_1, T_init, neg_init[0][1]) + x = [x_init_1] + [ + invert(U, b, U_n, T_init, f"{n} {suffix}") for U, b, n in neg_init[1:] + ] + if initialization_method == "voltage": + U_p = U_n + pybamm.InputParameter("V_init") + return x, [ + invert(U, b, U_p, T_init, f"{n} {suffix}") for U, b, n in pos_init + ] + # No cell-voltage relation: lithium conservation fixes the positive side. + remaining = Q_Li - sum(q * s for q, s in zip(Q_n, x, strict=True)) + if not pos_composite: + return x, [remaining / Q_p[0]] + U_p = pybamm._BrentUnknown(f"positive potential {suffix}") + lo, hi = bracket(pos_init, [0] * len(pos_init), T_init) + solved = pybamm._Brent( + sum( + q * invert(U, b, U_p, T_init, f"{n} {suffix} guess") + for q, (U, b, n) in zip(Q_p, pos_init, strict=True) ) - - V_init - ) - self.algebraic[y_init_1] = ( - _get_cyclable_lithium_equation(options, "init") - Q_Li + - remaining, + U_p, + (lo, hi), + abstol=_ABSTOL, + max_iter=_MAX_ITER, ) + return x, [ + invert(U, b, solved, T_init, f"{n} {suffix} solved") + for U, b, n in pos_init + ] + + unknown = pybamm._BrentUnknown("x_init_1") + x_guess, y_guess = init_states(unknown, "guess") + if initialization_method == "voltage": + closure = total_lithium(x_guess, y_guess) - Q_Li elif initialization_method == "SOC": - soc_init = pybamm.InputParameter("SOC_init") - negative_soc = x_init_1 * pybamm.InputParameter("Q_n_1") - if is_negative_composite: - x_init_2 = variables["x_init_2"] - negative_soc += x_init_2 * pybamm.InputParameter("Q_n_2") - - negative_0_soc = x_0_1 * pybamm.InputParameter("Q_n_1") - if is_negative_composite: - negative_0_soc += x_0_2 * pybamm.InputParameter("Q_n_2") - - negative_100_soc = x_100_1 * pybamm.InputParameter("Q_n_1") - if is_negative_composite: - negative_100_soc += x_100_2 * pybamm.InputParameter("Q_n_2") - self.algebraic[x_init_1] = ( - (negative_soc - negative_0_soc) / (negative_100_soc - negative_0_soc) - ) - soc_init - self.algebraic[y_init_1] = ( - _get_cyclable_lithium_equation(options, "init") - Q_Li - ) + + def charge(x): + return sum(q * s for q, s in zip(Q_n, x, strict=True)) + + closure = (charge(x_guess) - charge(x_0)) / ( + charge(x_100) - charge(x_0) + ) - pybamm.InputParameter("SOC_init") else: - raise ValueError("Invalid initialization method") - T = param.T_init if initialization_method == "voltage" else param.T_ref - if is_positive_composite: - y_init_2 = variables["y_init_2"] - self.algebraic[y_init_2] = param.p.prim.U( - y_init_1, - T, - get_lithiation_delithiation( - direction, "positive", options, phase="primary" - ), - ) - param.p.sec.U( - y_init_2, - T, - get_lithiation_delithiation( - direction, "positive", options, phase="secondary" - ), - ) - if is_negative_composite: - x_init_2 = variables["x_init_2"] - self.algebraic[x_init_2] = param.n.prim.U( - x_init_1, - T, - get_lithiation_delithiation( - direction, "negative", options, phase="primary" - ), - ) - param.n.sec.U( - x_init_2, - T, - get_lithiation_delithiation( - direction, "negative", options, phase="secondary" - ), + raise pybamm.OptionError( + f"Invalid initialization method '{initialization_method}', " + "expected 'voltage' or 'SOC'" ) - self.variables.update(variables) - if initialization_method == "SOC": - soc_init = pybamm.InputParameter("SOC_init") - else: - soc_init = (V_init - V_min) / (V_max - V_min) - self.initial_conditions.update(_get_initial_conditions(options, soc_init)) + # Solve for x_init_1 rather than the potential: the residual is sensitive to + # it, and a non-physical target still has an exact answer. + x_init_1 = pybamm._Brent( + closure, + unknown, + (_BRACKET_LO, _BRACKET_HI), + abstol=_ABSTOL, + max_iter=_MAX_ITER, + ) + x_init, y_init = init_states(x_init_1, "solved") + + # the stoichiometries are expressions, so the solver needs a state of its own; + # the placeholder is kept out of `variables` + placeholder = pybamm.Variable("ESOH placeholder") + self.algebraic = {placeholder: placeholder} + self.initial_conditions = {placeholder: pybamm.Scalar(0)} + + for index, value in enumerate(x_100): + self.variables[f"x_100_{index + 1}"] = value + for index, value in enumerate(y_100): + self.variables[f"y_100_{index + 1}"] = value + for index, value in enumerate(x_0): + self.variables[f"x_0_{index + 1}"] = value + for index, value in enumerate(y_0): + self.variables[f"y_0_{index + 1}"] = value + for index, value in enumerate(x_init): + self.variables[f"x_init_{index + 1}"] = value + for index, value in enumerate(y_init): + self.variables[f"y_init_{index + 1}"] = value + + # set by `_esoh_evaluator`, keyed on the parameter values + self._evaluator: tuple | None = None @property def default_solver(self): @@ -803,10 +788,11 @@ def solve_full( inputs : dict, optional Additional inputs initial_conditions : dict, optional - Dictionary of initial conditions for variables (e.g., from split solve) + Accepted and ignored. Each stoichiometry is found by a bracketed + rootfind, which needs a bracket rather than a starting guess. esoh_sim : :class:`pybamm.Simulation`, optional - A pre-built simulation wrapping an :class:`ElectrodeSOHComposite` model - to reuse across calls. If not provided, a new one is created. + A pre-built simulation wrapping an :class:`ElectrodeSOHComposite` model. + Passing one back reuses its compiled evaluator across calls. Returns ------- @@ -831,28 +817,29 @@ def solve_full( Q_Li = parameter_values.evaluate(param.Q_Li_particles_init, inputs=inputs) - if isinstance(initial_value, str) and initial_value.endswith("V"): + initialization_method = _initialization_method(initial_value) + if initialization_method == "voltage": V_init = float(initial_value[:-1]) - initialization_method = "voltage" - elif isinstance(initial_value, float): - initialization_method = "SOC" - if initial_value > 1: - warnings.warn( - message=f"Initial SoC {initial_value} is greater than 1", - category=UserWarning, - stacklevel=2, - ) - elif initial_value < 0: - warnings.warn( - message=f"Initial SoC {initial_value} is less than 0", - category=UserWarning, - stacklevel=2, - ) - else: - raise ValueError( - "Invalid initial value. Expected a float between 0 and 1 " - "(for SOC) or a string ending in 'V' (for voltage), got " - f"{initial_value!r} of type {type(initial_value).__name__}" + elif initial_value > 1: + warnings.warn( + message=f"Initial SoC {initial_value} is greater than 1", + category=UserWarning, + stacklevel=2, + ) + elif initial_value < 0: + warnings.warn( + message=f"Initial SoC {initial_value} is less than 0", + category=UserWarning, + stacklevel=2, + ) + + # The stoichiometries are one expression holding rootfinds nested several + # deep, which only the native plugin evaluates quickly enough; on Windows it + # cannot be reached, so the caller's split solve answers instead. + if pybamm.is_windows(): + raise pybamm.SolverError( + "the full composite electrode SOH solve needs the 'brent' CasADi " + "rootfinder plugin, which is not available on Windows" ) all_inputs = {**inputs, **Qs, "Q_Li": Q_Li} @@ -865,20 +852,114 @@ def solve_full( model = ElectrodeSOHComposite( options, direction, initialization_method=initialization_method ) - esoh_sim = pybamm.Simulation( - model, - parameter_values=parameter_values, - solver=get_esoh_default_solver(tol), + else: + model = esoh_sim.model + + names, input_names, function = _esoh_evaluator(model, parameter_values) + try: + values = np.asarray( + function(*[all_inputs[name] for name in input_names]) + ).reshape(-1) + except RuntimeError as error: + # the native rootfinder reports failure as RuntimeError; callers, and the + # fallback in get_initial_stoichiometries_composite, expect a SolverError + raise pybamm.SolverError( + f"Composite electrode SOH solve failed: {error}" + ) from error + if not np.all(np.isfinite(values)): + raise pybamm.SolverError( + "Composite electrode SOH solve returned a non-finite stoichiometry" ) + return dict(zip(names, values, strict=True)) - if initial_conditions is not None: - esoh_sim.build() - esoh_sim.built_model.set_initial_conditions_from( - initial_conditions, inputs=all_inputs - ) - sol = esoh_sim.solve([0], inputs=all_inputs) - return {var: sol[var].entries[0] for var in sol.all_models[0].variables} +def _unique_nodes(roots): + """Every symbol reachable from ``roots``, visited once. + + Unlike ``pre_order``, which re-yields a shared node once per path to it. + """ + seen: set = set() + nodes = [] + stack = list(roots) + while stack: + symbol = stack.pop() + if id(symbol) in seen: + continue + seen.add(id(symbol)) + nodes.append(symbol) + stack.extend(symbol.children) + return nodes + + +def _parameter_fingerprint(parameter_values, names): + """A comparable summary of the parameters in ``names``. + + Numbers compare by value and everything else by identity, so replacing an OCP + function counts as a change but re-reading the same one does not. + """ + return tuple( + (name, value) + if isinstance(value := parameter_values[name], (int, float)) + else (name, id(value)) + for name in sorted(names) + if name in parameter_values + ) + + +def _esoh_evaluator(model, parameter_values): + """Map capacities and target straight to stoichiometries. + + The model defines each stoichiometry as an expression, not a state, so one CasADi + function replaces a solve. Cached on ``model`` per parameter set. + + Returns + ------- + tuple + ``(names, input_names, function)``, where ``function`` maps the inputs in + ``input_names`` order to the stoichiometries in ``names`` order. + """ + cached = model._evaluator + if cached is not None: + baked, fingerprint, names, input_names, function = cached + if _parameter_fingerprint(parameter_values, baked) == fingerprint: + return names, input_names, function + + names = sorted(model.variables) + nodes = _unique_nodes(model.variables[name] for name in names) + input_names = sorted( + {s.name for s in nodes if isinstance(s, pybamm.InputParameter)} + ) + # Only the parameters the graph substitutes can invalidate it. The capacities + # arrive as InputParameter, so a caller ageing a cell does not rebuild anything. + baked = frozenset( + s.name + for s in nodes + if isinstance(s, pybamm.Parameter | pybamm.FunctionParameter) + ) + symbols = {name: casadi.MX.sym(name) for name in input_names} + time = casadi.MX.sym("t") + state = casadi.MX.sym("y", 1) + # one conversion cache: the variables are views of a single shared graph + converted: dict = {} + expressions = [ + parameter_values.process_symbol(model.variables[name]).to_casadi( + time, state, inputs=symbols, casadi_symbols=converted + ) + for name in names + ] + function = casadi.Function( + "electrode_soh_composite", + list(symbols.values()), + [casadi.vertcat(*expressions)], + ) + model._evaluator = ( + baked, + _parameter_fingerprint(parameter_values, baked), + names, + input_names, + function, + ) + return names, input_names, function def get_initial_stoichiometries_composite( @@ -953,6 +1034,10 @@ def get_initial_stoichiometries_composite( "Only `cyclable lithium capacity` is supported for composite electrodes" ) + # A value that is neither an SOC nor a voltage is the caller's mistake, not a + # solve that failed, so it must not reach the fallback. + _initialization_method(initial_value) + try: return ElectrodeSOHComposite.solve_full( initial_value, @@ -964,45 +1049,26 @@ def get_initial_stoichiometries_composite( inputs=inputs, esoh_sim=esoh_sim, ) - except (pybamm.SolverError, ValueError) as first_error: - if try_split_solve: - try: - split_results = ElectrodeSOHComposite.solve_split( - initial_value, - parameter_values, - direction=direction, - param=param, - options=options, - tol=tol, - inputs=inputs, - ) - - try: - return ElectrodeSOHComposite.solve_full( - initial_value, - parameter_values, - direction=direction, - param=param, - options=options, - tol=tol, - inputs=inputs, - initial_conditions=split_results, - ) - except (pybamm.SolverError, ValueError) as retry_error: - raise ValueError( - f"Failed to solve composite electrode SOH. " - f"Initial full solve error: {first_error}. " - f"Retry with split solve initial conditions also failed: " - f"{retry_error}" - ) from retry_error - - except (pybamm.SolverError, ValueError) as split_error: - raise ValueError( - f"Failed to solve composite electrode SOH. " - f"Full solve error: {first_error}. " - f"Split solve error: {split_error}" - ) from split_error - else: - raise ValueError( + except (pybamm.SolverError, ValueError, RuntimeError) as first_error: + if not try_split_solve: + raise pybamm.SolverError( f"Failed to solve composite electrode SOH: {first_error}" ) from first_error + # The split solve reaches the answer a different way, one electrode at a + # time, so its result is the fallback rather than a guess to re-solve from. + try: + return ElectrodeSOHComposite.solve_split( + initial_value, + parameter_values, + direction=direction, + param=param, + options=options, + tol=tol, + inputs=inputs, + ) + except (pybamm.SolverError, ValueError, RuntimeError) as split_error: + raise pybamm.SolverError( + f"Failed to solve composite electrode SOH. " + f"Full solve error: {first_error}. " + f"Split solve error: {split_error}" + ) from split_error diff --git a/packages/pybamm/src/pybamm/util.py b/packages/pybamm/src/pybamm/util.py index d7ecc7043b..46de9976f1 100644 --- a/packages/pybamm/src/pybamm/util.py +++ b/packages/pybamm/src/pybamm/util.py @@ -368,6 +368,13 @@ def is_macos_intel(): return sys.platform == "darwin" and platform.machine() == "x86_64" +def is_windows(): + """Check if running on Windows.""" + import sys + + return sys.platform == "win32" + + def raise_jax_not_found(): """Raise an appropriate error when JAX is not available.""" if is_macos_intel(): diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh_composite_sweeps.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh_composite_sweeps.py new file mode 100644 index 0000000000..24d0b25a29 --- /dev/null +++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh_composite_sweeps.py @@ -0,0 +1,258 @@ +# +# Sweep the composite electrode SOH solver densely over its targets +# +from __future__ import annotations + +import numpy as np +import pytest + +import pybamm +from pybamm.models.full_battery_models.lithium_ion import electrode_soh_composite as esc +from pybamm.models.full_battery_models.lithium_ion.util import ( + get_lithiation_delithiation, +) + +# Dense enough that a bracket that fails only on a sliver of the range is still hit. +SWEEP = 1000 + +OPTIONS = { + "particle phases": ("2", "1"), + "open-circuit potential": (("single", "current sigmoid"), "single"), +} + +# (Q_n_1, Q_n_2, Q_p_1, Q_Li) multipliers on the nominal capacities +WEAR = { + "nominal": (1.0, 1.0, 1.0, 1.0), + "very low lithium": (1.0, 1.0, 1.0, 0.60), + "lost secondary": (1.0, 0.50, 1.0, 0.90), + "worn all": (0.70, 1.0, 0.90, 0.75), +} + + +class Case: + """A composite SOH solver bound to one direction and initialisation method.""" + + def __init__(self, direction, method): + self.options = pybamm.BatteryModelOptions(OPTIONS) + self.parameter_values = pybamm.ParameterValues("Chen2020_composite") + self.param = pybamm.LithiumIonParameters(self.options) + self.method = method + model = pybamm.lithium_ion.ElectrodeSOHComposite( + self.options, direction, initialization_method=method + ) + self.names, self.input_names, self.function = esc._esoh_evaluator( + model, self.parameter_values + ) + + # mirror the temperature and branch choice in ElectrodeSOHComposite.__init__ + reference = self.parameter_values.evaluate(self.param.T_ref) + temperature = ( + self.parameter_values.evaluate(self.param.T_init) + if method == "voltage" + else reference + ) + + def potential(side, phase, electrode): + stoichiometry = pybamm.InputParameter("sto") + branch = get_lithiation_delithiation( + direction, electrode, self.options, phase=phase + ) + processed = self.parameter_values.process_symbol( + (side.prim if phase == "primary" else side.sec).U( + stoichiometry, temperature, branch + ) + ) + return lambda value: float( + np.asarray(processed.evaluate(inputs={"sto": value})).reshape(-1)[0] + ) + + self.U_n = potential(self.param.n, "primary", "negative") + self.U_n2 = potential(self.param.n, "secondary", "negative") + self.U_p = potential(self.param.p, "primary", "positive") + + def capacities(self, scales): + evaluate = self.parameter_values.evaluate + return { + "Q_n_1": evaluate(self.param.n.prim.Q_init) * scales[0], + "Q_n_2": evaluate(self.param.n.sec.Q_init) * scales[1], + "Q_p_1": evaluate(self.param.p.prim.Q_init) * scales[2], + "Q_Li": evaluate(self.param.Q_Li_particles_init) * scales[3], + } + + def solve(self, capacities, target): + key = "V_init" if self.method == "voltage" else "SOC_init" + inputs = {**capacities, key: float(target)} + values = np.asarray( + self.function(*[inputs[name] for name in self.input_names]) + ).reshape(-1) + return dict(zip(self.names, values, strict=True)) + + def lithium(self, capacities, state): + return ( + capacities["Q_n_1"] * state["x_init_1"] + + capacities["Q_n_2"] * state["x_init_2"] + + capacities["Q_p_1"] * state["y_init_1"] + ) + + def state_of_charge(self, capacities, state): + def charge(tag): + return ( + capacities["Q_n_1"] * state[f"x_{tag}_1"] + + capacities["Q_n_2"] * state[f"x_{tag}_2"] + ) + + return (charge("init") - charge("0")) / (charge("100") - charge("0")) + + +class TestCompositeElectrodeSOHSweeps: + """Dense sweeps checking physics rather than stored answers.""" + + @pytest.mark.parametrize("direction", ["discharge", "charge", None]) + @pytest.mark.parametrize("wear", list(WEAR)) + def test_a_thousand_states_of_charge(self, direction, wear): + case = Case(direction, "SOC") + capacities = case.capacities(WEAR[wear]) + for target in np.linspace(0.0, 1.0, SWEEP): + state = case.solve(capacities, target) + assert all(np.isfinite(v) for v in state.values()), target + assert case.lithium(capacities, state) == pytest.approx( + capacities["Q_Li"], rel=1e-9 + ), target + assert case.state_of_charge(capacities, state) == pytest.approx( + target, abs=1e-7 + ) + assert case.U_n(state["x_init_1"]) == pytest.approx( + case.U_n2(state["x_init_2"]), abs=1e-6 + ), target + + @pytest.mark.parametrize("direction", ["discharge", "charge", None]) + @pytest.mark.parametrize("wear", list(WEAR)) + def test_a_thousand_voltages(self, direction, wear): + case = Case(direction, "voltage") + capacities = case.capacities(WEAR[wear]) + for target in np.linspace(2.5, 4.2, SWEEP): + state = case.solve(capacities, target) + assert all(np.isfinite(v) for v in state.values()), target + assert case.lithium(capacities, state) == pytest.approx( + capacities["Q_Li"], rel=1e-9 + ), target + voltage = case.U_p(state["y_init_1"]) - case.U_n(state["x_init_1"]) + assert voltage == pytest.approx(target, abs=1e-7) + + def test_repeating_a_solve_returns_identical_bits(self): + # the rootfinder caches its last solve, which must not change an answer + case = Case("discharge", "SOC") + capacities = case.capacities(WEAR["nominal"]) + for target in np.linspace(0.0, 1.0, SWEEP): + first = case.solve(capacities, target) + assert case.solve(capacities, target) == first + + def test_the_answer_moves_smoothly_with_the_target(self): + # a bracket that flipped to another root would show up as a jump + case = Case("discharge", "SOC") + capacities = case.capacities(WEAR["nominal"]) + targets = np.linspace(0.0, 1.0, SWEEP) + x_init = np.array([case.solve(capacities, t)["x_init_1"] for t in targets]) + steps = np.diff(x_init) + assert np.all(steps > 0), "x_init_1 must increase with state of charge" + assert np.max(steps) < 50 * np.median(steps), "discontinuity in x_init_1" + + def test_a_state_of_charge_round_trips_through_its_voltage(self): + by_soc = Case("discharge", "SOC") + by_voltage = Case("discharge", "voltage") + soc_capacities = by_soc.capacities(WEAR["nominal"]) + voltage_capacities = by_voltage.capacities(WEAR["nominal"]) + for target in np.linspace(0.05, 0.95, SWEEP): + state = by_soc.solve(soc_capacities, target) + voltage = by_voltage.U_p(state["y_init_1"]) - by_voltage.U_n( + state["x_init_1"] + ) + back = by_voltage.solve(voltage_capacities, voltage) + assert back["x_init_1"] == pytest.approx(state["x_init_1"], abs=1e-6) + + @pytest.mark.parametrize( + ("method", "targets"), + [ + ("SOC", np.linspace(-5.0, 5.0, SWEEP)), + ("voltage", np.linspace(0.5, 6.0, SWEEP)), + ], + ) + def test_a_target_outside_the_window_still_solves_exactly(self, method, targets): + """A non-physical target has an exact answer, and must be given it.""" + case = Case("discharge", method) + capacities = case.capacities(WEAR["nominal"]) + for target in targets: + state = case.solve(capacities, target) + assert all(np.isfinite(v) for v in state.values()), target + assert case.lithium(capacities, state) == pytest.approx( + capacities["Q_Li"], rel=1e-9 + ), target + if method == "voltage": + voltage = case.U_p(state["y_init_1"]) - case.U_n(state["x_init_1"]) + assert voltage == pytest.approx(target, abs=1e-6), target + else: + assert case.state_of_charge(capacities, state) == pytest.approx( + target, abs=1e-7 + ), target + + def test_a_non_physical_target_returns_a_non_physical_stoichiometry(self): + case = Case("discharge", "SOC") + capacities = case.capacities(WEAR["nominal"]) + assert case.solve(capacities, -3.0)["x_init_1"] < 0 + assert case.solve(capacities, 3.0)["x_init_1"] > 1 + + +class TestCompositeElectrodeSOHReuse: + """Ageing a cell must not rebuild the evaluator. + + The capacities reach the compiled function as inputs, so only the parameters its + graph substitutes can invalidate it. + """ + + # both fall as a cell ages, and both feed the capacities the solve is given + FADE = ( + "Initial concentration in positive electrode [mol.m-3]", + "Primary: Negative electrode active material volume fraction", + ) + + @staticmethod + def _reused_simulation(parameter_values): + model = pybamm.lithium_ion.ElectrodeSOHComposite( + pybamm.BatteryModelOptions(OPTIONS), + "discharge", + initialization_method="SOC", + ) + return pybamm.Simulation(model, parameter_values=parameter_values) + + def _call(self, parameter_values, simulation): + return pybamm.lithium_ion.get_initial_stoichiometries_composite( + 0.5, + parameter_values, + direction="discharge", + options=OPTIONS, + esoh_sim=simulation, + ) + + def test_fading_capacities_reuse_the_compiled_evaluator(self): + parameter_values = pybamm.ParameterValues("Chen2020_composite") + simulation = self._reused_simulation(parameter_values) + self._call(parameter_values, simulation) + built = simulation.model._evaluator[-1] + + nominal = {key: parameter_values[key] for key in self.FADE} + for step in range(1, 11): + for key in self.FADE: + parameter_values[key] = nominal[key] * (1 - 0.002 * step) + state = self._call(parameter_values, simulation) + assert np.isfinite(state["x_init_1"]) + assert simulation.model._evaluator[-1] is built + + def test_a_new_open_circuit_potential_does_rebuild(self): + parameter_values = pybamm.ParameterValues("Chen2020_composite") + simulation = self._reused_simulation(parameter_values) + self._call(parameter_values, simulation) + built = simulation.model._evaluator[-1] + + parameter_values["Primary: Negative electrode OCP [V]"] = lambda sto: 1.5 - sto + self._call(parameter_values, simulation) + assert simulation.model._evaluator[-1] is not built diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py index 460a11179a..08b0a0f2e7 100644 --- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py +++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py @@ -4,9 +4,13 @@ import contextlib +import numpy as np import pytest import pybamm +from pybamm.models.full_battery_models.lithium_ion.util import ( + get_lithiation_delithiation, +) # Fixture for TestElectrodeSOHMSMR, TestCalculateTheoreticalEnergy and TestGetInitialOCPMSMR class. @@ -285,6 +289,12 @@ def ocv(x, y, negative_branch, positive_branch): assert ocv(x_0, y_0, None, None) == pytest.approx(V_min + 0.1, abs=1e-6) +needs_full_solve = pytest.mark.skipif( + pybamm.is_windows(), + reason="the full composite solve needs the brent plugin, which Windows cannot load", +) + + class TestElectrodeSOHComposite: @staticmethod def _check_phases_equal(results, xy, soc): @@ -373,6 +383,7 @@ def _get_params_and_options(composite_electrode): "positive", # positive-only composite ], ) + @needs_full_solve def test_half_cell_with_same_ocp_curves(self, composite_electrode, initial_value): pvals, options = self._get_params_and_options(composite_electrode) # Use composite ESOH helper to compute initial stoichiometries at a voltage @@ -431,6 +442,7 @@ def test_chen2020_composite_defaults(self, initial_value): - param.n.prim.U(results["x_init_1"], param.T_ref) ) == pytest.approx(V_target, abs=1e-05) + @needs_full_solve def test_chen2020_composite_default_solve(self): pvals = pybamm.ParameterValues("Chen2020_composite") options = {"particle phases": ("2", "1")} @@ -1056,15 +1068,16 @@ def test_error_warning(self): "5 A", parameter_values_composite, options=options_composite ) - with pytest.warns(UserWarning, match=r"is greater than 1"): - pybamm.lithium_ion.ElectrodeSOHComposite.solve_full( - 1.001, parameter_values_composite, options=options_composite - ) + if not pybamm.is_windows(): + with pytest.warns(UserWarning, match=r"is greater than 1"): + pybamm.lithium_ion.ElectrodeSOHComposite.solve_full( + 1.001, parameter_values_composite, options=options_composite + ) - with pytest.warns(UserWarning, match=r"is less than 0"): - pybamm.lithium_ion.ElectrodeSOHComposite.solve_full( - -0.001, parameter_values_composite, options=options_composite - ) + with pytest.warns(UserWarning, match=r"is less than 0"): + pybamm.lithium_ion.ElectrodeSOHComposite.solve_full( + -0.001, parameter_values_composite, options=options_composite + ) with pytest.warns(UserWarning, match=r"is greater than 1"): pybamm.lithium_ion.ElectrodeSOHComposite.solve_split( 1.001, parameter_values_composite, options=options_composite @@ -1163,3 +1176,198 @@ def test_min_max_ocp(self, options): ) assert Up_100 - Un_100 == pytest.approx(4.2) assert Up_0 - Un_0 == pytest.approx(2.8) + + +class TestElectrodeSOHCompositeHardCases: + """Composite-electrode states checked against the physics, not stored answers. + + The requested state must be the one that comes back, and lithium conserved. + """ + + OPTIONS = { + "particle phases": ("2", "1"), + "open-circuit potential": (("single", "current sigmoid"), "single"), + } + + # (Q_n_1, Q_n_2, Q_p_1, Q_Li) multipliers on the nominal capacities + NOMINAL = (1.0, 1.0, 1.0, 1.0) + LOW_LI = (1.0, 1.0, 1.0, 0.85) + VERY_LOW_LI = (1.0, 1.0, 1.0, 0.6) + WORN_NEGATIVE = (0.85, 0.85, 1.0, 0.85) + LOST_SECONDARY = (1.0, 0.5, 1.0, 0.9) + WORN_ALL = (0.7, 1.0, 0.9, 0.75) + + CASES = [ + ("voltage", NOMINAL, 3.35), + ("voltage", WORN_ALL, 2.6416666666666666), + ("SOC", NOMINAL, 0.5833333333333334), + ("SOC", LOW_LI, 0.0), + ("SOC", LOW_LI, 0.08333333333333333), + ("SOC", VERY_LOW_LI, 0.0), + ("SOC", VERY_LOW_LI, 0.08333333333333333), + ("SOC", WORN_NEGATIVE, 0.08333333333333333), + ("SOC", LOST_SECONDARY, 0.0), + ("SOC", LOST_SECONDARY, 0.08333333333333333), + ] + + @staticmethod + def _setup(scales): + options = pybamm.BatteryModelOptions(TestElectrodeSOHCompositeHardCases.OPTIONS) + parameter_values = pybamm.ParameterValues("Chen2020_composite") + param = pybamm.LithiumIonParameters(options) + scale_n1, scale_n2, scale_p1, scale_li = scales + capacities = { + "Q_Li": parameter_values.evaluate(param.Q_Li_particles_init) * scale_li, + "Q_n_1": parameter_values.evaluate(param.n.prim.Q_init) * scale_n1, + "Q_n_2": parameter_values.evaluate(param.n.sec.Q_init) * scale_n2, + "Q_p_1": parameter_values.evaluate(param.p.prim.Q_init) * scale_p1, + } + return options, parameter_values, param, capacities + + @staticmethod + def _ocp(parameter_values, potential, stoichiometry, temperature, branch): + sto = pybamm.InputParameter("sto") + processed = parameter_values.process_symbol(potential(sto, temperature, branch)) + return float( + np.asarray(processed.evaluate(inputs={"sto": stoichiometry})).reshape(-1)[0] + ) + + # solves the model itself, whose variables hold the rootfinds + @needs_full_solve + @pytest.mark.parametrize(("initialization_method", "scales", "target"), CASES) + def test_reaches_the_requested_state(self, initialization_method, scales, target): + options, parameter_values, param, capacities = self._setup(scales) + model = pybamm.lithium_ion.ElectrodeSOHComposite( + options, initialization_method=initialization_method + ) + key = "V_init" if initialization_method == "voltage" else "SOC_init" + inputs = {**capacities, key: target} + sim = pybamm.Simulation(model, parameter_values=parameter_values) + solution = sim.solve([0], inputs=inputs) + state = {name: float(solution[name](0)) for name in model.variables} + + def branch(electrode, phase): + return get_lithiation_delithiation(None, electrode, options, phase=phase) + + if initialization_method == "voltage": + voltage = self._ocp( + parameter_values, + param.p.prim.U, + state["y_init_1"], + param.T_init, + branch("positive", "primary"), + ) - self._ocp( + parameter_values, + param.n.prim.U, + state["x_init_1"], + param.T_init, + branch("negative", "primary"), + ) + assert voltage == pytest.approx(target, abs=1e-08) + else: + + def charge(tag): + return ( + capacities["Q_n_1"] * state[f"x_{tag}_1"] + + capacities["Q_n_2"] * state[f"x_{tag}_2"] + ) + + soc = (charge("init") - charge("0")) / (charge("100") - charge("0")) + assert soc == pytest.approx(target, abs=1e-08) + + lithium = ( + capacities["Q_n_1"] * state["x_init_1"] + + capacities["Q_n_2"] * state["x_init_2"] + + capacities["Q_p_1"] * state["y_init_1"] + ) + assert lithium == pytest.approx(capacities["Q_Li"], rel=1e-10) + + +class TestElectrodeSOHCompositeFallback: + """The split solve is the fallback when the bracketed solve cannot answer. + + Nothing in the parameter sweep fails, so these force the first attempt to. + """ + + OPTIONS = { + "particle phases": ("2", "1"), + "open-circuit potential": (("single", "current sigmoid"), "single"), + } + + def _call(self, **kwargs): + return pybamm.lithium_ion.get_initial_stoichiometries_composite( + 0.5, + pybamm.ParameterValues("Chen2020_composite"), + direction="discharge", + options=self.OPTIONS, + **kwargs, + ) + + @pytest.mark.parametrize( + "error", + [ + pybamm.SolverError("full solve failed"), + RuntimeError("rootfinder process failed"), + ValueError("bad bracket"), + ], + ) + def test_a_failed_full_solve_falls_back_to_the_split_solve( + self, monkeypatch, error + ): + def fail(*args, **kwargs): + raise error + + expected = pybamm.lithium_ion.ElectrodeSOHComposite.solve_split( + 0.5, + pybamm.ParameterValues("Chen2020_composite"), + direction="discharge", + options=self.OPTIONS, + ) + monkeypatch.setattr( + pybamm.lithium_ion.ElectrodeSOHComposite, "solve_full", fail + ) + from_split = self._call(try_split_solve=True) + assert from_split.keys() == expected.keys() + for name, value in expected.items(): + assert from_split[name] == pytest.approx(value, abs=1e-10), name + + def test_the_fallback_can_be_turned_off(self, monkeypatch): + def fail(*args, **kwargs): + raise RuntimeError("rootfinder process failed") + + monkeypatch.setattr( + pybamm.lithium_ion.ElectrodeSOHComposite, "solve_full", fail + ) + with pytest.raises(pybamm.SolverError, match="Failed to solve composite"): + self._call(try_split_solve=False) + + def test_both_failing_reports_both_errors(self, monkeypatch): + def fail_full(*args, **kwargs): + raise pybamm.SolverError("full is broken") + + def fail_split(*args, **kwargs): + raise pybamm.SolverError("split is broken") + + monkeypatch.setattr( + pybamm.lithium_ion.ElectrodeSOHComposite, "solve_full", fail_full + ) + monkeypatch.setattr( + pybamm.lithium_ion.ElectrodeSOHComposite, "solve_split", fail_split + ) + with pytest.raises( + pybamm.SolverError, match=r"full is broken.*split is broken" + ): + self._call(try_split_solve=True) + + def test_a_non_finite_stoichiometry_is_a_failure_not_a_result(self, monkeypatch): + # a solve that returns NaN must reach the fallback, not be handed back + def nan_evaluator(model, parameter_values): + names = sorted(model.variables) + return names, [], lambda *a: np.full(len(names), np.nan) + + monkeypatch.setattr( + pybamm.models.full_battery_models.lithium_ion.electrode_soh_composite, + "_esoh_evaluator", + nan_evaluator, + ) + assert np.isfinite(self._call(try_split_solve=True)["x_init_1"])