[Bug]: 23.9 - ShapeError from electrolyte equation when specifying "Ambient temperature [K]" as pybamm.Interpolant

[ejfdickinson]

PyBaMM Version

23.9

Python Version

3.9.13

Describe the bug

As of 23.9, using a pybamm.Interpolant to specify the "Ambient temperature [K]", to run a spatially isothermal model with a specified transient temperature profile, raises a ShapeError within Discretisation.process_dict() acting on the equation with key "Electrolyte current density [A.m-2]".

This parameter definition pattern was supported and functional in 23.5.

Steps to Reproduce

Works in 23.5. Fails in 23.9.

import pybamm
import numpy as np
import matplotlib.pyplot as plt

times = np.arange(0,1810,10)
tmax = max(times)
experiment = pybamm.Experiment([f"Discharge at 1C for {tmax} s"])

temp_drive_cycle = pybamm.Interpolant(
    times,
    298.15 + 20 * (times / tmax),
    pybamm.t
)

parameter_values = pybamm.ParameterValues("Chen2020")
parameter_values.update(
    {
        "Initial temperature [K]": 298.15,
        "Ambient temperature [K]": temp_drive_cycle
    }
)

sim = pybamm.Simulation(
    model=pybamm.lithium_ion.DFN(),
    parameter_values=parameter_values,
    experiment=experiment
)

sol = sim.solve()
sol.plot()

plt.plot(sol["Time [s]"].entries, sol["Volume-averaged cell temperature [K]"].entries - 273.15)
plt.xlabel("Time / s")
plt.ylabel("Cell temperature / degC");

Relevant log output

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ...\pybamm\discretisations\discretisation.py:728, in Discretisation.process_symbol(self, symbol)
    727 try:
--> 728     return self._discretised_symbols[symbol]
    729 except KeyError:

KeyError: Multiplication(0x31886afc2ee59eaf, *, children=['(concatenation(broadcast(1.6212500000000003e-11), broadcast(4.1791384878584733e-11), broadcast(2.5148216275389202e-11)) * ((maximum(concatenation(broadcast(4.0), broadcast(2.127659574468085), broadcast(2.9850746268656714)) * Porosity times concentration [mol.m-3](Negative electrode porosity times concentration [mol.m-3], Separator porosity times concentration [mol.m-3], Positive electrode porosity times concentration [mol.m-3]), 10.0)) ** 3.0)) - (concatenation(broadcast(9.92164615877829e-06), broadcast(2.557528655362438e-05), broadcast(1.5390081937322494e-05)) * ((maximum(concatenation(broadcast(4.0), broadcast(2.127659574468085), broadcast(2.9850746268656714)) * Porosity times concentration [mol.m-3](Negative electrode porosity times concentration [mol.m-3], Separator porosity times concentration [mol.m-3], Positive electrode porosity times concentration [mol.m-3]), 10.0)) ** 1.5)) + concatenation(broadcast(0.00041612500000000005), broadcast(0.0010726562857425487), broadcast(0.0006454773475772603)) * (maximum(concatenation(broadcast(4.0), broadcast(2.127659574468085), broadcast(2.9850746268656714)) * Porosity times concentration [mol.m-3](Negative electrode porosity times concentration [mol.m-3], Separator porosity times concentration [mol.m-3], Positive electrode porosity times concentration [mol.m-3]), 10.0))', '((0.00012763994028092067 * ing_functio(time) / (maximum(concatenation(broadcast(4.0), broadcast(2.127659574468085), broadcast(2.9850746268656714)) * Porosity times concentration [mol.m-3](Negative electrode porosity times concentration [mol.m-3], Separator porosity times concentration [mol.m-3], Positive electrode porosity times concentration [mol.m-3]), 0.01))) * grad(concatenation(broadcast(4.0), broadcast(2.127659574468085), broadcast(2.9850746268656714)) * Porosity times concentration [mol.m-3](Negative electrode porosity times concentration [mol.m-3], Separator porosity times concentration [mol.m-3], Positive electrode porosity times concentration [mol.m-3]))) - grad(Electrolyte potential [V](Negative electrolyte potential [V], Separator electrolyte potential [V], Positive electrolyte potential [V]))'], domains={'primary': ['negative electrode', 'separator', 'positive electrode'], 'secondary': ['current collector']})

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ...\pybamm\discretisations\discretisation.py:728, in Discretisation.process_symbol(self, symbol)
    727 try:
--> 728     return self._discretised_symbols[symbol]
    729 except KeyError:

KeyError: Subtraction(0x3831180c25b99da7, -, children=['(0.00012763994028092067 * ing_functio(time) / (maximum(concatenation(broadcast(4.0), broadcast(2.127659574468085), broadcast(2.9850746268656714)) * Porosity times concentration [mol.m-3](Negative electrode porosity times concentration [mol.m-3], Separator porosity times concentration [mol.m-3], Positive electrode porosity times concentration [mol.m-3]), 0.01))) * grad(concatenation(broadcast(4.0), broadcast(2.127659574468085), broadcast(2.9850746268656714)) * Porosity times concentration [mol.m-3](Negative electrode porosity times concentration [mol.m-3], Separator porosity times concentration [mol.m-3], Positive electrode porosity times concentration [mol.m-3]))', 'grad(Electrolyte potential [V](Negative electrolyte potential [V], Separator electrolyte potential [V], Positive electrolyte potential [V]))'], domains={'primary': ['negative electrode', 'separator', 'positive electrode'], 'secondary': ['current collector']})

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ...\pybamm\discretisations\discretisation.py:728, in Discretisation.process_symbol(self, symbol)
    727 try:
--> 728     return self._discretised_symbols[symbol]
    729 except KeyError:

KeyError: Multiplication(0x505505ba2e262de1, *, children=['0.00012763994028092067 * ing_functio(time) / (maximum(concatenation(broadcast(4.0), broadcast(2.127659574468085), broadcast(2.9850746268656714)) * Porosity times concentration [mol.m-3](Negative electrode porosity times concentration [mol.m-3], Separator porosity times concentration [mol.m-3], Positive electrode porosity times concentration [mol.m-3]), 0.01))', 'grad(concatenation(broadcast(4.0), broadcast(2.127659574468085), broadcast(2.9850746268656714)) * Porosity times concentration [mol.m-3](Negative electrode porosity times concentration [mol.m-3], Separator porosity times concentration [mol.m-3], Positive electrode porosity times concentration [mol.m-3]))'], domains={'primary': ['negative electrode', 'separator', 'positive electrode'], 'secondary': ['current collector']})

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ...\pybamm\discretisations\discretisation.py:728, in Discretisation.process_symbol(self, symbol)
    727 try:
--> 728     return self._discretised_symbols[symbol]
    729 except KeyError:

KeyError: Division(0xec2f03eefcb96a2, /, children=['0.00012763994028092067 * ing_functio(time)', 'maximum(concatenation(broadcast(4.0), broadcast(2.127659574468085), broadcast(2.9850746268656714)) * Porosity times concentration [mol.m-3](Negative electrode porosity times concentration [mol.m-3], Separator porosity times concentration [mol.m-3], Positive electrode porosity times concentration [mol.m-3]), 0.01)'], domains={'primary': ['negative electrode', 'separator', 'positive electrode'], 'secondary': ['current collector']})

During handling of the above exception, another exception occurred:

AttributeError                            Traceback (most recent call last)
File ...\pybamm\expression_tree\symbol.py:767, in Symbol.evaluate_for_shape(self)
    766 try:
--> 767     return self._saved_evaluate_for_shape
    768 except AttributeError:

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

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
File ...\pybamm\expression_tree\symbol.py:975, in Symbol.test_shape(self)
    974 try:
--> 975     self.shape_for_testing
    976 except ValueError as e:

File ...\pybamm\expression_tree\symbol.py:951, in Symbol.shape_for_testing(self)
    946 """
    947 Shape of an object for cases where it cannot be evaluated directly. If a symbol
    948 cannot be evaluated directly (e.g. it is a `Variable` or `Parameter`), it is
    949 instead given an arbitrary domain-dependent shape.
    950 """
--> 951 evaluated_self = self.evaluate_for_shape()
    952 if isinstance(evaluated_self, numbers.Number):

File ...\pybamm\expression_tree\symbol.py:769, in Symbol.evaluate_for_shape(self)
    768 except AttributeError:
--> 769     self._saved_evaluate_for_shape = self._evaluate_for_shape()
    770     return self._saved_evaluate_for_shape

File ...\pybamm\expression_tree\binary_operators.py:122, in BinaryOperator._evaluate_for_shape(self)
    121 right = self.children[1].evaluate_for_shape()
--> 122 return self._binary_evaluate(left, right)

File ...\pybamm\expression_tree\binary_operators.py:363, in Division._binary_evaluate(self, left, right)
    362 else:
--> 363     return left / right

ValueError: operands could not be broadcast together with shapes (41,1) (60,1) 

During handling of the above exception, another exception occurred:

ShapeError                                Traceback (most recent call last)
Cell In[1], line 29
     16 parameter_values.update(
     17     {
     18         "Initial temperature [K]": 298.15,
     19         "Ambient temperature [K]": temp_drive_cycle
     20     }
     21 )
     23 sim = pybamm.Simulation(
     24     model=pybamm.lithium_ion.DFN(),
     25     parameter_values=parameter_values,
     26     experiment=experiment
     27 )
---> 29 sol = sim.solve()
     30 sol.plot()
     32 plt.plot(sol["Time [s]"].entries, sol["Volume-averaged cell temperature [K]"].entries - 273.15)

File ...\pybamm\simulation.py:618, in Simulation.solve(self, t_eval, solver, check_model, save_at_cycles, calc_esoh, starting_solution, initial_soc, callbacks, showprogress, **kwargs)
    616 elif self.operating_mode == "with experiment":
    617     callbacks.on_experiment_start(logs)
--> 618     self.build_for_experiment(check_model=check_model, initial_soc=initial_soc)
    619     if t_eval is not None:
    620         pybamm.logger.warning(
    621             "Ignoring t_eval as solution times are specified by the experiment"
    622         )

File ...\pybamm\simulation.py:462, in Simulation.build_for_experiment(self, check_model, initial_soc)
    455 self.op_conds_to_built_solvers = {}
    456 for (
    457     op_cond,
    458     model_with_set_params,
    459 ) in self.experiment_unique_steps_to_model.items():
    460     # It's ok to modify the model with set parameters in place as it's
    461     # not returned anywhere
--> 462     built_model = self._disc.process_model(
    463         model_with_set_params, inplace=True, check_model=check_model
    464     )
    465     solver = self._solver.copy()
    466     self.op_conds_to_built_solvers[op_cond] = solver

File ...\pybamm\discretisations\discretisation.py:214, in Discretisation.process_model(self, model, inplace, check_model, remove_independent_variables_from_rhs)
    210 # Discretise variables (applying boundary conditions)
    211 # Note that we **do not** discretise the keys of model.rhs,
    212 # model.initial_conditions and model.boundary_conditions
    213 pybamm.logger.verbose("Discretise variables for {}".format(model.name))
--> 214 model_disc.variables = self.process_dict(model.variables)
    216 # Process parabolic and elliptic equations
    217 pybamm.logger.verbose("Discretise model equations for {}".format(model.name))

