Repository: pybamm-team/PyBaMM

Issue:
fix t_eval to data for drive cycles
when simulating drive cycles the solvers work if `t_eval` is the same as the times provided in the data. when solving we should check if the current function has "[current data]" and then force the solve times to be those in the data (with a warning message that the solver is fixing t_eval). see below example

```
#
# Simulate drive cycle loaded from csv file
#
import pybamm
import numpy as np
import pandas as pd
import os

# load model
pybamm.set_logging_level("INFO")
model = pybamm.lithium_ion.DFN()

# create geometry
geometry = model.default_geometry

# get drive cycle data
cycle_name = "US06"
drive_cycle = pd.read_csv(
    os.path.join(pybamm.root_dir(), "input", "drive_cycles", cycle_name+".csv"),
    comment="#",
    skip_blank_lines=True
)
time_dim = drive_cycle.values[:, 0]
current_dim = drive_cycle.values[:, 1]

# load parameter values and process model and geometry
param = model.default_parameter_values
param["Current function [A]"] = "[current data]" + cycle_name
param.process_model(model)
param.process_geometry(geometry)

# set mesh
var = pybamm.standard_spatial_vars
var_pts = {var.x_n: 11, var.x_s: 11, var.x_p: 11, var.r_n: 5, var.r_p: 5}
mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)

# discretise model
disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
disc.process_model(model)

# get dimensionless times for evaluation
tau = param.evaluate(pybamm.standard_parameters_lithium_ion.tau_discharge)
t_eval = time_dim / tau  # this works
# t_eval = np.linspace(0, 600 / tau, 600)  # this fails

# simulate US06 drive cycle
solver = pybamm.CasadiSolver(mode="fast")
solution = solver.solve(model, t_eval)

# plot
plot = pybamm.QuickPlot(solution)
plot.dynamic_plot()
```

