From eb0685e52250d7836366455eb3484fd2664774bd Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Wed, 6 Nov 2019 11:06:56 +0000
Subject: [PATCH 01/31] #692 added an get_external_variables method to
 submodels

---
 .../full_battery_models/base_battery_model.py | 14 +++++++
 pybamm/models/submodels/base_submodel.py      | 15 +++++++
 .../submodels/thermal/external/__init__.py    |  0
 .../models/submodels/thermal/external/full.py | 32 +++++++++++++++
 .../submodels/thermal/external/lumped.py      | 39 +++++++++++++++++++
 5 files changed, 100 insertions(+)
 create mode 100644 pybamm/models/submodels/thermal/external/__init__.py
 create mode 100644 pybamm/models/submodels/thermal/external/full.py
 create mode 100644 pybamm/models/submodels/thermal/external/lumped.py

diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index a0767a4573..6bd4a3b6d7 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -416,6 +416,19 @@ def build_model(self):
             )
             self.variables.update(submodel.get_fundamental_variables())
 
+        # Get any external variables
+        self.external_variables = []
+        for submodel_name, submodel in self.submodels.items():
+            pybamm.logger.debug(
+                "Getting external variables for {} submodel ({})".format(
+                    submodel_name, self.name
+                )
+            )
+            variables, external_variables = submodel.get_external_variables()
+
+            self.variables.update(variables)
+            self.external_variables += external_variables
+
         # Get coupled variables
         # Note: pybamm will try to get the coupled variables for the submodels in the
         # order they are set by the user. If this fails for a particular submodel,
@@ -454,6 +467,7 @@ def build_model(self):
             pybamm.logger.debug(
                 "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
             )