File ...\pybamm\discretisations\discretisation.py:698, in Discretisation.process_dict(self, var_eqn_dict, ics)
    695         eqn = pybamm.FullBroadcast(eqn, broadcast_domains=eqn_key.domains)
    697 pybamm.logger.debug("Discretise {!r}".format(eqn_key))
--> 698 processed_eqn = self.process_symbol(eqn)
    699 # Calculate scale if the key has a scale
    700 scale = getattr(eqn_key, "scale", 1)

File ...\pybamm\discretisations\discretisation.py:730, in Discretisation.process_symbol(self, symbol)
    728     return self._discretised_symbols[symbol]
    729 except KeyError:
--> 730     discretised_symbol = self._process_symbol(symbol)
    731     self._discretised_symbols[symbol] = discretised_symbol
    732     discretised_symbol.test_shape()

File ...\pybamm\discretisations\discretisation.py:766, in Discretisation._process_symbol(self, symbol)
    764 left, right = symbol.children
    765 disc_left = self.process_symbol(left)
--> 766 disc_right = self.process_symbol(right)
    767 if symbol.domain == []:
    768     return pybamm.simplify_if_constant(
    769         symbol._binary_new_copy(disc_left, disc_right)
    770     )

File ...\pybamm\discretisations\discretisation.py:730, in Discretisation.process_symbol(self, symbol)
    728     return self._discretised_symbols[symbol]
    729 except KeyError:
--> 730     discretised_symbol = self._process_symbol(symbol)
    731     self._discretised_symbols[symbol] = discretised_symbol
    732     discretised_symbol.test_shape()

File ...\pybamm\discretisations\discretisation.py:765, in Discretisation._process_symbol(self, symbol)
    762 if isinstance(symbol, pybamm.BinaryOperator):
    763     # Pre-process children
    764     left, right = symbol.children
--> 765     disc_left = self.process_symbol(left)
    766     disc_right = self.process_symbol(right)
    767     if symbol.domain == []:

File ...\pybamm\discretisations\discretisation.py:730, in Discretisation.process_symbol(self, symbol)
    728     return self._discretised_symbols[symbol]
    729 except KeyError:
--> 730     discretised_symbol = self._process_symbol(symbol)
    731     self._discretised_symbols[symbol] = discretised_symbol
    732     discretised_symbol.test_shape()

File ...\pybamm\discretisations\discretisation.py:765, in Discretisation._process_symbol(self, symbol)
    762 if isinstance(symbol, pybamm.BinaryOperator):
    763     # Pre-process children
    764     left, right = symbol.children
--> 765     disc_left = self.process_symbol(left)
    766     disc_right = self.process_symbol(right)
    767     if symbol.domain == []:

File ...\pybamm\discretisations\discretisation.py:732, in Discretisation.process_symbol(self, symbol)
    730 discretised_symbol = self._process_symbol(symbol)
    731 self._discretised_symbols[symbol] = discretised_symbol
--> 732 discretised_symbol.test_shape()
    734 # Assign mesh as an attribute to the processed variable
    735 if symbol.domain != []:

File ...\pybamm\expression_tree\symbol.py:977, in Symbol.test_shape(self)
    975     self.shape_for_testing
    976 except ValueError as e:
--> 977     raise pybamm.ShapeError("Cannot find shape (original error: {})".format(e))

ShapeError: Cannot find shape (original error: operands could not be broadcast together with shapes (41,1) (60,1) )

[ejfdickinson]

@rtimms Feels like this is a side effect of #3257 where an unidentified breaking change has crept in. Is it the case that in this parameter definition pattern, ambient temperature now has to be f(y,z,t), and can't just be f(t)?

[rtimms]

Correct, you need to wrap the Interpolant in a function of (y,z,t)

[rtimms]

This should have been documented as a breaking change. I don't know if there is some nice way to automatically detect if it is just a function of a single variable, in which case it is assumed to be a function of time only.

[ejfdickinson]

@rtimms If the thermal model is "isothermal" or "lumped", isn't every f(y,z,[vars]) by definition f([vars])?

