Fix Q_heat always zero for isothermal cells; add unit tests (#9)

* Fix Q_heat always zero for isothermal cells; add unit tests

- Enable 'calculate heat source for isothermal models' option for
  CellElectrical and CellCoSimElectrical via _thermal_extra_options
  class variable, wired into both _CellBase and _CoSimCellBase
- Fix _CoSimCellBase._initial_outputs: bm.y0 is None in PyBaMM 26.x
  before the first step; return zero-filled placeholder with SOC set
- Pass save=False to pybamm.Simulation.step() to prevent unbounded
  memory growth on long co-sim runs
- Add tests: Q_heat nonzero during discharge, temperature input
  effectiveness — for all four cell classes

* Fixes

Author: DavidMStraub

src/pathsim_batt/cells/pybamm_cell.py

@@ -62,6 +62,7 @@ class _CellBase(DynamicalSystem):
     """
 
     _thermal_option: str = ""
+    _thermal_extra_options: dict[str, str] = {}
     _pybamm_output_vars: list[str] = []
 
     def __init__(
@@ -74,7 +75,9 @@ def __init__(
         self._initial_soc = float(initial_soc)
 
         if model is None:
-            model = pybamm.lithium_ion.SPMe(options={"thermal": self._thermal_option})
+            model = pybamm.lithium_ion.SPMe(
+                options={"thermal": self._thermal_option, **self._thermal_extra_options}
+            )
 
         self._parameter_values = _prepare_parameter_values(parameter_values)
 
@@ -180,6 +183,7 @@ class _CoSimCellBase(Wrapper):
     """
 
     _thermal_option: str = ""
+    _thermal_extra_options: dict[str, str] = {}
     _pybamm_output_vars: list[str] = []
 
     def __init__(
@@ -196,7 +200,9 @@ def __init__(
             raise ValueError("dt must be positive")
 
         if model is None:
-            model = pybamm.lithium_ion.SPMe(options={"thermal": self._thermal_option})
+            model = pybamm.lithium_ion.SPMe(
+                options={"thermal": self._thermal_option, **self._thermal_extra_options}
+            )
 
         self._model = model
         self._parameter_values = _prepare_parameter_values(parameter_values)
@@ -278,7 +284,8 @@ class CellElectrical(_CellBase):
     Parameters
     ----------
     model : pybamm.BaseBatteryModel or None
-        PyBaMM lithium-ion model.  Defaults to ``SPMe(thermal="isothermal")``.
+        PyBaMM lithium-ion model.  Defaults to isothermal SPMe with heat
+        source calculation enabled.
     parameter_values : pybamm.ParameterValues or None
         PyBaMM parameter set.  Defaults to ``Chen2020``.
     initial_soc : float
@@ -300,6 +307,7 @@ class CellElectrical(_CellBase):
     """
 
     _thermal_option = "isothermal"
+    _thermal_extra_options = {"calculate heat source for isothermal models": "true"}
     _pybamm_output_vars = [
         "Terminal voltage [V]",
         "X-averaged total heating [W.m-3]",
@@ -371,7 +379,8 @@ class CellCoSimElectrical(_CoSimCellBase):
     Parameters
     ----------
     model : pybamm.BaseBatteryModel or None
-        PyBaMM lithium-ion model. Defaults to ``SPMe(thermal="isothermal")``.
+        PyBaMM lithium-ion model.  Defaults to isothermal SPMe with heat
+        source calculation enabled.
     parameter_values : pybamm.ParameterValues or None
         PyBaMM parameter set. Defaults to ``Chen2020``.
     initial_soc : float
@@ -384,6 +393,7 @@ class CellCoSimElectrical(_CoSimCellBase):
     """
 
     _thermal_option = "isothermal"
+    _thermal_extra_options = {"calculate heat source for isothermal models": "true"}
     _pybamm_output_vars = [
         "Terminal voltage [V]",
         "X-averaged total heating [W.m-3]",

tests/cells/test_pybamm_cell.py

@@ -190,6 +190,64 @@ def test_pathsim_state_advances(self):
         self.sim.run(2)
         self.assertFalse(np.allclose(self.cell.engine.state, state_before))
 
+    def test_q_heat_nonzero_during_discharge(self):
+        """Q_heat must be strictly positive when a discharge current flows.
+
+        With thermal='isothermal' PyBaMM does not compute heat source terms,
+        so Q_heat would be identically zero — this test guards against that.
+        """
+        cell = CellElectrical(initial_soc=1.0)
+        I_src = Constant(5.0)  # 1C-ish discharge
+        T_src = Constant(298.15)
+        sim = Simulation(
+            blocks=[I_src, T_src, cell],
+            connections=[
+                Connection(I_src, cell["I"]),
+                Connection(T_src, cell["T_cell"]),
+            ],
+            dt=10.0,
+            Solver=ESDIRK43,
+        )
+        sim.run(60)
+        self.assertGreater(
+            cell.outputs[1],
+            0.0,
+            "Q_heat is zero — thermal model may not compute heat sources",
+        )
+
+    def test_temperature_input_affects_voltage(self):
+        """T_cell must actually influence the electrochemistry.
+
+        Butler-Volmer kinetics are temperature-dependent, so discharging at
+        a significantly higher temperature must yield a measurably different
+        terminal voltage after the same duration.
+        """
+
+        def _run_and_get_voltage(T_cell):
+            cell = CellElectrical(initial_soc=1.0)
+            I_src = Constant(5.0)
+            T_src = Constant(T_cell)
+            sim = Simulation(
+                blocks=[I_src, T_src, cell],
+                connections=[
+                    Connection(I_src, cell["I"]),
+                    Connection(T_src, cell["T_cell"]),
+                ],
+                dt=10.0,
+                Solver=ESDIRK43,
+            )
+            sim.run(300)
+            return cell.outputs[0]  # terminal voltage [V]
+
+        V_cold = _run_and_get_voltage(278.15)  # 5 °C
+        V_hot = _run_and_get_voltage(318.15)  # 45 °C
+        self.assertNotAlmostEqual(
+            V_cold,
+            V_hot,
+            places=3,
+            msg="T_cell input has no effect on terminal voltage",
+        )
+
 
 class TestElectrothermal(unittest.TestCase):
     """Integration tests for CellElectrothermal — PathSim integrates the PyBaMM ODE."""
@@ -240,6 +298,58 @@ def test_pathsim_state_advances(self):
         self.sim.run(2)
         self.assertFalse(np.allclose(self.cell.engine.state, state_before))
 
+    def test_q_heat_nonzero_during_discharge(self):
+        """Q_heat must be strictly positive when a discharge current flows."""
+        cell = CellElectrothermal(initial_soc=1.0)
+        I_src = Constant(5.0)
+        T_src = Constant(298.15)
+        sim = Simulation(
+            blocks=[I_src, T_src, cell],
+            connections=[
+                Connection(I_src, cell["I"]),
+                Connection(T_src, cell["T_amb"]),
+            ],
+            dt=10.0,
+            Solver=ESDIRK43,
+        )
+        sim.run(60)
+        self.assertGreater(
+            cell.outputs[2],
+            0.0,
+            "Q_heat is zero — thermal model may not compute heat sources",
+        )
+
+    def test_tamb_input_affects_cell_temperature(self):
+        """T_amb must influence the output cell temperature.
+
+        With a lower ambient temperature the cell should run cooler after
+        the same discharge duration.
+        """
+
+        def _run_and_get_T_cell(T_amb):
+            cell = CellElectrothermal(initial_soc=1.0)
+            I_src = Constant(5.0)
+            T_src = Constant(T_amb)
+            sim = Simulation(
+                blocks=[I_src, T_src, cell],
+                connections=[
+                    Connection(I_src, cell["I"]),
+                    Connection(T_src, cell["T_amb"]),
+                ],
+                dt=10.0,
+                Solver=ESDIRK43,
+            )
+            sim.run(300)
+            return cell.outputs[1]  # cell temperature [K]
+
+        T_cell_cold_amb = _run_and_get_T_cell(278.15)  # 5 °C ambient
+        T_cell_hot_amb = _run_and_get_T_cell(318.15)  # 45 °C ambient
+        self.assertLess(
+            T_cell_cold_amb,
+            T_cell_hot_amb,
+            msg="T_amb input has no effect on output cell temperature",
+        )
+
 
 class TestCoSimulationElectrical(unittest.TestCase):
     """Integration tests for CellCoSimElectrical — PyBaMM performs the stepping."""
@@ -293,6 +403,55 @@ def test_dfn_step_outputs_physical(self):
         self.assertGreater(cell.outputs[2], 0.0)  # SOC
         self.assertLessEqual(cell.outputs[2], 1.0)
 
+    def test_q_heat_nonzero_during_discharge(self):
+        """Q_heat must be strictly positive when a discharge current flows."""
+        cell = CellCoSimElectrical(initial_soc=1.0, dt=10.0)
+        I_src = Constant(5.0)
+        T_src = Constant(298.15)
+        sim = Simulation(
+            blocks=[I_src, T_src, cell],
+            connections=[
+                Connection(I_src, cell["I"]),
+                Connection(T_src, cell["T_cell"]),
+            ],
+            dt=5.0,
+            Solver=ESDIRK43,
+        )
+        sim.run(60)
+        self.assertGreater(
+            cell.outputs[1],
+            0.0,
+            "Q_heat is zero — thermal model may not compute heat sources",
+        )
+
+    def test_temperature_input_affects_voltage(self):
+        """T_cell must actually influence the electrochemistry."""
+
+        def _run_and_get_voltage(T_cell):
+            cell = CellCoSimElectrical(initial_soc=1.0, dt=10.0)
+            I_src = Constant(5.0)
+            T_src = Constant(T_cell)
+            sim = Simulation(
+                blocks=[I_src, T_src, cell],
+                connections=[
+                    Connection(I_src, cell["I"]),
+                    Connection(T_src, cell["T_cell"]),
+                ],
+                dt=5.0,
+                Solver=ESDIRK43,
+            )
+            sim.run(300)
+            return cell.outputs[0]  # terminal voltage [V]
+
+        V_cold = _run_and_get_voltage(278.15)  # 5 °C
+        V_hot = _run_and_get_voltage(318.15)  # 45 °C
+        self.assertNotAlmostEqual(
+            float(V_cold),
+            float(V_hot),
+            places=3,
+            msg="T_cell input has no effect on terminal voltage",
+        )
+
 
 class TestCoSimulationElectrothermal(unittest.TestCase):
     """Integration tests for CellCoSimElectrothermal — PyBaMM performs the stepping."""
@@ -349,6 +508,52 @@ def test_dfn_step_outputs_physical(self):
         self.assertGreater(cell.outputs[3], 0.0)  # SOC
         self.assertLessEqual(cell.outputs[3], 1.0)
 
+    def test_q_heat_nonzero_during_discharge(self):
+        """Q_heat must be strictly positive when a discharge current flows."""
+        cell = CellCoSimElectrothermal(initial_soc=1.0, dt=10.0)
+        I_src = Constant(5.0)
+        T_src = Constant(298.15)
+        sim = Simulation(
+            blocks=[I_src, T_src, cell],
+            connections=[
+                Connection(I_src, cell["I"]),
+                Connection(T_src, cell["T_amb"]),
+            ],
+            dt=5.0,
+            Solver=ESDIRK43,
+        )
+        sim.run(60)
+        self.assertGreater(
+            cell.outputs[2],
+            0.0,
+            "Q_heat is zero — thermal model may not compute heat sources",
+        )
+
+    def test_tamb_input_affects_cell_temperature(self):
+        """T_amb must influence the output cell temperature."""
+
+        def _run_and_get_T_cell(T_amb):
+            cell = CellCoSimElectrothermal(initial_soc=1.0, dt=10.0)
+            I_src = Constant(5.0)
+            T_src = Constant(T_amb)
+            sim = Simulation(
+                blocks=[I_src, T_src, cell],
+                connections=[
+                    Connection(I_src, cell["I"]),
+                    Connection(T_src, cell["T_amb"]),
+                ],
+                dt=5.0,
+                Solver=ESDIRK43,
+            )
+            sim.run(300)
+            return cell.outputs[1]  # cell temperature [K]
+
+        T_cold = _run_and_get_T_cell(278.15)
+        T_hot = _run_and_get_T_cell(318.15)
+        self.assertLess(
+            T_cold, T_hot, msg="T_amb input has no effect on output cell temperature"
+        )
+
 
 if __name__ == "__main__":
     unittest.main()