+
             submodel.set_rhs(self.variables)
             pybamm.logger.debug(
                 "Setting algebraic for {} submodel ({})".format(
diff --git a/pybamm/models/submodels/base_submodel.py b/pybamm/models/submodels/base_submodel.py
index 26f41a2af9..25a1aaf7c1 100644
--- a/pybamm/models/submodels/base_submodel.py
+++ b/pybamm/models/submodels/base_submodel.py
@@ -101,6 +101,21 @@ def get_fundamental_variables(self):
         """
         return {}
 
+    def set_external_variables(self):
+        """
+        A public method that creates and returns the variables in a submodel which are
+        suppled external to the model.
+
+        Returns
+        -------
+        dict :
+            The variables created by the submodel which are independent of variables in
+            other submodels.
+        list :
+            A list of the external variables in the model.
+        """
+        return {}, []
+
     def get_coupled_variables(self, variables):
         """
         A public method that creates and returns the variables in a submodel which
diff --git a/pybamm/models/submodels/thermal/external/__init__.py b/pybamm/models/submodels/thermal/external/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/pybamm/models/submodels/thermal/external/full.py b/pybamm/models/submodels/thermal/external/full.py
new file mode 100644
index 0000000000..4c131ebbe3
--- /dev/null
+++ b/pybamm/models/submodels/thermal/external/full.py
@@ -0,0 +1,32 @@
+#
+# A class for full external thermal models
+#
+import pybamm
+
+from ..base_thermal import BaseThermal
+
+
+class Full(BaseThermal):
+    """Class to link full external thermal submodels.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+
+
+    **Extends:** :class:`pybamm.thermal.BaseModel`
+    """
+
+    def __init__(self, param):
+        super().__init__(param)
+
+    def get_external_variables(self):
+        T = pybamm.standard_variables.T
+        T_cn = pybamm.BoundaryValue(T, "left")
+        T_cp = pybamm.BoundaryValue(T, "right")
+
+        variables = self._get_standard_fundamental_variables(T, T_cn, T_cp)
+        external_variables = [T]
+
+        return variables, external_variables
diff --git a/pybamm/models/submodels/thermal/external/lumped.py b/pybamm/models/submodels/thermal/external/lumped.py
new file mode 100644
index 0000000000..7da09c5c2b
--- /dev/null
+++ b/pybamm/models/submodels/thermal/external/lumped.py
@@ -0,0 +1,39 @@
+#
+# A class for external lumped thermal models
+#
+import pybamm
+
+from ..base_thermal import BaseThermal
+
+
+class Full(BaseThermal):
+    """Class to link external lumped thermal submodels.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+
+
+    **Extends:** :class:`pybamm.thermal.BaseModel`
+    """
+
+    def __init__(self, param):
+        super().__init__(param)
+
+    def get_external_variables(self):
+        T_x_av = pybamm.standard_variables.T_av
+        T_n = pybamm.PrimaryBroadcast(T_x_av, "negative electrode")
+        T_s = pybamm.PrimaryBroadcast(T_x_av, "separator")
+        T_p = pybamm.PrimaryBroadcast(T_x_av, "positive electrode")
+        T = pybamm.Concatenation(T_n, T_s, T_p)
+
+        T_cn = T_x_av
+        T_cp = T_x_av
+
+        variables = self._get_standard_fundamental_variables(T, T_cn, T_cp)
+        variables = self._get_standard_fundamental_variables(T, T_cn, T_cp)
+
+        external_variables = [T_x_av]
+
+        return variables, external_variables

From faf0a59fb1baa22dfdb45903f49b3396a0670fab Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Fri, 8 Nov 2019 17:37:20 +0000
Subject: [PATCH 02/31] #692 changed to pass in external option

---
 pybamm/discretisations/discretisation.py      | 14 ++-
 pybamm/models/base_model.py                   |  7 +-
 .../full_battery_models/base_battery_model.py | 98 ++++++++++++-------
 pybamm/models/submodels/base_submodel.py      | 38 +++++--
 .../submodels/thermal/external/__init__.py    |  0
 .../models/submodels/thermal/external/full.py | 32 ------
 .../submodels/thermal/external/lumped.py      | 39 --------
 pybamm/simulation.py                          |  7 +-
 .../test_use_external_submodel.py             | 25 +++++
 .../standard_submodel_unit_tests.py           |  6 ++
 10 files changed, 143 insertions(+), 123 deletions(-)
 delete mode 100644 pybamm/models/submodels/thermal/external/__init__.py
 delete mode 100644 pybamm/models/submodels/thermal/external/full.py
 delete mode 100644 pybamm/models/submodels/thermal/external/lumped.py
 create mode 100644 tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 0c1ab0c6ab..632ae9d8c2 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -104,7 +104,11 @@ def process_model(self, model, inplace=True):
 
         # Prepare discretisation
         # set variables (we require the full variable not just id)
-        variables = list(model.rhs.keys()) + list(model.algebraic.keys())
+        variables = (
+            list(model.rhs.keys())
+            + list(model.algebraic.keys())
+            + model.external_variables
+        )
 
         # Set the y split for variables
         pybamm.logger.info("Set variable slices for {}".format(model.name))
@@ -132,6 +136,8 @@ def process_model(self, model, inplace=True):
 
         model_disc.bcs = self.bcs
 
+        self.external_variables = model.external_variables
+
         # Process initial condtions
         pybamm.logger.info("Discretise initial conditions for {}".format(model.name))
         ics, concat_ics = self.process_initial_conditions(model)
@@ -795,7 +801,11 @@ def _concatenate_in_order(self, var_eqn_dict, check_complete=False, sparse=False
         if check_complete:
             # Check keys from the given var_eqn_dict against self.y_slices
             ids = {v.id for v in unpacked_variables}
-            if ids != self.y_slices.keys():
+            external_id = {v.id for v in self.external_variables}
+            y_slices_with_external_removed = set(self.y_slices.keys()).difference(
+                external_id
+            )
+            if ids != y_slices_with_external_removed:
                 given_variable_names = [v.name for v in var_eqn_dict.keys()]
                 raise pybamm.ModelError(
                     "Initial conditions are insufficient. Only "
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 9d192af47a..ebaef339b6 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -371,7 +371,12 @@ def check_well_determined(self, post_discretisation):
         # If any variables in the equations don't appear in the keys then the model is
         # underdetermined
         vars_in_keys = vars_in_rhs_keys.union(vars_in_algebraic_keys)
-        extra_variables = vars_in_eqns.difference(vars_in_keys)
+        extra_variables_in_equations = vars_in_eqns.difference(vars_in_keys)
+
+        # get ids of external variables
+        external_ids = {var.id for var in self.external_variables}
+        extra_variables = extra_variables_in_equations.difference(external_ids)
+
         if extra_variables:
             raise pybamm.ModelError("model is underdetermined (too many variables)")
 
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 6bd4a3b6d7..a560dc5e3f 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -57,6 +57,11 @@ class BaseBatteryModel(pybamm.BaseModel):
                 only takes effect if "dimensionality" is 0. If "dimensionality"
                 is 1 or 2 current collector effects are always included. Must be 'False'
                 for lead-acid models.
+            * "external submodels" : list
+                A list of the submodels that you would like to supply an external
+                variable for instead of solving in PyBaMM. The entries of the lists
+                are strings that correspond to the submodel names in the keys
+                of `self.submodels`.
 
 
     **Extends:** :class:`pybamm.BaseModel`
@@ -68,6 +73,7 @@ def __init__(self, options=None, name="Unnamed battery model"):
         self.set_standard_output_variables()
         self.submodels = {}
         self._built = False
+        self._built_fundamental_and_external = False
 
     @property
     def default_parameter_values(self):
@@ -157,6 +163,7 @@ def options(self, extra_options):
             "particle": "Fickian diffusion",
             "thermal": "isothermal",
             "thermal current collector": False,
+            "external submodels": []
         }
         options = default_options
         # any extra options overwrite the default options
@@ -396,17 +403,7 @@ def set_standard_output_variables(self):
                 {"y": var.y, "y [m]": var.y * L_y, "z": var.z, "z [m]": var.z * L_z}
             )
 
-    def build_model(self):
-
-        # Check if already built
-        if self._built:
-            raise pybamm.ModelError(
-                """Model already built. If you are adding a new submodel, try using
-                `model.update` instead."""
-            )
-
-        pybamm.logger.info("Building {}".format(self.name))
-
+    def build_fundamental_and_external(self):
         # Get the fundamental variables
         for submodel_name, submodel in self.submodels.items():
             pybamm.logger.debug(
@@ -416,7 +413,11 @@ def build_model(self):
             )
             self.variables.update(submodel.get_fundamental_variables())
 
-        # Get any external variables
+        # set the submodels that are external
+        for sub in self.options["external submodels"]:
+            self.submodels[sub].external = True
+
+        # Set any external variables
         self.external_variables = []
         for submodel_name, submodel in self.submodels.items():
             pybamm.logger.debug(
@@ -424,12 +425,13 @@ def build_model(self):
                     submodel_name, self.name
                 )
             )
-            variables, external_variables = submodel.get_external_variables()
+            external_variables = submodel.get_external_variables()
 
-            self.variables.update(variables)
             self.external_variables += external_variables
 
-        # Get coupled variables
+        self._built_fundamental_and_external = True
+
+    def build_coupled_variables(self):
         # Note: pybamm will try to get the coupled variables for the submodels in the
         # order they are set by the user. If this fails for a particular submodel,
         # return to it later and try again. If setting coupled variables fails and
@@ -462,36 +464,56 @@ def build_model(self):
                             # try setting coupled variables on next loop through
                             pass
 
+    def build_model_equations(self):
         # Set model equations
         for submodel_name, submodel in self.submodels.items():
-            pybamm.logger.debug(
-                "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
-            )
+            if submodel.external is False:
+                pybamm.logger.debug(
+                    "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
+                )
 
-            submodel.set_rhs(self.variables)
-            pybamm.logger.debug(
-                "Setting algebraic for {} submodel ({})".format(
-                    submodel_name, self.name
+                submodel.set_rhs(self.variables)
+                pybamm.logger.debug(
+                    "Setting algebraic for {} submodel ({})".format(
+                        submodel_name, self.name
+                    )
                 )
-            )
-            submodel.set_algebraic(self.variables)
-            pybamm.logger.debug(
-                "Setting boundary conditions for {} submodel ({})".format(
-                    submodel_name, self.name
+                submodel.set_algebraic(self.variables)
+                pybamm.logger.debug(
+                    "Setting boundary conditions for {} submodel ({})".format(
+                        submodel_name, self.name
+                    )
                 )
-            )
-            submodel.set_boundary_conditions(self.variables)
-            pybamm.logger.debug(
-                "Setting initial conditions for {} submodel ({})".format(
-                    submodel_name, self.name
+                submodel.set_boundary_conditions(self.variables)
+                pybamm.logger.debug(
+                    "Setting initial conditions for {} submodel ({})".format(
+                        submodel_name, self.name
+                    )
                 )
+                submodel.set_initial_conditions(self.variables)
+                submodel.set_events(self.variables)
+                pybamm.logger.debug(
+                    "Updating {} submodel ({})".format(submodel_name, self.name)
+                )
+                self.update(submodel)
+
+    def build_model(self):
+
+        # Check if already built
+        if self._built:
+            raise pybamm.ModelError(
+                """Model already built. If you are adding a new submodel, try using
+                `model.update` instead."""
             )
-            submodel.set_initial_conditions(self.variables)
-            submodel.set_events(self.variables)
-            pybamm.logger.debug(
-                "Updating {} submodel ({})".format(submodel_name, self.name)
-            )
-            self.update(submodel)
+
+        pybamm.logger.info("Building {}".format(self.name))
+
+        if self._built_fundamental_and_external is False:
+            self.build_fundamental_and_external()
+
+        self.build_coupled_variables()
+
+        self.build_model_equations()
 
         pybamm.logger.debug("Setting voltage variables")
         self.set_voltage_variables()
diff --git a/pybamm/models/submodels/base_submodel.py b/pybamm/models/submodels/base_submodel.py
index 25a1aaf7c1..07df5c1e3d 100644
--- a/pybamm/models/submodels/base_submodel.py
+++ b/pybamm/models/submodels/base_submodel.py
@@ -46,7 +46,7 @@ class BaseSubModel:
         symbols.
     """
 
-    def __init__(self, param, domain=None, reactions=None):
+    def __init__(self, param, domain=None, reactions=None, external=False):
         super().__init__()
         self.param = param
         # Initialise empty variables (to avoid overwriting with 'None')
@@ -62,6 +62,8 @@ def __init__(self, param, domain=None, reactions=None):
         self.set_domain_for_broadcast()
         self.reactions = reactions
 
+        self.external = external
+
     @property
     def domain(self):
         return self._domain
@@ -101,20 +103,40 @@ def get_fundamental_variables(self):
         """
         return {}
 
-    def set_external_variables(self):
+    def get_external_variables(self):
         """
-        A public method that creates and returns the variables in a submodel which are
-        suppled external to the model.
+        A public method that returns the variables in a submodel which are
+        supplied by an external source.
 
         Returns
         -------
-        dict :
-            The variables created by the submodel which are independent of variables in
-            other submodels.
         list :
             A list of the external variables in the model.
         """
-        return {}, []
+
+        external_variables = []
+        list_of_vars = []
+
+        if self.external is True:
+            # look through all the variables in the submodel and get the
+            # variables which are state vectors
+            submodel_variables = self.get_fundamental_variables()
+            for var in submodel_variables.values():
+                if isinstance(var, pybamm.Variable):
+                    list_of_vars += [var]
+                elif isinstance(var, pybamm.Concatenation):
+                    for child in var.children:
+                        if isinstance(child, pybamm.Variable):
+                            list_of_vars += [child]
+
+            # remove duplicates
+            unique_ids = []
+            for var in list_of_vars:
+                if var.id not in unique_ids:
+                    external_variables += [var]
+                    unique_ids += [var.id]
+
+        return external_variables
 
     def get_coupled_variables(self, variables):
         """
diff --git a/pybamm/models/submodels/thermal/external/__init__.py b/pybamm/models/submodels/thermal/external/__init__.py
deleted file mode 100644
index e69de29bb2..0000000000
diff --git a/pybamm/models/submodels/thermal/external/full.py b/pybamm/models/submodels/thermal/external/full.py
deleted file mode 100644
index 4c131ebbe3..0000000000
--- a/pybamm/models/submodels/thermal/external/full.py
+++ /dev/null
@@ -1,32 +0,0 @@
-#
-# A class for full external thermal models
-#
-import pybamm
-
-from ..base_thermal import BaseThermal
-
-
-class Full(BaseThermal):
-    """Class to link full external thermal submodels.
-
-    Parameters
-    ----------
-    param : parameter class
-        The parameters to use for this submodel
-
-
-    **Extends:** :class:`pybamm.thermal.BaseModel`
-    """
-
-    def __init__(self, param):
-        super().__init__(param)
-
-    def get_external_variables(self):
-        T = pybamm.standard_variables.T
-        T_cn = pybamm.BoundaryValue(T, "left")
-        T_cp = pybamm.BoundaryValue(T, "right")
-
-        variables = self._get_standard_fundamental_variables(T, T_cn, T_cp)
-        external_variables = [T]
-
-        return variables, external_variables
diff --git a/pybamm/models/submodels/thermal/external/lumped.py b/pybamm/models/submodels/thermal/external/lumped.py
deleted file mode 100644
index 7da09c5c2b..0000000000
--- a/pybamm/models/submodels/thermal/external/lumped.py
+++ /dev/null
@@ -1,39 +0,0 @@
-#
-# A class for external lumped thermal models
-#
-import pybamm
-
-from ..base_thermal import BaseThermal
-
-
-class Full(BaseThermal):
-    """Class to link external lumped thermal submodels.
-
-    Parameters
-    ----------
-    param : parameter class
-        The parameters to use for this submodel
-
-
-    **Extends:** :class:`pybamm.thermal.BaseModel`
-    """
-
-    def __init__(self, param):
-        super().__init__(param)
-
-    def get_external_variables(self):
-        T_x_av = pybamm.standard_variables.T_av
-        T_n = pybamm.PrimaryBroadcast(T_x_av, "negative electrode")
-        T_s = pybamm.PrimaryBroadcast(T_x_av, "separator")
-        T_p = pybamm.PrimaryBroadcast(T_x_av, "positive electrode")
-        T = pybamm.Concatenation(T_n, T_s, T_p)
-
-        T_cn = T_x_av
-        T_cp = T_x_av
-
-        variables = self._get_standard_fundamental_variables(T, T_cn, T_cp)
-        variables = self._get_standard_fundamental_variables(T, T_cn, T_cp)
-
-        external_variables = [T_x_av]
-
-        return variables, external_variables
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 774a5750b8..5c6bf34a2a 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -33,7 +33,7 @@ def __init__(
         self._solver = solver or self._model.default_solver
         self._quick_plot_vars = quick_plot_vars
 
-        self.reset()
+        self.reset(update_model=False)
 
     def set_defaults(self):
         """
@@ -48,11 +48,12 @@ def set_defaults(self):
         self._solver = self._model.default_solver
         self._quick_plot_vars = None
 
-    def reset(self):
+    def reset(self, update_model=True):
         """
         A method to reset a simulation back to its unprocessed state.
         """
-        self.model = self._model_class(self._model_options)
+        if update_model:
+            self.model = self._model_class(self._model_options)
         self.geometry = copy.deepcopy(self._unprocessed_geometry)
         self._model_with_set_params = None
         self._built_model = None
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
new file mode 100644
index 0000000000..aa2e03a3b0
--- /dev/null
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
@@ -0,0 +1,25 @@
+#
+# Tests for external submodels
+#
+import pybamm
+import unittest
+
+
+# only works for statevector inputs
+class TestExternalSubmodel(unittest.TestCase):
+    def test_external_temperature(self):
+
+        model_options = {"thermal": "x-full", "external submodels": ["thermal"]}
+        model = pybamm.lithium_ion.SPMe(model_options)
+
+        sim = pybamm.Simulation(model)
+        sim.build()
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/unit/test_models/test_submodels/standard_submodel_unit_tests.py b/tests/unit/test_models/test_submodels/standard_submodel_unit_tests.py
index 033fb95d86..f0d45ae51f 100644
--- a/tests/unit/test_models/test_submodels/standard_submodel_unit_tests.py
+++ b/tests/unit/test_models/test_submodels/standard_submodel_unit_tests.py
@@ -14,10 +14,15 @@ def __init__(self, submodel, variables=None):
             variables = {}
         self.submodel = submodel
         self.variables = variables
+        self.external_variables = []
 
     def test_get_fundamental_variables(self):
         self.variables.update(self.submodel.get_fundamental_variables())
 
+    def test_get_external_variables(self):
+        external_variables = self.submodel.get_external_variables()
+        self.external_variables += external_variables
+
     def test_get_coupled_variables(self):
         self.variables.update(self.submodel.get_coupled_variables(self.variables))
 
@@ -38,6 +43,7 @@ def test_set_events(self):
 
     def test_all(self):
         self.test_get_fundamental_variables()
+        self.test_get_external_variables()
         self.test_get_coupled_variables()
         self.test_set_rhs()
         self.test_set_algebraic()

From 065d6e01678e09b6f0df67ce425c6c013168a966 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Sat, 9 Nov 2019 13:15:36 +0000
Subject: [PATCH 03/31] #692 build seems to work

---
 pybamm/discretisations/discretisation.py | 36 +++++++++++++++++++++++-
 pybamm/models/base_model.py              |  5 ++++
 pybamm/models/submodels/base_submodel.py | 19 ++++++++++---
 pybamm/spatial_methods/finite_volume.py  | 35 +++++++++++++++++++++++
 4 files changed, 90 insertions(+), 5 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 632ae9d8c2..559d6291d7 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -114,6 +114,10 @@ def process_model(self, model, inplace=True):
         pybamm.logger.info("Set variable slices for {}".format(model.name))
         self.set_variable_slices(variables)
 
+        # now add extrapolated external variables to the boundary conditions
+        # if required by the spatial method
+        self._preprocess_external_variables(model)
+
         # set boundary conditions (only need key ids for boundary_conditions)
         pybamm.logger.info("Discretise boundary conditions for {}".format(model.name))
         self.bcs = self.process_boundary_conditions(model)
@@ -223,6 +227,30 @@ def set_variable_slices(self, variables):
         # reset discretised_symbols
         self._discretised_symbols = {}
 
+    def _preprocess_external_variables(self, model):
+        """
+        A method to preprocess external variables so that they are
+        compatible with the spatial method. For example, in finite
+        volume, the user will supply a vector of values valid on the
+        cell centres. However, for model processing, we also require
+        the boundary edge fluxes. Therefore, we extrapolate and add
+        the boundary fluxes to the boundary conditions, which are
+        employed in generating the grad and div matrices.
+        The processing is delegated to spatial methods as
+        the preprocessing required for finite volume and finite
+        element will be different.
+        """
+
+        for var in model.external_variables:
+            if var.domain == []:
+                pass
+            else:
+                new_bcs = self.spatial_methods[
+                    var.domain[0]
+                ].preprocess_external_variables(var)
+
+                model.boundary_conditions.update(new_bcs)
+
     def set_internal_boundary_conditions(self, model):
         """
         A method to set the internal boundary conditions for the submodel.
@@ -615,7 +643,10 @@ def process_symbol(self, symbol):
         except KeyError:
             discretised_symbol = self._process_symbol(symbol)
             self._discretised_symbols[symbol.id] = discretised_symbol
-            discretised_symbol.test_shape()
+            try:
+                discretised_symbol.test_shape()
+            except:
+                discretised_symbol.test_shape()
             return discretised_symbol
 
     def _process_symbol(self, symbol):
@@ -802,6 +833,9 @@ def _concatenate_in_order(self, var_eqn_dict, check_complete=False, sparse=False
             # Check keys from the given var_eqn_dict against self.y_slices
             ids = {v.id for v in unpacked_variables}
             external_id = {v.id for v in self.external_variables}
+            for var in self.external_variables:
+                child_ids = {child.id for child in var.children}
+                external_id = external_id.union(child_ids)
             y_slices_with_external_removed = set(self.y_slices.keys()).difference(
                 external_id
             )
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index ebaef339b6..7968365e3e 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -375,6 +375,11 @@ def check_well_determined(self, post_discretisation):
 
         # get ids of external variables
         external_ids = {var.id for var in self.external_variables}
+        for var in self.external_variables:
+            if isinstance(var, pybamm.Concatenation):
+                child_ids = {child.id for child in var.children}
+                external_ids = external_ids.union(child_ids)
+
         extra_variables = extra_variables_in_equations.difference(external_ids)
 
         if extra_variables:
diff --git a/pybamm/models/submodels/base_submodel.py b/pybamm/models/submodels/base_submodel.py
index 07df5c1e3d..d90c864882 100644
--- a/pybamm/models/submodels/base_submodel.py
+++ b/pybamm/models/submodels/base_submodel.py
@@ -124,13 +124,24 @@ def get_external_variables(self):
             for var in submodel_variables.values():
                 if isinstance(var, pybamm.Variable):
                     list_of_vars += [var]
+
                 elif isinstance(var, pybamm.Concatenation):
-                    for child in var.children:
-                        if isinstance(child, pybamm.Variable):
-                            list_of_vars += [child]
+                    if all(
+                        isinstance(child, pybamm.Variable) for child in var.children
+                    ):
+                        list_of_vars += [var]
 
-            # remove duplicates
+            # first add only unique concatenations
             unique_ids = []
+            for var in list_of_vars:
+                if var.id not in unique_ids and isinstance(var, pybamm.Concatenation):
+                    external_variables += [var]
+                    unique_ids += [var.id]
+                    # also add the ids of the children to unique ids
+                    for child in var.children:
+                        unique_ids += [child.id]
+
+            # now add any unique variables that are not part of a concatentation
             for var in list_of_vars:
                 if var.id not in unique_ids:
                     external_variables += [var]
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 9698b4ffd2..530f77b0bb 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -83,6 +83,41 @@ def gradient(self, symbol, discretised_symbol, boundary_conditions):
         out = gradient_matrix @ discretised_symbol
         return out
 
+    def preprocess_external_variables(self, var):
+        """
+        For finite volumes, we need the boundary fluxes for discretising
+        properly. Here, we extrapolate and then add them to the boundary
+        conditions.
+
+        Parameters
+        ----------
+        var : :class:`pybamm.Variable` or :class:`pybamm.Concatenation`
+            The external variable that is to be processed
+
+        Returns
+        -------
+        new_bcs: dict
+            A dictionary containing the new boundary conditions
+        """
+
+        if isinstance(var, pybamm.Concatenation):
+            left_var = var.children[0]
+            right_var = var.children[-1]
+        else:
+            left_var = var
+            right_var = var
+
+        print("hello")
+
+        new_bcs = {
+            var: {
+                "left": (pybamm.BoundaryGradient(var, "left"), "Neumann"),
+                "right": (pybamm.BoundaryGradient(var, "right"), "Neumann"),
+            }
+        }
+
+        return new_bcs
+
     def gradient_matrix(self, domain):
         """
         Gradient matrix for finite volumes in the appropriate domain.

From 9fe3a3dedf07a6e07e27dc13daf6689ff90a6905 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Sat, 9 Nov 2019 17:05:15 +0000
Subject: [PATCH 04/31] #692 adding features to base solver to allow external
 variables

---
 pybamm/discretisations/discretisation.py      | 12 +++++++-
 pybamm/simulation.py                          | 21 ++++++++++++++
 pybamm/solvers/base_solver.py                 | 20 ++++++++++++-
 pybamm/spatial_methods/finite_volume.py       |  9 ------
 pybamm/spatial_methods/spatial_method.py      |  3 ++
 .../test_use_external_submodel.py             | 29 +++++++++++++++++--
 6 files changed, 81 insertions(+), 13 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 559d6291d7..55e0416e28 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -141,8 +141,17 @@ def process_model(self, model, inplace=True):
         model_disc.bcs = self.bcs
 
         self.external_variables = model.external_variables
+        # find where external variables begin in state vector
+        start_vals = []
+        for var in self.external_variables:
+            if isinstance(var, pybamm.Concatenation):
+                for child in var.children:
+                    start_vals += [self.y_slices[child.id][0].start]
+            elif isinstance(var, pybamm.Variable):
+                start_vals += [self.y_slices[var.id][0].start]
+
+        self.external_start = min(start_vals)
 
-        # Process initial condtions
         pybamm.logger.info("Discretise initial conditions for {}".format(model.name))
         ics, concat_ics = self.process_initial_conditions(model)
         model_disc.initial_conditions = ics
@@ -223,6 +232,7 @@ def set_variable_slices(self, variables):
                 start = end
 
         self.y_slices = y_slices
+        self.y_length = end
 
         # reset discretised_symbols
         self._discretised_symbols = {}
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 5c6bf34a2a..01674881f3 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -115,6 +115,27 @@ def solve(self, t_eval=None, solver=None):
 
         self._solution = solver.solve(self.built_model, t_eval)
 
+    def step(self, dt, solver=None, external_variables=None):
+        """
+        A method to step the model forward one timestep. This method will
+        automatically build and set the model parameters if not already done so.
+
+        Parameters
+        ----------
+        dt : numeric type
+            The timestep over which to step the solution
+        solver : :class:`pybamm.BaseSolver`
+            The solver to use to solve the model.
+        external_variables : dict
+            A dictionary of external variables and their corresponding
+            values at the current time
+        """
+
+        if solver is None:
+            solver = self.solver
+
+        self._solution = solver.step(self.built_model, dt, external_variables)
+
     def plot(self, quick_plot_vars=None):
         """
         A method to quickly plot the outputs of the simulation.
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 37f2fd819d..690086aca4 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -97,7 +97,7 @@ def solve(self, model, t_eval):
         )
         return solution
 
-    def step(self, model, dt, npts=2, log=True):
+    def step(self, model, dt, npts=2, log=True, external_variables=None):
         """
         Step the solution of the model forward by a given time increment. The
         first time this method is called it executes the necessary setup by
@@ -113,6 +113,9 @@ def step(self, model, dt, npts=2, log=True):
         npts : int, optional
             The number of points at which the solution will be returned during
             the step dt. default is 2 (returns the solution at t0 and t0 + dt).
+        external_variables : dict
+            A dictionary of external variables and their corresponding
+            values at the current time
 
         Raises
         ------
@@ -147,6 +150,21 @@ def step(self, model, dt, npts=2, log=True):
         else:
             set_up_time = None
 
+        # load external variables into a state vector
+        self.y_ext = np.zeros((model.y_length, 1))
+        for var_name, var_vals in external_variables.items():
+            var = model.variables[var_name]
+            if isinstance(var, pybamm.Concatenation):
+                start = model.y_slices[var.children[0].id][0].start
+                stop = model.y_slices[var.children[-1].id][-1].stop
+                y_slice = slice(start, stop)
+
+            elif isinstance(var, pybamm.Variable):
+                start = model.y_slices[var.id][0].start
+                stop = model.y_slices[var.id][-1].stop
+                y_slice = slice(start, stop)
+            self.y_ext[y_slice] = var_vals
+
         # Step
         t_eval = np.linspace(self.t, self.t + dt, npts)
         solution, solve_time, termination = self.compute_solution(model, t_eval)
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 530f77b0bb..cc0f1d93fd 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -100,15 +100,6 @@ def preprocess_external_variables(self, var):
             A dictionary containing the new boundary conditions
         """
 
-        if isinstance(var, pybamm.Concatenation):
-            left_var = var.children[0]
-            right_var = var.children[-1]
-        else:
-            left_var = var
-            right_var = var
-
-        print("hello")
-
         new_bcs = {
             var: {
                 "left": (pybamm.BoundaryGradient(var, "left"), "Neumann"),
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 1a8fbd7a24..eba3e77c24 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -278,6 +278,9 @@ def internal_neumann_condition(
 
         raise NotImplementedError
 
+    def preprocess_external_variables(self, var):
+        return {}
+
     def boundary_value_or_flux(self, symbol, discretised_child):
         """
         Returns the boundary value or flux using the approriate expression for the
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
index aa2e03a3b0..f1e1976891 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
@@ -3,18 +3,43 @@
 #
 import pybamm
 import unittest
+import numpy as np
 
 
 # only works for statevector inputs
 class TestExternalSubmodel(unittest.TestCase):
     def test_external_temperature(self):
 
-        model_options = {"thermal": "x-full", "external submodels": ["thermal"]}
+        model_options = {
+            "thermal": "x-full",
+            "external submodels": ["thermal"],
+        }
+
         model = pybamm.lithium_ion.SPMe(model_options)
 
-        sim = pybamm.Simulation(model)
+        neg_pts = 5
+        sep_pts = 3
+        pos_pts = 5
+        tot_pts = neg_pts + sep_pts + pos_pts
+
+        var_pts = {
+            pybamm.standard_spatial_vars.x_n: neg_pts,
+            pybamm.standard_spatial_vars.x_s: sep_pts,
+            pybamm.standard_spatial_vars.x_p: pos_pts,
+            pybamm.standard_spatial_vars.r_n: 5,
+            pybamm.standard_spatial_vars.r_p: 5,
+        }
+
+        sim = pybamm.Simulation(model, var_pts=var_pts)
         sim.build()
 
+        t_eval = np.linspace(0, 0.17, 100)
+
+        for t in t_eval:
+            T = np.zeros((tot_pts, 1))
+            external_variables = {"Cell temperature": T}
+            sim.step(external_variables=external_variables)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

From 97f9a67d2b0aeeb04e212f64498bb0e4ab6506fe Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Sat, 9 Nov 2019 18:10:58 +0000
Subject: [PATCH 05/31] #692 structure for external variables in place but need
 to debug

---
 pybamm/discretisations/discretisation.py      |  5 ++++-
 pybamm/simulation.py                          |  4 +++-
 pybamm/solvers/base_solver.py                 | 19 +++++++++++++++++++
 pybamm/solvers/dae_solver.py                  |  2 ++
 pybamm/solvers/ode_solver.py                  |  6 ++++++
 .../test_use_external_submodel.py             |  5 +++--
 6 files changed, 37 insertions(+), 4 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 55e0416e28..57e5726c95 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -150,7 +150,10 @@ def process_model(self, model, inplace=True):
             elif isinstance(var, pybamm.Variable):
                 start_vals += [self.y_slices[var.id][0].start]
 
-        self.external_start = min(start_vals)
+        model_disc.external_variables = model.external_variables
+        model_disc.external_start = min(start_vals)
+        model_disc.y_length = self.y_length
+        model_disc.y_slices = self.y_slices
 
         pybamm.logger.info("Discretise initial conditions for {}".format(model.name))
         ics, concat_ics = self.process_initial_conditions(model)
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 01674881f3..148c83976b 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -134,7 +134,9 @@ def step(self, dt, solver=None, external_variables=None):
         if solver is None:
             solver = self.solver
 
-        self._solution = solver.step(self.built_model, dt, external_variables)
+        self._solution = solver.step(
+            self.built_model, dt, external_variables=external_variables
+        )
 
     def plot(self, quick_plot_vars=None):
         """
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 690086aca4..17e9f38a28 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -22,6 +22,9 @@ def __init__(self, method=None, rtol=1e-6, atol=1e-6):
         self._atol = atol
         self.name = "Base solver"
 
+        self.y_pad = None
+        self.y_ext = None
+
     @property
     def method(self):
         return self._method
@@ -147,9 +150,16 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
                 self.set_up(model)
             self.t = 0.0
             set_up_time = timer.time()
+
+            # create a y_pad vector of the correct size:
+            self.y_pad = np.zeros((model.external_start))
+
         else:
             set_up_time = None
 
+        if external_variables is None:
+            external_variables = {}
+
         # load external variables into a state vector
         self.y_ext = np.zeros((model.y_length, 1))
         for var_name, var_vals in external_variables.items():
@@ -193,6 +203,15 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
             )
         return solution
 
+    def add_external(self, y):
+        """
+        Pad the state vector and then add the external variables so that
+        it is of the correct shape for evaluate
+        """
+        if self.y_pad and self.y_ext:
+            y = np.concatenate(y, self.y_pad) + self.y_ext
+        return y
+
     def compute_solution(self, model, t_eval):
         """Calculate the solution of the model at specified times. Note: this
         does *not* execute the solver setup.
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 74001b7b9c..1839329e45 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -188,6 +188,7 @@ def residuals(t, y, ydot):
             pybamm.logger.debug(
                 "Evaluating residuals for {} at t={}".format(model.name, t)
             )
+            y = self.add_external(y)
             y = y[:, np.newaxis]
             rhs_eval, known_evals = concatenated_rhs.evaluate(t, y, known_evals={})
             # reuse known_evals
@@ -202,6 +203,7 @@ def residuals(t, y, ydot):
         # Create event-dependent function to evaluate events
         def event_fun(event):
             def eval_event(t, y):
+                y = self.add_external(y)
                 return event.evaluate(t, y)
 
             return eval_event
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index a1655d745f..bc779b24ec 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -121,6 +121,7 @@ def set_up(self, model):
         # Create function to evaluate rhs
         def dydt(t, y):
             pybamm.logger.debug("Evaluating RHS for {} at t={}".format(model.name, t))
+            y = self.add_external(y)
             y = y[:, np.newaxis]
             dy = concatenated_rhs.evaluate(t, y, known_evals={})[0]
             return dy[:, 0]
@@ -128,6 +129,7 @@ def dydt(t, y):
         # Create event-dependent function to evaluate events
         def event_fun(event):
             def eval_event(t, y):
+                y = self.add_external(y)
                 return event.evaluate(t, y)
 
             return eval_event
@@ -138,6 +140,7 @@ def eval_event(t, y):
         if jac_rhs is not None:
 
             def jacobian(t, y):
+                y = self.add_external(y)
                 return jac_rhs.evaluate(t, y, known_evals={})[0]
 
         else:
@@ -193,6 +196,7 @@ def set_up_casadi(self, model):
 
         def dydt(t, y):
             pybamm.logger.debug("Evaluating RHS for {} at t={}".format(model.name, t))
+            y = self.add_external(y)
             dy = concatenated_rhs_fn(t, y).full()
             return dy[:, 0]
 
@@ -201,6 +205,7 @@ def event_fun(event):
             casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi], [event])
 
             def eval_event(t, y):
+                y = self.add_external(y)
                 return casadi_event_fn(t, y)
 
             return eval_event
@@ -216,6 +221,7 @@ def eval_event(t, y):
             )
 
             def jacobian(t, y):
+                y = self.add_external(y)
                 return casadi_jac_fn(t, y)
 
         else:
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
index f1e1976891..fee8cd3582 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
@@ -35,10 +35,11 @@ def test_external_temperature(self):
 
         t_eval = np.linspace(0, 0.17, 100)
 
-        for t in t_eval:
+        for i, t in enumerate(t_eval):
+            dt = t_eval[i + 1] - t_eval[i]
             T = np.zeros((tot_pts, 1))
             external_variables = {"Cell temperature": T}
-            sim.step(external_variables=external_variables)
+            sim.step(dt, external_variables=external_variables)
 
 
 if __name__ == "__main__":

From 99aa4ea6bc0ccee7d8501ba3cdee81cd7d64e2a0 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Sat, 9 Nov 2019 18:33:10 +0000
Subject: [PATCH 06/31] #692 fixed yslice issue

---
 pybamm/solvers/base_solver.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 17e9f38a28..6d7eef4846 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -165,13 +165,13 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
         for var_name, var_vals in external_variables.items():
             var = model.variables[var_name]
             if isinstance(var, pybamm.Concatenation):
-                start = model.y_slices[var.children[0].id][0].start
-                stop = model.y_slices[var.children[-1].id][-1].stop
+                start = var.children[0].y_slices[0].start
+                stop = var.children[-1].y_slices[-1].stop
                 y_slice = slice(start, stop)
 
             elif isinstance(var, pybamm.Variable):
-                start = model.y_slices[var.id][0].start
-                stop = model.y_slices[var.id][-1].stop
+                start = var.y_slices[0].start
+                stop = var.y_slices[-1].stop
                 y_slice = slice(start, stop)
             self.y_ext[y_slice] = var_vals
 
@@ -208,8 +208,8 @@ def add_external(self, y):
         Pad the state vector and then add the external variables so that
         it is of the correct shape for evaluate
         """
-        if self.y_pad and self.y_ext:
-            y = np.concatenate(y, self.y_pad) + self.y_ext
+        if self.y_pad is not None and self.y_ext is not None:
+            y = np.concatenate([y, self.y_pad]) + self.y_ext
         return y
 
     def compute_solution(self, model, t_eval):

From a2834426528eeb39e3c95c38f1aa613c8e0c79da Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Sat, 9 Nov 2019 19:23:44 +0000
Subject: [PATCH 07/31] #692 simulation runs

---
 pybamm/solvers/base_solver.py                             | 5 +++--
 pybamm/solvers/ode_solver.py                              | 8 +++++++-
 .../test_lithium_ion/test_use_external_submodel.py        | 4 +++-
 3 files changed, 13 insertions(+), 4 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 6d7eef4846..78f26fa00c 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -152,7 +152,7 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
             set_up_time = timer.time()
 
             # create a y_pad vector of the correct size:
-            self.y_pad = np.zeros((model.external_start))
+            self.y_pad = np.zeros((model.y_length - model.external_start, 1))
 
         else:
             set_up_time = None
@@ -272,8 +272,9 @@ def get_termination_reason(self, solution, events):
             # Get final event value
             final_event_values = {}
             for name, event in events.items():
+                y_event = self.add_external(solution.y_event)
                 final_event_values[name] = abs(
-                    event.evaluate(solution.t_event, solution.y_event)
+                    event.evaluate(solution.t_event, y_event)
                 )
             termination_event = min(final_event_values, key=final_event_values.get)
             # Add the event to the solution object
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index bc779b24ec..1cd00681ca 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -45,6 +45,7 @@ def compute_solution(self, model, t_eval):
             mass_matrix=model.mass_matrix.entries,
             jacobian=self.jacobian,
         )
+
         solve_time = timer.time() - solve_start_time
 
         # Identify the event that caused termination
@@ -121,14 +122,15 @@ def set_up(self, model):
         # Create function to evaluate rhs
         def dydt(t, y):
             pybamm.logger.debug("Evaluating RHS for {} at t={}".format(model.name, t))
-            y = self.add_external(y)
             y = y[:, np.newaxis]
+            y = self.add_external(y)
             dy = concatenated_rhs.evaluate(t, y, known_evals={})[0]
             return dy[:, 0]
 
         # Create event-dependent function to evaluate events
         def event_fun(event):
             def eval_event(t, y):
+                y = y[:, np.newaxis]
                 y = self.add_external(y)
                 return event.evaluate(t, y)
 
@@ -140,6 +142,7 @@ def eval_event(t, y):
         if jac_rhs is not None:
 
             def jacobian(t, y):
+                y = y[:, np.newaxis]
                 y = self.add_external(y)
                 return jac_rhs.evaluate(t, y, known_evals={})[0]
 
@@ -196,6 +199,7 @@ def set_up_casadi(self, model):
 
         def dydt(t, y):
             pybamm.logger.debug("Evaluating RHS for {} at t={}".format(model.name, t))
+            y = y[:, np.newaxis]
             y = self.add_external(y)
             dy = concatenated_rhs_fn(t, y).full()
             return dy[:, 0]
@@ -205,6 +209,7 @@ def event_fun(event):
             casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi], [event])
 
             def eval_event(t, y):
+                y = y[:, np.newaxis]
                 y = self.add_external(y)
                 return casadi_event_fn(t, y)
 
@@ -221,6 +226,7 @@ def eval_event(t, y):
             )
 
             def jacobian(t, y):
+                y = y[:, np.newaxis]
                 y = self.add_external(y)
                 return casadi_jac_fn(t, y)
 
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
index fee8cd3582..0e731b8506 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
@@ -17,6 +17,8 @@ def test_external_temperature(self):
 
         model = pybamm.lithium_ion.SPMe(model_options)
 
+        # model.convert_to_format = False
+
         neg_pts = 5
         sep_pts = 3
         pos_pts = 5
@@ -35,7 +37,7 @@ def test_external_temperature(self):
 
         t_eval = np.linspace(0, 0.17, 100)
 
-        for i, t in enumerate(t_eval):
+        for i in np.arange(1, len(t_eval) - 1):
             dt = t_eval[i + 1] - t_eval[i]
             T = np.zeros((tot_pts, 1))
             external_variables = {"Cell temperature": T}

From 3e8ec1e995001a0ea4c4e5c60977c8a9cf7ace0e Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Sat, 9 Nov 2019 20:18:29 +0000
Subject: [PATCH 08/31] #692 simulation runs and can save a reasonable looking
 solution however with plotting remain

---
 pybamm/simulation.py                          | 28 +++++++++++++++++--
 pybamm/solvers/base_solver.py                 |  7 ++---
 pybamm/solvers/ode_solver.py                  |  8 ++++++
 .../test_use_external_submodel.py             |  2 ++
 4 files changed, 39 insertions(+), 6 deletions(-)

diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 148c83976b..44918ad0b0 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -33,6 +33,8 @@ def __init__(
         self._solver = solver or self._model.default_solver
         self._quick_plot_vars = quick_plot_vars
 
+        self._made_step = False
+
         self.reset(update_model=False)
 
     def set_defaults(self):
@@ -60,6 +62,7 @@ def reset(self, update_model=True):
         self._mesh = None
         self._disc = None
         self._solution = None
+        self._made_step = False
 
     def set_parameters(self):
         """
@@ -115,7 +118,7 @@ def solve(self, t_eval=None, solver=None):
 
         self._solution = solver.solve(self.built_model, t_eval)
 
-    def step(self, dt, solver=None, external_variables=None):
+    def step(self, dt, solver=None, external_variables=None, save=True):
         """
         A method to step the model forward one timestep. This method will
         automatically build and set the model parameters if not already done so.
@@ -129,15 +132,36 @@ def step(self, dt, solver=None, external_variables=None):
         external_variables : dict
             A dictionary of external variables and their corresponding
             values at the current time
+        save : bool
+            Turn on to store the solution of all previous timesteps
         """
 
         if solver is None:
             solver = self.solver
 
-        self._solution = solver.step(
+        solution = solver.step(
             self.built_model, dt, external_variables=external_variables
         )
 
+        if save is False or self._made_step is False:
+            self._solution = solution
+        else:
+            self.update_solution(solution)
+
+        self._made_step = True
+
+    def update_solution(self, solution):
+
+        self._solution.set_up_time += solution.set_up_time
+        self._solution.solve_time += solution.solve_time
+        self._solution.t = np.append(self._solution.t, solution.t[-1])
+        self._solution.t_event = solution.t_event
+        self._solution.termination = solution.termination
+        self._solution.y = np.concatenate(
+            [self._solution.y, solution.y[:, -1][:, np.newaxis]], axis=1
+        )
+        self._solution.y_event = solution.y_event
+
     def plot(self, quick_plot_vars=None):
         """
         A method to quickly plot the outputs of the simulation.
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 78f26fa00c..3233e2b8cd 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -155,7 +155,7 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
             self.y_pad = np.zeros((model.y_length - model.external_start, 1))
 
         else:
-            set_up_time = None
+            set_up_time = 0
 
         if external_variables is None:
             external_variables = {}
@@ -181,12 +181,11 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
 
         # Assign times
         solution.solve_time = solve_time
-        if set_up_time:
-            solution.set_up_time = set_up_time
+        solution.set_up_time = set_up_time
 
         # Set self.t and self.y0 to their values at the final step
         self.t = solution.t[-1]
-        self.y0 = solution.y[:, -1]
+        self.y0 = solution.y[: model.external_start, -1]
 
         pybamm.logger.debug("Finish stepping {} ({})".format(model.name, termination))
         if set_up_time:
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 1cd00681ca..c777a0091b 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -46,6 +46,14 @@ def compute_solution(self, model, t_eval):
             jacobian=self.jacobian,
         )
 
+        # add the external points onto the solution vector
+        full_y = np.zeros((model.y_length, solution.y.shape[1]))
+        for i in np.arange(solution.y.shape[1]):
+            sol_y = solution.y[:, i]
+            sol_y = sol_y[:, np.newaxis]
+            full_y[:, i] = self.add_external(sol_y)[:, 0]
+        solution.y = full_y
+
         solve_time = timer.time() - solve_start_time
 
         # Identify the event that caused termination
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
index 0e731b8506..0a56a053aa 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
@@ -43,6 +43,8 @@ def test_external_temperature(self):
             external_variables = {"Cell temperature": T}
             sim.step(dt, external_variables=external_variables)
 
+        sim.plot(["Negative particle surface concentration"])
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

From 08c9a4ab39cd960093c71cc5793cecd05c97293b Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Sat, 9 Nov 2019 20:56:23 +0000
Subject: [PATCH 09/31] #692 trying to locate strange quick plot error

---
 pybamm/discretisations/discretisation.py                      | 3 ++-
 pybamm/simulation.py                                          | 4 ++--
 .../test_lithium_ion/test_use_external_submodel.py            | 3 ++-
 3 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 57e5726c95..f58c8cafcb 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -151,7 +151,8 @@ def process_model(self, model, inplace=True):
                 start_vals += [self.y_slices[var.id][0].start]
 
         model_disc.external_variables = model.external_variables
-        model_disc.external_start = min(start_vals)
+        if start_vals:
+            model_disc.external_start = min(start_vals)
         model_disc.y_length = self.y_length
         model_disc.y_slices = self.y_slices
 
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 44918ad0b0..258a4f21b2 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -146,11 +146,11 @@ def step(self, dt, solver=None, external_variables=None, save=True):
         if save is False or self._made_step is False:
             self._solution = solution
         else:
-            self.update_solution(solution)
+            self._update_solution(solution)
 
         self._made_step = True
 
-    def update_solution(self, solution):
+    def _update_solution(self, solution):
 
         self._solution.set_up_time += solution.set_up_time
         self._solution.solve_time += solution.solve_time
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
index 0a56a053aa..cd12ea7c64 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
@@ -43,7 +43,8 @@ def test_external_temperature(self):
             external_variables = {"Cell temperature": T}
             sim.step(dt, external_variables=external_variables)
 
-        sim.plot(["Negative particle surface concentration"])
+        # sim.plot(["Negative particle surface concentration"])
+        sim.plot(["Electrolyte concentration"])
 
 
 if __name__ == "__main__":

From ed6a08334362fdbfd99ed0689bfb267861b89450 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Mon, 11 Nov 2019 12:08:48 +0000
Subject: [PATCH 10/31] #692 prevented repeated timesteps being saved to
 solution

---
 pybamm/simulation.py                               | 14 ++++++++++----
 .../test_lithium_ion/test_use_external_submodel.py |  3 +--
 2 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 258a4f21b2..f07b9fb2d8 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -33,7 +33,7 @@ def __init__(
         self._solver = solver or self._model.default_solver
         self._quick_plot_vars = quick_plot_vars
 
-        self._made_step = False
+        self._made_first_step = False
 
         self.reset(update_model=False)
 
@@ -62,7 +62,7 @@ def reset(self, update_model=True):
         self._mesh = None
         self._disc = None
         self._solution = None
-        self._made_step = False
+        self._made_first_step = False
 
     def set_parameters(self):
         """
@@ -143,12 +143,14 @@ def step(self, dt, solver=None, external_variables=None, save=True):
             self.built_model, dt, external_variables=external_variables
         )
 
-        if save is False or self._made_step is False:
+        if save is False or self._made_first_step is False:
             self._solution = solution
+        elif self._solution.t[-1] == solution.t[-1]:
+            pass
         else:
             self._update_solution(solution)
 
-        self._made_step = True
+        self._made_first_step = True
 
     def _update_solution(self, solution):
 
@@ -239,6 +241,10 @@ def var_pts(self):
     def spatial_methods(self):
         return self._spatial_methods
 
+    @property
+    def mesh(self):
+        return self._mesh
+
     @property
     def solver(self):
         return self._solver
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
index cd12ea7c64..6c8a7e97ec 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
@@ -43,8 +43,7 @@ def test_external_temperature(self):
             external_variables = {"Cell temperature": T}
             sim.step(dt, external_variables=external_variables)
 
-        # sim.plot(["Negative particle surface concentration"])
-        sim.plot(["Electrolyte concentration"])
+        sim.plot()
 
 
 if __name__ == "__main__":

From 121f4793a77a72417af8f77778256befda42925e Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Mon, 11 Nov 2019 12:42:39 +0000
Subject: [PATCH 11/31] #692 working on ode models

---
 pybamm/__init__.py                            |  1 +
 pybamm/expression_tree/exceptions.py          |  7 +++
 pybamm/simulation.py                          |  1 +
 pybamm/solvers/base_solver.py                 |  9 +++-
 .../test_use_external_submodel.py             | 47 +++++++++++++++++--
 5 files changed, 61 insertions(+), 4 deletions(-)

diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 69cda501f7..4de9763e49 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -149,6 +149,7 @@ def version(formatted=False):
     ModelWarning,
     UndefinedOperationError,
     GeometryError,
+    InputError,
 )
 
 # Operations
diff --git a/pybamm/expression_tree/exceptions.py b/pybamm/expression_tree/exceptions.py
index f007201340..a647b8ab1b 100644
--- a/pybamm/expression_tree/exceptions.py
+++ b/pybamm/expression_tree/exceptions.py
@@ -59,3 +59,10 @@ class UndefinedOperationError(Exception):
     """
 
     pass
+
+class InputError(Exception):
+    """
+    An external variable has been input incorrectly into PyBaMM
+    """
+
+    pass
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index f07b9fb2d8..96be32b733 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -135,6 +135,7 @@ def step(self, dt, solver=None, external_variables=None, save=True):
         save : bool
             Turn on to store the solution of all previous timesteps
         """
+        self.build()
 
         if solver is None:
             solver = self.solver
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 3233e2b8cd..e9677c99c9 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -169,10 +169,17 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
                 stop = var.children[-1].y_slices[-1].stop
                 y_slice = slice(start, stop)
 
-            elif isinstance(var, pybamm.Variable):
+            elif isinstance(var, pybamm.StateVector):
                 start = var.y_slices[0].start
                 stop = var.y_slices[-1].stop
                 y_slice = slice(start, stop)
+            else:
+                raise pybamm.InputError(
+                    """The variable you have inputted is not a StateVector or Concatenation
+            of StateVectors. Please check the submodel you have made "external" and
+            ensure that the variable you
+            are passing in is the variable that is solved for in that submodel"""
+                )
             self.y_ext[y_slice] = var_vals
 
         # Step
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
index 6c8a7e97ec..68e57a781d 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
@@ -8,6 +8,35 @@
 
 # only works for statevector inputs
 class TestExternalSubmodel(unittest.TestCase):
+    def test_external_lumped_temperature(self):
+        model_options = {
+            "thermal": "x-lumped",
+            "external submodels": ["thermal"],
+        }
+        model = pybamm.lithium_ion.SPMe(model_options)
+        sim = pybamm.Simulation(model)
+
+        t_eval = np.linspace(0, 0.17, 100)
+
+        T_av = 0
+
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+            external_variables = {"X-averaged cell temperature": T_av}
+            T_av += 1
+            sim.step(dt, external_variables=external_variables)
+
+        sim.plot(
+            [
+                "Negative particle surface concentration",
+                "Electrolyte concentration",
+                "Positive particle surface concentration",
+                "Cell temperature [K]",
+                "Terminal voltage [V]",
+                "Total heating [W.m-3]",
+            ]
+        )
+
     def test_external_temperature(self):
 
         model_options = {
@@ -15,7 +44,7 @@ def test_external_temperature(self):
             "external submodels": ["thermal"],
         }
 
-        model = pybamm.lithium_ion.SPMe(model_options)
+        model = pybamm.lithium_ion.SPM(model_options)
 
         # model.convert_to_format = False
 
@@ -36,14 +65,26 @@ def test_external_temperature(self):
         sim.build()
 
         t_eval = np.linspace(0, 0.17, 100)
+        x = np.linspace(0, 1, tot_pts)
 
         for i in np.arange(1, len(t_eval) - 1):
             dt = t_eval[i + 1] - t_eval[i]
-            T = np.zeros((tot_pts, 1))
+            T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
+                :, np.newaxis
+            ]
             external_variables = {"Cell temperature": T}
             sim.step(dt, external_variables=external_variables)
 
