Use T_init in initial_state (#5189)

* use T_init in initial_state

* changelog

* Update src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_composite.py

Co-authored-by: Brady Planden <[email protected]>

---------

Co-authored-by: Brady Planden <[email protected]>

Author: rtimms

CHANGELOG.md

@@ -9,6 +9,7 @@
 
 ## Bug fixes
 
+- Fix a bug in setting initial stoichiometries where the reference temperature was used instead of the initial temperature. ([#5189](https://github.com/pybamm-team/PyBaMM/pull/5189))
 - Fix a bug in the calculation of "Bulk" OCP terms in hysteresis models ([#5169](https://github.com/pybamm-team/PyBaMM/pull/5169))
 - Fixed a bug where the final duration of a drive cycle would not be inferred correctly. ([#5153](https://github.com/pybamm-team/PyBaMM/pull/5153))
 

src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py

@@ -711,21 +711,21 @@ def get_initial_stoichiometries(self, initial_value, tol=1e-6, inputs=None):
             soc = pybamm.Variable("soc")
             x = x_0 + soc * (x_100 - x_0)
             y = y_0 - soc * (y_0 - y_100)
-            T_ref = parameter_values["Reference temperature [K]"]
+            T_init = parameter_values["Initial temperature [K]"]
             if self.options["open-circuit potential"] == "MSMR":
                 xn = self.param.n.prim.x
                 xp = self.param.p.prim.x
                 Up = pybamm.Variable("Up")
                 Un = pybamm.Variable("Un")
-                soc_model.algebraic[Up] = x - xn(Un, T_ref)
-                soc_model.algebraic[Un] = y - xp(Up, T_ref)
+                soc_model.algebraic[Up] = x - xn(Un, T_init)
+                soc_model.algebraic[Un] = y - xp(Up, T_init)
                 soc_model.initial_conditions[Un] = 0
                 soc_model.initial_conditions[Up] = V_max
                 soc_model.algebraic[soc] = Up - Un - V_init
             else:
                 Up = self.param.p.prim.U
                 Un = self.param.n.prim.U
-                soc_model.algebraic[soc] = Up(y, T_ref) - Un(x, T_ref) - V_init
+                soc_model.algebraic[soc] = Up(y, T_init) - Un(x, T_init) - V_init
             # initial guess for soc linearly interpolates between 0 and 1
             # based on V linearly interpolating between V_max and V_min
             soc_model.initial_conditions[soc] = (V_init - V_min) / (V_max - V_min)
@@ -821,9 +821,13 @@ def get_initial_ocps(self, initial_value, tol=1e-6, inputs=None):
             Un = sol["Un"].data[0]
             Up = sol["Up"].data[0]
         else:
-            T_ref = parameter_values["Reference temperature [K]"]
-            Un = parameter_values.evaluate(self.param.n.prim.U(x, T_ref), inputs=inputs)
-            Up = parameter_values.evaluate(self.param.p.prim.U(y, T_ref), inputs=inputs)
+            T_init = parameter_values["Initial temperature [K]"]
+            Un = parameter_values.evaluate(
+                self.param.n.prim.U(x, T_init), inputs=inputs
+            )
+            Up = parameter_values.evaluate(
+                self.param.p.prim.U(y, T_init), inputs=inputs
+            )
         return Un, Up
 
     def get_min_max_ocps(self, inputs=None):

src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_composite.py

@@ -130,6 +130,8 @@ def __init__(
         y_0_1 = variables["y_0_1"]
         V_max = param.voltage_high_cut
         V_min = param.voltage_low_cut
+        # Here we use T_ref as the stoichiometry limits are defined using the reference
+        # state
         if is_negative_composite:
             x_100_2 = variables["x_100_2"]
             x_0_2 = variables["x_0_2"]
@@ -167,10 +169,12 @@ def __init__(
         x_init_1 = variables["x_init_1"]
         y_init_1 = variables["y_init_1"]
         if initialization_method == "voltage":
+            # Here we use T_init so that the initial voltage is correct, including the
+            # contribution from the entropic change
             V_init = pybamm.InputParameter("V_init")
             self.algebraic[x_init_1] = (
-                param.p.prim.U(y_init_1, param.T_ref)
-                - param.n.prim.U(x_init_1, param.T_ref)
+                param.p.prim.U(y_init_1, param.T_init)
+                - param.n.prim.U(x_init_1, param.T_init)
                 - V_init
             )
             self.algebraic[y_init_1] = (
@@ -198,17 +202,21 @@ def __init__(
             )
         else:
             raise ValueError("Invalid initialization method")
-        # Add voltage equations for secondary phases (init)
+        # Add voltage equations for secondary phases (init), we use T_ref if setting
+        # based on SOC since the stoichiometry limits are defined using the reference
+        # state, and T_init if setting based on voltage since the entropic change is
+        # included in the voltage equation
+        T = param.T_init if initialization_method == "voltage" else param.T_ref
         if is_positive_composite:
             y_init_2 = variables["y_init_2"]
-            self.algebraic[y_init_2] = param.p.prim.U(
-                y_init_1, param.T_ref
-            ) - param.p.sec.U(y_init_2, param.T_ref)
+            self.algebraic[y_init_2] = param.p.prim.U(y_init_1, T) - param.p.sec.U(
+                y_init_2, T
+            )
         if is_negative_composite:
             x_init_2 = variables["x_init_2"]
-            self.algebraic[x_init_2] = param.n.prim.U(
-                x_init_1, param.T_ref
-            ) - param.n.sec.U(x_init_2, param.T_ref)
+            self.algebraic[x_init_2] = param.n.prim.U(x_init_1, T) - param.n.sec.U(
+                x_init_2, T
+            )
 
         self.variables.update(variables)
         if initialization_method == "SOC":

src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_half_cell.py

@@ -145,17 +145,17 @@ def get_initial_stoichiometry_half_cell(
         model = pybamm.BaseModel()
         x = pybamm.Variable("x")
         Up = param.p.prim.U
-        T_ref = parameter_values["Reference temperature [K]"]
+        T_init = parameter_values["Initial temperature [K]"]
         model.variables["x"] = x
-        model.algebraic[x] = Up(x, T_ref) - V_init
+        model.algebraic[x] = Up(x, T_init) - V_init
         # initial guess for x linearly interpolates between 0 and 1
         # based on V linearly interpolating between V_max and V_min
         soc_initial_guess = (V_init - V_min) / (V_max - V_min)
         model.initial_conditions[x] = 1 - soc_initial_guess
         if is_composite:
             Up_2 = param.p.sec.U
             x_2 = pybamm.Variable("x_2")
-            model.algebraic[x_2] = Up_2(x_2, T_ref) - V_init
+            model.algebraic[x_2] = Up_2(x_2, T_init) - V_init
             model.variables["x_2"] = x_2
             model.initial_conditions[x_2] = 1 - soc_initial_guess
 
@@ -177,6 +177,7 @@ def get_initial_stoichiometry_half_cell(
             x_2 = pybamm.Variable("x_2")
             U_p = param.p.prim.U
             U_p_2 = param.p.sec.U
+            # here we use T_ref and SOC 0 and 1 are defined using the reference state
             T_ref = parameter_values["Reference temperature [K]"]
             model.algebraic[x] = U_p(x, T_ref) - U_p_2(x_2, T_ref)
             model.initial_conditions[x] = x_0 + initial_value * (x_100 - x_0)

tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py

@@ -179,6 +179,12 @@ def _get_params_and_options(composite_electrode):
                     "Secondary: Negative electrode OCP [V]": params[
                         "Negative electrode OCP [V]"
                     ],
+                    "Primary: Negative electrode OCP entropic change [V.K-1]": params[
+                        "Negative electrode OCP entropic change [V.K-1]"
+                    ],
+                    "Secondary: Negative electrode OCP entropic change [V.K-1]": params[
+                        "Negative electrode OCP entropic change [V.K-1]"
+                    ],
                     "Primary: Negative electrode active material volume fraction": 0.5,
                     "Secondary: Negative electrode active material volume fraction": 0.5,
                     "Primary: Maximum concentration in negative electrode [mol.m-3]": params[
@@ -206,6 +212,12 @@ def _get_params_and_options(composite_electrode):
                     "Secondary: Positive electrode OCP [V]": params[
                         "Positive electrode OCP [V]"
                     ],
+                    "Primary: Positive electrode OCP entropic change [V.K-1]": params[
+                        "Positive electrode OCP entropic change [V.K-1]"
+                    ],
+                    "Secondary: Positive electrode OCP entropic change [V.K-1]": params[
+                        "Positive electrode OCP entropic change [V.K-1]"
+                    ],
                     "Primary: Positive electrode active material volume fraction": 0.5,
                     "Secondary: Positive electrode active material volume fraction": 0.5,
                     "Primary: Maximum concentration in positive electrode [mol.m-3]": params[