Repository: pybamm-team/PyBaMM (python)

Issue:
`ParameterValues.create_from_bpx` - support BPX v0.4
### Description

Add support to `ParameterValues.create_from_bpx()` for new [BPX 0.4](https://bpxstandard.com/bpx-standard/) features.

- Import blended electrode parameters according to specification.
- Consider whether PyBaMM could define a standard for importing user-defined parameters for model extensions. As an example, this might include a read from BPX for the necessary branches for `current sigmoid` hysteresis from user-defined BPX `FloatFunctionTable` features.

### Motivation

_No response_

### Possible Implementation

A simple PyBaMM-internal standard might be that the contents of the BPX user-defined parameters node is read directly into the `ParameterValues` object. Then users could add any PyBaMM-compatible parameter names to the BPX, under this node.

Since the content of a BPX file implies the type of model that should be used (DFN vs SPM, blended electrode or not, hysteresis or not), could `create_from_BPX` optionally return a `dict` of appropriate model options? Or else another utility function could do this?

### Additional context

_No response_

Pull request description:
[WIP] updates for BPX v0.4.0 compatibility 
# Description
Add support for blended electrodes and user-defined parameters in BPX. Create parameters from BPX file in a loop using `partial`.

PyBaMM now requires `bpx>=0.4.0`

Note: BPX set initial concentrations based on SOC, which is not yet supported in pybamm for blended electrodes (see  https://github.com/pybamm-team/PyBaMM/issues/2682). If a user supplies a "blended electrode" BPX file then they must manually provide initial concentrations after the `ParameterValues` class has been created (or supply them in the "user-defined" section).

Fixes #3679 

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [ ] New feature (non-breaking change which adds functionality)
- [ ] Optimization (back-end change that speeds up the code)
- [ ] Bug fix (non-breaking change which fixes an issue)

# Key checklist:

- [ ] No style issues: `$ pre-commit run` (or `$ nox -s pre-commit`) (see [CONTRIBUTING.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CONTRIBUTING.md#installing-and-using-pre-commit) for how to set this up to run automatically when committing locally, in just two lines of code)
- [ ] All tests pass: `$ python run-tests.py --all` (or `$ nox -s tests`)
- [ ] The documentation builds: `$ python run-tests.py --doctest` (or `$ nox -s doctests`)

You can run integration tests, unit tests, and doctests together at once, using `$ python run-tests.py --quick` (or `$ nox -s quick`).

## Further checks:

- [ ] Code is commented, particularly in hard-to-understand areas
- [ ] Tests added that prove fix is effective or that feature works


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7d36a0c4b..2e0ba2178 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 a97288b06..e5ed5bb72 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 constants
 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 @@ electrolyte = Domain(name="electrolyte", pre_name="Electrolyte ", short_pre_name
 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 _bpx_to_param_dict(bpx: BPX) -> dict:
         {"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 _bpx_to_param_dict(bpx: BPX) -> dict:
     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 _bpx_to_param_dict(bpx: BPX) -> dict:
     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 _bpx_to_param_dict(bpx: BPX) -> dict:
     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 4256561c9..85c8f02b5 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -94,6 +94,7 @@ class ParameterValues:
             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 @@ class ParameterValues:
             # 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 a16f7819c..7b7416870 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 = [


Test patch:
diff --git a/tests/unit/test_parameters/test_bpx.py b/tests/unit/test_parameters/test_bpx.py
index e131e906c..694020f15 100644
--- a/tests/unit/test_parameters/test_bpx.py
+++ b/tests/unit/test_parameters/test_bpx.py
@@ -299,6 +299,165 @@ class TestBPX(TestCase):
         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")