[Bug]: create_from_bpx(): reaction rate activation energies are always zero

[ejfdickinson]

PyBaMM Version

23.4

Python Version

3.9.16

Describe the bug

Issue

The effective values of reaction rate constant activation energies specified in a BPX JSON file are always zero after processing by create_from_bpx().

Note: identified in 23.4 but no relevant code changes in latest version. See also #3051 which reports other out-of-specification behaviour of the same function.

Cause:

This is because BPX specifies that these have the name Reaction rate constant activation energy [J.mol-1] which is rendered by _bpx_to_domain_param_dict() into

Negative electrode reaction rate constant activation energy [J.mol-1]
Positive electrode reaction rate constant activation energy [J.mol-1]

However, _bpx_to_param_dict() checks as follows when recreating exchange current density functions:

E_a_n = pybamm_dict.get(
    negative_electrode.pre_name + "reaction rate activation energy [J.mol-1]", 0.0
)
[...]
E_a_p = pybamm_dict.get(
    positive_electrode.pre_name + "reaction rate activation energy [J.mol-1]", 0.0
)

Because the word "constant" is missing from the lookup key, the dict entries are not found, and the imported reaction rate constant activation energies are always replaced with the value field 0.0 in the resulting exchange current density function in the ParameterValues object.

Proposed fix:

Either:

• Correct the name lookups in _bpx_to_param_dict(), or
• Delete the substring "constant " when processing these fields in _bpx_to_domain_param_dict().

Steps to Reproduce

Using nmc_pouch_cell_BPX.json from the BPX example download:

import pybamm

from pybamm import constants, exp, log

# From BPX example download
params = pybamm.ParameterValues.create_from_bpx("nmc_pouch_cell_BPX.json")

cs_max = params["Maximum concentration in positive electrode [mol.m-3]"]
cl_ref = params["Initial concentration in electrolyte [mol.m-3]"]

Ea_kin_pos = params["Positive electrode reaction rate constant activation energy [J.mol-1]"]

# Evaluate exchange current density at Tref
T_ref  = params["Reference temperature [K]"]
j0_pos_ref = params["Positive electrode exchange-current density [A.m-2]"](cl_ref, cs_max/2, cs_max, T_ref)

# Evaluate exchange current density at Ttest =/= Tref
T_test       = 290
j0_pos_test = params["Positive electrode exchange-current density [A.m-2]"](cl_ref, cs_max/2, cs_max, T_test)

# Evaluate expected exchange current density at Ttest
j0_pos_test_expected = j0_pos_test * exp(-Ea_kin_pos / constants.R * (1/T_test - 1/T_ref))

# Evaluate apparent activation energy
Ea_kin_pos_apparent = -constants.R / (1/T_test - 1/T_ref) * log(j0_pos_test/j0_pos_ref)

print("Positive electrode reaction rate constant activation energy (test temperature, expected):", j0_pos_test_expected)
print("Positive electrode reaction rate constant activation energy (test temperature, actual):", j0_pos_test)

print("Positive electrode reaction rate constant activation energy activation energy (stated):", Ea_kin_pos)
print("Positive electrode reaction rate constant activation energy activation energy (apparent):", Ea_kin_pos_apparent)

### Relevant log output

```shell
Positive electrode reaction rate constant activation energy (test temperature, expected): 0.7477868882542091
Positive electrode reaction rate constant activation energy (test temperature, actual): 1.111989625
Positive electrode reaction rate constant activation energy activation energy (stated): 35000.0
Positive electrode reaction rate constant activation energy activation energy (apparent): -0.0

[rtimms]

Thanks for pointing this out. I think the best solution is to correct the name lookups in _bpx_to_param_dict() to match the schema.

For #3051 we should just use distinct names, as suggested.

We clearly need a better test for create_from_bpx(). At the moment we just test that you can load a BPX and run a simulation or that (some of) the variables end up being the correct types. We should do a comparison between a ParameterValues created from BPX and the correct ParameterValues.

@ejfdickinson it seems like you may have some code snippets like the one above for testing specific parameters (especially those embedded within functions, like activation energies). Do you think you could have a go at providing a fix plus updating the test?

[ejfdickinson]

@rtimms I think that @darryl-ad may be able to take on this issue and #3051, not least since we have a self-interest in the fixes being ready for the next release.

For a test, I think probably the most appropriate method would be to use the ParameterValues object to simulate under a couple of known conditions and check for equality of results. 'Equivalence' is otherwise quite nebulous as one can always write functions in slightly different ways under the hood.

[rtimms]

Thanks, sounds good. I guess for "equivalence" I was thinking let's at least check all the scalar values we can read directly or calculate (e.g. by evaluating functions at some fixed values, like in your above example) are correct.

@darryl-ad let me know if you need any help.

[rtimms]

@all-contributors please add @ejfdickinson for ideas and bug reports

[allcontributors]

@rtimms

I've put up a pull request to add @ejfdickinson! 🎉

[darryl-ad]

No problem I'll take this and #3051

[darryl-ad]

@rtimms @ejfdickinson I've fixed both of the discussed issues in PR #3242.

I tackled both issues in a single PR to prevent merge conflicts.

@rtimms certain parameters are both scalar values in the parameter dictionary, and also hard-coded into other functional parameters during their creation in _bpx_to_param_dict(). These include activation energies, k_ref/k_norm etc.

Is it possible to replace the hard-coded values inside the functional parameters, with references using pybamm.Parameter()? This would mean that updating the activation energy entries in the parameter dictionary will also update their values in the functional parameters. If this is not possible it might be better to delete these values where possible from the parameter dictionary, as it's ambiguous which parameters are used by PyBaMM and which aren't.

[rtimms]

Yeah, I agree, for BPX the values that already appear as scalars in ParameterValues should be referenced using pybamm.Parameter() in the functions. Thanks for the suggestion.

[rtimms]

@all-contributors please add @darryl-ad for ideas

[allcontributors]

@rtimms

I've put up a pull request to add @darryl-ad! 🎉