-        sim.plot()
+        sim.plot(
+            [
+                "Negative particle surface concentration",
+                "Electrolyte concentration",
+                "Positive particle surface concentration",
+                "Cell temperature [K]",
+                "Terminal voltage [V]",
+                "Total heating [W.m-3]",
+            ]
+        )
 
 
 if __name__ == "__main__":

From 8979e1978f8f76f06be91c4fce113d59dea7532d Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Mon, 11 Nov 2019 16:37:39 +0000
Subject: [PATCH 12/31] #692 starting to put in dae support

---
 pybamm/solvers/base_solver.py                 |  60 ++++-----
 pybamm/solvers/dae_solver.py                  |  11 +-
 .../test_use_external_submodel.py             | 114 ++++++++++++------
 3 files changed, 122 insertions(+), 63 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index e9677c99c9..d56e8b1a62 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -133,6 +133,12 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
         # Set timer
         timer = pybamm.Timer()
 
+        if not hasattr(self, "y0"):
+            # create a y_pad vector of the correct size:
+            self.y_pad = np.zeros((model.y_length - model.external_start, 1))
+
+        self.set_external_variables(model, external_variables)
+
         # Run set up on first step
         if not hasattr(self, "y0"):
             pybamm.logger.info(
@@ -151,37 +157,9 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
             self.t = 0.0
             set_up_time = timer.time()
 
-            # create a y_pad vector of the correct size:
-            self.y_pad = np.zeros((model.y_length - model.external_start, 1))
-
         else:
             set_up_time = 0
 
-        if external_variables is None:
-            external_variables = {}
-
-        # load external variables into a state vector
-        self.y_ext = np.zeros((model.y_length, 1))
-        for var_name, var_vals in external_variables.items():
-            var = model.variables[var_name]
-            if isinstance(var, pybamm.Concatenation):
-                start = var.children[0].y_slices[0].start
-                stop = var.children[-1].y_slices[-1].stop
-                y_slice = slice(start, stop)
-
-            elif isinstance(var, pybamm.StateVector):
-                start = var.y_slices[0].start
-                stop = var.y_slices[-1].stop
-                y_slice = slice(start, stop)
-            else:
-                raise pybamm.InputError(
-                    """The variable you have inputted is not a StateVector or Concatenation
-            of StateVectors. Please check the submodel you have made "external" and
-            ensure that the variable you
-            are passing in is the variable that is solved for in that submodel"""
-                )
-            self.y_ext[y_slice] = var_vals
-
         # Step
         t_eval = np.linspace(self.t, self.t + dt, npts)
         solution, solve_time, termination = self.compute_solution(model, t_eval)
@@ -209,6 +187,32 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
             )
         return solution
 
