explicitly create interpolants from functions

CHANGELOG.md

@@ -4,6 +4,7 @@
 
 ## Features
 
+- Updated parameter sets so that interpolants are created explicitly in the parameter set python file. This does not change functionality but allows finer control, e.g. specifying a "cubic" interpolator instead of the default "linear" ([#2510](https://github.com/pybamm-team/PyBaMM/pull/2510))
 - Equivalent circuit models ([#2478](https://github.com/pybamm-team/PyBaMM/pull/2478))
 - New Idaklu solver options for jacobian type and linear solver, support Sundials v6 ([#2444](https://github.com/pybamm-team/PyBaMM/pull/2444))
 - Added `scale` and `reference` attributes to `Variable` objects, which can be use to make the ODE/DAE solver better conditioned ([#2440](https://github.com/pybamm-team/PyBaMM/pull/2440))

pybamm/input/parameters/ecm/example_set.py

@@ -5,13 +5,38 @@
 
 path, _ = os.path.split(os.path.abspath(__file__))
 
-ocv = pybamm.parameters.process_1D_data("ecm_example_ocv.csv", path=path)
+ocv_data = pybamm.parameters.process_1D_data("ecm_example_ocv.csv", path=path)
 
-r0 = pybamm.parameters.process_3D_data_csv("ecm_example_r0.csv", path=path)
-r1 = pybamm.parameters.process_3D_data_csv("ecm_example_r1.csv", path=path)
-c1 = pybamm.parameters.process_3D_data_csv("ecm_example_c1.csv", path=path)
+r0_data = pybamm.parameters.process_3D_data_csv("ecm_example_r0.csv", path=path)
+r1_data = pybamm.parameters.process_3D_data_csv("ecm_example_r1.csv", path=path)
+c1_data = pybamm.parameters.process_3D_data_csv("ecm_example_c1.csv", path=path)
 
-dUdT = pybamm.parameters.process_2D_data_csv("ecm_example_dudt.csv", path=path)
+dUdT_data = pybamm.parameters.process_2D_data_csv("ecm_example_dudt.csv", path=path)
+
+
+def ocv(sto):
+    name, (x, y) = ocv_data
+    return pybamm.Interpolant(x, y, sto, name)
+
+
+def r0(T_cell, current, soc):
+    name, (x, y) = r0_data
+    return pybamm.Interpolant(x, y, [T_cell, current, soc], name)
+
+
+def r1(T_cell, current, soc):
+    name, (x, y) = r1_data
+    return pybamm.Interpolant(x, y, [T_cell, current, soc], name)
+
+
+def c1(T_cell, current, soc):
+    name, (x, y) = c1_data
+    return pybamm.Interpolant(x, y, [T_cell, current, soc], name)
+
+
+def dUdT(ocv, T_cell):
+    name, (x, y) = dUdT_data
+    return pybamm.Interpolant(x, y, [ocv, T_cell], name)
 
 
 def get_parameter_values():

pybamm/input/parameters/lithium_ion/Ai2020.py

@@ -498,10 +498,22 @@ def electrolyte_conductivity_Ai2020(c_e, T):
 
 # Load data in the appropriate format
 path, _ = os.path.split(os.path.abspath(__file__))
-graphite_ocp_Enertech_Ai2020 = pybamm.parameters.process_1D_data(
+graphite_ocp_Enertech_Ai2020_data = pybamm.parameters.process_1D_data(
     "graphite_ocp_Enertech_Ai2020.csv", path=path
 )
-lico2_ocp_Ai2020 = pybamm.parameters.process_1D_data("lico2_ocp_Ai2020.csv", path=path)
+lico2_ocp_Ai2020_data = pybamm.parameters.process_1D_data(
+    "lico2_ocp_Ai2020.csv", path=path
+)
+
+
+def graphite_ocp_Enertech_Ai2020(sto):
+    name, (x, y) = graphite_ocp_Enertech_Ai2020_data
+    return pybamm.Interpolant(x, y, sto, name=name, interpolator="cubic")
+
+
+def lico2_ocp_Ai2020(sto):
+    name, (x, y) = lico2_ocp_Ai2020_data
+    return pybamm.Interpolant(x, y, sto, name=name, interpolator="cubic")
 
 
 # Call dict via a function to avoid errors when editing in place

pybamm/input/parameters/lithium_ion/Chen2020_composite.py

@@ -308,11 +308,16 @@ def electrolyte_conductivity_Nyman2008(c_e, T):
 
 # Load data in the appropriate format
 path, _ = os.path.split(os.path.abspath(__file__))
-graphite_ocp_Enertech_Ai2020 = pybamm.parameters.process_1D_data(
+graphite_ocp_Enertech_Ai2020_data = pybamm.parameters.process_1D_data(
     "graphite_ocp_Enertech_Ai2020.csv", path=path
 )
 
 
+def graphite_ocp_Enertech_Ai2020(sto):
+    name, (x, y) = graphite_ocp_Enertech_Ai2020_data
+    return pybamm.Interpolant(x, y, sto, name=name, interpolator="cubic")
+
+
 # Call dict via a function to avoid errors when editing in place
 def get_parameter_values():
     """

pybamm/input/parameters/lithium_ion/Ecker2015.py

@@ -1,5 +1,4 @@
 import pybamm
-import os
 
 
 def graphite_diffusivity_Ecker2015(sto, T):
@@ -40,7 +39,7 @@ def graphite_diffusivity_Ecker2015(sto, T):
     return D_ref * arrhenius
 
 
-def graphite_ocp_Ecker2015_function(sto):
+def graphite_ocp_Ecker2015(sto):
     """
     Graphite OCP as a function of stochiometry [1, 2, 3].
 
@@ -190,7 +189,7 @@ def nco_diffusivity_Ecker2015(sto, T):
     return D_ref * arrhenius
 
 
-def nco_ocp_Ecker2015_function(sto):
+def nco_ocp_Ecker2015(sto):
     """
     NCO OCP as a function of stochiometry [1, 2, 3].
 
@@ -383,22 +382,6 @@ def electrolyte_conductivity_Ecker2015(c_e, T):
     return sigma_e
 
 
-# Load data in the appropriate format
-path, _ = os.path.split(os.path.abspath(__file__))
-measured_graphite_diffusivity_Ecker2015 = pybamm.parameters.process_1D_data(
-    "measured_graphite_diffusivity_Ecker2015.csv", path=path
-)
-graphite_ocp_Ecker2015 = pybamm.parameters.process_1D_data(
-    "graphite_ocp_Ecker2015.csv", path=path
-)
-measured_nco_diffusivity_Ecker2015 = pybamm.parameters.process_1D_data(
-    "measured_nco_diffusivity_Ecker2015.csv", path=path
-)
-nco_ocp_Ecker2015 = pybamm.parameters.process_1D_data(
-    "nco_ocp_Ecker2015.csv", path=path
-)
-
-
 # Call dict via a function to avoid errors when editing in place
 def get_parameter_values():
     """
@@ -518,11 +501,8 @@ def get_parameter_values():
         # negative electrode
         "Negative electrode conductivity [S.m-1]": 14.0,
         "Maximum concentration in negative electrode [mol.m-3]": 31920.0,
-        "Measured negative electrode diffusivity [m2.s-1]"
-        "": measured_graphite_diffusivity_Ecker2015,
         "Negative electrode diffusivity [m2.s-1]": graphite_diffusivity_Ecker2015,
-        "Measured negative electrode OCP [V]": graphite_ocp_Ecker2015,
-        "Negative electrode OCP [V]": graphite_ocp_Ecker2015_function,
+        "Negative electrode OCP [V]": graphite_ocp_Ecker2015,
         "Negative electrode porosity": 0.329,
         "Negative electrode active material volume fraction": 0.372403,
         "Negative particle radius [m]": 1.37e-05,
@@ -539,11 +519,8 @@ def get_parameter_values():
         # positive electrode
         "Positive electrode conductivity [S.m-1]": 68.1,
         "Maximum concentration in positive electrode [mol.m-3]": 48580.0,
-        "Measured positive electrode diffusivity [m2.s-1]"
-        "": measured_nco_diffusivity_Ecker2015,
         "Positive electrode diffusivity [m2.s-1]": nco_diffusivity_Ecker2015,
-        "Measured positive electrode OCP [V]": nco_ocp_Ecker2015,
-        "Positive electrode OCP [V]": nco_ocp_Ecker2015_function,
+        "Positive electrode OCP [V]": nco_ocp_Ecker2015,
         "Positive electrode porosity": 0.296,
         "Positive electrode active material volume fraction": 0.40832,
         "Positive particle radius [m]": 6.5e-06,

pybamm/input/parameters/lithium_ion/NCA_Kim2011.py

@@ -267,6 +267,11 @@ def electrolyte_conductivity_Kim2011(c_e, T):
 )
 
 
+def nca_ocp_Kim2011(sto):
+    name, (x, y) = nca_ocp_Kim2011_data
+    return pybamm.Interpolant(x, y, sto, name=name, interpolator="linear")
+
+
 # Call dict via a function to avoid errors when editing in place
 def get_parameter_values():
     """
@@ -388,7 +393,7 @@ def get_parameter_values():
         "Positive electrode conductivity [S.m-1]": 10.0,
         "Maximum concentration in positive electrode [mol.m-3]": 49000.0,
         "Positive electrode diffusivity [m2.s-1]": nca_diffusivity_Kim2011,
-        "Positive electrode OCP [V]": nca_ocp_Kim2011_data,
+        "Positive electrode OCP [V]": nca_ocp_Kim2011,
         "Positive electrode porosity": 0.4,
         "Positive electrode active material volume fraction": 0.41,
         "Positive particle radius [m]": 1.633e-06,

pybamm/input/parameters/lithium_ion/OKane2022.py

@@ -487,11 +487,16 @@ def electrolyte_conductivity_Nyman2008_arrhenius(c_e, T):
 
 # Load data in the appropriate format
 path, _ = os.path.split(os.path.abspath(__file__))
-graphite_LGM50_ocp_Chen2020 = pybamm.parameters.process_1D_data(
+graphite_LGM50_ocp_Chen2020_data = pybamm.parameters.process_1D_data(
     "graphite_LGM50_ocp_Chen2020.csv", path=path
 )
 
 
+def graphite_LGM50_ocp_Chen2020(sto):
+    name, (x, y) = graphite_LGM50_ocp_Chen2020_data
+    return pybamm.Interpolant(x, y, sto, name=name, interpolator="cubic")
+
+
 # Call dict via a function to avoid errors when editing in place
 def get_parameter_values():
     """

pybamm/models/full_battery_models/equivalent_circuit/ecm_model_options.py

@@ -2,8 +2,8 @@
 
 
 class NaturalNumberOption:
-    def __init__(self, defualt_value):
-        self.value = defualt_value
+    def __init__(self, default_value):
+        self.value = default_value
 
     def __contains__(self, value):
         is_an_integer = isinstance(value, int)

pybamm/models/full_battery_models/lithium_ion/electrode_soh.py

@@ -277,19 +277,24 @@ def solve(self, inputs):
         return sol
 
     def _set_up_solve(self, inputs):
+        # Try with full sim
         sim = self._get_electrode_soh_sims_full()
-        x0_min, x100_max, _, _ = self._get_lims(inputs)
-
-        x100_init = x100_max
-        x0_init = x0_min
         if sim.solution is not None:
-            # Update the initial conditions if they are valid
-            x100_init_sol = sim.solution["x_100"].data[0]
-            if x0_min < x100_init_sol < x100_max:
-                x100_init = x100_init_sol
-            x0_init_sol = sim.solution["x_0"].data[0]
-            if x0_min < x0_init_sol < x100_max:
-                x0_init = x0_init_sol
+            x100_sol = sim.solution["x_100"].data
+            x0_sol = sim.solution["x_0"].data
+            return {"x_100": x100_sol, "x_0": x0_sol}
+
+        # Try with split sims
+        x100_sim, x0_sim = self._get_electrode_soh_sims_split()
+        if x100_sim.solution is not None and x0_sim.solution is not None:
+            x100_sol = x100_sim.solution["x_100"].data
+            x0_sol = x0_sim.solution["x_0"].data
+            return {"x_100": x100_sol, "x_0": x0_sol}
+
+        # Fall back to initial conditions calculated from limits
+        x0_min, x100_max, _, _ = self._get_lims(inputs)
+        x100_init = min(x100_max, 0.8)
+        x0_init = max(x0_min, 0.2)
         return {"x_100": np.array(x100_init), "x_0": np.array(x0_init)}
 
     def _solve_full(self, inputs, ics):

pybamm/parameters/parameter_values.py

@@ -644,7 +644,10 @@ def _process_symbol(self, symbol):
                     # For parameters provided as data we use a cubic interpolant
                     # Note: the cubic interpolant can be differentiated
                     function = pybamm.Interpolant(
-                        input_data[0], input_data[-1], new_children, name=name
+                        input_data[0],
+                        input_data[-1],
+                        new_children,
+                        name=name,
                     )
 
                 else:  # pragma: no cover

tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_initial_soc.py

@@ -8,7 +8,16 @@
 class TestInitialSOC(unittest.TestCase):
     def test_interpolant_parameter_sets(self):
         model = pybamm.lithium_ion.SPM()
-        for param in ["OKane2022", "Ai2020"]:
+        params = [
+            "Ai2020",
+            "Chen2020",
+            "Ecker2015",
+            "Marquis2019",
+            "Mohtat2020",
+            "OKane2022",
+            "ORegan2022",
+        ]
+        for param in params:
             with self.subTest(param=param):
                 parameter_values = pybamm.ParameterValues(param)
                 sim = pybamm.Simulation(model=model, parameter_values=parameter_values)

tests/unit/test_models/test_full_battery_models/test_equivalent_circuit/test_thevenin.py

@@ -102,6 +102,7 @@ def test_get_default_parameters(self):
         model = pybamm.equivalent_circuit.Thevenin()
         values = model.default_parameter_values
         self.assertIn("Initial SoC", list(values.keys()))
+        values.process_model(model)
 
     def test_get_default_quick_plot_variables(self):
         model = pybamm.equivalent_circuit.Thevenin()

tests/unit/test_solvers/test_casadi_solver.py

@@ -509,7 +509,7 @@ def test_interpolant_extrapolate(self):
         model = pybamm.lithium_ion.DFN()
         param = pybamm.ParameterValues("NCA_Kim2011")
         experiment = pybamm.Experiment(
-            ["Charge at 1C until 4.6 V"], period="10 seconds"
+            ["Charge at 1C until 4.2 V"], period="10 seconds"
         )
 
         param["Upper voltage cut-off [V]"] = 4.8
@@ -528,18 +528,6 @@ def test_interpolant_extrapolate(self):
         with self.assertRaisesRegex(pybamm.SolverError, "interpolation bounds"):
             sim.solve()
 
-        ci = param["Initial concentration in positive electrode [mol.m-3]"]
-        param["Initial concentration in positive electrode [mol.m-3]"] = 0.8 * ci
-
-        sim = pybamm.Simulation(
-            model,
-            parameter_values=param,
-            experiment=experiment,
-            solver=pybamm.CasadiSolver(mode="safe", dt_max=0.05),
-        )
-        with self.assertRaisesRegex(pybamm.SolverError, "interpolation bounds"):
-            sim.solve()
-
     def test_casadi_safe_no_termination(self):
         model = pybamm.BaseModel()
         v = pybamm.Variable("v")