Repository: pybamm-team/PyBaMM

Issue:
add `DiscreteTimeSum` subclass of `pybamm.Symbol`
### Description

This would add a new unary operator in a pybamm expression tree that would represent a discrete sum over time. This would be a discrete version of the already existing `pybamm.ExplicitTimeIntegral` class.

### Motivation

This would give model developers and users the ability to add "sum of squares" type variables to a pybamm model that would calculate the difference between, for example, a solution variable and a user-provided dataset. This would be useful for implementing parameter inference using pybamm



### Possible Implementation

Similar to the `ExplicitTimeIntegral`, this would only be evaluated in the `Solution` class (see `_update_variable` function). I would propose that the sum is done over the list of time points in the solution (rather than have a separate list of points in the expression tree node), that way a user can specify the time-points via `t_interp`.

I'm a bit unsure how to allow users to provide their data to the expression. Perhaps they would wrap it in an Interpolant:

```python
data = pybamm.Interpolant(data_t, data_y, pybamm.t)
model.variable["data_comparison"] = pybamm.DiscreteTimeSum((model.variable["Voltage"] - data)**2)
```

### Additional context

see https://github.com/pybop-team/PyBOP/issues/513

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4501?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) Report
Attention: Patch coverage is `96.92308%` with `2 lines` in your changes missing coverage. Please review.
> Project coverage is 99.41%. Comparing base [(`9e62b66`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/9e62b66c53f423ceb5bb337d4fd893f71358c8d4?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) to head [(`4cb5360`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/4cb53603e4378c09c1bf48171dd937a01d2c1adc?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).

| [Files with missing lines](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4501?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) | Patch % | Lines |
|---|---|---|
| [src/pybamm/solvers/processed\_variable.py](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4501?src=pr&el=tree&filepath=src%2Fpybamm%2Fsolvers%2Fprocessed_variable.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team#diff-c3JjL3B5YmFtbS9zb2x2ZXJzL3Byb2Nlc3NlZF92YXJpYWJsZS5weQ==) | 92.30% | [1 Missing :warning: ](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4501?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) |
| [...pybamm/solvers/processed\_variable\_time\_integral.py](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4501?src=pr&el=tree&filepath=src%2Fpybamm%2Fsolvers%2Fprocessed_variable_time_integral.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team#diff-c3JjL3B5YmFtbS9zb2x2ZXJzL3Byb2Nlc3NlZF92YXJpYWJsZV90aW1lX2ludGVncmFsLnB5) | 93.75% | [1 Missing :warning: ](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4501?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) |

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##           develop    #4501      +/-   ##
===========================================
- Coverage    99.42%   99.41%   -0.01%     
===========================================
  Files          297      299       +2     
  Lines        22685    22739      +54     
===========================================
+ Hits         22554    22606      +52     
- Misses         131      133       +2     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4501?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8ffc210a7..fd69942fd 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@
 - Added phase-dependent particle options to LAM
   ([#4369](https://github.com/pybamm-team/PyBaMM/pull/4369))
 - Added a lithium ion equivalent circuit model with split open circuit voltages for each electrode (`SplitOCVR`). ([#4330](https://github.com/pybamm-team/PyBaMM/pull/4330))
+- Added the `pybamm.DiscreteTimeSum` expression node to sum an expression over a sequence of data times, and accompanying `pybamm.DiscreteTimeData` class to store the data times and values ([#4501](https://github.com/pybamm-team/PyBaMM/pull/4501))
 
 ## Optimizations
 
diff --git a/src/pybamm/__init__.py b/src/pybamm/__init__.py
index f51ba05d0..2deb30f30 100644
--- a/src/pybamm/__init__.py
+++ b/src/pybamm/__init__.py
@@ -36,6 +36,7 @@ from .expression_tree.averages import _BaseAverage
 from .expression_tree.broadcasts import *
 from .expression_tree.functions import *
 from .expression_tree.interpolant import Interpolant
+from .expression_tree.discrete_time_sum import *
 from .expression_tree.input_parameter import InputParameter
 from .expression_tree.parameter import Parameter, FunctionParameter
 from .expression_tree.scalar import Scalar
@@ -158,6 +159,7 @@ from .spatial_methods.scikit_finite_element import ScikitFiniteElement
 
 # Solver classes
 from .solvers.solution import Solution, EmptySolution, make_cycle_solution
+from .solvers.processed_variable_time_integral import ProcessedVariableTimeIntegral
 from .solvers.processed_variable import ProcessedVariable, process_variable
 from .solvers.processed_variable_computed import ProcessedVariableComputed
 from .solvers.base_solver import BaseSolver
diff --git a/src/pybamm/expression_tree/__init__.py b/src/pybamm/expression_tree/__init__.py
index 0b06746e6..7ac80e535 100644
--- a/src/pybamm/expression_tree/__init__.py
+++ b/src/pybamm/expression_tree/__init__.py
@@ -2,4 +2,4 @@ __all__ = ['array', 'averages', 'binary_operators', 'broadcasts',
            'concatenations', 'exceptions', 'functions', 'independent_variable',
            'input_parameter', 'interpolant', 'matrix', 'operations',
            'parameter', 'printing', 'scalar', 'state_vector', 'symbol',
-           'unary_operators', 'variable', 'vector']
+           'unary_operators', 'variable', 'vector', 'discrete_time_sum' ]
diff --git a/src/pybamm/expression_tree/discrete_time_sum.py b/src/pybamm/expression_tree/discrete_time_sum.py
new file mode 100644
index 000000000..41cd14960
--- /dev/null
+++ b/src/pybamm/expression_tree/discrete_time_sum.py
@@ -0,0 +1,88 @@
+import pybamm
+import numpy as np
+
+
+class DiscreteTimeData(pybamm.Interpolant):
+    """
+    A class for representing data that is only defined at discrete points in time.
+    This is implemented as a 1D interpolant with the time points as the nodes.
+
+    Parameters
+    ----------
+
+    time_points : :class:`numpy.ndarray`
+        The time points at which the data is defined
+    data : :class:`numpy.ndarray`
+        The data to be interpolated
+    name : str
+        The name of the data
+
+    """
+
+    def __init__(self, time_points: np.ndarray, data: np.ndarray, name: str):
+        super().__init__(time_points, data, pybamm.t, name)
+
+    def create_copy(self, new_children=None, perform_simplifications=True):
+        """See :meth:`pybamm.Symbol.new_copy()`."""
+        return pybamm.DiscreteTimeData(self.x[0], self.y, self.name)
+
+
+class DiscreteTimeSum(pybamm.UnaryOperator):
+    """
+    A node in the expression tree representing a discrete time sum operator.
+
+    .. math::
+        \\sum_{i=0}^{N} f(y(t_i), t_i)
+
+    where f is the expression given by the child, and the sum is over the discrete
+    time points t_i. The set of time points is given by the :class:`pybamm.DiscreteTimeData` node,
+    which must be somewhere in the expression tree given by the child. If the child
+    does not contain a :class:`pybamm.DiscreteTimeData` node, then an error will be raised when
+    the node is created. If the child contains multiple :class:`pybamm.DiscreteTimeData` nodes,
+    an error will be raised when the node is created.
+
+
+    Parameters
+    ----------
+    child: :class:`pybamm.Symbol`
+        The symbol to be summed
+
+    Attributes
+    ----------
+    data: :class:`pybamm.DiscreteTimeData`
+        The discrete time data node in the child
+
+    Raises
+    ------
+    :class:`pybamm.ModelError`
+        If the child does not contain a :class:`pybamm.DiscreteTimeData` node, or if the child
+        contains multiple :class:`pybamm.DiscreteTimeData` nodes.
+    """
+
+    def __init__(self, child: pybamm.Symbol):
+        self.data = None
+        for node in child.pre_order():
+            if isinstance(node, DiscreteTimeData):
+                # Check that there is exactly one DiscreteTimeData node in the child
+                if self.data is not None:
+                    raise pybamm.ModelError(
+                        "DiscreteTimeSum can only have one DiscreteTimeData node in the child"
+                    )
+                self.data = node
+        if self.data is None:
+            raise pybamm.ModelError(
+                "DiscreteTimeSum must contain a DiscreteTimeData node"
+            )
+        super().__init__("discrete time sum", child)
+
+    @property
+    def sum_values(self):
+        return self.data.y
+
+    @property
+    def sum_times(self):
+        return self.data.x[0]
+
+    def _unary_evaluate(self, child):
+        # return result of evaluating the child, we'll only implement the sum once the model is solved (in pybamm.ProcessedVariable)
+        return child
diff --git a/src/pybamm/solvers/__init__.py b/src/pybamm/solvers/__init__.py
index fc8be7e2f..e9d9d306f 100644
--- a/src/pybamm/solvers/__init__.py
+++ b/src/pybamm/solvers/__init__.py
@@ -2,4 +2,4 @@ __all__ = ['algebraic_solver', 'base_solver', 'c_solvers',
            'casadi_algebraic_solver', 'casadi_solver', 'dummy_solver',
            'idaklu_jax', 'idaklu_solver', 'jax_bdf_solver', 'jax_solver',
            'lrudict', 'processed_variable', 'processed_variable_computed',
-           'scipy_solver', 'solution']
+           'scipy_solver', 'solution', 'processed_variable_time_integral']
diff --git a/src/pybamm/solvers/processed_variable.py b/src/pybamm/solvers/processed_variable.py
index 5cf928ca7..12cf140b3 100644
--- a/src/pybamm/solvers/processed_variable.py
+++ b/src/pybamm/solvers/processed_variable.py
@@ -1,6 +1,7 @@
 #
 # Processed Variable class
 #
+from typing import Optional
 import casadi
 import numpy as np
 import pybamm
@@ -29,6 +30,8 @@ class ProcessedVariable:
         `base_Variable.evaluate` (but more efficiently).
     solution : :class:`pybamm.Solution`
         The solution object to be used to create the processed variables
+    time_integral : :class:`pybamm.ProcessedVariableTimeIntegral`, optional
+        Not none if the variable is to be time-integrated (default is None)
     """
 
     def __init__(
@@ -36,7 +39,7 @@ class ProcessedVariable:
         base_variables,
         base_variables_casadi,
         solution,
-        cumtrapz_ic=None,
+        time_integral: Optional[pybamm.ProcessedVariableTimeIntegral] = None,
     ):
         self.base_variables = base_variables
         self.base_variables_casadi = base_variables_casadi
@@ -50,7 +53,7 @@ class ProcessedVariable:
         self.mesh = base_variables[0].mesh
         self.domain = base_variables[0].domain
         self.domains = base_variables[0].domains
-        self.cumtrapz_ic = cumtrapz_ic
+        self.time_integral = time_integral
 
         # Process spatial variables
         geometry = solution.all_models[0].geometry
@@ -271,18 +274,21 @@ class ProcessedVariable:
                     self._coords_raw,
                 )
 
-            processed_entries = self._xr_interpolate(
-                entries_for_interp,
-                coords,
-                observe_raw,
-                t,
-                x,
-                r,
-                y,
-                z,
-                R,
-                fill_value,
-            )
+            if self.time_integral is None:
+                processed_entries = self._xr_interpolate(
+                    entries_for_interp,
+                    coords,
+                    observe_raw,
+                    t,
+                    x,
+                    r,
+                    y,
+                    z,
+                    R,
+                    fill_value,
+                )
+            else:
+                processed_entries = entries_for_interp
         else:
             processed_entries = entries
 
@@ -343,6 +349,16 @@ class ProcessedVariable:
             t_observe (np.ndarray): time points to observe
             observe_raw (bool): True if observing the raw data
         """
+        # if this is a time integral variable, t must be None and we observe either the
+        # data times (for a discrete sum) or the solution times (for a continuous sum)
+        if self.time_integral is not None:
+            if self.time_integral.method == "discrete":
+                # discrete sum should be observed at the discrete times
+                t = self.time_integral.discrete_times
+            else:
+                # assume we can do a sufficiently accurate trapezoidal integration at t_pts
+                t = self.t_pts
+
         observe_raw = (t is None) or (
             np.asarray(t).size == len(self.t_pts) and np.all(t == self.t_pts)
         )
@@ -483,14 +499,14 @@ class ProcessedVariable0D(ProcessedVariable):
         base_variables,
         base_variables_casadi,
         solution,
-        cumtrapz_ic=None,
+        time_integral: Optional[pybamm.ProcessedVariableTimeIntegral] = None,
     ):
         self.dimensions = 0
         super().__init__(
             base_variables,
             base_variables_casadi,
             solution,
-            cumtrapz_ic=cumtrapz_ic,
+            time_integral=time_integral,
         )
 
     def _observe_raw_python(self):
@@ -510,13 +526,19 @@ class ProcessedVariable0D(ProcessedVariable):
                 idx += 1
         return entries
 
-    def _observe_postfix(self, entries, _):
-        if self.cumtrapz_ic is None:
+    def _observe_postfix(self, entries, t):
+        if self.time_integral is None:
             return entries
-
-        return cumulative_trapezoid(
-            entries, self.t_pts, initial=float(self.cumtrapz_ic)
-        )
+        if self.time_integral.method == "discrete":
+            return np.sum(entries, axis=0, initial=self.time_integral.initial_condition)
+        elif self.time_integral.method == "continuous":
+            return cumulative_trapezoid(
+                entries, self.t_pts, initial=float(self.time_integral.initial_condition)
+            )
+        else:
+            raise ValueError(
+                "time_integral method must be 'discrete' or 'continuous'"
+            )  # pragma: no cover
 
     def _interp_setup(self, entries, t):
         # save attributes for interpolation
@@ -556,14 +578,14 @@ class ProcessedVariable1D(ProcessedVariable):
         base_variables,
         base_variables_casadi,
         solution,
-        cumtrapz_ic=None,
+        time_integral: Optional[pybamm.ProcessedVariableTimeIntegral] = None,
     ):
         self.dimensions = 1
         super().__init__(
             base_variables,
             base_variables_casadi,
             solution,
-            cumtrapz_ic=cumtrapz_ic,
+            time_integral=time_integral,
         )
 
     def _observe_raw_python(self):
@@ -653,14 +675,14 @@ class ProcessedVariable2D(ProcessedVariable):
         base_variables,
         base_variables_casadi,
         solution,
-        cumtrapz_ic=None,
+        time_integral: Optional[pybamm.ProcessedVariableTimeIntegral] = None,
     ):
         self.dimensions = 2
         super().__init__(
             base_variables,
             base_variables_casadi,
             solution,
-            cumtrapz_ic=cumtrapz_ic,
+            time_integral=time_integral,
         )
         first_dim_nodes = self.mesh.nodes
         first_dim_edges = self.mesh.edges
@@ -819,14 +841,14 @@ class ProcessedVariable2DSciKitFEM(ProcessedVariable2D):
         base_variables,
         base_variables_casadi,
         solution,
-        cumtrapz_ic=None,
+        time_integral: Optional[pybamm.ProcessedVariableTimeIntegral] = None,
     ):
         self.dimensions = 2
         super(ProcessedVariable2D, self).__init__(
             base_variables,
             base_variables_casadi,
             solution,
-            cumtrapz_ic=cumtrapz_ic,
+            time_integral=time_integral,
         )
         y_sol = self.mesh.edges["y"]
         z_sol = self.mesh.edges["z"]
diff --git a/src/pybamm/solvers/processed_variable_time_integral.py b/src/pybamm/solvers/processed_variable_time_integral.py
new file mode 100644
index 000000000..4fcdfb56b
--- /dev/null
+++ b/src/pybamm/solvers/processed_variable_time_integral.py
@@ -0,0 +1,28 @@
+from dataclasses import dataclass
+from typing import Literal, Optional, Union
+import numpy as np
+import pybamm
+
+
+@dataclass
+class ProcessedVariableTimeIntegral:
+    method: Literal["discrete", "continuous"]
+    initial_condition: np.ndarray
+    discrete_times: Optional[np.ndarray]
+
+    @staticmethod
+    def from_pybamm_var(
+        var: Union[pybamm.DiscreteTimeSum, pybamm.ExplicitTimeIntegral],
+    ) -> "ProcessedVariableTimeIntegral":
+        if isinstance(var, pybamm.DiscreteTimeSum):
+            return ProcessedVariableTimeIntegral(
+                method="discrete", initial_condition=0.0, discrete_times=var.sum_times
+            )
+        elif isinstance(var, pybamm.ExplicitTimeIntegral):
+            return ProcessedVariableTimeIntegral(
+                method="continuous",
+                initial_condition=var.initial_condition.evaluate(),
+                discrete_times=None,
+            )
+        else:
+            raise ValueError("Unsupported variable type")  # pragma: no cover
diff --git a/src/pybamm/solvers/solution.py b/src/pybamm/solvers/solution.py
index c884e79e3..1aa540ab3 100644
--- a/src/pybamm/solvers/solution.py
+++ b/src/pybamm/solvers/solution.py
@@ -571,7 +571,7 @@ class Solution:
             self._update_variable(variable)
 
     def _update_variable(self, variable):
-        cumtrapz_ic = None
+        time_integral = None
         pybamm.logger.debug(f"Post-processing {variable}")
         vars_pybamm = [
             model.variables_and_events[variable] for model in self.all_models
@@ -591,16 +591,22 @@ class Solution:
                     "solve. Please re-run the solve with `output_variables` set to "
                     "include this variable."
                 )
-            elif isinstance(var_pybamm, pybamm.ExplicitTimeIntegral):
-                cumtrapz_ic = var_pybamm.initial_condition
-                cumtrapz_ic = cumtrapz_ic.evaluate()
-                var_pybamm = var_pybamm.child
-                var_casadi = self.process_casadi_var(
-                    var_pybamm,
-                    inputs,
-                    ys.shape,
+            elif isinstance(
+                var_pybamm, (pybamm.ExplicitTimeIntegral, pybamm.DiscreteTimeSum)
+            ):
+                time_integral = pybamm.ProcessedVariableTimeIntegral.from_pybamm_var(
+                    var_pybamm
                 )
-                model._variables_casadi[variable] = var_casadi
+                var_pybamm = var_pybamm.child
+                if variable in model._variables_casadi:
+                    var_casadi = model._variables_casadi[variable]
+                else:
+                    var_casadi = self.process_casadi_var(
+                        var_pybamm,
+                        inputs,
+                        ys.shape,
+                    )
+                    model._variables_casadi[variable] = var_casadi
                 vars_pybamm[i] = var_pybamm
             elif variable in model._variables_casadi:
                 var_casadi = model._variables_casadi[variable]
@@ -613,7 +619,7 @@ class Solution:
                 model._variables_casadi[variable] = var_casadi
             vars_casadi.append(var_casadi)
         var = pybamm.process_variable(
-            vars_pybamm, vars_casadi, self, cumtrapz_ic=cumtrapz_ic
+            vars_pybamm, vars_casadi, self, time_integral=time_integral
         )
 
         self._variables[variable] = var


Test patch:
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 812cbd8f6..0dbafa38c 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -724,6 +724,36 @@ class TestUnaryOperators:
         assert expr.create_copy() == expr
         assert not expr.is_constant()
 
+    def test_discrete_time_sum(self):
+        times = np.array([1, 2, 3, 4, 5])
+        values = np.array([2, 2, 3, 3, 1])
+        data = pybamm.DiscreteTimeData(times, values, "test")
+        assert data.name == "test"
+        np.testing.assert_array_equal(data.x[0], times)
+        np.testing.assert_array_equal(data.y, values)
+
+        y = pybamm.StateVector(slice(0, 1))
+
+        # check that raises error if data is not present
+        with pytest.raises(pybamm.ModelError, match="must contain a DiscreteTimeData"):
+            pybamm.DiscreteTimeSum(2 * y)
+
+        # check that raises error if two data are present
+        data2 = pybamm.DiscreteTimeData(values, times, "test2")
+        with pytest.raises(pybamm.ModelError, match="only have one DiscreteTimeData"):
+            pybamm.DiscreteTimeSum(data + data2)
+
+        sum = pybamm.DiscreteTimeSum(2 * data - y)
+        np.testing.assert_array_equal(sum.sum_times, times)
+        np.testing.assert_array_equal(sum.sum_values, values)
+        y = np.array([1])
+
+        # evaluate should return the values to sum up
+        for i in range(len(times)):
+            eval = sum.evaluate(y=y, t=times[i])
+            expect = 2 * values[i] - y
+            np.testing.assert_array_equal(eval[0], expect)
+
     def test_to_from_json(self, mocker):
         # UnaryOperator
         a = pybamm.Symbol("a", domain=["test"])
diff --git a/tests/unit/test_solvers/test_processed_variable.py b/tests/unit/test_solvers/test_processed_variable.py
index 04de88963..b3ab85060 100644
--- a/tests/unit/test_solvers/test_processed_variable.py
+++ b/tests/unit/test_solvers/test_processed_variable.py
@@ -160,6 +160,67 @@ class TestProcessedVariable:
                 processed_var._observe_raw_cpp(), processed_var._observe_raw_python()
             )
 
+    @pytest.mark.parametrize("hermite_interp", _hermite_args)
+    def test_processed_variable_0D_discrete_data(self, hermite_interp):
+        y = pybamm.StateVector(slice(0, 1))
+        t_sol = np.linspace(0, 1)
+        y_sol = np.array([np.linspace(0, 5)])
+        data_const = 3.6
+        if hermite_interp:
+            yp_sol = 5 * np.ones_like(y_sol)
+        else:
+            yp_sol = None
+
+        # hermite interpolation can do order 2 interpolation, otherwise make sure result is linear
+        order = 2 if hermite_interp else 1
+
+        # data is same timepoints as solution
+        data_t = t_sol
+        data_v = -data_const * data_t
+        data = pybamm.DiscreteTimeData(data_t, data_v, "test_data")
+        var = (y - data) ** order
+        expected_entries = (y_sol - data_v) ** order
+        var.mesh = None
+        model = pybamm.BaseModel()
+        var_casadi = to_casadi(var, y_sol)
+        processed_var = pybamm.process_variable(
+            [var],
+            [var_casadi],
+            self._sol_default(t_sol, y_sol, yp_sol, model),
+        )
+        np.testing.assert_array_equal(processed_var.entries, expected_entries.flatten())
+        np.testing.assert_array_equal(processed_var(data_t), expected_entries.flatten())
+
+        # data is different timepoints as solution
+        data_t = np.linspace(0, 1, 7)
+        data_v = -data_const * data_t
+        y_sol_interp = (np.interp(data_t, t_sol, y_sol[0]),)
+        data_v_interp = np.interp(t_sol, data_t, data_v)
+        data = pybamm.DiscreteTimeData(data_t, data_v, "test_data")
+
+        # check data interp
+        np.testing.assert_array_almost_equal(
+            data.evaluate(t=t_sol).flatten(), data_v_interp
+        )
+
+        var = (y - data) ** order
+        expected = (y_sol_interp - data_v) ** order
+        expected_entries = (y_sol - data_v_interp) ** order
+        var.mesh = None
+        model = pybamm.BaseModel()
+        var_casadi = to_casadi(var, y_sol)
+        processed_var = pybamm.process_variable(
+            [var],
+            [var_casadi],
+            self._sol_default(t_sol, y_sol, yp_sol, model),
+        )
+        np.testing.assert_array_almost_equal(
+            processed_var.entries, expected_entries.flatten(), decimal=10
+        )
+        np.testing.assert_array_almost_equal(
+            processed_var(t=data_t), expected.flatten(), decimal=10
+        )
+
     @pytest.mark.parametrize("hermite_interp", _hermite_args)
     def test_processed_variable_0D_no_sensitivity(self, hermite_interp):
         # without space
diff --git a/tests/unit/test_solvers/test_solution.py b/tests/unit/test_solvers/test_solution.py
index 4ac931253..2fb25f79a 100644
--- a/tests/unit/test_solvers/test_solution.py
+++ b/tests/unit/test_solvers/test_solution.py
@@ -420,3 +420,43 @@ class TestSolution:
         sim.solve(t_eval=np.linspace(0, 10, 10), inputs=inputs)
         time = sim.solution["Time [h]"](sim.solution.t)
         assert len(time) == 10
+
+    _solver_classes = [pybamm.CasadiSolver]
+    if pybamm.has_idaklu():
+        _solver_classes.append(pybamm.IDAKLUSolver)
+
+    @pytest.mark.parametrize("solver_class", _solver_classes)
+    def test_discrete_data_sum(self, solver_class):
+        model = pybamm.BaseModel(name="test_model")
+        c = pybamm.Variable("c")
+        model.rhs = {c: -c}
+        model.initial_conditions = {c: 1}
+        model.variables["c"] = c
+
+        data_times = np.linspace(0, 1, 10)
+        if solver_class == pybamm.IDAKLUSolver:
+            t_eval = [data_times[0], data_times[-1]]
+            t_interp = data_times
+        else:
+            t_eval = data_times
+            t_interp = None
+        solver = solver_class()
+        data_values = solver.solve(model, t_eval=t_eval, t_interp=t_interp)["c"].entries
+
+        data = pybamm.DiscreteTimeData(data_times, data_values, "test_data")
+        data_comparison = pybamm.DiscreteTimeSum((c - data) ** 2)
+
+        model = pybamm.BaseModel(name="test_model2")
+        a = pybamm.InputParameter("a")
+        model.rhs = {c: -a * c}
+        model.initial_conditions = {c: 1}
+        model.variables["data_comparison"] = data_comparison
+
+        solver = solver_class()
+        for a in [0.5, 1.0, 2.0]:
+            sol = solver.solve(model, t_eval=t_eval, inputs={"a": a})
+            y_sol = np.exp(-a * data_times)
+            expected = np.sum((y_sol - data_values) ** 2)
+            np.testing.assert_array_almost_equal(
+                sol["data_comparison"](), expected, decimal=2
+            )