+    def set_external_variables(self, model, external_variables):
+        if external_variables is None:
+            external_variables = {}
+
+        # load external variables into a state vector
+        self.y_ext = np.zeros((model.y_length, 1))
+        for var_name, var_vals in external_variables.items():
+            var = model.variables[var_name]
+            if isinstance(var, pybamm.Concatenation):
+                start = var.children[0].y_slices[0].start
+                stop = var.children[-1].y_slices[-1].stop
+                y_slice = slice(start, stop)
+
+            elif isinstance(var, pybamm.StateVector):
+                start = var.y_slices[0].start
+                stop = var.y_slices[-1].stop
+                y_slice = slice(start, stop)
+            else:
+                raise pybamm.InputError(
+                    """The variable you have inputted is not a StateVector or Concatenation
+            of StateVectors. Please check the submodel you have made "external" and
+            ensure that the variable you
+            are passing in is the variable that is solved for in that submodel"""
+                )
+            self.y_ext[y_slice] = var_vals
+
     def add_external(self, y):
         """
         Pad the state vector and then add the external variables so that
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 1839329e45..5268b7f76e 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -146,9 +146,13 @@ def set_up(self, model):
                 jac = pybamm.EvaluatorPython(jac)
 
             def jacobian(t, y):
+                y = y[:, np.newaxis]
+                y = self.add_external(y)
                 return jac.evaluate(t, y, known_evals={})[0]
 
             def jacobian_alg(t, y):
+                y = y[:, np.newaxis]
+                y = self.add_external(y)
                 return jac_algebraic.evaluate(t, y, known_evals={})[0]
 
         else:
@@ -167,9 +171,13 @@ def jacobian_alg(t, y):
 
         # Calculate consistent initial conditions for the algebraic equations
         def rhs(t, y):
+            y = y[:, np.newaxis]
+            y = self.add_external(y)
             return concatenated_rhs.evaluate(t, y, known_evals={})[0][:, 0]
 
         def algebraic(t, y):
+            y = y[:, np.newaxis]
+            y = self.add_external(y)
             return concatenated_algebraic.evaluate(t, y, known_evals={})[0][:, 0]
 
         if len(model.algebraic) > 0:
@@ -188,8 +196,8 @@ def residuals(t, y, ydot):
             pybamm.logger.debug(
                 "Evaluating residuals for {} at t={}".format(model.name, t)
             )
-            y = self.add_external(y)
             y = y[:, np.newaxis]
+            y = self.add_external(y)
             rhs_eval, known_evals = concatenated_rhs.evaluate(t, y, known_evals={})
             # reuse known_evals
             alg_eval = concatenated_algebraic.evaluate(t, y, known_evals=known_evals)[0]
@@ -203,6 +211,7 @@ def residuals(t, y, ydot):
         # Create event-dependent function to evaluate events
         def event_fun(event):
             def eval_event(t, y):
+                y = y[:, np.newaxis]
                 y = self.add_external(y)
                 return event.evaluate(t, y)
 
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
index 68e57a781d..5765404baf 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
@@ -6,47 +6,93 @@
 import numpy as np
 
 
-# only works for statevector inputs
 class TestExternalSubmodel(unittest.TestCase):
-    def test_external_lumped_temperature(self):
-        model_options = {
-            "thermal": "x-lumped",
-            "external submodels": ["thermal"],
-        }
-        model = pybamm.lithium_ion.SPMe(model_options)
-        sim = pybamm.Simulation(model)
-
-        t_eval = np.linspace(0, 0.17, 100)
-
-        T_av = 0
-
-        for i in np.arange(1, len(t_eval) - 1):
-            dt = t_eval[i + 1] - t_eval[i]
-            external_variables = {"X-averaged cell temperature": T_av}
-            T_av += 1
-            sim.step(dt, external_variables=external_variables)
-
-        sim.plot(
-            [
-                "Negative particle surface concentration",
-                "Electrolyte concentration",
-                "Positive particle surface concentration",
-                "Cell temperature [K]",
-                "Terminal voltage [V]",
-                "Total heating [W.m-3]",
-            ]
-        )
-
-    def test_external_temperature(self):
+    # def test_external_lumped_temperature(self):
+    #     model_options = {
+    #         "thermal": "x-lumped",
+    #         "external submodels": ["thermal"],
+    #     }
+    #     model = pybamm.lithium_ion.SPMe(model_options)
+    #     sim = pybamm.Simulation(model)
+
+    #     t_eval = np.linspace(0, 0.17, 100)
+
+    #     T_av = 0
+
+    #     for i in np.arange(1, len(t_eval) - 1):
+    #         dt = t_eval[i + 1] - t_eval[i]
+    #         external_variables = {"X-averaged cell temperature": T_av}
+    #         T_av += 1
+    #         sim.step(dt, external_variables=external_variables)
+
+    #     sim.plot(
+    #         [
+    #             "Negative particle surface concentration",
+    #             "Electrolyte concentration",
+    #             "Positive particle surface concentration",
+    #             "Cell temperature [K]",
+    #             "Terminal voltage [V]",
+    #             "Total heating [W.m-3]",
+    #         ]
+    #     )
+
+    # def test_external_temperature(self):
+
+    #     model_options = {
+    #         "thermal": "x-full",
+    #         "external submodels": ["thermal"],
+    #     }
+
+    #     model = pybamm.lithium_ion.SPM(model_options)
+
+    #     # model.convert_to_format = False
+
+    #     neg_pts = 5
+    #     sep_pts = 3
+    #     pos_pts = 5
+    #     tot_pts = neg_pts + sep_pts + pos_pts
+
+    #     var_pts = {
+    #         pybamm.standard_spatial_vars.x_n: neg_pts,
+    #         pybamm.standard_spatial_vars.x_s: sep_pts,
+    #         pybamm.standard_spatial_vars.x_p: pos_pts,
+    #         pybamm.standard_spatial_vars.r_n: 5,
+    #         pybamm.standard_spatial_vars.r_p: 5,
+    #     }
+
+    #     sim = pybamm.Simulation(model, var_pts=var_pts)
+    #     sim.build()
+
+    #     t_eval = np.linspace(0, 0.17, 100)
+    #     x = np.linspace(0, 1, tot_pts)
+
+    #     for i in np.arange(1, len(t_eval) - 1):
+    #         dt = t_eval[i + 1] - t_eval[i]
+    #         T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
+    #             :, np.newaxis
+    #         ]
+    #         external_variables = {"Cell temperature": T}
+    #         sim.step(dt, external_variables=external_variables)
+
+    #     sim.plot(
+    #         [
+    #             "Negative particle surface concentration",
+    #             "Electrolyte concentration",
+    #             "Positive particle surface concentration",
+    #             "Cell temperature [K]",
+    #             "Terminal voltage [V]",
+    #             "Total heating [W.m-3]",
+    #         ]
+    #     )
+
+    def test_dae_external_temperature(self):
 
         model_options = {
             "thermal": "x-full",
             "external submodels": ["thermal"],
         }
 
-        model = pybamm.lithium_ion.SPM(model_options)
-
-        # model.convert_to_format = False
+        model = pybamm.lithium_ion.DFN(model_options)
 
         neg_pts = 5
         sep_pts = 3

From f9370adfb5a1b22f8ff2b83f9fb4b388480555df Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Mon, 11 Nov 2019 17:47:42 +0000
Subject: [PATCH 13/31] #692 external works for IDAKLU

---
 pybamm/solvers/base_solver.py                    | 16 ++++++++++++----
 pybamm/solvers/dae_solver.py                     |  2 ++
 pybamm/solvers/idaklu_solver.py                  |  2 +-
 pybamm/solvers/ode_solver.py                     |  8 --------
 .../test_use_external_submodel.py                |  7 +++++--
 5 files changed, 20 insertions(+), 15 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index d56e8b1a62..fff622829a 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -164,14 +164,22 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
         t_eval = np.linspace(self.t, self.t + dt, npts)
         solution, solve_time, termination = self.compute_solution(model, t_eval)
 
+        # Set self.t and self.y0 to their values at the final step
+        self.t = solution.t[-1]
+        self.y0 = solution.y[:, -1]
+
+        # add the external points onto the solution
+        full_y = np.zeros((model.y_length, solution.y.shape[1]))
+        for i in np.arange(solution.y.shape[1]):
+            sol_y = solution.y[:, i]
+            sol_y = sol_y[:, np.newaxis]
+            full_y[:, i] = self.add_external(sol_y)[:, 0]
+        solution.y = full_y
+
         # Assign times
         solution.solve_time = solve_time
         solution.set_up_time = set_up_time
 
-        # Set self.t and self.y0 to their values at the final step
-        self.t = solution.t[-1]
-        self.y0 = solution.y[: model.external_start, -1]
-
         pybamm.logger.debug("Finish stepping {} ({})".format(model.name, termination))
         if set_up_time:
             pybamm.logger.debug(
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 5268b7f76e..a79c6d83a6 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -89,6 +89,7 @@ def compute_solution(self, model, t_eval):
             mass_matrix=model.mass_matrix.entries,
             jacobian=self.jacobian,
         )
+
         solve_time = timer.time() - solve_start_time
 
         # Identify the event that caused termination
@@ -240,6 +241,7 @@ def set_up_casadi(self, model):
         # Convert model attributes to casadi
         t_casadi = casadi.SX.sym("t")
         y0 = model.concatenated_initial_conditions
+
         y_diff = casadi.SX.sym("y_diff", len(model.concatenated_rhs.evaluate(0, y0)))
         y_alg = casadi.SX.sym(
             "y_alg", len(model.concatenated_algebraic.evaluate(0, y0))
diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py
index 716e20b70c..b32efe1565 100644
--- a/pybamm/solvers/idaklu_solver.py
+++ b/pybamm/solvers/idaklu_solver.py
@@ -247,7 +247,7 @@ def rootfn(t, y):
             elif sol.flag == 2:
                 termination = "event"
             return pybamm.Solution(
-                sol.t, np.transpose(y_out), t[-1], np.transpose(y_out[-1]), termination
+                sol.t, np.transpose(y_out), t[-1], np.transpose(y_out[-1])[:, np.newaxis], termination
             )
         else:
             raise pybamm.SolverError(sol.message)
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index c777a0091b..1cd00681ca 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -46,14 +46,6 @@ def compute_solution(self, model, t_eval):
             jacobian=self.jacobian,
         )
 
-        # add the external points onto the solution vector
-        full_y = np.zeros((model.y_length, solution.y.shape[1]))
-        for i in np.arange(solution.y.shape[1]):
-            sol_y = solution.y[:, i]
-            sol_y = sol_y[:, np.newaxis]
-            full_y[:, i] = self.add_external(sol_y)[:, 0]
-        solution.y = full_y
-
         solve_time = timer.time() - solve_start_time
 
         # Identify the event that caused termination
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
index 5765404baf..1c3173c530 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
@@ -107,10 +107,11 @@ def test_dae_external_temperature(self):
             pybamm.standard_spatial_vars.r_p: 5,
         }
 
-        sim = pybamm.Simulation(model, var_pts=var_pts)
+        solver = pybamm.IDAKLUSolver()
+        sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver)
         sim.build()
 
-        t_eval = np.linspace(0, 0.17, 100)
+        t_eval = np.linspace(0, 0.17, 20)
         x = np.linspace(0, 1, tot_pts)
 
         for i in np.arange(1, len(t_eval) - 1):
@@ -118,8 +119,10 @@ def test_dae_external_temperature(self):
             T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
                 :, np.newaxis
             ]
+            # T = np.ones((tot_pts, 1))
             external_variables = {"Cell temperature": T}
             sim.step(dt, external_variables=external_variables)
+            print(i)
 
         sim.plot(
             [

From d64511a51552ccdb71681ec7c1c2603ca687484c Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Mon, 11 Nov 2019 17:53:03 +0000
Subject: [PATCH 14/31] #692 moved integration tests

---
 .../test_external/__init__.py                 |   0
 .../test_external_current_collector.py        |  15 ++
 .../test_external_temperature.py              | 145 ++++++++++++++++++
 .../test_use_external_submodel.py             | 145 ------------------
 4 files changed, 160 insertions(+), 145 deletions(-)
 create mode 100644 tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/__init__.py
 create mode 100644 tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
 create mode 100644 tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
 delete mode 100644 tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py

diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/__init__.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
new file mode 100644
index 0000000000..73e347b470
--- /dev/null
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
@@ -0,0 +1,15 @@
+#
+# Tests for external submodels
+#
+import pybamm
+import unittest
+import numpy as np
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
new file mode 100644
index 0000000000..2d56e1d42d
--- /dev/null
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
@@ -0,0 +1,145 @@
+#
+# Tests for external submodels
+#
+import pybamm
+import unittest
+import numpy as np
+
+
+class TestExternalThermalModels(unittest.TestCase):
+    def test_external_lumped_temperature(self):
+        model_options = {
+            "thermal": "x-lumped",
+            "external submodels": ["thermal"],
+        }
+        model = pybamm.lithium_ion.SPMe(model_options)
+        sim = pybamm.Simulation(model)
+
+        t_eval = np.linspace(0, 0.17, 100)
+
+        T_av = 0
+
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+            external_variables = {"X-averaged cell temperature": T_av}
+            T_av += 1
+            sim.step(dt, external_variables=external_variables)
+
+        sim.plot(
+            [
+                "Negative particle surface concentration",
+                "Electrolyte concentration",
+                "Positive particle surface concentration",
+                "Cell temperature [K]",
+                "Terminal voltage [V]",
+                "Total heating [W.m-3]",
+            ]
+        )
+
+    def test_external_temperature(self):
+
+        model_options = {
+            "thermal": "x-full",
+            "external submodels": ["thermal"],
+        }
+
+        model = pybamm.lithium_ion.SPM(model_options)
+
+        # model.convert_to_format = False
+
+        neg_pts = 5
+        sep_pts = 3
+        pos_pts = 5
+        tot_pts = neg_pts + sep_pts + pos_pts
+
+        var_pts = {
+            pybamm.standard_spatial_vars.x_n: neg_pts,
+            pybamm.standard_spatial_vars.x_s: sep_pts,
+            pybamm.standard_spatial_vars.x_p: pos_pts,
+            pybamm.standard_spatial_vars.r_n: 5,
+            pybamm.standard_spatial_vars.r_p: 5,
+        }
+
+        sim = pybamm.Simulation(model, var_pts=var_pts)
+        sim.build()
+
+        t_eval = np.linspace(0, 0.17, 100)
+        x = np.linspace(0, 1, tot_pts)
+
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+            T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
+                :, np.newaxis
+            ]
+            external_variables = {"Cell temperature": T}
+            sim.step(dt, external_variables=external_variables)
+
+        sim.plot(
+            [
+                "Negative particle surface concentration",
+                "Electrolyte concentration",
+                "Positive particle surface concentration",
+                "Cell temperature [K]",
+                "Terminal voltage [V]",
+                "Total heating [W.m-3]",
+            ]
+        )
+
+    def test_dae_external_temperature(self):
+
+        model_options = {
+            "thermal": "x-full",
+            "external submodels": ["thermal"],
+        }
+
+        model = pybamm.lithium_ion.DFN(model_options)
+
+        neg_pts = 5
+        sep_pts = 3
+        pos_pts = 5
+        tot_pts = neg_pts + sep_pts + pos_pts
+
+        var_pts = {
+            pybamm.standard_spatial_vars.x_n: neg_pts,
+            pybamm.standard_spatial_vars.x_s: sep_pts,
+            pybamm.standard_spatial_vars.x_p: pos_pts,
+            pybamm.standard_spatial_vars.r_n: 5,
+            pybamm.standard_spatial_vars.r_p: 5,
+        }
+
+        solver = pybamm.IDAKLUSolver()
+        sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver)
+        sim.build()
+
+        t_eval = np.linspace(0, 0.17, 20)
+        x = np.linspace(0, 1, tot_pts)
+
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+            T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
+                :, np.newaxis
+            ]
+            # T = np.ones((tot_pts, 1))
+            external_variables = {"Cell temperature": T}
+            sim.step(dt, external_variables=external_variables)
+            print(i)
+
+        sim.plot(
+            [
+                "Negative particle surface concentration",
+                "Electrolyte concentration",
+                "Positive particle surface concentration",
+                "Cell temperature [K]",
+                "Terminal voltage [V]",
+                "Total heating [W.m-3]",
+            ]
+        )
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
deleted file mode 100644
index 1c3173c530..0000000000
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_use_external_submodel.py
+++ /dev/null
@@ -1,145 +0,0 @@
-#
-# Tests for external submodels
-#
-import pybamm
-import unittest
-import numpy as np
-
-
-class TestExternalSubmodel(unittest.TestCase):
-    # def test_external_lumped_temperature(self):
-    #     model_options = {
-    #         "thermal": "x-lumped",
-    #         "external submodels": ["thermal"],
-    #     }
-    #     model = pybamm.lithium_ion.SPMe(model_options)
-    #     sim = pybamm.Simulation(model)
-
-    #     t_eval = np.linspace(0, 0.17, 100)
-
-    #     T_av = 0
-
-    #     for i in np.arange(1, len(t_eval) - 1):
-    #         dt = t_eval[i + 1] - t_eval[i]
-    #         external_variables = {"X-averaged cell temperature": T_av}
-    #         T_av += 1
-    #         sim.step(dt, external_variables=external_variables)
-
-    #     sim.plot(
-    #         [
-    #             "Negative particle surface concentration",
-    #             "Electrolyte concentration",
-    #             "Positive particle surface concentration",
-    #             "Cell temperature [K]",
-    #             "Terminal voltage [V]",
-    #             "Total heating [W.m-3]",
-    #         ]
-    #     )
-
-    # def test_external_temperature(self):
-
-    #     model_options = {
-    #         "thermal": "x-full",
-    #         "external submodels": ["thermal"],
-    #     }
-
-    #     model = pybamm.lithium_ion.SPM(model_options)
-
-    #     # model.convert_to_format = False
-
-    #     neg_pts = 5
-    #     sep_pts = 3
-    #     pos_pts = 5
-    #     tot_pts = neg_pts + sep_pts + pos_pts
-
-    #     var_pts = {
-    #         pybamm.standard_spatial_vars.x_n: neg_pts,
-    #         pybamm.standard_spatial_vars.x_s: sep_pts,
-    #         pybamm.standard_spatial_vars.x_p: pos_pts,
-    #         pybamm.standard_spatial_vars.r_n: 5,
-    #         pybamm.standard_spatial_vars.r_p: 5,
-    #     }
-
-    #     sim = pybamm.Simulation(model, var_pts=var_pts)
-    #     sim.build()
-
-    #     t_eval = np.linspace(0, 0.17, 100)
-    #     x = np.linspace(0, 1, tot_pts)
-
-    #     for i in np.arange(1, len(t_eval) - 1):
-    #         dt = t_eval[i + 1] - t_eval[i]
-    #         T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
-    #             :, np.newaxis
-    #         ]
-    #         external_variables = {"Cell temperature": T}
-    #         sim.step(dt, external_variables=external_variables)
-
-    #     sim.plot(
-    #         [
-    #             "Negative particle surface concentration",
-    #             "Electrolyte concentration",
-    #             "Positive particle surface concentration",
-    #             "Cell temperature [K]",
-    #             "Terminal voltage [V]",
-    #             "Total heating [W.m-3]",
-    #         ]
-    #     )
-
-    def test_dae_external_temperature(self):
-
-        model_options = {
-            "thermal": "x-full",
-            "external submodels": ["thermal"],
-        }
-
-        model = pybamm.lithium_ion.DFN(model_options)
-
-        neg_pts = 5
-        sep_pts = 3
-        pos_pts = 5
-        tot_pts = neg_pts + sep_pts + pos_pts
-
-        var_pts = {
-            pybamm.standard_spatial_vars.x_n: neg_pts,
-            pybamm.standard_spatial_vars.x_s: sep_pts,
-            pybamm.standard_spatial_vars.x_p: pos_pts,
-            pybamm.standard_spatial_vars.r_n: 5,
-            pybamm.standard_spatial_vars.r_p: 5,
-        }
-
-        solver = pybamm.IDAKLUSolver()
-        sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver)
-        sim.build()
-
-        t_eval = np.linspace(0, 0.17, 20)
-        x = np.linspace(0, 1, tot_pts)
-
-        for i in np.arange(1, len(t_eval) - 1):
-            dt = t_eval[i + 1] - t_eval[i]
-            T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
-                :, np.newaxis
-            ]
-            # T = np.ones((tot_pts, 1))
-            external_variables = {"Cell temperature": T}
-            sim.step(dt, external_variables=external_variables)
-            print(i)
-
-        sim.plot(
-            [
-                "Negative particle surface concentration",
-                "Electrolyte concentration",
-                "Positive particle surface concentration",
-                "Cell temperature [K]",
-                "Terminal voltage [V]",
-                "Total heating [W.m-3]",
-            ]
-        )
-
-
-if __name__ == "__main__":
-    print("Add -v for more debug output")
-    import sys
-
-    if "-v" in sys.argv:
-        debug = True
-    unittest.main()

From 6b084be3f46da3f23e4327b1f75bd6cfc3edc505 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Tue, 12 Nov 2019 10:43:12 +0000
Subject: [PATCH 15/31] #692 begining to add cc tests

---
 .../test_external_current_collector.py          | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
index 73e347b470..2d8d991029 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
@@ -6,6 +6,23 @@
 import numpy as np
 
 
+class TestExternalCC(unittest.TestCase):
+
+    model_options = {
+        "current collector": "potential pair",
+        "dimensionality": 1,
+        "external submodels": ["current collector"],
+    }
+    model = pybamm.lithium_ion.SPM(model_options)
+    var_pts = model.default_var_pts
+    var_pts.update({pybamm.standard_spatial_vars.y: 5})
+    sim = pybamm.Simulation(model)
+
+    # provide phi_s_n and i_cc
+
+    # return phi_s_p
+
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys

From 1b9d527fa0d723a3af148929982eaac222c68a8c Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Tue, 12 Nov 2019 14:57:22 +0000
Subject: [PATCH 16/31] #692 external works for current collector models

---
 pybamm/simulation.py                          | 29 ++++++++++
 .../test_external_current_collector.py        | 53 ++++++++++++++-----
 .../test_external_temperature.py              |  2 -
 3 files changed, 68 insertions(+), 16 deletions(-)

diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 96be32b733..a5c7704265 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -165,6 +165,35 @@ def _update_solution(self, solution):
         )
         self._solution.y_event = solution.y_event
 
+    def get_variable_array(self, *variables):
+        """
+        A helper function to easily obtain a dictionary of arrays of values
+        for a list of variables at the latest timestep.
+
+        Parameters
+        ----------
+        variable: str
+            The name of the variable/variables you wish to obtain the arrays for.
+
+        Returns
+        -------
+        variable_arrays: dict
+            A dictionary of the variable names and their corresponding
+            arrays.
+        """
+
+        variable_arrays = [
+            self.built_model.variables[var].evaluate(
+                self.solution.t[-1], self.solution.y[:, -1]
+            )
+            for var in variables
+        ]
+
+        if len(variable_arrays) == 1:
+            return variable_arrays[0]
+        else:
+            return tuple(variable_arrays)
+
     def plot(self, quick_plot_vars=None):
         """
         A method to quickly plot the outputs of the simulation.
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
index 2d8d991029..57cc570061 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
@@ -7,20 +7,45 @@
 
 
 class TestExternalCC(unittest.TestCase):
