diff --git a/CHANGELOG.md b/CHANGELOG.md index 7d36a0c4bc..2e0ba2178a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,14 +1,19 @@ # [Unreleased](https://github.com/pybamm-team/PyBaMM/) +## Features + +- Add support for BPX version 0.4.0 which allows for blended electrodes and user-defined parameters in BPX([#3414](https://github.com/pybamm-team/PyBaMM/pull/3414)) +- Added the ability to specify a custom solver tolerance in `get_initial_stoichiometries` and related functions ([#3714](https://github.com/pybamm-team/PyBaMM/pull/3714)) + ## Bug Fixes - Updated `_steps_util.py` to throw a specific exception when drive cycle starts at t>0 ([#3756](https://github.com/pybamm-team/PyBaMM/pull/3756)) - Updated `plot_voltage_components.py` to support both `Simulation` and `Solution` objects. Added new methods in both `Simulation` and `Solution` classes for allow the syntax `simulation.plot_voltage_components` and `solution.plot_voltage_components`. Updated `test_plot_voltage_components.py` to reflect these changes ([#3723](https://github.com/pybamm-team/PyBaMM/pull/3723)). - The SEI thickness decreased at some intervals when the 'electron-migration limited' model was used. It has been corrected ([#3622](https://github.com/pybamm-team/PyBaMM/pull/3622)) -## Features +## Breaking changes -- Added the ability to specify a custom solver tolerance in `get_initial_stoichiometries` and related functions ([#3714](https://github.com/pybamm-team/PyBaMM/pull/3714)) +- Dropped support for BPX version 0.3.0 and below ([#3414](https://github.com/pybamm-team/PyBaMM/pull/3414)) # [v24.1rc2](https://github.com/pybamm-team/PyBaMM/tree/v24.1rc2) - 2024-01-24 diff --git a/pybamm/parameters/bpx.py b/pybamm/parameters/bpx.py index a97288b062..e5ed5bb724 100644 --- a/pybamm/parameters/bpx.py +++ b/pybamm/parameters/bpx.py @@ -1,4 +1,5 @@ from bpx import BPX, Function, InterpolatedTable +from bpx.schema import ElectrodeBlended, ElectrodeBlendedSPM import pybamm import math from dataclasses import dataclass @@ -7,22 +8,39 @@ from pybamm import exp -import types -import functools +from functools import partial -def _copy_func(f): - """Based on http://stackoverflow.com/a/6528148/190597 (Glenn Maynard)""" - g = types.FunctionType( - f.__code__, - f.__globals__, - name=f.__name__, - argdefs=f.__defaults__, - closure=f.__closure__, - ) - g = functools.update_wrapper(g, f) - g.__kwdefaults__ = f.__kwdefaults__ - return g +def _callable_func(var, fun): + return fun(var) + + +def _interpolant_func(var, name, x, y): + return pybamm.Interpolant(x, y, var, name=name, interpolator="linear") + + +preamble = "from pybamm import exp, tanh, cosh\n\n" + + +def process_float_function_table(value, name): + """ + Process BPX FloatFunctionTable to a float, python function or data for a pybamm + Interpolant. + """ + if isinstance(value, Function): + value = value.to_python_function(preamble=preamble) + elif isinstance(value, InterpolatedTable): + # return (name, (x, y)) to match the output of + # `pybamm.parameters.process_1D_data` we will create an interpolant on a + # case-by-case basis to get the correct argument for each parameter + x = np.array(value.x) + y = np.array(value.y) + # sort the arrays as CasADi requires x to be in ascending order + sort_idx = np.argsort(x) + x = x[sort_idx] + y = y[sort_idx] + value = (name, (x, y)) + return value @dataclass @@ -59,18 +77,49 @@ class Domain: separator = Domain(name="separator", pre_name="Separator ", short_pre_name="") experiment = Domain(name="experiment", pre_name="", short_pre_name="") +PHASE_NAMES = ["Primary: ", "Secondary: "] + + +def _get_phase_names(domain): + """ + Return a list of the phase names in a given domain + """ + if isinstance(domain, (ElectrodeBlended, ElectrodeBlendedSPM)): + phases = len(getattr(domain, "particle").keys()) + else: + phases = 1 + if phases == 1: + return [""] + elif phases == 2: + return ["Primary: ", "Secondary: "] + else: + raise NotImplementedError( + "PyBaMM does not support more than two " + "particle phases in blended electrodes" + ) + def _bpx_to_param_dict(bpx: BPX) -> dict: + """ + Turns a BPX object in to a dictionary of parameters for PyBaMM + """ + domain_phases = { + "negative electrode": _get_phase_names(bpx.parameterisation.negative_electrode), + "positive electrode": _get_phase_names(bpx.parameterisation.positive_electrode), + } + + # Loop over each component of BPX and add to pybamm dictionary pybamm_dict = {} pybamm_dict = _bpx_to_domain_param_dict( - bpx.parameterisation.cell, pybamm_dict, cell + bpx.parameterisation.positive_electrode, pybamm_dict, positive_electrode ) pybamm_dict = _bpx_to_domain_param_dict( - bpx.parameterisation.negative_electrode, pybamm_dict, negative_electrode + bpx.parameterisation.cell, pybamm_dict, cell ) pybamm_dict = _bpx_to_domain_param_dict( - bpx.parameterisation.positive_electrode, pybamm_dict, positive_electrode + bpx.parameterisation.negative_electrode, pybamm_dict, negative_electrode ) + pybamm_dict = _bpx_to_domain_param_dict( bpx.parameterisation.electrolyte, pybamm_dict, electrolyte ) @@ -88,7 +137,7 @@ def _bpx_to_param_dict(bpx: BPX) -> dict: # activity pybamm_dict["Thermodynamic factor"] = 1.0 - # assume Bruggeman relation for effection electrolyte properties + # assume Bruggeman relation for effective electrolyte properties for domain in [negative_electrode, separator, positive_electrode]: pybamm_dict[domain.pre_name + "Bruggeman coefficient (electrolyte)"] = 1.5 @@ -119,9 +168,6 @@ def _bpx_to_param_dict(bpx: BPX) -> dict: # reference temperature T_ref = pybamm_dict["Reference temperature [K]"] - def arrhenius(Ea, T): - return exp(Ea / constants.R * (1 / T_ref - 1 / T)) - # lumped parameters for name in [ "Specific heat capacity [J.K-1.kg-1]", @@ -163,92 +209,23 @@ def arrhenius(Ea, T): {"Total heat transfer coefficient [W.m-2.K-1]": 0}, check_already_exists=False ) - # BET surface area - for domain in [negative_electrode, positive_electrode]: - pybamm_dict[domain.pre_name + "active material volume fraction"] = ( - pybamm_dict[domain.pre_name + "surface area per unit volume [m-1]"] - * pybamm_dict[domain.short_pre_name + "particle radius [m]"] - ) / 3.0 - # transport efficiency for domain in [negative_electrode, separator, positive_electrode]: pybamm_dict[domain.pre_name + "porosity"] = pybamm_dict[ domain.pre_name + "transport efficiency" ] ** (1.0 / 1.5) - # TODO: allow setting function parameters in a loop over domains - - # ocp - - # negative electrode (only need to check for data, other cases pass through) - U_n = pybamm_dict[negative_electrode.pre_name + "OCP [V]"] - if isinstance(U_n, tuple): - - def _negative_electrode_ocp(sto): - name, (x, y) = U_n - return pybamm.Interpolant(x, y, sto, name=name, interpolator="linear") - - pybamm_dict[negative_electrode.pre_name + "OCP [V]"] = _negative_electrode_ocp - - # positive electrode (only need to check for data, other cases pass through) - U_p = pybamm_dict[positive_electrode.pre_name + "OCP [V]"] - if isinstance(U_p, tuple): - - def _positive_electrode_ocp(sto): - name, (x, y) = U_p - return pybamm.Interpolant(x, y, sto, name=name, interpolator="linear") + # define functional forms for pybamm parameters that depend on more than one + # variable - pybamm_dict[positive_electrode.pre_name + "OCP [V]"] = _positive_electrode_ocp - - # entropic change - - # negative electrode - dUdT_n = pybamm_dict[ - negative_electrode.pre_name + "entropic change coefficient [V.K-1]" - ] - if callable(dUdT_n): - - def _negative_electrode_entropic_change(sto, c_s_max): - return dUdT_n(sto) - - elif isinstance(dUdT_n, tuple): - - def _negative_electrode_entropic_change(sto, c_s_max): - name, (x, y) = dUdT_n - return pybamm.Interpolant(x, y, sto, name=name, interpolator="linear") - - else: - - def _negative_electrode_entropic_change(sto, c_s_max): - return dUdT_n - - pybamm_dict[ - negative_electrode.pre_name + "OCP entropic change [V.K-1]" - ] = _negative_electrode_entropic_change - - # positive electrode - dUdT_p = pybamm_dict[ - positive_electrode.pre_name + "entropic change coefficient [V.K-1]" - ] - if callable(dUdT_p): - - def _positive_electrode_entropic_change(sto, c_s_max): - return dUdT_p(sto) - - elif isinstance(dUdT_p, tuple): - - def _positive_electrode_entropic_change(sto, c_s_max): - name, (x, y) = dUdT_p - return pybamm.Interpolant(x, y, sto, name=name, interpolator="linear") - - else: - - def _positive_electrode_entropic_change(sto, c_s_max): - return dUdT_p + def _arrhenius(Ea, T): + return exp(Ea / constants.R * (1 / T_ref - 1 / T)) - pybamm_dict[ - positive_electrode.pre_name + "OCP entropic change [V.K-1]" - ] = _positive_electrode_entropic_change + def _entropic_change(sto, c_s_max, dUdT, constant=False): + if constant: + return dUdT + else: + return dUdT(sto) # reaction rates in pybamm exchange current is defined j0 = k * sqrt(ce * cs * # (cs-cs_max)) in BPX exchange current is defined j0 = F * k_norm * sqrt((ce/ce0) * @@ -256,125 +233,117 @@ def _positive_electrode_entropic_change(sto, c_s_max): c_e = pybamm_dict["Initial concentration in electrolyte [mol.m-3]"] F = 96485 - # negative electrode - c_n_max = pybamm_dict[ - "Maximum concentration in " + negative_electrode.pre_name.lower() + "[mol.m-3]" - ] - k_n_norm = pybamm_dict[ - negative_electrode.pre_name + "reaction rate constant [mol.m-2.s-1]" - ] - Ea_k_n = pybamm_dict.get( - negative_electrode.pre_name - + "reaction rate constant activation energy [J.mol-1]", - 0.0, - ) - # Note that in BPX j = 2*F*k_norm*sqrt((ce/ce0)*(c/c_max)*(1-c/c_max))*sinh(...), - # and in PyBaMM j = 2*k*sqrt(ce*c*(c_max - c))*sinh(...) - k_n = k_n_norm * F / (c_n_max * c_e**0.5) - - def _negative_electrode_exchange_current_density(c_e, c_s_surf, c_s_max, T): - k_ref = k_n # (A/m2)(m3/mol)**1.5 - includes ref concentrations - + def _exchange_current_density(c_e, c_s_surf, c_s_max, T, k_ref, Ea): return ( k_ref - * arrhenius(Ea_k_n, T) + * _arrhenius(Ea, T) * c_e**0.5 * c_s_surf**0.5 * (c_s_max - c_s_surf) ** 0.5 ) - pybamm_dict[ - negative_electrode.pre_name + "exchange-current density [A.m-2]" - ] = _copy_func(_negative_electrode_exchange_current_density) + def _diffusivity(sto, T, D_ref, Ea, constant=False): + if constant: + return _arrhenius(Ea, T) * D_ref + else: + return _arrhenius(Ea, T) * D_ref(sto) - # positive electrode - c_p_max = pybamm_dict[ - "Maximum concentration in " + positive_electrode.pre_name.lower() + "[mol.m-3]" - ] - k_p_norm = pybamm_dict[ - positive_electrode.pre_name + "reaction rate constant [mol.m-2.s-1]" - ] - Ea_k_p = pybamm_dict.get( - positive_electrode.pre_name - + "reaction rate constant activation energy [J.mol-1]", - 0.0, - ) - # Note that in BPX j = 2*F*k_norm*sqrt((ce/ce0)*(c/c_max)*(1-c/c_max))*sinh(...), - # and in PyBaMM j = 2*k*sqrt(ce*c*(c_max - c))*sinh(...) - k_p = k_p_norm * F / (c_p_max * c_e**0.5) + def _conductivity(c_e, T, Ea, sigma_ref, constant=False): + if constant: + return _arrhenius(Ea, T) * sigma_ref + else: + return _arrhenius(Ea, T) * sigma_ref(c_e) - def _positive_electrode_exchange_current_density(c_e, c_s_surf, c_s_max, T): - k_ref = k_p # (A/m2)(m3/mol)**1.5 - includes ref concentrations - - return ( - k_ref - * arrhenius(Ea_k_p, T) - * c_e**0.5 - * c_s_surf**0.5 - * (c_s_max - c_s_surf) ** 0.5 - ) - - pybamm_dict[domain.pre_name + "exchange-current density [A.m-2]"] = _copy_func( - _positive_electrode_exchange_current_density - ) - - # diffusivity - - # negative electrode - Ea_D_n = pybamm_dict.get( - negative_electrode.pre_name + "diffusivity activation energy [J.mol-1]", 0.0 - ) - D_n_ref = pybamm_dict[negative_electrode.pre_name + "diffusivity [m2.s-1]"] - - if callable(D_n_ref): - - def _negative_electrode_diffusivity(sto, T): - return arrhenius(Ea_D_n, T) * D_n_ref(sto) - - elif isinstance(D_n_ref, tuple): - - def _negative_electrode_diffusivity(sto, T): - name, (x, y) = D_n_ref - return arrhenius(Ea_D_n, T) * pybamm.Interpolant( - x, y, sto, name=name, interpolator="linear" + # Loop over electrodes and construct derived parameters + for domain in [negative_electrode, positive_electrode]: + for phase_pre_name in domain_phases[domain.name]: + phase_domain_pre_name = phase_pre_name + domain.pre_name + + # BET surface area + pybamm_dict[phase_domain_pre_name + "active material volume fraction"] = ( + pybamm_dict[ + phase_domain_pre_name + "surface area per unit volume [m-1]" + ] + * pybamm_dict[ + phase_pre_name + domain.short_pre_name + "particle radius [m]" + ] + ) / 3.0 + + # ocp + U = pybamm_dict[phase_domain_pre_name + "OCP [V]"] + if isinstance(U, tuple): + pybamm_dict[phase_domain_pre_name + "OCP [V]"] = partial( + _interpolant_func, name=U[0], x=U[1][0], y=U[1][1] + ) + + # entropic change + dUdT = pybamm_dict[ + phase_domain_pre_name + "entropic change coefficient [V.K-1]" + ] + if callable(dUdT): + pybamm_dict[ + phase_domain_pre_name + "OCP entropic change [V.K-1]" + ] = partial(_entropic_change, dUdT=dUdT) + elif isinstance(dUdT, tuple): + pybamm_dict[ + phase_domain_pre_name + "OCP entropic change [V.K-1]" + ] = partial( + _entropic_change, + dUdT=partial( + _interpolant_func, name=dUdT[0], x=dUdT[1][0], y=dUdT[1][1] + ), + ) + else: + pybamm_dict[ + phase_domain_pre_name + "OCP entropic change [V.K-1]" + ] = partial(_entropic_change, dUdT=dUdT, constant=True) + + # reaction rate + c_max = pybamm_dict[ + phase_pre_name + + "Maximum concentration in " + + domain.pre_name.lower() + + "[mol.m-3]" + ] + k_norm = pybamm_dict[ + phase_domain_pre_name + "reaction rate constant [mol.m-2.s-1]" + ] + Ea_k = pybamm_dict.get( + phase_domain_pre_name + + "reaction rate constant activation energy [J.mol-1]", + 0.0, ) - - else: - - def _negative_electrode_diffusivity(sto, T): - return arrhenius(Ea_D_n, T) * D_n_ref - - pybamm_dict[negative_electrode.pre_name + "diffusivity [m2.s-1]"] = _copy_func( - _negative_electrode_diffusivity - ) - - # positive electrode - Ea_D_p = pybamm_dict.get( - positive_electrode.pre_name + "diffusivity activation energy [J.mol-1]", 0.0 - ) - D_p_ref = pybamm_dict[positive_electrode.pre_name + "diffusivity [m2.s-1]"] - - if callable(D_p_ref): - - def _positive_electrode_diffusivity(sto, T): - return arrhenius(Ea_D_p, T) * D_p_ref(sto) - - elif isinstance(D_p_ref, tuple): - - def _positive_electrode_diffusivity(sto, T): - name, (x, y) = D_p_ref - return arrhenius(Ea_D_p, T) * pybamm.Interpolant( - x, y, sto, name=name, interpolator="linear" + # Note that in BPX j = 2*F*k_norm*sqrt((ce/ce0)*(c/c_max)*(1-c/c_max))... + # *sinh(), + # and in PyBaMM j = 2*k*sqrt(ce*c*(c_max - c))*sinh() + k = k_norm * F / (c_max * c_e**0.5) + pybamm_dict[ + phase_domain_pre_name + "exchange-current density [A.m-2]" + ] = partial(_exchange_current_density, k_ref=k, Ea=Ea_k) + + # diffusivity + Ea_D = pybamm_dict.get( + phase_domain_pre_name + "diffusivity activation energy [J.mol-1]", + 0.0, ) - - else: - - def _positive_electrode_diffusivity(sto, T): - return arrhenius(Ea_D_p, T) * D_p_ref - - pybamm_dict[positive_electrode.pre_name + "diffusivity [m2.s-1]"] = _copy_func( - _positive_electrode_diffusivity - ) + D_ref = pybamm_dict[phase_domain_pre_name + "diffusivity [m2.s-1]"] + + if callable(D_ref): + pybamm_dict[phase_domain_pre_name + "diffusivity [m2.s-1]"] = partial( + _diffusivity, D_ref=D_ref, Ea=Ea_D + ) + elif isinstance(D_ref, tuple): + pybamm_dict[phase_domain_pre_name + "diffusivity [m2.s-1]"] = partial( + _diffusivity, + D_ref=partial( + _interpolant_func, name=D_ref[0], x=D_ref[1][0], y=D_ref[1][1] + ), + Ea=Ea_D, + ) + else: + pybamm_dict[phase_domain_pre_name + "diffusivity [m2.s-1]"] = partial( + _diffusivity, D_ref=D_ref, Ea=Ea_D, constant=True + ) # electrolyte Ea_D_e = pybamm_dict.get( @@ -383,26 +352,21 @@ def _positive_electrode_diffusivity(sto, T): D_e_ref = pybamm_dict[electrolyte.pre_name + "diffusivity [m2.s-1]"] if callable(D_e_ref): - - def _electrolyte_diffusivity(sto, T): - return arrhenius(Ea_D_e, T) * D_e_ref(sto) - + pybamm_dict[electrolyte.pre_name + "diffusivity [m2.s-1]"] = partial( + _diffusivity, D_ref=D_e_ref, Ea=Ea_D_e + ) elif isinstance(D_e_ref, tuple): - - def _electrolyte_diffusivity(sto, T): - name, (x, y) = D_e_ref - return arrhenius(Ea_D_e, T) * pybamm.Interpolant( - x, y, sto, name=name, interpolator="linear" - ) - + pybamm_dict[electrolyte.pre_name + "diffusivity [m2.s-1]"] = partial( + _diffusivity, + D_ref=partial( + _interpolant_func, name=D_e_ref[0], x=D_e_ref[1][0], y=D_e_ref[1][1] + ), + Ea=Ea_D_e, + ) else: - - def _electrolyte_diffusivity(sto, T): - return arrhenius(Ea_D_e, T) * D_e_ref - - pybamm_dict[electrolyte.pre_name + "diffusivity [m2.s-1]"] = _copy_func( - _electrolyte_diffusivity - ) + pybamm_dict[electrolyte.pre_name + "diffusivity [m2.s-1]"] = partial( + _diffusivity, D_ref=D_e_ref, Ea=Ea_D_e, constant=True + ) # conductivity Ea_sigma_e = pybamm_dict.get( @@ -411,68 +375,96 @@ def _electrolyte_diffusivity(sto, T): sigma_e_ref = pybamm_dict[electrolyte.pre_name + "conductivity [S.m-1]"] if callable(sigma_e_ref): - - def _conductivity(c_e, T): - return arrhenius(Ea_sigma_e, T) * sigma_e_ref(c_e) - + pybamm_dict[electrolyte.pre_name + "conductivity [S.m-1]"] = partial( + _conductivity, sigma_ref=sigma_e_ref, Ea=Ea_sigma_e + ) elif isinstance(sigma_e_ref, tuple): - - def _conductivity(c_e, T): - name, (x, y) = sigma_e_ref - return arrhenius(Ea_sigma_e, T) * pybamm.Interpolant( - x, y, c_e, name=name, interpolator="linear" - ) - + pybamm_dict[electrolyte.pre_name + "conductivity [S.m-1]"] = partial( + _conductivity, + sigma_ref=partial( + _interpolant_func, + name=sigma_e_ref[0], + x=sigma_e_ref[1][0], + y=sigma_e_ref[1][1], + ), + Ea=Ea_sigma_e, + ) else: + pybamm_dict[electrolyte.pre_name + "conductivity [S.m-1]"] = partial( + _conductivity, sigma_ref=sigma_e_ref, Ea=Ea_sigma_e, constant=True + ) - def _conductivity(c_e, T): - return arrhenius(Ea_sigma_e, T) * sigma_e_ref - - pybamm_dict[electrolyte.pre_name + "conductivity [S.m-1]"] = _copy_func( - _conductivity - ) - + # Add user-defined parameters, if any + user_defined = bpx.parameterisation.user_defined + if user_defined: + for name in user_defined.__dict__.keys(): + value = getattr(user_defined, name) + value = process_float_function_table(value, name) + if callable(value): + pybamm_dict[name] = partial(_callable_func, fun=value) + elif isinstance(value, tuple): + pybamm_dict[name] = partial( + _interpolant_func, name=value[0], x=value[1][0], y=value[1][1] + ) + else: + pybamm_dict[name] = value return pybamm_dict -preamble = "from pybamm import exp, tanh, cosh\n\n" +def _get_pybamm_name(pybamm_name, domain): + """ + Process pybamm name to include domain name and handle special cases + """ + pybamm_name_lower = pybamm_name[:1].lower() + pybamm_name[1:] + if pybamm_name.startswith("Initial concentration") or pybamm_name.startswith( + "Maximum concentration" + ): + init_len = len("Initial concentration ") + pybamm_name = ( + pybamm_name[:init_len] + + "in " + + domain.pre_name.lower() + + pybamm_name[init_len:] + ) + elif pybamm_name.startswith("Particle radius"): + pybamm_name = domain.short_pre_name + pybamm_name_lower + elif pybamm_name.startswith("OCP"): + pybamm_name = domain.pre_name + pybamm_name + elif pybamm_name.startswith("Cation transference number"): + pybamm_name = pybamm_name + elif domain.pre_name != "": + pybamm_name = domain.pre_name + pybamm_name_lower + return pybamm_name def _bpx_to_domain_param_dict(instance: BPX, pybamm_dict: dict, domain: Domain) -> dict: + """ + Turns a BPX instance in to a dictionary of parameters for PyBaMM for a given domain + """ + # Loop over fields in BPX instance and add to pybamm dictionary for name, field in instance.__fields__.items(): value = getattr(instance, name) - if value is None: - continue - elif isinstance(value, Function): - value = value.to_python_function(preamble=preamble) - elif isinstance(value, InterpolatedTable): - # return (name, (x, y)) to match the output of - # `pybamm.parameters.process_1D_data` we will create an interpolant on a - # case-by-case basis to get the correct argument for each parameter - x = np.array(value.x) - y = np.array(value.y) - value = (name, (x, y)) - - pybamm_name = field.field_info.alias - pybamm_name_lower = pybamm_name[:1].lower() + pybamm_name[1:] - if pybamm_name.startswith("Initial concentration") or pybamm_name.startswith( - "Maximum concentration" + # Handle blended electrodes, where the field is now an instance of + # ElectrodeBlended or ElectrodeBlendedSPM + if ( + isinstance(instance, (ElectrodeBlended, ElectrodeBlendedSPM)) + and name == "particle" ): - init_len = len("Initial concentration ") - pybamm_name = ( - pybamm_name[:init_len] - + "in " - + domain.pre_name.lower() - + pybamm_name[init_len:] - ) - elif pybamm_name.startswith("Particle radius"): - pybamm_name = domain.short_pre_name + pybamm_name_lower - elif pybamm_name.startswith("OCP"): - pybamm_name = domain.pre_name + pybamm_name - elif pybamm_name.startswith("Cation transference number"): - pybamm_name = pybamm_name - elif domain.pre_name != "": - pybamm_name = domain.pre_name + pybamm_name_lower - - pybamm_dict[pybamm_name] = value + particle_instance = getattr(instance, "particle") + # Loop over phases + for i, phase_name in enumerate(particle_instance.keys()): + phase_instance = particle_instance[phase_name] + # Loop over fields in phase instance and add to pybamm dictionary + for name, field in phase_instance.__fields__.items(): + value = getattr(phase_instance, name) + pybamm_name = PHASE_NAMES[i] + _get_pybamm_name( + field.field_info.alias, domain + ) + value = process_float_function_table(value, name) + pybamm_dict[pybamm_name] = value + # Handle other fields, which correspond directly to parameters + else: + pybamm_name = _get_pybamm_name(field.field_info.alias, domain) + value = process_float_function_table(value, name) + pybamm_dict[pybamm_name] = value return pybamm_dict diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py index 4256561c96..85c8f02b53 100644 --- a/pybamm/parameters/parameter_values.py +++ b/pybamm/parameters/parameter_values.py @@ -94,6 +94,7 @@ def create_from_bpx(filename, target_soc=1): raise ValueError("Target SOC should be between 0 and 1") from bpx import parse_bpx_file, get_electrode_concentrations + from bpx.schema import ElectrodeBlended, ElectrodeBlendedSPM from .bpx import _bpx_to_param_dict # parse bpx @@ -111,9 +112,25 @@ def create_from_bpx(filename, target_soc=1): # ahead with the low voltage limit. # get initial concentrations based on SOC - c_n_init, c_p_init = get_electrode_concentrations(target_soc, bpx) - pybamm_dict["Initial concentration in negative electrode [mol.m-3]"] = c_n_init - pybamm_dict["Initial concentration in positive electrode [mol.m-3]"] = c_p_init + # Note: we cannot set SOC for blended electrodes, + # see https://github.com/pybamm-team/PyBaMM/issues/2682 + bpx_neg = bpx.parameterisation.negative_electrode + bpx_pos = bpx.parameterisation.positive_electrode + if isinstance(bpx_neg, (ElectrodeBlended, ElectrodeBlendedSPM)) or isinstance( + bpx_pos, (ElectrodeBlended, ElectrodeBlendedSPM) + ): + pybamm.logger.warning( + "Initial concentrations cannot be set using stoichiometry limits for " + "blend electrodes. Please set the initial concentrations manually." + ) + else: + c_n_init, c_p_init = get_electrode_concentrations(target_soc, bpx) + pybamm_dict[ + "Initial concentration in negative electrode [mol.m-3]" + ] = c_n_init + pybamm_dict[ + "Initial concentration in positive electrode [mol.m-3]" + ] = c_p_init return pybamm.ParameterValues(pybamm_dict) diff --git a/pyproject.toml b/pyproject.toml index a16f7819c6..7b74168703 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -92,7 +92,7 @@ latexify = [ ] # Battery Parameter eXchange format bpx = [ - "bpx", + "bpx>=0.4.0", ] # Low-overhead progress bars tqdm = [ diff --git a/tests/unit/test_parameters/test_bpx.py b/tests/unit/test_parameters/test_bpx.py index e131e906c4..694020f157 100644 --- a/tests/unit/test_parameters/test_bpx.py +++ b/tests/unit/test_parameters/test_bpx.py @@ -299,6 +299,165 @@ def arrhenius_assertion(pv, param_key, Ea_key): for param_key, Ea_key in zip(param_keys, Ea_keys): arrhenius_assertion(pv, param_key, Ea_key) + def test_bpx_blended(self): + bpx_obj = copy.copy(self.base) + bpx_obj["Parameterisation"]["Positive electrode"] = { + "Thickness [m]": 5.23e-05, + "Conductivity [S.m-1]": 0.789, + "Porosity": 0.277493, + "Transport efficiency": 0.1462, + "Particle": { + "Large Particles": { + "Diffusivity [m2.s-1]": 3.2e-14, + "Particle radius [m]": 8e-06, + "OCP [V]": "-3.04420906 * x + 10.04892207 - 0.65637536 * tanh(-4.02134095 * (x - 0.80063948)) + 4.24678547 * tanh(12.17805062 * (x - 7.57659337)) - 0.3757068 * tanh(59.33067782 * (x - 0.99784492))", + "Entropic change coefficient [V.K-1]": -1e-4, + "Surface area per unit volume [m-1]": 186331, + "Reaction rate constant [mol.m-2.s-1]": 2.305e-05, + "Minimum stoichiometry": 0.42424, + "Maximum stoichiometry": 0.96210, + "Maximum concentration [mol.m-3]": 46200, + "Diffusivity activation energy [J.mol-1]": 15000, + "Reaction rate constant activation energy [J.mol-1]": 3500, + }, + "Small Particles": { + "Diffusivity [m2.s-1]": 3.2e-14, + "Particle radius [m]": 1e-06, + "OCP [V]": "-3.04420906 * x + 10.04892207 - 0.65637536 * tanh(-4.02134095 * (x - 0.80063948)) + 4.24678547 * tanh(12.17805062 * (x - 7.57659337)) - 0.3757068 * tanh(59.33067782 * (x - 0.99784492))", + "Entropic change coefficient [V.K-1]": -1e-4, + "Surface area per unit volume [m-1]": 496883, + "Reaction rate constant [mol.m-2.s-1]": 2.305e-05, + "Minimum stoichiometry": 0.42424, + "Maximum stoichiometry": 0.96210, + "Maximum concentration [mol.m-3]": 46200, + "Diffusivity activation energy [J.mol-1]": 15000, + "Reaction rate constant activation energy [J.mol-1]": 3500, + }, + }, + } + + filename = "tmp.json" + with tempfile.NamedTemporaryFile( + suffix=filename, delete=False, mode="w" + ) as tmp: + # write to a tempory file so we can + # get the source later on using inspect.getsource + # (as long as the file still exists) + json.dump(bpx_obj, tmp) + tmp.flush() + + pv = pybamm.ParameterValues.create_from_bpx(tmp.name) + # initial concentration must be set manually for blended models (for now) + pv.update( + { + "Initial concentration in negative electrode [mol.m-3]": 22000, + "Primary: Initial concentration in positive electrode [mol.m-3]": 19404, + "Secondary: Initial concentration in positive electrode [mol.m-3]": 19404, + }, + check_already_exists=False, + ) + model = pybamm.lithium_ion.SPM({"particle phases": ("1", "2")}) + experiment = pybamm.Experiment( + [ + "Discharge at C/5 for 1 hour", + ] + ) + sim = pybamm.Simulation(model, parameter_values=pv, experiment=experiment) + sim.solve(calc_esoh=False) + + def test_bpx_blended_error(self): + bpx_obj = copy.copy(self.base) + bpx_obj["Parameterisation"]["Positive electrode"] = { + "Thickness [m]": 5.23e-05, + "Conductivity [S.m-1]": 0.789, + "Porosity": 0.277493, + "Transport efficiency": 0.1462, + "Particle": { + "Large Particles": { + "Diffusivity [m2.s-1]": 3.2e-14, + "Particle radius [m]": 8e-06, + "OCP [V]": "-3.04420906 * x + 10.04892207 - 0.65637536 * tanh(-4.02134095 * (x - 0.80063948)) + 4.24678547 * tanh(12.17805062 * (x - 7.57659337)) - 0.3757068 * tanh(59.33067782 * (x - 0.99784492))", + "Entropic change coefficient [V.K-1]": -1e-4, + "Surface area per unit volume [m-1]": 186331, + "Reaction rate constant [mol.m-2.s-1]": 2.305e-05, + "Minimum stoichiometry": 0.42424, + "Maximum stoichiometry": 0.96210, + "Maximum concentration [mol.m-3]": 46200, + "Diffusivity activation energy [J.mol-1]": 15000, + "Reaction rate constant activation energy [J.mol-1]": 3500, + }, + "Medium Particles": { + "Diffusivity [m2.s-1]": 3.2e-14, + "Particle radius [m]": 4e-06, + "OCP [V]": "-3.04420906 * x + 10.04892207 - 0.65637536 * tanh(-4.02134095 * (x - 0.80063948)) + 4.24678547 * tanh(12.17805062 * (x - 7.57659337)) - 0.3757068 * tanh(59.33067782 * (x - 0.99784492))", + "Entropic change coefficient [V.K-1]": -1e-4, + "Surface area per unit volume [m-1]": 186331, + "Reaction rate constant [mol.m-2.s-1]": 2.305e-05, + "Minimum stoichiometry": 0.42424, + "Maximum stoichiometry": 0.96210, + "Maximum concentration [mol.m-3]": 46200, + "Diffusivity activation energy [J.mol-1]": 15000, + "Reaction rate constant activation energy [J.mol-1]": 3500, + }, + "Small Particles": { + "Diffusivity [m2.s-1]": 3.2e-14, + "Particle radius [m]": 1e-06, + "OCP [V]": "-3.04420906 * x + 10.04892207 - 0.65637536 * tanh(-4.02134095 * (x - 0.80063948)) + 4.24678547 * tanh(12.17805062 * (x - 7.57659337)) - 0.3757068 * tanh(59.33067782 * (x - 0.99784492))", + "Entropic change coefficient [V.K-1]": -1e-4, + "Surface area per unit volume [m-1]": 186331, + "Reaction rate constant [mol.m-2.s-1]": 2.305e-05, + "Minimum stoichiometry": 0.42424, + "Maximum stoichiometry": 0.96210, + "Maximum concentration [mol.m-3]": 46200, + "Diffusivity activation energy [J.mol-1]": 15000, + "Reaction rate constant activation energy [J.mol-1]": 3500, + }, + }, + } + + filename = "tmp.json" + with tempfile.NamedTemporaryFile( + suffix=filename, delete=False, mode="w" + ) as tmp: + # write to a tempory file so we can + # get the source later on using inspect.getsource + # (as long as the file still exists) + json.dump(bpx_obj, tmp) + tmp.flush() + + with self.assertRaisesRegex(NotImplementedError, "PyBaMM does not support"): + pybamm.ParameterValues.create_from_bpx(tmp.name) + + def test_bpx_user_defined(self): + bpx_obj = copy.copy(self.base) + data = {"x": [0, 1], "y": [0, 1]} + bpx_obj["Parameterisation"]["User-defined"] = { + "User-defined scalar parameter": 1.0, + "User-defined parameter data": data, + "User-defined parameter data function": "x**2", + } + + filename = "tmp.json" + with tempfile.NamedTemporaryFile( + suffix=filename, delete=False, mode="w" + ) as tmp: + # write to a tempory file so we can + # get the source later on using inspect.getsource + # (as long as the file still exists) + json.dump(bpx_obj, tmp) + tmp.flush() + + param = pybamm.ParameterValues.create_from_bpx(tmp.name) + + self.assertEqual(param["User-defined scalar parameter"], 1.0) + var = pybamm.Variable("var") + self.assertIsInstance( + param["User-defined parameter data"](var), pybamm.Interpolant + ) + self.assertIsInstance( + param["User-defined parameter data function"](var), pybamm.Power + ) + if __name__ == "__main__": print("Add -v for more debug output")