Repository: pybamm-team/PyBaMM
Pull request #5092 (2025-07-01 19:23:32)

bug: explicit-time-integral returns integral
# Description

an addendum to #5063, which fixes explicit time integral so that the return value is the integral rather than the cumulative integral

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CHANGELOG.md) to document the change (include PR #)

# Important checks:

Please confirm the following before marking the PR as ready for review:
- No style issues: `nox -s pre-commit`
- All tests pass: `nox -s tests`
- The documentation builds: `nox -s doctests`
- Code is commented for hard-to-understand areas
- Tests added that prove fix is effective or that feature works



Patch:
diff --git a/src/pybamm/solvers/processed_variable_time_integral.py b/src/pybamm/solvers/processed_variable_time_integral.py
index 4148ba453..22a7aba3a 100644
--- a/src/pybamm/solvers/processed_variable_time_integral.py
+++ b/src/pybamm/solvers/processed_variable_time_integral.py
@@ -6,7 +6,7 @@ from typing import Literal
 import casadi
 import numpy as np
 import numpy.typing as npt
-from scipy.integrate import cumulative_trapezoid
+from scipy.integrate import trapezoid
 
 import pybamm
 
@@ -26,8 +26,8 @@ class ProcessedVariableTimeIntegral:
                 entries, axis=0, initial=self.initial_condition, keepdims=True
             )
         else:
-            return cumulative_trapezoid(
-                entries, t_pts, initial=float(self.initial_condition), axis=0
+            return np.array(
+                [trapezoid(entries, t_pts, axis=0) + float(self.initial_condition)]
             )
 
     def postfix(self, entries, t_pts, inputs) -> np.ndarray:
@@ -150,11 +150,15 @@ class ProcessedVariableTimeIntegral:
                 discrete_times=sum_node.sum_times,
             )
         elif isinstance(sum_node, pybamm.ExplicitTimeIntegral):
+            if isinstance(sum_node.initial_condition, pybamm.Symbol):
+                initial_condition = sum_node.initial_condition.evaluate()
+            else:
+                initial_condition = sum_node.initial_condition
             return ProcessedVariableTimeIntegral(
                 method="continuous",
                 post_sum_node=post_sum_node,
                 sum_node=sum_node,
-                initial_condition=sum_node.initial_condition.evaluate(),
+                initial_condition=initial_condition,
                 discrete_times=None,
             )
         else:


Test patch:
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index 6f6646a62..275c912b5 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -3,7 +3,7 @@ import os
 import numpy as np
 import pandas as pd
 import pytest
-from scipy.integrate import cumulative_trapezoid
+from scipy.integrate import trapezoid
 
 import pybamm
 from tests import no_internet_connection
@@ -94,9 +94,7 @@ class TestSimulation:
         t = sol["Time [s]"].data
         I = sol["Current [A]"].data
         q = sol["Discharge capacity [A.h]"].data
-        np.testing.assert_allclose(
-            q, cumulative_trapezoid(I, t, initial=0) / 3600, rtol=1e-7, atol=1e-6
-        )
+        np.testing.assert_allclose(q, trapezoid(I, t) / 3600, rtol=1e-7, atol=1e-6)
 
     def test_solve_non_battery_model(self):
         model = pybamm.BaseModel()
diff --git a/tests/unit/test_solvers/test_solution.py b/tests/unit/test_solvers/test_solution.py
index 348eab5b4..2057e5092 100644
--- a/tests/unit/test_solvers/test_solution.py
+++ b/tests/unit/test_solvers/test_solution.py
@@ -8,6 +8,7 @@ import logging
 import numpy as np
 import pandas as pd
 import pytest
+import scipy
 from scipy.io import loadmat
 
 import pybamm
@@ -653,3 +654,87 @@ class TestSolution:
                         inputs={"a": a, "b": b},
                         calculate_sensitivities=True,
                     )["data_comparison"].sensitivities["a"]
+
+    @pytest.mark.parametrize(
+        "solver_class,use_post_sum,use_output_var", _solver_classes
+    )
+    def test_explicit_time_integral(self, solver_class, use_post_sum, use_output_var):
+        times = np.linspace(0, 1, 10)
+        c = pybamm.Variable("c")
+        if solver_class == pybamm.IDAKLUSolver:
+            t_eval = [times[0], times[-1]]
+            t_interp = times
+        else:
+            t_eval = times
+            t_interp = None
+
+        if use_post_sum:
+            integral = pybamm.ExplicitTimeIntegral(c, 0) ** 2
+        else:
+            integral = pybamm.ExplicitTimeIntegral(c, 0)
+
+        model = pybamm.BaseModel(name="test_model")
+        a = pybamm.InputParameter("a")
+        b = pybamm.InputParameter("b")
+        c2 = pybamm.Variable("c2")
+        model.rhs = {c: b * -a * c, c2: -2 * c2}
+        model.initial_conditions = {c: 1, c2: 1}
+        model.variables["integral"] = integral
+        model.variables["c"] = c
+
+        if use_output_var:
+            output_variables = ["integral", "c"]
+            solver = solver_class(output_variables=output_variables)
+        else:
+            solver = solver_class()
+        range = [0.5, 1.0, 2.0]
+        range2 = np.ones(3)
+        for a, b in zip(range, range2, strict=False):
+            sol = solver.solve(
+                model, t_eval=t_eval, t_interp=t_interp, inputs={"a": a, "b": b}
+            )
+            y_sol = np.exp(b * -a * times)
+            expected = -(1.0 / b / a) * (
+                np.exp(b * -a * times[-1]) - np.exp(b * -a * times[0])
+            )
+            if use_post_sum:
+                expected = expected**2
+            np.testing.assert_allclose(
+                sol["integral"](), expected, rtol=1e-3, atol=1e-2
+            )
+            assert isinstance(sol["integral"].data, np.ndarray)
+
+            # sensitivity calculation only supported for IDAKLUSolver
+            if solver_class == pybamm.IDAKLUSolver:
+                sol = solver.solve(
+                    model,
+                    t_eval=t_eval,
+                    t_interp=t_interp,
+                    inputs={"a": a, "b": b},
+                    calculate_sensitivities=True,
+                )
+                y_sol = np.exp(b * -a * times)
+                dy_sol_da = -b * times * y_sol
+                expected_sens = scipy.integrate.trapezoid(dy_sol_da, times)
+                if use_post_sum:
+                    expected_sens = 2 * expected * expected_sens
+
+                np.testing.assert_allclose(
+                    sol["c"].data,
+                    y_sol,
+                    rtol=1e-3,
+                    atol=1e-2,
+                )
+                np.testing.assert_allclose(
+                    sol["c"].sensitivities["a"].flatten(),
+                    dy_sol_da,
+                    rtol=1e-3,
+                    atol=1e-2,
+                )
+                np.testing.assert_allclose(
+                    sol["integral"].sensitivities["a"],
+                    expected_sens,
+                    rtol=1e-3,
+                    atol=1e-2,
+                )
+                assert isinstance(sol["integral"].sensitivities["a"], np.ndarray)