-
-    model_options = {
-        "current collector": "potential pair",
-        "dimensionality": 1,
-        "external submodels": ["current collector"],
-    }
-    model = pybamm.lithium_ion.SPM(model_options)
-    var_pts = model.default_var_pts
-    var_pts.update({pybamm.standard_spatial_vars.y: 5})
-    sim = pybamm.Simulation(model)
-
-    # provide phi_s_n and i_cc
-
-    # return phi_s_p
+    def test_2p1d(self):
+        model_options = {
+            "current collector": "potential pair",
+            "dimensionality": 2,
+            "external submodels": ["current collector"],
+        }
+        model = pybamm.lithium_ion.DFN(model_options)
+        yz_pts = 3
+        var_pts = {
+            pybamm.standard_spatial_vars.x_n: 4,
+            pybamm.standard_spatial_vars.x_s: 4,
+            pybamm.standard_spatial_vars.x_p: 4,
+            pybamm.standard_spatial_vars.r_n: 4,
+            pybamm.standard_spatial_vars.r_p: 4,
+            pybamm.standard_spatial_vars.y: yz_pts,
+            pybamm.standard_spatial_vars.z: yz_pts,
+        }
+        sim = pybamm.Simulation(model, var_pts=var_pts)
+
+        t_eval = np.linspace(0, 0.08, 3)
+
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+            print(t_eval[i])
+
+            # provide phi_s_n and i_cc
+            phi_s_n = np.zeros((yz_pts**2, 1))
+            i_boundary_cc = np.ones((yz_pts**2, 1))
+            external_variables = {
+                "Negative current collector potential": phi_s_n,
+                "Current collector current density": i_boundary_cc,
+            }
+
+            sim.step(dt, external_variables=external_variables)
+
+            # obtain phi_s_n from the pybamm solution at the current time
+            phi_s_p = sim.get_variable_array("Positive current collector potential")
+
+        self.assertTrue(phi_s_p.shape, (yz_pts**2, 1))
 
 
 if __name__ == "__main__":
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
index 2d56e1d42d..ba11935215 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
@@ -45,8 +45,6 @@ def test_external_temperature(self):
 
         model = pybamm.lithium_ion.SPM(model_options)
 
-        # model.convert_to_format = False
-
         neg_pts = 5
         sep_pts = 3
         pos_pts = 5

From 607df1f039c7b13ec92052977f0afdc1a455de89 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Tue, 12 Nov 2019 15:06:01 +0000
Subject: [PATCH 17/31] #692 added unit test for external variables in
 discretization

---
 tests/unit/test_discretisations/test_discretisation.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 0a61b05078..957ec214c6 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -66,6 +66,10 @@ def test_add_internal_boundary_conditions(self):
         for child in c_e.children:
             self.assertTrue(child.id in disc.bcs.keys())
 
+    def test_adding_external_variables(self):
+        print("hello")
+
+
     def test_discretise_slicing(self):
         # create discretisation
         mesh = get_mesh_for_testing()

From 270d207f41157e04087b7dd6560b41d3c78451cf Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Tue, 12 Nov 2019 16:20:55 +0000
Subject: [PATCH 18/31] #692 can calculate jacobian using casadi on external
 variables

---
 pybamm/discretisations/discretisation.py      |   1 +
 pybamm/models/base_model.py                   |   7 +-
 pybamm/solvers/dae_solver.py                  |  34 ++--
 pybamm/solvers/ode_solver.py                  |  16 +-
 .../test_external_temperature.py              | 153 +++++++++---------
 .../test_discretisation.py                    |  15 +-
 6 files changed, 134 insertions(+), 92 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index c57979bece..8e4f9d43fc 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -52,6 +52,7 @@ def __init__(self, mesh=None, spatial_methods=None):
         self.bcs = {}
         self.y_slices = {}
         self._discretised_symbols = {}
+        self.external_variables = []
 
     @property
     def mesh(self):
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 4fab9f7752..c6a1818ff5 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -99,6 +99,7 @@ def __init__(self, name="Unnamed model"):
         self._mass_matrix = None
         self._jacobian = None
         self._jacobian_algebraic = None
+        self.external_variables = []
 
         # Default behaviour is to use the jacobian and simplify
         self.use_jacobian = True
@@ -472,19 +473,21 @@ def check_variables(self):
                 {x.id: x for x in eqn.pre_order() if isinstance(x, pybamm.Variable)}
             )
         var_ids_in_keys = set()
+
         for var in {**self.rhs, **self.algebraic}.keys():
             if isinstance(var, pybamm.Variable):
                 var_ids_in_keys.add(var.id)
             # Key can be a concatenation
             elif isinstance(var, pybamm.Concatenation):
                 var_ids_in_keys.update([child.id for child in var.children])
+
         for var_id, var in all_vars.items():
             if var_id not in var_ids_in_keys:
                 raise pybamm.ModelError(
                     """
                     No key set for variable '{}'. Make sure it is included in either
-                    model.rhs or model.algebraic in an unmodified form (e.g. not
-                    Broadcasted)
+                    model.rhs, model.algebraic, or model.external_variables in an
+                    unmodified form (e.g. not Broadcasted)
                     """.format(
                         var
                     )
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index aac2ebd589..ad106777e3 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -242,49 +242,57 @@ def set_up_casadi(self, model):
         # Convert model attributes to casadi
         t_casadi = casadi.MX.sym("t")
         y0 = model.concatenated_initial_conditions
+        y0 = self.add_external(y0)
+
         y_diff = casadi.MX.sym("y_diff", len(model.concatenated_rhs.evaluate(0, y0)))
         y_alg = casadi.MX.sym(
             "y_alg", len(model.concatenated_algebraic.evaluate(0, y0))
         )
+        y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
         y_casadi = casadi.vertcat(y_diff, y_alg)
+        y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
         pybamm.logger.info("Converting RHS to CasADi")
-        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi)
+        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi_w_ext)
         pybamm.logger.info("Converting algebraic to CasADi")
         concatenated_algebraic = model.concatenated_algebraic.to_casadi(
-            t_casadi, y_casadi
+            t_casadi, y_casadi_w_ext
         )
         all_states = casadi.vertcat(concatenated_rhs, concatenated_algebraic)
         pybamm.logger.info("Converting events to CasADi")
         casadi_events = {
-            name: event.to_casadi(t_casadi, y_casadi)
+            name: event.to_casadi(t_casadi, y_casadi_w_ext)
             for name, event in model.events.items()
         }
 
         # Create functions to evaluate rhs and algebraic
         concatenated_rhs_fn = casadi.Function(
-            "rhs", [t_casadi, y_casadi], [concatenated_rhs]
+            "rhs", [t_casadi, y_casadi_w_ext], [concatenated_rhs]
         )
         concatenated_algebraic_fn = casadi.Function(
-            "algebraic", [t_casadi, y_casadi], [concatenated_algebraic]
+            "algebraic", [t_casadi, y_casadi_w_ext], [concatenated_algebraic]
         )
-        all_states_fn = casadi.Function("all", [t_casadi, y_casadi], [all_states])
+        all_states_fn = casadi.Function("all", [t_casadi, y_casadi_w_ext], [all_states])
 
         if model.use_jacobian:
 
             pybamm.logger.info("Calculating jacobian")
             casadi_jac = casadi.jacobian(all_states, y_casadi)
             casadi_jac_fn = casadi.Function(
-                "jacobian", [t_casadi, y_casadi], [casadi_jac]
+                "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac]
             )
             casadi_jac_alg = casadi.jacobian(concatenated_algebraic, y_casadi)
             casadi_jac_alg_fn = casadi.Function(
-                "jacobian", [t_casadi, y_casadi], [casadi_jac_alg]
+                "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac_alg]
             )
 
             def jacobian(t, y):
+                y = y[:, np.newaxis]
+                y = self.add_external(y)
                 return casadi_jac_fn(t, y)
 
             def jacobian_alg(t, y):
+                y = y[:, np.newaxis]
+                y = self.add_external(y)
                 return casadi_jac_alg_fn(t, y)
 
         else:
@@ -293,9 +301,13 @@ def jacobian_alg(t, y):
 
         # Calculate consistent initial conditions for the algebraic equations
         def rhs(t, y):
+            y = y[:, np.newaxis]
+            y = self.add_external(y)
             return concatenated_rhs_fn(t, y).full()[:, 0]
 
         def algebraic(t, y):
+            y = y[:, np.newaxis]
+            y = self.add_external(y)
             return concatenated_algebraic_fn(t, y).full()[:, 0]
 
         if len(model.algebraic) > 0:
@@ -316,14 +328,18 @@ def residuals(t, y, ydot):
             pybamm.logger.debug(
                 "Evaluating residuals for {} at t={}".format(model.name, t)
             )
+            y = y[:, np.newaxis]
+            y = self.add_external(y)
             states_eval = all_states_fn(t, y).full()[:, 0]
             return states_eval - model.mass_matrix.entries @ ydot
 
         # Create event-dependent function to evaluate events
         def event_fun(event):
-            casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi], [event])
+            casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi_w_ext], [event])
 
             def eval_event(t, y):
+                y = y[:, np.newaxis]
+                y = self.add_external(y)
                 return casadi_event_fn(t, y)
 
             return eval_event
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 9b04d68356..6aaab8d092 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -184,17 +184,21 @@ def set_up_casadi(self, model):
 
         t_casadi = casadi.MX.sym("t")
         y_casadi = casadi.MX.sym("y", len(y0))
+
+        y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
+        y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
+
         pybamm.logger.info("Converting RHS to CasADi")
-        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi)
+        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi_w_ext)
         pybamm.logger.info("Converting events to CasADi")
         casadi_events = {
-            name: event.to_casadi(t_casadi, y_casadi)
+            name: event.to_casadi(t_casadi, y_casadi_w_ext)
             for name, event in model.events.items()
         }
 
         # Create function to evaluate rhs
         concatenated_rhs_fn = casadi.Function(
-            "rhs", [t_casadi, y_casadi], [concatenated_rhs]
+            "rhs", [t_casadi, y_casadi_w_ext], [concatenated_rhs]
         )
 
         def dydt(t, y):
@@ -206,7 +210,9 @@ def dydt(t, y):
 
         # Create event-dependent function to evaluate events
         def event_fun(event):
-            casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi], [event])
+            casadi_event_fn = casadi.Function(
+                "event", [t_casadi, y_casadi_w_ext], [event]
+            )
 
             def eval_event(t, y):
                 y = y[:, np.newaxis]
@@ -222,7 +228,7 @@ def eval_event(t, y):
             pybamm.logger.info("Calculating jacobian")
             casadi_jac = casadi.jacobian(concatenated_rhs, y_casadi)
             casadi_jac_fn = casadi.Function(
-                "jacobian", [t_casadi, y_casadi], [casadi_jac]
+                "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac]
             )
 
             def jacobian(t, y):
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
index ba11935215..b136b0bcf3 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
@@ -7,81 +7,84 @@
 
 
 class TestExternalThermalModels(unittest.TestCase):
-    def test_external_lumped_temperature(self):
-        model_options = {
-            "thermal": "x-lumped",
-            "external submodels": ["thermal"],
-        }
-        model = pybamm.lithium_ion.SPMe(model_options)
-        sim = pybamm.Simulation(model)
-
-        t_eval = np.linspace(0, 0.17, 100)
-
-        T_av = 0
-
-        for i in np.arange(1, len(t_eval) - 1):
-            dt = t_eval[i + 1] - t_eval[i]
-            external_variables = {"X-averaged cell temperature": T_av}
-            T_av += 1
-            sim.step(dt, external_variables=external_variables)
-
-        sim.plot(
-            [
-                "Negative particle surface concentration",
-                "Electrolyte concentration",
-                "Positive particle surface concentration",
-                "Cell temperature [K]",
-                "Terminal voltage [V]",
-                "Total heating [W.m-3]",
-            ]
-        )
-
-    def test_external_temperature(self):
-
-        model_options = {
-            "thermal": "x-full",
-            "external submodels": ["thermal"],
-        }
-
-        model = pybamm.lithium_ion.SPM(model_options)
-
-        neg_pts = 5
-        sep_pts = 3
-        pos_pts = 5
-        tot_pts = neg_pts + sep_pts + pos_pts
-
-        var_pts = {
-            pybamm.standard_spatial_vars.x_n: neg_pts,
-            pybamm.standard_spatial_vars.x_s: sep_pts,
-            pybamm.standard_spatial_vars.x_p: pos_pts,
-            pybamm.standard_spatial_vars.r_n: 5,
-            pybamm.standard_spatial_vars.r_p: 5,
-        }
-
-        sim = pybamm.Simulation(model, var_pts=var_pts)
-        sim.build()
-
-        t_eval = np.linspace(0, 0.17, 100)
-        x = np.linspace(0, 1, tot_pts)
-
-        for i in np.arange(1, len(t_eval) - 1):
-            dt = t_eval[i + 1] - t_eval[i]
-            T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
-                :, np.newaxis
-            ]
-            external_variables = {"Cell temperature": T}
-            sim.step(dt, external_variables=external_variables)
-
-        sim.plot(
-            [
-                "Negative particle surface concentration",
-                "Electrolyte concentration",
-                "Positive particle surface concentration",
-                "Cell temperature [K]",
-                "Terminal voltage [V]",
-                "Total heating [W.m-3]",
-            ]
-        )
+    # def test_external_lumped_temperature(self):
+    #     model_options = {
+    #         "thermal": "x-lumped",
+    #         "external submodels": ["thermal"],
+    #     }
+    #     model = pybamm.lithium_ion.SPMe(model_options)
+    #     sim = pybamm.Simulation(model)
+
+    #     t_eval = np.linspace(0, 0.17, 100)
+
+    #     T_av = 0
+
+    #     for i in np.arange(1, len(t_eval) - 1):
+    #         dt = t_eval[i + 1] - t_eval[i]
+    #         external_variables = {"X-averaged cell temperature": T_av}
+    #         T_av += 1
+    #         sim.step(dt, external_variables=external_variables)
+
+    #     sim.solver = pybamm.IDAKLUSolver()
+
+    #     sim.plot(
+    #         [
+    #             "Negative particle surface concentration",
+    #             "Electrolyte concentration",
+    #             "Positive particle surface concentration",
+    #             "Cell temperature [K]",
+    #             "Terminal voltage [V]",
+    #             "Total heating [W.m-3]",
+    #         ]
+    #     )
+
+    # def test_external_temperature(self):
+
+    #     model_options = {
+    #         "thermal": "x-full",
+    #         "external submodels": ["thermal"],
+    #     }
+
+    #     model = pybamm.lithium_ion.SPM(model_options)
+
+    #     neg_pts = 5
+    #     sep_pts = 3
+    #     pos_pts = 5
+    #     tot_pts = neg_pts + sep_pts + pos_pts
+
+    #     var_pts = {
+    #         pybamm.standard_spatial_vars.x_n: neg_pts,
+    #         pybamm.standard_spatial_vars.x_s: sep_pts,
+    #         pybamm.standard_spatial_vars.x_p: pos_pts,
+    #         pybamm.standard_spatial_vars.r_n: 5,
+    #         pybamm.standard_spatial_vars.r_p: 5,
+    #     }
+
+    #     sim = pybamm.Simulation(model, var_pts=var_pts)
+
+    #     sim.solver = pybamm.IDAKLUSolver()
+
+    #     t_eval = np.linspace(0, 0.17, 100)
+    #     x = np.linspace(0, 1, tot_pts)
+
+    #     for i in np.arange(1, len(t_eval) - 1):
+    #         dt = t_eval[i + 1] - t_eval[i]
+    #         T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
+    #             :, np.newaxis
+    #         ]
+    #         external_variables = {"Cell temperature": T}
+    #         sim.step(dt, external_variables=external_variables)
+
+    #     sim.plot(
+    #         [
+    #             "Negative particle surface concentration",
+    #             "Electrolyte concentration",
+    #             "Positive particle surface concentration",
+    #             "Cell temperature [K]",
+    #             "Terminal voltage [V]",
+    #             "Total heating [W.m-3]",
+    #         ]
+    #     )
 
     def test_dae_external_temperature(self):
 
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 957ec214c6..8877f10b0c 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -67,7 +67,20 @@ def test_add_internal_boundary_conditions(self):
             self.assertTrue(child.id in disc.bcs.keys())
 
     def test_adding_external_variables(self):
-        print("hello")
+        model = pybamm.BaseModel()
+        a = pybamm.Variable("a")
+        b = pybamm.Variable("b")
+
+        model.rhs = {a: a * b}
+        model.initial_conditions = {a: 0}
+        model.external_variables = [b]
+        model.variables = {"a": a, "b": b, "c": a * b}
+
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        model
+
 
 
     def test_discretise_slicing(self):