It feels like this drives required usage in the opposite direction from the geometry-agnosticism defined for these thermal models, also in #3257.

[rtimms]

Yes, that is a fair point. At the moment, every FunctionParameter is defined as a function of the variables it could be a function of, in any model/geometry, so we just define the function once in the parameters class. This is ok for variables that depend on other variables, e.g. D(c,T), but I agree shouldn't be the case for variables that depend on space and time. We should define these differently based on the model geometry so that we can throw an explicit and informative error.

[rtimms]

Historically, most of the thermal stuff in PyBaMM was developed for single-layer pouch cell models, so there are lots of hangovers from design choices based on that.

[ejfdickinson]

I'd feel it's most Pythonic to put optional arguments last, and explicitly allow them to be optional. In general we don't want to have to write functions whose return values don't depend on their arguments, right?

What's happened here is that you've enabled complicated behaviour by requiring simple behaviour to have a complicated calling pattern!

[ejfdickinson]

Correct, you need to wrap the Interpolant in a function of (y,z,t)

Sorry, can you offer example code? I tried substituting the interpolant with a callable as follows, but got the same ShapeError:

def temp_drive_cycle(y,z,t):
    return pybamm.Interpolant(
        times,
        298.15 + 20 * (times / tmax),
        t
    )

and got the same error.

[rtimms]

#
# Example showing how to solve the DFN with a varying ambient temperature
#

import pybamm
import numpy as np

pybamm.set_logging_level("DEBUG")


# load model
options = {"thermal": "lumped"}
model = pybamm.lithium_ion.DFN(options)

# create geometry
geometry = model.default_geometry

# load parameter values and process model and geometry
times = np.arange(0, 1810, 10)
tmax = max(times)


def ambient_temperature(y, z, t):
    return pybamm.Interpolant(times, 298.15 + 20 * (times / tmax), pybamm.t)


param = model.default_parameter_values
param.update(
    {"Ambient temperature [K]": ambient_temperature}, check_already_exists=False
)
param.process_model(model)
param.process_geometry(geometry)

# set mesh
var_pts = {"x_n": 30, "x_s": 30, "x_p": 30, "r_n": 10, "r_p": 10}
mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)

# discretise model
disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
disc.process_model(model)

# solve model
t_eval = np.linspace(0, 3600 / 2, 100)
solver = pybamm.CasadiSolver(mode="fast", atol=1e-6, rtol=1e-3)
solution = solver.solve(model, t_eval)

# plot
plot = pybamm.QuickPlot(
    solution, ["X-averaged cell temperature [K]", "Ambient temperature [K]"]
)
plot.dynamic_plot()

[rtimms]

That works, but I get an error with your example. Will look into it.

[ejfdickinson]

@rtimms Any further insight into why the fix is not effective for the test case I suggested?

[rtimms]

@ejfdickinson I'll make it a priority to look this week

[rtimms]

@ejfdickinson I've identified the bug and am working on a fix. For now, a very hacky workaround is to make sure y or z is in the returned function so that the discretisation identifies the shape correctly. E.g.

import pybamm
import numpy as np

times = np.arange(0, 1810, 10)
tmax = max(times)
experiment = pybamm.Experiment([f"Discharge at 1C for {tmax} s"])


def temp_drive_cycle(y, z, t):
    return 1e-32 * y + pybamm.Interpolant(times, 298.15 + 20 * (times / tmax), t)


parameter_values = pybamm.ParameterValues("Chen2020")
parameter_values.update(
    {"Initial temperature [K]": 298.15, "Ambient temperature [K]": temp_drive_cycle}
)

sim = pybamm.Simulation(
    model=pybamm.lithium_ion.DFN(),
    parameter_values=parameter_values,
    experiment=experiment,
)

sol = sim.solve()
sol.plot(["Voltage [V]", "Volume-averaged cell temperature [K]"])