Defining diffusivity as a function of lithiation / SOC

[SimMLPhys]

I have been trying to figure out how to define diffusivity as a function of lithiation/SOC. As a rudimentary example, consider I am using a predefined parameter set:

parameters_1 = pybamm.ParameterValues("Chen2020")
parameters_1['Negative electrode diffusivity [m2.s-1]'] = lambda c_s_n, T: neg_diff(c_s_n, T, my_graphite_diff_parameter) 

where I have a previously defined function as follows:

      def neg_diff(c_s_n,T,diff_param):

        
        max_val = 28746 #maximum concentration in negative electrode
        
        ratio = c_s_n/max_val

        return pybamm.maximum(ratio*diff_param*3e-13, 5e-15)

In this case, even if diff_param = 1 and the ratio is higher than 0.1 (so the return statement with the maximum should return the first argument instead of 5e-15), the value for the diffusivity that gets returned is 5e-15. I think this has to do with the fact that the "ratio" term isn't getting treated like a number. When I print ratio, this is the output:

Graphite ratio: (maximum(minimum(3.478744868851319e-05 * X-averaged negative particle concentration [mol.m-3], 0.9999999999), 1e-10)) / max(maximum(minimum(3.478744868851319e-05 * X-averaged negative particle concentration [mol.m-3], 0.9999999999), 1e-10))

How can I get this working so that the ratio time some number gets returned?

[SimMLPhys]

I am adding more details here because I did not want to obscure the original post with too much extraneous information:

Notably, when I have a function where it's just returning ratio*3e-13, then it does return numbers near 3e-13. However, my ratio at this time should have been around 0.5, as I had defined
parameters_1['Initial concentration in negative electrode [mol.m-3]'] = 13400 parameters_1['Maximum concentration in negative electrode [mol.m-3]'] = 28746

So the returned diffusivity should have been 1.5e-13, not around 3e-13 (the initial part of my experiment is just a rest for 30 seconds).
experiment = pybamm.Experiment( [("Rest for 30 seconds", f"Charge at C/2 for 600 seconds or until 4.2 V")])

This is how I am checking my returned diffusivity values -

negative_diffusivity = solution["Negative particle effective diffusivity [m2.s-1]"]

time = solution["Time [s]"].entries
for t in time:
    diffusivity_value_at_t = negative_diffusivity(t)
    # Extract the surface diffusivity value, which is the last element
    surface_diffusivity = diffusivity_value_at_t[-1, -1]
    print(f"At time {t} s, negative electrode surface diffusivity = {surface_diffusivity} m2/s")

`
I am using a SPMe model so surface diffusivity = diffusivity throughout the particle.

[valentinsulzer]

The input for the diffusivity function is the stoichiometry (already divided by max_val), not the concentration.

[SimMLPhys]

Thanks, that worked out well. That's also why my ratios were so low - I was dividing by an extra term with magnitude 1e-4.

What's the best way for me to check these definitions? I figured it out by going through the particle submodel here https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/models/submodels/particle/base_particle.py, but is there a way for to do it within my pybamm code while debugging?

[valentinsulzer]

It's not perfect, but the inputs of the functions are generally documented in the parameter sets https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/input/parameters/lithium_ion/Chen2020.py and in https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/parameters/lithium_ion_parameters.py