From 832940a07603a2278708635531b75ab43a92087e Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Tue, 12 Nov 2019 17:31:53 +0000
Subject: [PATCH 19/31] #692 added tests for discretization

---
 pybamm/models/base_model.py                   |  8 +++-
 .../test_discretisation.py                    | 43 ++++++++++++++++++-
 2 files changed, 48 insertions(+), 3 deletions(-)

diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index c6a1818ff5..efec42e12f 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -474,7 +474,13 @@ def check_variables(self):
             )
         var_ids_in_keys = set()
 
-        for var in {**self.rhs, **self.algebraic}.keys():
+        model_and_external_variables = (
+            list(self.rhs.keys())
+            + list(self.algebraic.keys())
+            + self.external_variables
+        )
+
+        for var in model_and_external_variables:
             if isinstance(var, pybamm.Variable):
                 var_ids_in_keys.add(var.id)
             # Key can be a concatenation
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 8877f10b0c..39797750cc 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -66,7 +66,7 @@ def test_add_internal_boundary_conditions(self):
         for child in c_e.children:
             self.assertTrue(child.id in disc.bcs.keys())
 
-    def test_adding_external_variables(self):
+    def test_adding_0D_external_variable(self):
         model = pybamm.BaseModel()
         a = pybamm.Variable("a")
         b = pybamm.Variable("b")
@@ -79,9 +79,48 @@ def test_adding_external_variables(self):
         disc = pybamm.Discretisation()
         disc.process_model(model)
 
-        model
+        self.assertEqual(len(model.y_slices), 2)
+        self.assertEqual(model.y_slices[b.id][0], slice(1, 2, None))
 
+    def test_adding_1D_external_variable(self):
+        model = pybamm.BaseModel()
+
+        a = pybamm.Variable("a", domain=["test"])
+        b = pybamm.Variable("b", domain=["test"])
+
+        model.rhs = {a: a * b}
+        model.boundary_conditions = {
+            a: {"left": (0, "Dirichlet"), "right": (0, "Dirichlet")}
+        }
+        model.initial_conditions = {a: 0}
+        model.external_variables = [b]
+        model.variables = {"a": a, "b": b, "c": a * b, "grad b": pybamm.grad(b), "div grad b": pybamm.div(pybamm.grad(b))}
+
+        x = pybamm.SpatialVariable("x", domain="test", coord_sys="cartesian")
+        geometry = {
+            "test": {"primary": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}}
+        }
+
+        submesh_types = {"test": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh)}
+        var_pts = {x: 10}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        spatial_methods = {"test": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+
+        self.assertEqual(len(model.y_slices), 2)
+        self.assertEqual(model.y_slices[a.id][0], slice(0, 10, None))
+        self.assertEqual(model.y_slices[b.id][0], slice(10, 20, None))
+
+        # check that b is added to the boundary conditions
+        model.bcs[b.id]["left"]
+        model.bcs[b.id]["right"]
 
+        # check that grad and div(grad ) produce the correct shapes
+        self.assertEqual(model.variables["b"].shape_for_testing, (10, 1))
+        self.assertEqual(model.variables["grad b"].shape_for_testing, (11, 1))
+        self.assertEqual(model.variables["div grad b"].shape_for_testing, (10, 1))
 
     def test_discretise_slicing(self):
         # create discretisation

From d36d5cc46858921928de89deabf713bf8a1ca666 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Wed, 13 Nov 2019 10:00:35 +0000
Subject: [PATCH 20/31] 692 added tests for simulation

---
 pybamm/discretisations/discretisation.py      |   6 +-
 pybamm/solvers/ode_solver.py                  |   7 +-
 .../test_external_temperature.py              | 134 +++++++-----------
 tests/unit/test_simulation.py                 |  53 +++++++
 4 files changed, 117 insertions(+), 83 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 8e4f9d43fc..e1bb40c5c7 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -167,10 +167,12 @@ def process_model(self, model, inplace=True):
                 start_vals += [self.y_slices[var.id][0].start]
 
         model_disc.external_variables = model.external_variables
-        if start_vals:
-            model_disc.external_start = min(start_vals)
         model_disc.y_length = self.y_length
         model_disc.y_slices = self.y_slices
+        if start_vals:
+            model_disc.external_start = min(start_vals)
+        else:
+            model_disc.external_start = self.y_length
 
         pybamm.logger.info("Discretise initial conditions for {}".format(model.name))
         ics, concat_ics = self.process_initial_conditions(model)
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 6aaab8d092..42e460fbf4 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -185,8 +185,11 @@ def set_up_casadi(self, model):
         t_casadi = casadi.MX.sym("t")
         y_casadi = casadi.MX.sym("y", len(y0))
 
-        y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
-        y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
+        if self.y_pad is not None:
+            y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
+            y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
+        else:
+            y_casadi_w_ext = y_casadi
 
         pybamm.logger.info("Converting RHS to CasADi")
         concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi_w_ext)
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
index b136b0bcf3..00bca26527 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
@@ -7,84 +7,60 @@
 
 
 class TestExternalThermalModels(unittest.TestCase):
-    # def test_external_lumped_temperature(self):
-    #     model_options = {
-    #         "thermal": "x-lumped",
-    #         "external submodels": ["thermal"],
-    #     }
-    #     model = pybamm.lithium_ion.SPMe(model_options)
-    #     sim = pybamm.Simulation(model)
-
-    #     t_eval = np.linspace(0, 0.17, 100)
-
-    #     T_av = 0
-
-    #     for i in np.arange(1, len(t_eval) - 1):
-    #         dt = t_eval[i + 1] - t_eval[i]
-    #         external_variables = {"X-averaged cell temperature": T_av}
-    #         T_av += 1
-    #         sim.step(dt, external_variables=external_variables)
-
-    #     sim.solver = pybamm.IDAKLUSolver()
-
-    #     sim.plot(
-    #         [
-    #             "Negative particle surface concentration",
-    #             "Electrolyte concentration",
-    #             "Positive particle surface concentration",
-    #             "Cell temperature [K]",
-    #             "Terminal voltage [V]",
-    #             "Total heating [W.m-3]",
-    #         ]
-    #     )
-
-    # def test_external_temperature(self):
-
-    #     model_options = {
-    #         "thermal": "x-full",
-    #         "external submodels": ["thermal"],
-    #     }
-
-    #     model = pybamm.lithium_ion.SPM(model_options)
-
-    #     neg_pts = 5
-    #     sep_pts = 3
-    #     pos_pts = 5
-    #     tot_pts = neg_pts + sep_pts + pos_pts
-
-    #     var_pts = {
-    #         pybamm.standard_spatial_vars.x_n: neg_pts,
-    #         pybamm.standard_spatial_vars.x_s: sep_pts,
-    #         pybamm.standard_spatial_vars.x_p: pos_pts,
-    #         pybamm.standard_spatial_vars.r_n: 5,
-    #         pybamm.standard_spatial_vars.r_p: 5,
-    #     }
-
-    #     sim = pybamm.Simulation(model, var_pts=var_pts)
-
-    #     sim.solver = pybamm.IDAKLUSolver()
-
-    #     t_eval = np.linspace(0, 0.17, 100)
-    #     x = np.linspace(0, 1, tot_pts)
-
-    #     for i in np.arange(1, len(t_eval) - 1):
-    #         dt = t_eval[i + 1] - t_eval[i]
-    #         T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
-    #             :, np.newaxis
-    #         ]
-    #         external_variables = {"Cell temperature": T}
-    #         sim.step(dt, external_variables=external_variables)
-
-    #     sim.plot(
-    #         [
-    #             "Negative particle surface concentration",
-    #             "Electrolyte concentration",
-    #             "Positive particle surface concentration",
-    #             "Cell temperature [K]",
-    #             "Terminal voltage [V]",
-    #             "Total heating [W.m-3]",
-    #         ]
-    #     )
+    def test_external_lumped_temperature(self):
+        model_options = {
+            "thermal": "x-lumped",
+            "external submodels": ["thermal"],
+        }
+        model = pybamm.lithium_ion.SPMe(model_options)
+        sim = pybamm.Simulation(model)
+
+        t_eval = np.linspace(0, 0.01, 3)
+
+        T_av = 0
+
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+            external_variables = {"X-averaged cell temperature": T_av}
+            T_av += 1
+            sim.step(dt, external_variables=external_variables)
+
+    def test_external_temperature(self):
+
+        model_options = {
+            "thermal": "x-full",
+            "external submodels": ["thermal"],
+        }
+
+        model = pybamm.lithium_ion.SPM(model_options)
+
+        neg_pts = 5
+        sep_pts = 3
+        pos_pts = 5
+        tot_pts = neg_pts + sep_pts + pos_pts
+
+        var_pts = {
+            pybamm.standard_spatial_vars.x_n: neg_pts,
+            pybamm.standard_spatial_vars.x_s: sep_pts,
+            pybamm.standard_spatial_vars.x_p: pos_pts,
+            pybamm.standard_spatial_vars.r_n: 5,
+            pybamm.standard_spatial_vars.r_p: 5,
+        }
+
+        sim = pybamm.Simulation(model, var_pts=var_pts)
+
+        sim.solver = pybamm.IDAKLUSolver()
+
+        t_eval = np.linspace(0, 0.01, 3)
+        x = np.linspace(0, 1, tot_pts)
+
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+            T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
+                :, np.newaxis
+            ]
+            external_variables = {"Cell temperature": T}
+            sim.step(dt, external_variables=external_variables)
 
     def test_dae_external_temperature(self):
 
@@ -112,7 +88,7 @@ def test_dae_external_temperature(self):
         sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver)
         sim.build()
 
-        t_eval = np.linspace(0, 0.17, 20)
+        t_eval = np.linspace(0, 0.01, 3)
         x = np.linspace(0, 1, tot_pts)
 
         for i in np.arange(1, len(t_eval) - 1):
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index 8dd3aa1b36..594a3f093d 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -1,4 +1,5 @@
 import pybamm
+import numpy as np
 import unittest
 
 
@@ -129,6 +130,58 @@ def test_specs(self):
         sim.specs(spatial_methods=spatial_methods)
         sim.build()
 
+    def test_get_variable_array(self):
+
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
+        sim.solve()
+
+        phi_s_n = sim.get_variable_array("Negative electrode potential")
+
+        self.assertIsInstance(phi_s_n, np.ndarray)
+
+        c_s_n_surf, c_e = sim.get_variable_array(
+            "Negative particle surface concentration", "Electrolyte concentration"
+        )
+
+        self.assertIsInstance(c_s_n_surf, np.ndarray)
+        self.assertIsInstance(c_e, np.ndarray)
+
+    def test_set_external_variable(self):
+        model_options = {
+            "thermal": "x-lumped",
+            "external submodels": ["thermal"],
+        }
+        model = pybamm.lithium_ion.SPMe(model_options)
+        sim = pybamm.Simulation(model)
+
+        T_av = 0
+
+        dt = 0.001
+
+        external_variables = {"X-averaged cell temperature": T_av}
+        sim.step(dt, external_variables=external_variables)
+
+    def test_step(self):
+
+        dt = 0.001
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
+        sim.step(dt)  # 1 step stores first two points
+        self.assertEqual(sim.solution.t.size, 2)
+        self.assertEqual(sim.solution.y[0, :].size, 2)
+        self.assertEqual(sim.solution.t[0], 0)
+        self.assertEqual(sim.solution.t[1], dt)
+        sim.step(dt)  # automatically append the next step
+        self.assertEqual(sim.solution.t.size, 3)
+        self.assertEqual(sim.solution.y[0, :].size, 3)
+        self.assertEqual(sim.solution.t[0], 0)
+        self.assertEqual(sim.solution.t[1], dt)
+        self.assertEqual(sim.solution.t[2], 2 * dt)
+        sim.step(dt, save=False)  # now only store the two end step points
+        self.assertEqual(sim.solution.t.size, 2)
+        self.assertEqual(sim.solution.y[0, :].size, 2)
+        self.assertEqual(sim.solution.t[0], 2 * dt)
+        self.assertEqual(sim.solution.t[1], 3 * dt)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

From e57de2c03585cc2fcd18ba96c9a4f271e2082165 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Wed, 13 Nov 2019 12:32:27 +0000
Subject: [PATCH 21/31] 692 fixed flake8 and added solver tests

---
 pybamm/discretisations/discretisation.py      |  5 +---
 pybamm/expression_tree/exceptions.py          |  1 +
 pybamm/simulation.py                          |  4 ---
 pybamm/solvers/dae_solver.py                  | 12 ++++++---
 pybamm/solvers/idaklu_solver.py               |  7 +++++-
 .../test_discretisation.py                    |  8 +++++-
 tests/unit/test_simulation.py                 | 25 +++++++++++++++++++
 tests/unit/test_solvers/test_base_solver.py   | 12 +++++++++
 8 files changed, 61 insertions(+), 13 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index e1bb40c5c7..46a7aa7462 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -682,10 +682,7 @@ def process_symbol(self, symbol):
         except KeyError:
             discretised_symbol = self._process_symbol(symbol)
             self._discretised_symbols[symbol.id] = discretised_symbol
-            try:
-                discretised_symbol.test_shape()
-            except:
-                discretised_symbol.test_shape()
+            discretised_symbol.test_shape()
             return discretised_symbol
 
     def _process_symbol(self, symbol):
diff --git a/pybamm/expression_tree/exceptions.py b/pybamm/expression_tree/exceptions.py
index a647b8ab1b..cbeca72b5e 100644
--- a/pybamm/expression_tree/exceptions.py
+++ b/pybamm/expression_tree/exceptions.py
@@ -60,6 +60,7 @@ class UndefinedOperationError(Exception):
 
     pass
 
+
 class InputError(Exception):
     """
     An external variable has been input incorrectly into PyBaMM
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index b7832ed52b..10d980ab61 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -275,10 +275,6 @@ def var_pts(self):
     def spatial_methods(self):
         return self._spatial_methods
 
-    @property
-    def mesh(self):
-        return self._mesh
-
     @property
     def solver(self):
         return self._solver
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index ad106777e3..580bdaf3de 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -248,9 +248,13 @@ def set_up_casadi(self, model):
         y_alg = casadi.MX.sym(
             "y_alg", len(model.concatenated_algebraic.evaluate(0, y0))
         )
-        y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
         y_casadi = casadi.vertcat(y_diff, y_alg)
-        y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
+        if self.y_pad is not None:
+            y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
+            y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
+        else:
+            y_casadi_w_ext = y_casadi
+
         pybamm.logger.info("Converting RHS to CasADi")
         concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi_w_ext)
         pybamm.logger.info("Converting algebraic to CasADi")
@@ -335,7 +339,9 @@ def residuals(t, y, ydot):
 
         # Create event-dependent function to evaluate events
         def event_fun(event):
-            casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi_w_ext], [event])
+            casadi_event_fn = casadi.Function(
+                "event", [t_casadi, y_casadi_w_ext], [event]
+            )
 
             def eval_event(t, y):
                 y = y[:, np.newaxis]
diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py
index 76fd6a631c..6c16cf0082 100644
--- a/pybamm/solvers/idaklu_solver.py
+++ b/pybamm/solvers/idaklu_solver.py
@@ -260,8 +260,13 @@ def rootfn(t, y):
             # 2 = found root(s)
             elif sol.flag == 2:
                 termination = "event"
+
             return pybamm.Solution(
-                sol.t, np.transpose(y_out), t[-1], np.transpose(y_out[-1])[:, np.newaxis], termination
+                sol.t,
+                np.transpose(y_out),
+                t[-1],
+                np.transpose(y_out[-1])[:, np.newaxis],
+                termination,
             )
         else:
             raise pybamm.SolverError(sol.message)
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 39797750cc..b664f1e434 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -94,7 +94,13 @@ def test_adding_1D_external_variable(self):
         }
         model.initial_conditions = {a: 0}
         model.external_variables = [b]
-        model.variables = {"a": a, "b": b, "c": a * b, "grad b": pybamm.grad(b), "div grad b": pybamm.div(pybamm.grad(b))}
+        model.variables = {
+            "a": a,
+            "b": b,
+            "c": a * b,
+            "grad b": pybamm.grad(b),
+            "div grad b": pybamm.div(pybamm.grad(b)),
+        }
 
         x = pybamm.SpatialVariable("x", domain="test", coord_sys="cartesian")
         geometry = {
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index 594a3f093d..d2a1566262 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -130,6 +130,31 @@ def test_specs(self):
         sim.specs(spatial_methods=spatial_methods)
         sim.build()
 
+    def test_set_defaults(self):
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
+
+        model_options = {"thermal": "x-full"}
+        submesh_types = {
+            "Negative particle": pybamm.MeshGenerator(pybamm.Exponential1DSubMesh)
+        }
+        solver = pybamm.IDAKLUSolver()
+        quick_plot_vars = ["Negative particle surface concentration"]
+        sim.specs(
+            model_options=model_options,
+            submesh_types=submesh_types,
+            solver=solver,
+            quick_plot_vars=quick_plot_vars,
+        )
+
+        sim.set_defaults()
+
+        self.assertEqual(sim.model_options["thermal"], "x-full")
+        self.assertEqual(
+            sim.submesh_types["negative particle"].submesh_type, pybamm.Uniform1DSubMesh
+        )
+        self.assertEqual(sim.quick_plot_vars, None)
+        self.assertIsInstance(sim.solver, pybamm.ScipySolver)
+
     def test_get_variable_array(self):
 
         sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index 9d0075e3ef..01fb85bcc9 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -2,6 +2,7 @@
 # Tests for the Base Solver class
 #
 import pybamm
+import numpy as np
 
 import unittest
 
@@ -30,6 +31,17 @@ def test_step_or_solve_empty_model(self):
         with self.assertRaisesRegex(pybamm.ModelError, "Cannot solve empty model"):
             solver.solve(model, None)
 
+    def test_set_external_variables(self):
+        options = {"thermal": "x-full", "external submodels": ["thermal"]}
+        model = pybamm.lithium_ion.SPM(options)
+        sim = pybamm.Simulation(model)
+        sim.build()
+        solver = pybamm.BaseSolver()
+
+        T = np.ones((60, 1))
+        external_variables = {"Cell temperature": T}
+        solver.set_external_variables(sim.built_model, external_variables)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

From dfd292610a1abf0c643ef7a1b00646a7c821973c Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Wed, 13 Nov 2019 12:43:37 +0000
Subject: [PATCH 22/31] 692 updated changelog

---
 CHANGELOG.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 51b38d3992..65b407176c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,9 @@
 
 ## Features
 
+-   Generalized importing of external variables ([#728](https://github.com/pybamm-team/PyBaMM/pull/728))
+-   New extrapolation methods ([#707](https://github.com/pybamm-team/PyBaMM/pull/707))
+-   Allow abs tolerance to be set by variable for IDA KLU solver ([#700](https://github.com/pybamm-team/PyBaMM/pull/700))
 -   Added Simulation class ([#693](https://github.com/pybamm-team/PyBaMM/pull/693))
 -   Added interface to CasADi solver ([#687](https://github.com/pybamm-team/PyBaMM/pull/687), [#691](https://github.com/pybamm-team/PyBaMM/pull/691), [#714](https://github.com/pybamm-team/PyBaMM/pull/714)). This makes the SUNDIALS DAE solvers (Scikits and KLU) truly optional (though IDA KLU is recommended for solving the DFN).
 -   Added option to use CasADi's Algorithmic Differentiation framework to calculate Jacobians ([#687](https://github.com/pybamm-team/PyBaMM/pull/687))

From 751500b1d0e9bb438e0a9c6c6094eff41a976545 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Wed, 13 Nov 2019 12:45:42 +0000
Subject: [PATCH 23/31] 692 updated documentation in simulation slightly

---
 pybamm/simulation.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 10d980ab61..68ce30d09f 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -131,7 +131,9 @@ def step(self, dt, solver=None, external_variables=None, save=True):
             The solver to use to solve the model.
         external_variables : dict
             A dictionary of external variables and their corresponding
-            values at the current time
+            values at the current time. The variables must correspond to
+            the variables that would normally be found by solving the
+            submodels that have been made external.
         save : bool
             Turn on to store the solution of all previous timesteps
         """

