Issue 1310 interpolant

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added new functionality for `Interpolant` ([#1312](https://github.com/pybamm-team/PyBaMM/pull/1312))
 -   Added option to express experiments (and extract solutions) in terms of cycles of operating condition ([#1309](https://github.com/pybamm-team/PyBaMM/pull/1309))
 -   Reformatted the `BasicDFNHalfCell` to be consistent with the other models ([#1282](https://github.com/pybamm-team/PyBaMM/pull/1282))
 -   Added option to make the total interfacial current density a state ([#1280](https://github.com/pybamm-team/PyBaMM/pull/1280))
@@ -19,6 +20,7 @@
 
 ## Breaking changes
 
+-   `Interpolant` now takes `x` and `y` instead of a single `data` entry ([#1312](https://github.com/pybamm-team/PyBaMM/pull/1312))
 -   Boolean model options ('sei porosity change', 'convection') must now be given in string format ('true' or 'false' instead of True or False) ([#1280](https://github.com/pybamm-team/PyBaMM/pull/1280))
 -   Operations such as `1*x` and `0+x` now directly return `x`. This can be bypassed by explicitly creating the binary operators, e.g. `pybamm.Multiplication(1, x)` ([#1252](https://github.com/pybamm-team/PyBaMM/pull/1252))
 

examples/notebooks/Getting Started/Tutorial 4 - Setting parameter values.ipynb

@@ -887,7 +887,7 @@
     "\n",
     "# Create interpolant\n",
     "timescale = parameter_values.evaluate(model.timescale)\n",
-    "current_interpolant = pybamm.Interpolant(drive_cycle, timescale * pybamm.t)\n",
+    "current_interpolant = pybamm.Interpolant(drive_cycle[:, 0], drive_cycle[:, 1], timescale * pybamm.t)\n",
     "\n",
     "# Set drive cycle\n",
     "parameter_values[\"Current function [A]\"] = current_interpolant"

examples/notebooks/change-input-current.ipynb

@@ -158,7 +158,7 @@
     "\n",
     "# create interpolant - must be a function of *dimensional* time\n",
     "timescale = param.evaluate(model.timescale)\n",
-    "current_interpolant = pybamm.Interpolant(drive_cycle, timescale * pybamm.t)\n",
+    "current_interpolant = pybamm.Interpolant(drive_cycle[:, 0], drive_cycle[:, 1], timescale * pybamm.t)\n",
     "\n",
     "# set drive cycle\n",
     "param[\"Current function [A]\"] = current_interpolant\n",

examples/notebooks/initialize-model-with-solution.ipynb

@@ -82,7 +82,7 @@
     "# create interpolant\n",
     "param = model.default_parameter_values\n",
     "timescale = param.evaluate(model.timescale)\n",
-    "current_interpolant = pybamm.Interpolant(drive_cycle, timescale * pybamm.t)\n",
+    "current_interpolant = pybamm.Interpolant(drive_cycle[:, 0], drive_cycle[:, 1], timescale * pybamm.t)\n",
     "# set drive cycle\n",
     "param[\"Current function [A]\"] = current_interpolant"
    ]

examples/scripts/compare_comsol/compare_comsol_DFN.py

@@ -1,5 +1,4 @@
 import pybamm
-import numpy as np
 import os
 import pickle
 import scipy.interpolate as interp
@@ -78,29 +77,17 @@ def get_interp_fun(variable_name, domain):
         comsol_x = comsol_variables["x_p"]
     elif domain == whole_cell:
         comsol_x = comsol_variables["x"]
+
     # Make sure to use dimensional space
     pybamm_x = mesh.combine_submeshes(*domain).nodes * L_x
     variable = interp.interp1d(comsol_x, variable, axis=0)(pybamm_x)
 
-    def myinterp(t):
-        try:
-            return interp.interp1d(
-                comsol_t, variable, fill_value="extrapolate", bounds_error=False
-            )(t)[:, np.newaxis]
-        except ValueError as err:
-            raise ValueError(
-                (
-                    "Failed to interpolate '{}' with time range [{}, {}] at time {}."
-                    + "Original error: {}"
-                ).format(variable_name, comsol_t[0], comsol_t[-1], t, err)
-            )
-
-    # Make sure to use dimensional time
-    fun = pybamm.Function(
-        myinterp,
+    fun = pybamm.Interpolant(
+        comsol_t,
+        variable.T,
         pybamm.t * pybamm_model.timescale.evaluate(),
-        name=variable_name + "_comsol",
     )
+
     fun.domain = domain
     fun.mesh = mesh.combine_submeshes(*domain)
     fun.secondary_mesh = None
@@ -113,15 +100,12 @@ def myinterp(t):
 comsol_phi_n = get_interp_fun("phi_n", ["negative electrode"])
 comsol_phi_e = get_interp_fun("phi_e", whole_cell)
 comsol_phi_p = get_interp_fun("phi_p", ["positive electrode"])
-comsol_voltage = pybamm.Function(
-    interp.interp1d(
-        comsol_t,
-        comsol_variables["voltage"],
-        fill_value="extrapolate",
-        bounds_error=False,
-    ),
+comsol_voltage = pybamm.Interpolant(
+    comsol_t,
+    comsol_variables["voltage"],
     pybamm.t * pybamm_model.timescale.evaluate(),
 )
+
 comsol_voltage.mesh = None
 comsol_voltage.secondary_mesh = None
 

examples/scripts/drive_cycle.py

@@ -21,7 +21,9 @@
 
 # create interpolant
 timescale = param.evaluate(model.timescale)
-current_interpolant = pybamm.Interpolant(drive_cycle, timescale * pybamm.t)
+current_interpolant = pybamm.Interpolant(
+    drive_cycle[:, 0], drive_cycle[:, 1], timescale * pybamm.t
+)
 
 # set drive cycle
 param["Current function [A]"] = current_interpolant

pybamm/discretisations/discretisation.py

@@ -787,8 +787,8 @@ def process_dict(self, var_eqn_dict):
         """
         new_var_eqn_dict = {}
         for eqn_key, eqn in var_eqn_dict.items():
-            # Broadcast if the equation evaluates to a number(e.g. Scalar)
-            if eqn.evaluates_to_number() and not isinstance(eqn_key, str):
+            # Broadcast if the equation evaluates to a number (e.g. Scalar)
+            if np.prod(eqn.shape_for_testing) == 1 and not isinstance(eqn_key, str):
                 eqn = pybamm.FullBroadcast(
                     eqn, eqn_key.domain, eqn_key.auxiliary_domains
                 )

pybamm/expression_tree/interpolant.py

@@ -2,6 +2,7 @@
 # Interpolating class
 #
 import pybamm
+import numpy as np
 from scipy import interpolate
 
 
@@ -11,11 +12,13 @@ class Interpolant(pybamm.Function):
 
     Parameters
     ----------
-    data : :class:`numpy.ndarray`
-        Numpy array of data to use for interpolation. Must have exactly two columns (x
-        and y data)
-    child : :class:`pybamm.Symbol`
-        Node to use when evaluating the interpolant
+    x : iterable of :class:`numpy.ndarray`
+        1-D array(s) of real values defining the data point coordinates.
+    y : :class:`numpy.ndarray`
+        The values of the function to interpolate at the data points.
+    children : iterable of :class:`pybamm.Symbol`
+        Node(s) to use when evaluating the interpolant. Each child corresponds to an
+        entry of x
     name : str, optional
         Name of the interpolant. Default is None, in which case the name "interpolating
         function" is given.
@@ -32,28 +35,73 @@ class Interpolant(pybamm.Function):
 
     def __init__(
         self,
-        data,
-        child,
+        x,
+        y,
+        children,
         name=None,
-        interpolator="cubic spline",
+        interpolator=None,
         extrapolate=True,
         entries_string=None,
     ):
-        if data.ndim != 2 or data.shape[1] != 2:
-            raise ValueError(
-                """
-                data should have exactly two columns (x and y) but has shape {}
-                """.format(
-                    data.shape
+        if isinstance(x, (tuple, list)) and len(x) == 2:
+            interpolator = interpolator or "linear"
+            if interpolator != "linear":
+                raise ValueError(
+                    "interpolator should be 'linear' if x is two-dimensional"
                 )
+            x1, x2 = x
+            if y.ndim != 2:
+                raise ValueError("y should be two-dimensional if len(x)=2")
+        else:
+            interpolator = interpolator or "cubic spline"
+            if isinstance(x, (tuple, list)):
+                x1 = x[0]
+            else:
+                x1 = x
+                x = [x]
+            x2 = None
+        if x1.shape[0] != y.shape[0]:
+            raise ValueError(
+                "len(x1) should equal y=shape[0], "
+                "but x1.shape={} and y.shape={}".format(x1.shape, y.shape)
+            )
+        if x2 is not None and x2.shape[0] != y.shape[1]:
+            raise ValueError(
+                "len(x2) should equal y=shape[1], "
+                "but x2.shape={} and y.shape={}".format(x2.shape, y.shape)
             )
+        if isinstance(children, pybamm.Symbol):
+            children = [children]
+        # Either a single x is provided and there is one child
+        # or x is a 2-tuple and there are two children
+        if len(x) != len(children):
+            raise ValueError("len(x) should equal len(children)")
+        # if there is only one x, y can be 2-dimensional but the child must have
+        # length 1
+        if len(x) == 1 and y.ndim == 2 and children[0].size != 1:
+            raise ValueError(
+                "child should have size 1 if y is two-dimensional and len(x)==1"
+            )
+
+        if interpolator == "linear":
+            if len(x) == 1:
+                if extrapolate is False:
+                    interpolating_function = interpolate.interp1d(
+                        x1, y.T, bounds_error=False, fill_value=np.nan
+                    )
+                elif extrapolate is True:
+                    interpolating_function = interpolate.interp1d(
+                        x1, y.T, bounds_error=False, fill_value="extrapolate"
+                    )
+            elif len(x) == 2:
+                interpolating_function = interpolate.interp2d(x1, x2, y)
         elif interpolator == "pchip":
             interpolating_function = interpolate.PchipInterpolator(
-                data[:, 0], data[:, 1], extrapolate=extrapolate
+                x1, y, extrapolate=extrapolate
             )
         elif interpolator == "cubic spline":
             interpolating_function = interpolate.CubicSpline(
-                data[:, 0], data[:, 1], extrapolate=extrapolate
+                x1, y, extrapolate=extrapolate
             )
         else:
             raise ValueError("interpolator '{}' not recognised".format(interpolator))
@@ -62,14 +110,13 @@ def __init__(
             name = "interpolating function ({})".format(name)
         else:
             name = "interpolating function"
-        self.data = data
+        self.x = x
+        self.y = y
         self.entries_string = entries_string
         super().__init__(
-            interpolating_function, child, name=name, derivative="derivative"
+            interpolating_function, *children, name=name, derivative="derivative"
         )
         # Store information as attributes
-        self.x = data[:, 0]
-        self.y = data[:, 1]
         self.interpolator = interpolator
         self.extrapolate = extrapolate
 
@@ -85,8 +132,10 @@ def entries_string(self, value):
         if value is not None:
             self._entries_string = value
         else:
-            entries = self.data
-            self._entries_string = entries.tobytes()
+            self._entries_string = ""
+            for i, x in enumerate(self.x):
+                self._entries_string += "x" + str(i) + "_" + str(x.tobytes())
+            self._entries_string += "y_" + str(self.y.tobytes())
 
     def set_id(self):
         """ See :meth:`pybamm.Symbol.set_id()`. """
@@ -97,10 +146,21 @@ def set_id(self):
     def _function_new_copy(self, children):
         """ See :meth:`Function._function_new_copy()` """
         return pybamm.Interpolant(
-            self.data,
-            *children,
+            self.x,
+            self.y,
+            children,
             name=self.name,
             interpolator=self.interpolator,
             extrapolate=self.extrapolate,
-            entries_string=self.entries_string
+            entries_string=self.entries_string,
         )
+
+    def _function_evaluate(self, evaluated_children):
+        children_eval_flat = []
+        for child in evaluated_children:
+            if isinstance(child, np.ndarray):
+                children_eval_flat.append(child.flatten())
+            else:
+                children_eval_flat.append(child)
+
+        return self.function(*children_eval_flat).flatten()[:, np.newaxis]

pybamm/expression_tree/operations/convert_to_casadi.py

@@ -4,7 +4,6 @@
 import pybamm
 import casadi
 import numpy as np
-from scipy.interpolate import PchipInterpolator, CubicSpline
 from scipy import special
 
 
@@ -132,10 +131,10 @@ def _convert(self, symbol, t, y, y_dot, inputs):
                 return casadi.sign(*converted_children)
             elif symbol.function == special.erf:
                 return casadi.erf(*converted_children)
-            elif isinstance(symbol.function, (PchipInterpolator, CubicSpline)):
-                return casadi.interpolant("LUT", "bspline", [symbol.x], symbol.y)(
-                    *converted_children
-                )
+            elif isinstance(symbol, pybamm.Interpolant):
+                return casadi.interpolant(
+                    "LUT", "bspline", symbol.x, symbol.y.flatten()
+                )(*converted_children)
             elif symbol.function.__name__.startswith("elementwise_grad_of_"):
                 differentiating_child_idx = int(symbol.function.__name__[-1])
                 # Create dummy symbolic variables in order to differentiate using CasADi

pybamm/input/parameters/lithium-ion/anodes/graphite_Ai2020/graphite_entropy_Enertech_Ai2020_function.py

@@ -1,26 +1,25 @@
-def graphite_entropy_Enertech_Ai2020_function(sto, T):
+def graphite_entropy_Enertech_Ai2020_function(sto):
     """
-        Lithium Cobalt Oxide (LiCO2) entropic change in open circuit potential (OCP) at
-        a temperature of 298.15K as a function of the stochiometry. The fit is taken
-        from Ref [1], which is only accurate
-       for 0.43 < sto < 0.9936.
+    Lithium Cobalt Oxide (LiCO2) entropic change in open circuit potential (OCP) at
+    a temperature of 298.15K as a function of the stochiometry. The fit is taken
+    from Ref [1], which is only accurate
+    for 0.43 < sto < 0.9936.
 
-        References
-        ----------
-        .. [1] Ai, W., Kraft, L., Sturm, J., Jossen, A., & Wu, B. (2020).
-        Electrochemical Thermal-Mechanical Modelling of Stress Inhomogeneity in Lithium-Ion Pouch Cells. # noqa
-        Journal of The Electrochemical Society, 167(1), 013512. DOI: 10.1149/2.0122001JES # noqa
+    References
+    ----------
+    .. [1] Ai, W., Kraft, L., Sturm, J., Jossen, A., & Wu, B. (2020).
+    Electrochemical Thermal-Mechanical Modelling of Stress Inhomogeneity in Lithium-Ion Pouch Cells. # noqa
+    Journal of The Electrochemical Society, 167(1), 013512. DOI: 10.1149/2.0122001JES # noqa
 
     Parameters
     ----------
     sto: double
        Stochiometry of material (li-fraction)
-    T : :class:`pybamm.Symbol`
-        Temperature [K]
+
     Returns
     -------
     :class:`pybamm.Symbol`
-        Entropic change [v.K-1]
+        Entropic change [V.K-1]
     """
 
     du_dT = (
@@ -48,5 +47,5 @@ def graphite_entropy_Enertech_Ai2020_function(sto, T):
             + 165705.8597 * sto ** 8
         )
     )
-    # show no temperature dependence
+
     return du_dT

pybamm/input/parameters/lithium-ion/anodes/graphite_Ramadass2004/graphite_entropic_change_Moura2016.py

@@ -1,7 +1,7 @@
-from pybamm import exp, cosh
+from pybamm import exp, cosh, Parameter
 
 
-def graphite_entropic_change_Moura2016(sto, c_n_max):
+def graphite_entropic_change_Moura2016(sto):
     """
     Graphite entropic change in open circuit potential (OCP) at a temperature of
     298.15K as a function of the stochiometry taken from Scott Moura's FastDFN code
@@ -11,12 +11,13 @@ def graphite_entropic_change_Moura2016(sto, c_n_max):
     ----------
     .. [1] https://github.com/scott-moura/fastDFN
 
-      Parameters
-      ----------
-      sto : :class:`pybamm.Symbol`
-           Stochiometry of material (li-fraction)
+    Parameters
+    ----------
+    sto : :class:`pybamm.Symbol`
+        Stochiometry of material (li-fraction)
 
     """
+    c_n_max = Parameter("Maximum concentration in negative electrode [mol.m-3]")
 
     du_dT = (
         -1.5 * (120.0 / c_n_max) * exp(-120 * sto)