2D Interpolation

[chmabaur]

I struggle to get the 2D interpolation working (1D interpolation works just fine). However, I would like to have the exchange current density dependent on the stoichiometry and temperature, using lookup tables and 2d interpolation. For that i use the function:

def exchange_current_density(c_e, c_s_surf, c_s_max, T):
    temperatures = np.array([0, 25, 45])
    sto = np.linspace(0, 1, 100)
    x = (sto, temperatures)
    y = np.ones((100, 3)) # just an example to illustrate the interpolation issue
       
    c_s_rel = c_s_surf/c_s_max
    z = [c_s_rel, T]
    
    i_0 = pybamm.Interpolant(x, y, z)
    
    alpha = 0.5
    c_el_initial = 1000.0
    return (c_e/c_el_initial)**alpha * i_0

But when trying to use this in a DFN model, i receive the error:

During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\expression_tree\symbol.py:758 in evaluate_for_shape
    return self._saved_evaluate_for_shape

AttributeError: 'Subtraction' object has no attribute '_saved_evaluate_for_shape'


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\expression_tree\symbol.py:758 in evaluate_for_shape
    return self._saved_evaluate_for_shape

AttributeError: 'Concatenation' object has no attribute '_saved_evaluate_for_shape'


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File c:\repositories\batmodel\pybamm simulations\dfn_isothermal_gtw.py:66
    solution = sim.solve()

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\simulation.py:578 in solve
    self.build(check_model=check_model, initial_soc=initial_soc)

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\simulation.py:463 in build
    self.set_parameters()

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\simulation.py:413 in set_parameters
    self._model_with_set_params = self._parameter_values.process_model(

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\parameters\parameter_values.py:447 in process_model
    new_algebraic[new_variable] = self.process_symbol(equation)

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\parameters\parameter_values.py:605 in process_symbol
    processed_symbol = self._process_symbol(symbol)

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\parameters\parameter_values.py:719 in _process_symbol
    new_symbol = symbol._binary_new_copy(new_left, new_right)

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\expression_tree\binary_operators.py:102 in _binary_new_copy
    return self._binary_evaluate(left, right)

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\expression_tree\binary_operators.py:273 in _binary_evaluate
    return left * right

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\expression_tree\symbol.py:559 in __mul__
    return pybamm.simplified_multiplication(self, other)

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\expression_tree\binary_operators.py:967 in simplified_multiplication
    if left.shape_for_testing == right.shape_for_testing and all(

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\expression_tree\symbol.py:942 in shape_for_testing
    evaluated_self = self.evaluate_for_shape()

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\expression_tree\symbol.py:760 in evaluate_for_shape
    self._saved_evaluate_for_shape = self._evaluate_for_shape()

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\expression_tree\binary_operators.py:113 in _evaluate_for_shape
    right = self.children[1].evaluate_for_shape()

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\expression_tree\symbol.py:760 in evaluate_for_shape
    self._saved_evaluate_for_shape = self._evaluate_for_shape()

  File ~\AppData\Local\Programs\Python\Environments\gen-env\lib\site-packages\pybamm\expression_tree\concatenations.py:128 in _evaluate_for_shape
    return concatenation_function(

  File <__array_function__ internals>:180 in concatenate

ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 1 and the array at index 2 has size 51

[brosaplanella]

Note that for 2D Interpolants, you need to pass a tuple with x and y as first argument, and z as second (see the docs for interpolant).

[chmabaur]

@brosaplanella Isn't this what I do already? My x is a tuple of two 1D arrays, y the 2d array of function values to interpolate, and z the coordinates where the interpolation should be evaluated at (Interpolant(x, y, z))

[brosaplanella]

Sorry, I read too quickly! Your syntax is correct. Can you share a MWE so I can debug it at my end?

[chmabaur]

Thanks for looking into this, very much appreciated! PyBamm Version 23.5.

import pybamm
import numpy as np

def exchange_current_density(c_e, c_s_surf, c_s_max, T):
    temperatures = np.array([0, 25, 45])
    sto = np.linspace(0, 1, 100)
    x = (sto, temperatures)
    y = np.ones((100, 3)) # just an example to illustrate the interpolation issue
       
    c_s_rel = c_s_surf/c_s_max
    z = [c_s_rel, T]
    
    i_0 = pybamm.Interpolant(x, y, z)
    
    alpha = 0.5
    c_el_initial = 1000.0
    return (c_e/c_el_initial)**alpha * i_0


model = pybamm.lithium_ion.DFN()
geometry = model.default_geometry
param = pybamm.ParameterValues("Chen2020")
param.update({
    "Negative electrode exchange-current density [A.m-2]": exchange_current_density
    })

experiment = pybamm.Experiment([
    ("Discharge at 1C until 1.5V",
     "Rest for 1 hour"
    )
    ], period='1 second')   

solver = pybamm.CasadiSolver(mode="safe", atol=1e-6, rtol=1e-6, dt_max=10)

sim = pybamm.Simulation(model, experiment=experiment, parameter_values=param,
                        solver=solver)
solution = sim.solve(calc_esoh=False)

[ChenZHAO-factorial]

hello, is there any updates on this? I am having the same problem. Thanks