From 3932c2300c167d3c58120adb70b2d3f9a161a2ae Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Wed, 13 Nov 2019 15:50:09 +0000
Subject: [PATCH 24/31] 692 fixed external temperature tests

---
 .../test_external/test_external_temperature.py    | 15 ---------------
 1 file changed, 15 deletions(-)

diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
index 00bca26527..a4999d31d0 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
@@ -49,8 +49,6 @@ def test_external_temperature(self):
 
         sim = pybamm.Simulation(model, var_pts=var_pts)
 
-        sim.solver = pybamm.IDAKLUSolver()
-
         t_eval = np.linspace(0, 0.01, 3)
         x = np.linspace(0, 1, tot_pts)
 
@@ -96,21 +94,8 @@ def test_dae_external_temperature(self):
             T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
                 :, np.newaxis
             ]
-            # T = np.ones((tot_pts, 1))
             external_variables = {"Cell temperature": T}
             sim.step(dt, external_variables=external_variables)
-            print(i)
-
-        sim.plot(
-            [
-                "Negative particle surface concentration",
-                "Electrolyte concentration",
-                "Positive particle surface concentration",
-                "Cell temperature [K]",
-                "Terminal voltage [V]",
-                "Total heating [W.m-3]",
-            ]
-        )
 
 
 if __name__ == "__main__":

From b4e8ace3559c7f634d7eb4397846c455018967d1 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Wed, 13 Nov 2019 17:36:05 +0000
Subject: [PATCH 25/31] 692 skip tests that require idaklu

---
 .../test_external/test_external_current_collector.py   | 10 ++++++----
 .../test_external/test_external_temperature.py         |  1 +
 tests/unit/test_simulation.py                          |  2 +-
 3 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
index 57cc570061..0787d851cf 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
@@ -7,6 +7,7 @@
 
 
 class TestExternalCC(unittest.TestCase):
+    @unittest.skipIf(not pybamm.have_idaklu(), "idaklu solver is not installed")
     def test_2p1d(self):
         model_options = {
             "current collector": "potential pair",
@@ -24,7 +25,8 @@ def test_2p1d(self):
             pybamm.standard_spatial_vars.y: yz_pts,
             pybamm.standard_spatial_vars.z: yz_pts,
         }
-        sim = pybamm.Simulation(model, var_pts=var_pts)
+        solver = pybamm.IDAKLUSolver()
+        sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver)
 
         t_eval = np.linspace(0, 0.08, 3)
 
@@ -33,8 +35,8 @@ def test_2p1d(self):
             print(t_eval[i])
 
             # provide phi_s_n and i_cc