Note that this works with Casadi in "fast" mode!

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3d68502a4..74a4bec60 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Features
 
 -   Added NCA parameter set ([#824](https://github.com/pybamm-team/PyBaMM/pull/824))
+-   Added functionality to `Solution` that automatically gets `t_eval` from the data when simulating drive cycles and performs checks to ensure the output has the required resolution to accurately capture the input current ([#819](https://github.com/pybamm-team/PyBaMM/pull/819))
 -   Added options to export a solution to matlab or csv ([#811](https://github.com/pybamm-team/PyBaMM/pull/811))
 -   Allow porosity to vary in space ([#809](https://github.com/pybamm-team/PyBaMM/pull/809))
 -   Added functionality to solve DAE models with non-smooth current inputs ([#808](https://github.com/pybamm-team/PyBaMM/pull/808))
@@ -52,6 +53,7 @@
 
 ## Bug fixes
 
+-   Fixed a bug where the first line of the data wasn't loaded when parameters are loaded from data ([#819](https://github.com/pybamm-team/PyBaMM/pull/819))
 -   Made `graphviz` an optional dependency ([#810](https://github.com/pybamm-team/PyBaMM/pull/810))
 -   Fixed examples to run with basic pip installation ([#800](https://github.com/pybamm-team/PyBaMM/pull/800))
 -   Added events for CasADi solver when stepping ([#800](https://github.com/pybamm-team/PyBaMM/pull/800))
diff --git a/examples/scripts/drive_cycle.py b/examples/scripts/drive_cycle.py
new file mode 100644
index 000000000..1e9c708da
--- /dev/null
+++ b/examples/scripts/drive_cycle.py
@@ -0,0 +1,17 @@
+#
+# Simulate drive cycle loaded from csv file
+#
+import pybamm
+
+# load model and update parameters so the input current is the US06 drive cycle
+model = pybamm.lithium_ion.DFN()
+param = model.default_parameter_values
+param["Current function [A]"] = "[current data]US06"
+
+# create and run simulation using the CasadiSolver in "fast" mode, remembering to
+# pass in the updated parameters
+sim = pybamm.Simulation(
+    model, parameter_values=param, solver=pybamm.CasadiSolver(mode="fast")
+)
+sim.solve()
+sim.plot()
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 5f48a2648..022ac6437 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -120,6 +120,7 @@ from .expression_tree.exceptions import (
     OptionError,
     ModelError,
     SolverError,
+    SolverWarning,
     ShapeError,
     ModelWarning,
     UndefinedOperationError,
diff --git a/pybamm/expression_tree/exceptions.py b/pybamm/expression_tree/exceptions.py
index 7097c475a..a71172cc4 100644
--- a/pybamm/expression_tree/exceptions.py
+++ b/pybamm/expression_tree/exceptions.py
@@ -37,6 +37,14 @@ class SolverError(Exception):
     pass
 
 
+class SolverWarning(UserWarning):
+    """
+    Solver warning: the chosen solver settings may not give the desired output
+    """
+
+    pass
+
+
 class ShapeError(Exception):
     """
     Shape error: cannot evaluate an object to find its shape
@@ -67,4 +75,3 @@ class InputError(Exception):
     """
 
     pass
-
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index a5a41afbf..aeb0fe7fa 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -227,7 +227,7 @@ class ParameterValues:
                         filename = os.path.join(path, value[6:] + ".csv")
                         function_name = value[6:]
                     data = pd.read_csv(
-                        filename, comment="#", skip_blank_lines=True
+                        filename, comment="#", skip_blank_lines=True, header=None
                     ).to_numpy()
                     # Save name and data
                     self._dict_items[name] = (function_name, data)
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index b16c83f46..9f3f02069 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -5,6 +5,8 @@ import pickle
 import pybamm
 import numpy as np
 import copy
+import warnings
+import sys
 
 
 def isnotebook():
@@ -317,10 +319,11 @@ class Simulation:
         Parameters
         ----------
         t_eval : numeric type, optional
-            The times at which to compute the solution. If None the model will
-            be solved for a full discharge (1 hour / C_rate) if the discharge
-            timescale is provided. Otherwise the model will be solved up to a
-            non-dimensional time of 1.
+            The times at which to compute the solution. If None and the parameter
+            "Current function [A]" is not read from data the model will
+            be solved for a full discharge (1 hour / C_rate). If None and the
+            parameter "Current function [A]" is read from data the model will be
+            solved at the times provided in the data.
         solver : :class:`pybamm.BaseSolver`
             The solver to use to solve the model.
         external_variables : dict
@@ -340,19 +343,68 @@ class Simulation:
             solver = self.solver
 
         if self.operating_mode == "without experiment":
-            # Solve the normal way, with a single solve
-            if t_eval is None:
-                try:
-                    # Try to compute discharge time
-                    tau = self._parameter_values.evaluate(self.model.param.timescale)
-                    C_rate = self._parameter_values["C-rate"]
-                    t_end = 3600 / tau / C_rate
-                    t_eval = np.linspace(0, t_end, 100)
-                except AttributeError:
-                    t_eval = np.linspace(0, 1, 100)
+            # For drive cycles (current provided as data) we perform additional tests
+            # on t_eval (if provided) to ensure the returned solution captures the
+            # input. If the current is provided as data then the "Current function [A]"
+            # is the tuple (filename, data).
+            if isinstance(self._parameter_values["Current function [A]"], tuple):
+                filename = self._parameter_values["Current function [A]"][0]
+                tau = self._parameter_values.evaluate(self.model.param.timescale)
+                time_data = (
+                    self._parameter_values["Current function [A]"][1][:, 0] / tau
+                )
+                # If no t_eval is provided, we use the times provided in the data.
+                if t_eval is None:
+                    pybamm.logger.info(
+                        "Setting t_eval as specified by the data '{}'".format(filename)
+                    )
+                    t_eval = time_data
+                # If t_eval is provided we first check if it contains all of the times
+                # in the data to within 10-12. If it doesn't, we then check
+                # that the largest gap in t_eval is smaller than the smallest gap in the
+                # time data (to ensure the resolution of t_eval is fine enough).
+                # We only raise a warning here as users may genuinely only want
+                # the solution returned at some specified points.
+                elif (
+                    set(np.round(time_data, 12)).issubset(set(np.round(t_eval, 12)))
+                ) is False:
+                    warnings.warn(
+                        """
+                        t_eval does not contain all of the time points in the data
+                        '{}'. Note: passing t_eval = None automatically sets t_eval
+                        to be the points in the data.
+                        """.format(
+                            filename
+                        ),
+                        pybamm.SolverWarning,
+                    )
+                    dt_data_min = np.min(np.diff(time_data))
+                    dt_eval_max = np.max(np.diff(t_eval))
+                    if dt_eval_max > dt_data_min + sys.float_info.epsilon:
+                        warnings.warn(
+                            """
+                            The largest timestep in t_eval ({}) is larger than
+                            the smallest timestep in the data ({}). The returned
+                            solution may not have the correct resolution to accurately
+                            capture the input. Try refining t_eval. Alternatively,
+                            passing t_eval = None automatically sets t_eval to be the
+                            points in the data.
+                            """.format(
+                                dt_eval_max, dt_data_min
+                            ),
+                            pybamm.SolverWarning,
+                        )
+            # If not using a drive cycle and t_eval is not provided, set t_eval
+            # to correspond to a single discharge
+            elif t_eval is None:
+                tau = self._parameter_values.evaluate(self.model.param.timescale)
+                C_rate = self._parameter_values["C-rate"]
+                t_end = 3600 / tau / C_rate
+                t_eval = np.linspace(0, t_end, 100)
 
             self.t_eval = t_eval
             self._solution = solver.solve(self.built_model, t_eval, inputs=inputs)
+
         elif self.operating_mode == "with experiment":
             if t_eval is not None:
                 pybamm.logger.warning(


Test patch:
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index 5e3a090a4..c2aa884c2 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -1,5 +1,7 @@
 import pybamm
 import numpy as np
+import pandas as pd
+import os
 import unittest
 
 
@@ -347,6 +349,40 @@ class TestSimulation(unittest.TestCase):
         sim.solve(t_eval=t_eval)
         sim.plot(testing=True)
 
+    def test_drive_cycle_data(self):
+        model = pybamm.lithium_ion.SPM()
+        param = model.default_parameter_values
+        param["Current function [A]"] = "[current data]US06"
+
+        drive_cycle = pd.read_csv(
+            os.path.join(pybamm.root_dir(), "input", "drive_cycles", "US06.csv"),
+            comment="#",
+            skip_blank_lines=True,
+            header=None,
+        )
+        time_data = drive_cycle.values[:, 0]
+        tau = param.evaluate(pybamm.standard_parameters_lithium_ion.tau_discharge)
+
+        sim = pybamm.Simulation(model, parameter_values=param)
+
+        # check solution is returned at the times in the data (only almost equal due
+        # to multiplication by tau)
+        sim.solve()
+        np.testing.assert_array_almost_equal(sim.solution.t * tau, time_data)
+
+        # check warning raised if the largest gap in t_eval is bigger than the
+        # smallest gap in the data
+        sim.reset()
+        with self.assertWarns(pybamm.SolverWarning):
+            sim.solve(t_eval=np.linspace(0, 1, 100))
+
+        # check warning raised if t_eval doesnt conatin time_data , but has a finer
+        # resolution (can still solve, but good for users to know they dont have
+        # the solution returned at the data points)
+        sim.reset()
+        with self.assertWarns(pybamm.SolverWarning):
+            sim.solve(t_eval=np.linspace(0, time_data[-1] / tau, 800))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")