[Bug]: Interpolant getting `nans` as input

[brosaplanella]

PyBaMM Version

develop

Python Version

3.11

Describe the bug

When defining the diffusivity as an interpolant, the stoichiometry passed is as nan. This means that the function gets evaluated at the midpoint of the interpolation range rather than at its actual value.

Steps to Reproduce

• Set debug breakpoint in line 270 of interpolant.py (i.e. inside _function_evaluate)
• Run minimal_example_of_lookup_tables.py
• Check that evaluated_children contains nans in the second argument (stoichiometry). The first time I believe it is for shape evaluation only, so it's fine to have nans, but the subsequent evaluations are in the solver.

Relevant log output

[DrSOKane]

Does this happen for any quantity defined as an interpolant, or only diffusivity?

[ahdezm]

Running the same minimal_example_of_lookup_tables.py calculation but interpolating the "Negative electrode OCP [V]" or "Electrolyte conductivity [S.m-1]" parameters gives similar results, with the concentration or stoichiometry being NaNs.

I was the person who initially noticed this bug and don't know the PyBaMM internals very well, but it seems like some of the parameters which sto depends on are not numbers (i.e. "Negative particle concentration [mol.m-3]") for some reason. Setting sto = pybamm.full_like([sto], 0.5) in D_s_n before creating the pybamm.Interpolant makes the issue disappear.

I have also noticed that running minimal_interp3d_example.py with an additional print statement for evaluated_children in line 270 of interpolant.py also gives NaNs on every variable, which I'm not sure how to explain.

[ahdezm]

I'm not entirely sure if this is expected, but the NaNs seem to be coming from the y StateVector, which always has that value when evaluated on an Interpolant. This can be confirmed by adding the following snippet before line 130 in functions.py (in the evaluate function) and running the example:

if type(self).__name__ == "Interpolant":
    print("all_nan?", np.isnan(y).all())

[martinjrobins]

if _function_evaluate() is being called it is not part of the solve. Before the solve happens the expression tree gets converted to a casadi function, and this is the only thing that is called as part of a solve, not the python code. When an expression tree is evaluated to determine the shape of a symbol, nans are used for the state vector values, so this might be what you are seeing.

Do you have a test case where the output of a simulation is affected?

[ahdezm]

I assume that at some point in this conversion the pybamm.Interpolant should be evaluated with non-NaN values in all children (see line 316 in interpolant.py) so that the resulting interpolated values can be used in the solve, but this doesn't seem to be happening.

If I understand correctly, any interpolated function which depends on the state (which will always be NaN during conversion) will always be evaluated at the same point (the midpoint), giving incorrect results.

[martinjrobins]

no, the solver will never evaluate any python code, by that time the interpolant has been converted to a casadi function. The only time the code in pybamm.Interpolant will be called is before or after the solve (unless you are using the scipy solver with convert_to=python)

[ahdezm]

Thank you for explaining, that makes sense - the conversion for interpolants seems to happen in convert_to_casadi.py L143.

I suspect that this issue can be closed now, but just for my peace of mind is there an easy way of extracting the interpolation inputs and outputs for diffusivity from the CasADi solution (I'm not sure how to match the input temperature and concentration with parts of the Solution object) or should I just run the calculation with the scipy solver?

[martinjrobins]

just add the inputs and outputs as additional variables of the model (i.e. add to the model.variables attribute)

[ahdezm]

Thank you again for your help, the values for linear interpolation seem to match.

When I use the cubic method instead, however, the resulting diffusivity doesn't match what I expect (it's smaller and negative). I might be doing something wrong but please find attached my testing script below.

pybamm_interpolant_test.txt

[brosaplanella]

Thanks @martinjrobins for the input on this! Closing it as the issue seems resolved.

[brosaplanella]

Sorry, I closed this without noting the comment on the cubic interpolant. @martinjrobins, do you have any insight on this? My best guess is that it is to do with the shape of the 2D geometry.

Maybe using the SPM model, for which the concentration would only be a function of time and radius, might make it easier to debug.