-            phi_s_n = np.zeros((yz_pts**2, 1))
-            i_boundary_cc = np.ones((yz_pts**2, 1))
+            phi_s_n = np.zeros((yz_pts ** 2, 1))
+            i_boundary_cc = np.ones((yz_pts ** 2, 1))
             external_variables = {
                 "Negative current collector potential": phi_s_n,
                 "Current collector current density": i_boundary_cc,
@@ -45,7 +47,7 @@ def test_2p1d(self):
             # obtain phi_s_n from the pybamm solution at the current time
             phi_s_p = sim.get_variable_array("Positive current collector potential")
 
-        self.assertTrue(phi_s_p.shape, (yz_pts**2, 1))
+        self.assertTrue(phi_s_p.shape, (yz_pts ** 2, 1))
 
 
 if __name__ == "__main__":
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
index a4999d31d0..29141aa8ba 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
@@ -25,6 +25,7 @@ def test_external_lumped_temperature(self):
             T_av += 1
             sim.step(dt, external_variables=external_variables)
 
+    @unittest.skipIf(not pybamm.have_idaklu(), "idaklu solver is not installed")
     def test_external_temperature(self):
 
         model_options = {
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index d2a1566262..092ef2ebe2 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -137,7 +137,7 @@ def test_set_defaults(self):
         submesh_types = {
             "Negative particle": pybamm.MeshGenerator(pybamm.Exponential1DSubMesh)
         }
-        solver = pybamm.IDAKLUSolver()
+        solver = pybamm.BaseSolver()
         quick_plot_vars = ["Negative particle surface concentration"]
         sim.specs(
             model_options=model_options,

From 18e90bb3df5f840707610a50e7846c063627ab4b Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Thu, 14 Nov 2019 13:51:23 +0000
Subject: [PATCH 26/31] 692 adressed tinos comments

---
 pybamm/discretisations/discretisation.py                    | 6 +++---
 .../test_external/test_external_temperature.py              | 1 +
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index e21c4efe47..51b5f6750d 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -166,6 +166,8 @@ def process_model(self, model, inplace=True):
             elif isinstance(var, pybamm.Variable):
                 start_vals += [self.y_slices[var.id][0].start]
 
+        # attach properties of the state vector so that it
+        # can be divided correctly during the solving stage
         model_disc.external_variables = model.external_variables
         model_disc.y_length = self.y_length
         model_disc.y_slices = self.y_slices
@@ -274,9 +276,7 @@ def _preprocess_external_variables(self, model):
         """
 
         for var in model.external_variables:
-            if var.domain == []:
-                pass
-            else:
+            if var.domain != []:
                 new_bcs = self.spatial_methods[
                     var.domain[0]
                 ].preprocess_external_variables(var)
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
index 29141aa8ba..26600fc463 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
@@ -61,6 +61,7 @@ def test_external_temperature(self):
             external_variables = {"Cell temperature": T}
             sim.step(dt, external_variables=external_variables)
 
+    @unittest.skipIf(not pybamm.have_idaklu(), "idaklu solver is not installed")
     def test_dae_external_temperature(self):
 
         model_options = {

From 4528403b98494c86fd8a5d15307d4b68775d74e7 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Tue, 19 Nov 2019 11:04:37 +0000
Subject: [PATCH 27/31] 692 addressed roberts comments

---
 pybamm/discretisations/discretisation.py      |   2 +
 .../full_battery_models/base_battery_model.py |   7 +
 pybamm/simulation.py                          |  37 +++-
 .../2plus1D/set_temperature_spm_1plus1D.py    | 205 +++---------------
 tests/unit/test_simulation.py                 |   6 +
 5 files changed, 80 insertions(+), 177 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 51b5f6750d..5b9a88beb2 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -158,6 +158,8 @@ def process_model(self, model, inplace=True):
 
         self.external_variables = model.external_variables
         # find where external variables begin in state vector
+        # we always append external variables to the end, so
+        # it is sufficient to only know the starting location
         start_vals = []
         for var in self.external_variables:
             if isinstance(var, pybamm.Concatenation):
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 88fceec674..b16e88a637 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -726,6 +726,11 @@ def set_voltage_variables(self):
             V = pybamm.BoundaryValue(phi_s_cp, "positive tab")
             V_dim = pybamm.BoundaryValue(phi_s_cp_dim, "positive tab")
 
+        phi_s_cn = self.variables["Negative current collector potential"]
+        phi_s_cn_dim = self.variables["Negative current collector potential [V]"]
+        V_local = phi_s_cp - phi_s_cn
+        V_local_dim = phi_s_cp_dim - phi_s_cn_dim
+
         # TODO: add current collector losses to the voltage in 3D
 
         self.variables.update(
@@ -738,6 +743,8 @@ def set_voltage_variables(self):
                 "X-averaged reaction overpotential [V]": eta_r_av_dim,
                 "X-averaged solid phase ohmic losses": delta_phi_s_av,
                 "X-averaged solid phase ohmic losses [V]": delta_phi_s_av_dim,
+                "Local voltage": V_local,
+                "Local voltage [V]": V_local_dim,
                 "Terminal voltage": V,
                 "Terminal voltage [V]": V_dim,
             }
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 68ce30d09f..187af10688 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -10,6 +10,26 @@ class Simulation:
     ----------
     model : :class:`pybamm.BaseModel`
         The model to be simulated
+    geometry: :class:`pybamm.Geometry` (optional)
+        The geometry upon which to solve the model
+    parameter_values: dict (optional)
+        A dictionary of parameters and their corresponding numerical
+        values
+    submesh_types: dict (optional)
+        A dictionary of the types of submesh to use on each subdomain
+    var_pts: dict (optional)
+        A dictionary of the number of points used by each spatial
+        variable
+    spatial_methods: dict (optional)
+        A dictionary of the types of spatial method to use on each
+        domain (e.g. pybamm.FiniteVolume)
+    solver: :class:`pybamm.BaseSolver` (optional)
+        The solver to use to solve the model.
+    quick_plot_vars: list (optional)
+        A list of variables to plot automatically
+    C_rate: float (optional)
+        The C_rate at which you would like to run a constant current
+        experiment at.
     """
 
     def __init__(
@@ -22,6 +42,7 @@ def __init__(
         spatial_methods=None,
         solver=None,
         quick_plot_vars=None,
+        C_rate=None,
     ):
         self.model = model
 
@@ -33,6 +54,10 @@ def __init__(
         self._solver = solver or self._model.default_solver
         self._quick_plot_vars = quick_plot_vars
 
+        self.C_rate = C_rate
+        if self.C_rate:
+            self._parameter_values.update({"C-rate": self.C_rate})
+
         self._made_first_step = False
 
         self.reset(update_model=False)
@@ -307,6 +332,7 @@ def specs(
         spatial_methods=None,
         solver=None,
         quick_plot_vars=None,
+        C_rate=None,
     ):
         """
         A method to set the various specs of the simulation. This method
@@ -329,10 +355,13 @@ def specs(
         spatial_methods: dict (optional)
             A dictionary of the types of spatial method to use on each
             domain (e.g. pybamm.FiniteVolume)
-        solver: :class:`pybamm.BaseSolver`
+        solver: :class:`pybamm.BaseSolver` (optional)
             The solver to use to solve the model.
-        quick_plot_vars: list
+        quick_plot_vars: list (optional)
             A list of variables to plot automatically
+        C_rate: float (optional)
+            The C_rate at which you would like to run a constant current
+            experiment at.
         """
 
         if model_options:
@@ -354,6 +383,10 @@ def specs(
         if quick_plot_vars:
             self._quick_plot_vars = quick_plot_vars
 
+        if C_rate:
+            self.C_rate = C_rate
+            self._parameter_values.update({"C-rate": self.C_rate})
+
         if (
             model_options
             or geometry
diff --git a/results/2plus1D/set_temperature_spm_1plus1D.py b/results/2plus1D/set_temperature_spm_1plus1D.py
index f6101a8824..90c2bd475d 100644
--- a/results/2plus1D/set_temperature_spm_1plus1D.py
+++ b/results/2plus1D/set_temperature_spm_1plus1D.py
@@ -4,189 +4,44 @@
 
 import pybamm
 import numpy as np
-import matplotlib.pyplot as plt
-import sys
 
 # set logging level
 pybamm.set_logging_level("INFO")
 
-# load (1+1D) SPM model
-options = {
+model_options = {
     "current collector": "potential pair",
     "dimensionality": 1,
-    "thermal": "set external temperature",
+    "thermal": "x-lumped",
+    "external submodels": ["thermal"],
 }
-model = pybamm.lithium_ion.SPM(options)
+model = pybamm.lithium_ion.SPMe(model_options)
 
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-param.process_model(model)
-param.process_geometry(geometry)
-
-# set mesh
-nbat = 5
 var = pybamm.standard_spatial_vars
-var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 5, var.r_p: 5, var.z: nbat}
-# depending on number of points in y-z plane may need to increase recursion depth...
-sys.setrecursionlimit(10000)
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-
-# define a method which updates statevector
-def update_statevector(variables, statevector):
-    "takes in a dict of variable name and vector of updated state"
-    for name, new_vector in variables.items():
-        var_slice = model.variables[name].y_slices
-        statevector[var_slice] = new_vector
-    return statevector
-
-
-# define a method to non-dimensionalise a temperature
-def non_dim_temperature(temperature):
-    "takes in a temperature and returns the non-dimensional version"
-    Delta_T = param.process_symbol(model.submodels["thermal"].param.Delta_T).evaluate()
-    T_ref = param.process_symbol(model.submodels["thermal"].param.T_ref).evaluate()
-    return (temperature - T_ref) / Delta_T
-
-
-# step model in time and process variables for later plotting
-solver = model.default_solver
-dt = 0.1  # timestep to take
-npts = 20  # number of points to store the solution at during this step
-solution1 = solver.step(model, dt, npts=npts)
-# create dict of variables to post process
-output_variables = [
-    "Negative current collector potential [V]",
-    "Positive current collector potential [V]",
-    "Current [A]",
-    "X-averaged total heating [W.m-3]",
-    "X-averaged positive particle surface concentration [mol.m-3]",
-    "X-averaged cell temperature [K]",
-]
-output_variables_dict = {}
-for var in output_variables:
-    output_variables_dict[var] = model.variables[var]
-processed_vars_step1 = pybamm.post_process_variables(
-    output_variables_dict, solution1.t, solution1.y, mesh
+z_pts = 5
+var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 5, var.r_p: 5, var.z: z_pts}
+
+sim = pybamm.Simulation(model, var_pts=var_pts, C_rate=2)
+
+# Set the temperature (in dimensionless form)
+T_av = np.linspace(0, 1, z_pts)[:, np.newaxis]
+
+t_eval = np.linspace(0, 0.13, 50)
+# step through the solver, setting the temperature at each timestep
+for i in np.arange(1, len(t_eval) - 1):
+    dt = t_eval[i + 1] - t_eval[i]
+    external_variables = {"X-averaged cell temperature": T_av}
+    sim.step(dt, external_variables=external_variables)
+
+sim.plot(
+    [
+        "Terminal voltage [V]",
+        "X-averaged total heating [W.m-3]",
+        "X-averaged cell temperature [K]",
+        "X-averaged negative particle surface concentration [mol.m-3]",
+        "X-averaged positive particle surface concentration [mol.m-3]",
+        "Negative current collector potential [V]",
+        "Positive current collector potential [V]",
+        "Local voltage [V]",
+    ]
 )
 
-# get the current state and temperature
-current_state = solution1.y[:, -1]
-temp_ave = current_state[model.variables["X-averaged cell temperature"].y_slices]
-
-# update the temperature
-T_ref = param.process_symbol(model.submodels["thermal"].param.T_ref).evaluate()
-t_external = np.linspace(T_ref, T_ref + 6.0, nbat)
-non_dim_t_external = non_dim_temperature(t_external)
-variables = {"X-averaged cell temperature": non_dim_t_external}
-new_state = update_statevector(variables, current_state)
-
-# step in time again and process variables for later plotting
-# use new state as initial condition. Note: need to to recompute consistent initial
-# values for the algebraic part of the model. Since the (dummy) equation for the
-# temperature is an ODE, the imposed change in temperature is unaffected by this
-# process
-solver.y0 = solver.calculate_consistent_initial_conditions(
-    solver.rhs, solver.algebraic, new_state
-)
-solution2 = solver.step(model, dt, npts=npts)
-processed_vars_step2 = pybamm.post_process_variables(
-    output_variables_dict, solution2.t, solution2.y, mesh
-)
-
-# plots
-t_sec = param.evaluate(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_hour = t_sec / (3600)
-z = np.linspace(0, 1, nbat)
-
-# local voltage
-plt.figure()
-for bat_id in range(nbat):
-    plt.plot(
-        solution1.t * t_hour,
-        processed_vars_step1["Positive current collector potential [V]"](
-            solution1.t, z=z
-        )[bat_id, :]
-        - processed_vars_step1["Negative current collector potential [V]"](
-            solution1.t, z=z
-        )[bat_id, :],
-        solution2.t * t_hour,
-        processed_vars_step2["Positive current collector potential [V]"](
-            solution2.t, z=z
-        )[bat_id, :]
-        - processed_vars_step1["Negative current collector potential [V]"](
-            solution2.t, z=z
-        )[bat_id, :],
-    )
-plt.xlabel("t [hrs]")
-plt.ylabel("Local voltage [V]")
-
-# applied current
-plt.figure()
-plt.plot(
-    solution1.t,
-    processed_vars_step1["Current [A]"](solution1.t),
-    solution2.t,
-    processed_vars_step2["Current [A]"](solution2.t),
-)
-plt.xlabel("t")
-plt.ylabel("Current [A]")
-
-# local heating
-plt.figure()
-z = np.linspace(0, 1, nbat)
-for bat_id in range(nbat):
-    plt.plot(
-        solution1.t * t_hour,
-        processed_vars_step1["X-averaged total heating [W.m-3]"](solution1.t, z=z)[
-            bat_id, :
-        ],
-        solution2.t * t_hour,
-        processed_vars_step2["X-averaged total heating [W.m-3]"](solution2.t, z=z)[
-            bat_id, :
-        ],
-    )
-plt.xlabel("t [hrs]")
-plt.ylabel("X-averaged total heating [W.m-3]")
-plt.yscale("log")
-
-# local concentration
-plt.figure()
-for bat_id in range(nbat):
-    plt.plot(
-        solution1.t * t_hour,
-        processed_vars_step1[
-            "X-averaged positive particle surface concentration [mol.m-3]"
-        ](solution1.t, z=z)[bat_id, :],
-        solution2.t * t_hour,
-        processed_vars_step2[
-            "X-averaged positive particle surface concentration [mol.m-3]"
-        ](solution2.t, z=z)[bat_id, :],
-    )
-plt.xlabel("t [hrs]")
-plt.ylabel("X-averaged positive particle surface concentration [mol.m-3]")
-
-# local temperature
-plt.figure()
-for bat_id in range(nbat):
-    plt.plot(
-        solution1.t * t_hour,
-        processed_vars_step1["X-averaged cell temperature [K]"](solution1.t, z=z)[
-            bat_id, :
-        ],
-        solution2.t * t_hour,
-        processed_vars_step2["X-averaged cell temperature [K]"](solution2.t, z=z)[
-            bat_id, :
-        ],
-    )
-plt.xlabel("t [hrs]")
-plt.ylabel("X-averaged cell temperature [K]")
-
-plt.show()
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index 092ef2ebe2..1c1cfb5e49 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -130,6 +130,12 @@ def test_specs(self):
         sim.specs(spatial_methods=spatial_methods)
         sim.build()
 
+    def test_set_crate(self):
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM(), C_rate=2)
+        self.assertEqual(sim.parameter_values["C-rate"], 2)
+        sim.specs(C_rate=3)
+        self.assertEqual(sim.parameter_values["C-rate"], 3)
+
     def test_set_defaults(self):
         sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
 

From faca42f0fe10b2aa06871a52adcfa1315973e898 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Tue, 19 Nov 2019 11:56:14 +0000
Subject: [PATCH 28/31] 692 updated example

---
 results/2plus1D/set_temperature_spm_1plus1D.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/results/2plus1D/set_temperature_spm_1plus1D.py b/results/2plus1D/set_temperature_spm_1plus1D.py
index 90c2bd475d..c469d76079 100644
--- a/results/2plus1D/set_temperature_spm_1plus1D.py
+++ b/results/2plus1D/set_temperature_spm_1plus1D.py
@@ -17,18 +17,22 @@
 model = pybamm.lithium_ion.SPMe(model_options)
 
 var = pybamm.standard_spatial_vars
-z_pts = 5
+z_pts = 20
 var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 5, var.r_p: 5, var.z: z_pts}
 
 sim = pybamm.Simulation(model, var_pts=var_pts, C_rate=2)
 
 # Set the temperature (in dimensionless form)
-T_av = np.linspace(0, 1, z_pts)[:, np.newaxis]
+# T_av = np.linspace(0, 1, z_pts)[:, np.newaxis]
 
+z = np.linspace(0, 1, z_pts)
 t_eval = np.linspace(0, 0.13, 50)
 # step through the solver, setting the temperature at each timestep
 for i in np.arange(1, len(t_eval) - 1):
     dt = t_eval[i + 1] - t_eval[i]
+    T_av = (np.sin(2 * np.pi * z) * np.sin(2 * np.pi * 100 * t_eval[i]))[
+        :, np.newaxis
+    ]
     external_variables = {"X-averaged cell temperature": T_av}
     sim.step(dt, external_variables=external_variables)
 

From 41694ca94bd12453a54b6a57d0ab42d02b90cdc6 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Wed, 20 Nov 2019 12:14:35 +0000
Subject: [PATCH 29/31] 692 passes external variables test

---
 pybamm/solvers/base_solver.py                 |  23 ++--
 pybamm/solvers/dae_solver.py                  | 100 ++++++------------
 pybamm/solvers/ode_solver.py                  |  57 ++++------
 .../test_external_current_collector.py        |   1 -
 4 files changed, 64 insertions(+), 117 deletions(-)

diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index fff622829a..5ef78b690e 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -173,7 +173,7 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
         for i in np.arange(solution.y.shape[1]):
             sol_y = solution.y[:, i]
             sol_y = sol_y[:, np.newaxis]
-            full_y[:, i] = self.add_external(sol_y)[:, 0]
+            full_y[:, i] = add_external(sol_y, self.y_pad, self.y_ext)[:, 0]
         solution.y = full_y
 
         # Assign times
@@ -221,15 +221,6 @@ def set_external_variables(self, model, external_variables):
                 )
             self.y_ext[y_slice] = var_vals
 
-    def add_external(self, y):
-        """
-        Pad the state vector and then add the external variables so that
-        it is of the correct shape for evaluate
-        """
-        if self.y_pad is not None and self.y_ext is not None:
-            y = np.concatenate([y, self.y_pad]) + self.y_ext
-        return y
-
     def compute_solution(self, model, t_eval):
         """Calculate the solution of the model at specified times. Note: this
         does *not* execute the solver setup.
@@ -290,7 +281,7 @@ def get_termination_reason(self, solution, events):
             # Get final event value
             final_event_values = {}
             for name, event in events.items():
-                y_event = self.add_external(solution.y_event)
+                y_event = add_external(solution.y_event, self.y_pad, self.y_ext)
                 final_event_values[name] = abs(
                     event.evaluate(solution.t_event, y_event)
                 )
@@ -298,3 +289,13 @@ def get_termination_reason(self, solution, events):
             # Add the event to the solution object
             solution.termination = "event: {}".format(termination_event)
             return "the termination event '{}' occurred".format(termination_event)
+
+
+def add_external(y, y_pad, y_ext):
+    """
+    Pad the state vector and then add the external variables so that
+    it is of the correct shape for evaluate
+    """
+    if y_pad is not None and y_ext is not None:
+        y = np.concatenate([y, y_pad]) + y_ext
+    return y
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 554eb263d0..0e073e6284 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -7,6 +7,8 @@
 from scipy import optimize
 from scipy.sparse import issparse
 
+from .base_solver import add_external
+
 
 class DaeSolver(pybamm.BaseSolver):
     """Solve a discretised model.
@@ -79,6 +81,14 @@ def compute_solution(self, model, t_eval):
         """
         timer = pybamm.Timer()
 
+        # update y_pad and y_ext
+        self.rhs.set_pad_ext(self.y_pad, self.y_ext)
+        self.algebraic.set_pad_ext(self.y_pad, self.y_ext)
+        self.residuals.set_pad_ext(self.y_pad, self.y_ext)
+        for evnt in self.event_funs:
+            evnt.set_pad_ext(self.y_pad, self.y_ext)
+        self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+
         solve_start_time = timer.time()
         pybamm.logger.info("Calling DAE solver")
         solution = self.integrate(
@@ -147,11 +157,11 @@ def set_up(self, model):
                 jac_algebraic = pybamm.EvaluatorPython(jac_algebraic)
                 jac = pybamm.EvaluatorPython(jac)
 
-            jacobian = Jacobian(jac.evaluate, self.y_pad, self.y_ext)
+            jacobian = Jacobian(jac.evaluate)
 
             def jacobian_alg(t, y):
                 y = y[:, np.newaxis]
-                y = self.add_external(y)
+                y = add_external(y, self.y_pad, self.y_ext)
                 return jac_algebraic.evaluate(t, y, known_evals={})[0]
 
         else:
@@ -169,8 +179,8 @@ def jacobian_alg(t, y):
             }
 
         # Calculate consistent initial conditions for the algebraic equations
-        rhs = Rhs(concatenated_rhs.evaluate, self.y_pad, self.y_ext)
-        algebraic = Algebraic(concatenated_algebraic.evaluate, self.y_pad, self.y_ext)
+        rhs = Rhs(concatenated_rhs.evaluate)
+        algebraic = Algebraic(concatenated_algebraic.evaluate)
 
         if len(model.algebraic) > 0:
             y0 = self.calculate_consistent_initial_conditions(
@@ -185,7 +195,7 @@ def jacobian_alg(t, y):
 
         # Create event-dependent function to evaluate events
         def get_event_class(event):
-            return EvalEvent(event.evaluate, self.y_pad, self.y_ext)
+            return EvalEvent(event.evaluate)
 
         # Add the solver attributes
         self.y0 = y0
@@ -210,7 +220,7 @@ def set_up_casadi(self, model):
         # Convert model attributes to casadi
         t_casadi = casadi.MX.sym("t")
         y0 = model.concatenated_initial_conditions
-        y0 = self.add_external(y0)
+        y0 = add_external(y0, self.y_pad, self.y_ext)
 
         y_diff = casadi.MX.sym("y_diff", len(model.concatenated_rhs.evaluate(0, y0)))
         y_alg = casadi.MX.sym(
@@ -257,19 +267,22 @@ def set_up_casadi(self, model):
                 "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac_alg]
             )
 
-            jacobian = JacobianCasadi(casadi_jac_fn, self.y_pad, self.y_ext)
+            jacobian = JacobianCasadi(casadi_jac_fn)
 
             def jacobian_alg(t, y):
                 y = y[:, np.newaxis]
-                y = self.add_external(y)
+                y = add_external(y, self.y_pad, self.y_ext)
                 return casadi_jac_alg_fn(t, y)
 
         else:
             jacobian = None
             jacobian_alg = None
 
-        rhs = RhsCasadi(concatenated_rhs_fn, self.y_pad, self.y_ext)
-        algebraic = AlgebraicCasadi(concatenated_algebraic_fn, self.y_pad, self.y_ext)
+        rhs = RhsCasadi(concatenated_rhs_fn)
+        algebraic = AlgebraicCasadi(concatenated_algebraic_fn)
+
+        rhs.set_pad_ext(self.y_pad, self.y_ext)
+        algebraic.set_pad_ext(self.y_pad, self.y_ext)
 
         if len(model.algebraic) > 0:
 
@@ -285,8 +298,10 @@ def jacobian_alg(t, y):
 
         # Create event-dependent function to evaluate events
         def get_event_class(event):
-            casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi], [event])
-            return EvalEvent(casadi_event_fn, self.y_pad, self.y_ext)
+            casadi_event_fn = casadi.Function(
+                "event", [t_casadi, y_casadi_w_ext], [event]
+            )
+            return EvalEvent(casadi_event_fn)
 
         # Add the solver attributes
         # Note: these are the (possibly) converted to python version rhs, algebraic
@@ -426,10 +441,8 @@ def integrate(
 class Rhs:
     "Returns information about rhs at time t and state y"
 
-    def __init__(self, concatenated_rhs_fn, y_pad, y_ext):
+    def __init__(self, concatenated_rhs_fn):
         self.concatenated_rhs_fn = concatenated_rhs_fn
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -444,14 +457,6 @@ def __call__(self, t, y):
 class RhsCasadi(Rhs):
     "Returns information about rhs at time t and state y, with CasADi"
 
-    def __init__(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
@@ -461,10 +466,8 @@ def __call__(self, t, y):
 class Algebraic:
     "Returns information about algebraic equations at time t and state y"
 
-    def __init__(self, concatenated_algebraic_fn, y_pad, y_ext):
+    def __init__(self, concatenated_algebraic_fn):
         self.concatenated_algebraic_fn = concatenated_algebraic_fn
-        self.y_pad
-        self.y_ext
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -479,14 +482,6 @@ def __call__(self, t, y):
 class AlgebraicCasadi(Algebraic):
     "Returns information about algebraic equations at time t and state y, with CasADi"
 
-    def __init__(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
@@ -496,13 +491,11 @@ def __call__(self, t, y):
 class Residuals:
     "Returns information about residuals at time t and state y"
 
-    def __init__(self, model, concatenated_rhs_fn, concatenated_algebraic_fn, y_pad, y_ext):
+    def __init__(self, model, concatenated_rhs_fn, concatenated_algebraic_fn):
         self.model = model
         self.concatenated_rhs_fn = concatenated_rhs_fn
         self.concatenated_algebraic_fn = concatenated_algebraic_fn
         self.mass_matrix = model.mass_matrix.entries
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -526,16 +519,10 @@ def __call__(self, t, y, ydot):
 class ResidualsCasadi(Residuals):
     "Returns information about residuals at time t and state y, with CasADi"
 
-    def __init__(self, model, all_states_fn, y_pad, y_ext):
+    def __init__(self, model, all_states_fn):
         self.model = model
         self.all_states_fn = all_states_fn
         self.mass_matrix = model.mass_matrix.entries
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def __call__(self, t, y, ydot):
         pybamm.logger.debug(
@@ -550,10 +537,8 @@ def __call__(self, t, y, ydot):
 class EvalEvent:
     "Returns information about events at time t and state y"
 
-    def __init__(self, event_fn, y_pad, y_ext):
+    def __init__(self, event_fn):
         self.event_fn = event_fn
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -568,10 +553,8 @@ def __call__(self, t, y):
 class Jacobian:
     "Returns information about the jacobian at time t and state y"
 
-    def __init__(self, jac_fn, y_pad, y_ext):
+    def __init__(self, jac_fn):
         self.jac_fn = jac_fn
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -586,25 +569,8 @@ def __call__(self, t, y):
 class JacobianCasadi(Jacobian):
     "Returns information about the jacobian at time t and state y, with CasADi"
 
-    def __init__(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
         return self.jac_fn(t, y)
 
-
-def add_external(y, y_pad, y_ext):
-    """
-    Pad the state vector and then add the external variables so that
-    it is of the correct shape for evaluate
-    """
-    if y_pad is not None and y_ext is not None:
-        y = np.concatenate([y, y_pad]) + y_ext
-    return y
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 4b11bcc8bc..1f1440919c 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -5,6 +5,8 @@
 import pybamm
 import numpy as np
 
+from .base_solver import add_external
+
 
 class OdeSolver(pybamm.BaseSolver):
     """Solve a discretised model.
@@ -35,6 +37,11 @@ def compute_solution(self, model, t_eval):
         """
         timer = pybamm.Timer()
 
+        self.dydt.set_pad_ext(self.y_pad, self.y_ext)
+        for evnt in self.event_funs:
+            evnt.set_pad_ext(self.y_pad, self.y_ext)
+        self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+
         solve_start_time = timer.time()
         pybamm.logger.info("Calling ODE solver")
         solution = self.integrate(
@@ -121,11 +128,11 @@ def set_up(self, model):
 
         # Create event-dependent function to evaluate events
         def get_event_class(event):
-            return EvalEvent(event.evaluate, self.y_pad, self.y_ext)
+            return EvalEvent(event.evaluate)
 
         # Create function to evaluate jacobian
         if jac_rhs is not None:
-            jacobian = Jacobian(jac_rhs.evaluate, self.y_pad, self.y_ext)
+            jacobian = Jacobian(jac_rhs.evaluate)
         else:
             jacobian = None
 
@@ -186,8 +193,10 @@ def set_up_casadi(self, model):
 
         # Create event-dependent function to evaluate events
         def get_event_class(event):
-            casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi], [event])
-            return EvalEvent(casadi_event_fn, self.y_pad, self.y_ext)
+            casadi_event_fn = casadi.Function(
+                "event", [t_casadi, y_casadi_w_ext], [event]
+            )
+            return EvalEvent(casadi_event_fn)
 
         # Create function to evaluate jacobian
         if model.use_jacobian:
@@ -197,7 +206,7 @@ def get_event_class(event):
                 "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac]
             )
 
-            jacobian = JacobianCasadi(casadi_jac_fn, self.y_pad, self.y_ext)
+            jacobian = JacobianCasadi(casadi_jac_fn)
 
         else:
             jacobian = None
@@ -238,9 +247,11 @@ def integrate(
 class Dydt:
     "Returns information about time derivatives at time t and state y"
 
-    def __init__(self, model, concatenated_rhs_fn, y_pad, y_ext):
+    def __init__(self, model, concatenated_rhs_fn):
         self.model = model
         self.concatenated_rhs_fn = concatenated_rhs_fn
+
+    def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
         self.y_ext = y_ext
 
@@ -255,14 +266,6 @@ def __call__(self, t, y):
 class DydtCasadi(Dydt):
     "Returns information about time derivatives at time t and state y, with CasADi"
 
-    def __init__(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
     def __call__(self, t, y):
         pybamm.logger.debug("Evaluating RHS for {} at t={}".format(self.model.name, t))
         y = y[:, np.newaxis]
@@ -274,10 +277,8 @@ def __call__(self, t, y):
 class EvalEvent:
     "Returns information about events at time t and state y"
 
-    def __init__(self, event_fn, y_pad, y_ext):
+    def __init__(self, event_fn):
         self.event_fn = event_fn
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -292,10 +293,8 @@ def __call__(self, t, y):
 class Jacobian:
     "Returns information about the jacobian at time t and state y"
 
-    def __init__(self, jac_fn, y_pad, y_ext):
+    def __init__(self, jac_fn):
         self.jac_fn = jac_fn
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -310,25 +309,7 @@ def __call__(self, t, y):
 class JacobianCasadi(Jacobian):
     "Returns information about the jacobian at time t and state y, with CasADi"
 
-    def __init__(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
-
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
         return self.jac_fn(t, y)
-
-
-def add_external(y, y_pad, y_ext):
-    """
-    Pad the state vector and then add the external variables so that
-    it is of the correct shape for evaluate
-    """
-    if y_pad is not None and y_ext is not None:
-        y = np.concatenate([y, y_pad]) + y_ext
-    return y
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
index 0787d851cf..7f0a661fda 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
@@ -32,7 +32,6 @@ def test_2p1d(self):
 
         for i in np.arange(1, len(t_eval) - 1):
             dt = t_eval[i + 1] - t_eval[i]
-            print(t_eval[i])
 
             # provide phi_s_n and i_cc
             phi_s_n = np.zeros((yz_pts ** 2, 1))

From 86788391c3c271b961f79bcb52b0c0d9ba15bc48 Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Wed, 20 Nov 2019 12:28:45 +0000
Subject: [PATCH 30/31] 692 fixed some tests

---
 pybamm/solvers/dae_solver.py  | 15 ++++++++++++++-
 pybamm/solvers/ode_solver.py  |  9 ++++++++-
 tests/unit/test_simulation.py |  3 +--
 3 files changed, 23 insertions(+), 4 deletions(-)

diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 0e073e6284..33d312de49 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -87,7 +87,8 @@ def compute_solution(self, model, t_eval):
         self.residuals.set_pad_ext(self.y_pad, self.y_ext)
         for evnt in self.event_funs:
             evnt.set_pad_ext(self.y_pad, self.y_ext)
-        self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+        if self.jacobian:
+            self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
 
         solve_start_time = timer.time()
         pybamm.logger.info("Calling DAE solver")
@@ -443,6 +444,8 @@ class Rhs:
 
     def __init__(self, concatenated_rhs_fn):
         self.concatenated_rhs_fn = concatenated_rhs_fn
+        self.y_pad = None
+        self.y_ext = None
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -468,6 +471,8 @@ class Algebraic:
 
     def __init__(self, concatenated_algebraic_fn):
         self.concatenated_algebraic_fn = concatenated_algebraic_fn
+        self.y_pad = None
+        self.y_ext = None
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -496,6 +501,8 @@ def __init__(self, model, concatenated_rhs_fn, concatenated_algebraic_fn):
         self.concatenated_rhs_fn = concatenated_rhs_fn
         self.concatenated_algebraic_fn = concatenated_algebraic_fn
         self.mass_matrix = model.mass_matrix.entries
+        self.y_pad = None
+        self.y_ext = None
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -523,6 +530,8 @@ def __init__(self, model, all_states_fn):
         self.model = model
         self.all_states_fn = all_states_fn
         self.mass_matrix = model.mass_matrix.entries
+        self.y_pad = None
+        self.y_ext = None
 
     def __call__(self, t, y, ydot):
         pybamm.logger.debug(
@@ -539,6 +548,8 @@ class EvalEvent:
 
     def __init__(self, event_fn):
         self.event_fn = event_fn
+        self.y_pad = None
+        self.y_ext = None
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -555,6 +566,8 @@ class Jacobian:
 
     def __init__(self, jac_fn):
         self.jac_fn = jac_fn
+        self.y_pad = None
+        self.y_ext = None
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 1f1440919c..934e521ef5 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -40,7 +40,8 @@ def compute_solution(self, model, t_eval):
         self.dydt.set_pad_ext(self.y_pad, self.y_ext)
         for evnt in self.event_funs:
             evnt.set_pad_ext(self.y_pad, self.y_ext)
-        self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+        if self.jacobian:
+            self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
 
         solve_start_time = timer.time()
         pybamm.logger.info("Calling ODE solver")
@@ -250,6 +251,8 @@ class Dydt:
     def __init__(self, model, concatenated_rhs_fn):
         self.model = model
         self.concatenated_rhs_fn = concatenated_rhs_fn
+        self.y_pad = None
+        self.y_ext = None
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -279,6 +282,8 @@ class EvalEvent:
 
     def __init__(self, event_fn):
         self.event_fn = event_fn
+        self.y_pad = None
+        self.y_ext = None
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
@@ -295,6 +300,8 @@ class Jacobian:
 
     def __init__(self, jac_fn):
         self.jac_fn = jac_fn
+        self.y_pad = None
+        self.y_ext = None
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index 5cb7dcdc7e..21755dbd7f 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -1,7 +1,6 @@
 import pybamm
 import numpy as np
 import unittest
-import numpy as np
 
 
 class TestSimulation(unittest.TestCase):
@@ -277,7 +276,7 @@ def test_save_load_klu(self):
         sim_load = pybamm.load_sim("test.pickle")
         self.assertEqual(sim.model.name, sim_load.model.name)
 
-    def test_set_defaults(self):
+    def test_set_defaults2(self):
         model = pybamm.lithium_ion.SPM()
 
         # make simulation with silly options (should this be allowed?)

From dbdd10c75f32ff43fa6d6a74ceaced7d052a4b4d Mon Sep 17 00:00:00 2001
From: Scott Marquis <marquis@maths.ox.ac.uk>
Date: Wed, 20 Nov 2019 14:00:21 +0000
Subject: [PATCH 31/31] 692 added line to restart tests

---
 .../test_full_battery_models/test_lithium_ion/test_spm.py        | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
index 7e56b605e3..ec064b75ca 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
@@ -16,6 +16,7 @@ def test_basic_processing(self):
 
     def test_basic_processing_1plus1D(self):
         options = {"current collector": "potential pair", "dimensionality": 1}
+
         model = pybamm.lithium_ion.SPM(options)
         var = pybamm.standard_spatial_vars
         var_pts = {