From 9ee941a792f38e6b8ed1a7dc3575786c9320683f Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Tue, 29 Nov 2022 14:02:06 +0000
Subject: [PATCH 01/15] remove external_variables

---
 pybamm/discretisations/discretisation.py      | 163 +++-------------
 pybamm/models/base_model.py                   |  51 +----
 .../full_battery_models/base_battery_model.py |  12 --
 pybamm/models/submodels/base_submodel.py      |  46 -----
 pybamm/parameters/parameter_values.py         |   5 -
 pybamm/solvers/base_solver.py                 |  51 ++---
 pybamm/spatial_methods/finite_volume.py       |  26 ---
 pybamm/spatial_methods/spatial_method.py      |   3 -
 pybamm/spatial_methods/spectral_volume.py     |   2 +-
 .../test_external_current_collector.py        |   1 +
 .../test_external_temperature.py              |   1 +
 .../test_solvers/test_external_variables.py   |  59 ------
 .../integration/test_solvers/test_solution.py |  58 ------
 .../test_discretisation.py                    | 178 ------------------
 .../test_operations/test_convert_to_casadi.py |  25 ---
 .../test_expression_tree/test_variable.py     |  36 ----
 .../test_lithium_ion/test_spme.py             |   1 +
 tests/unit/test_simulation.py                 |  28 ---
 tests/unit/test_solvers/test_casadi_solver.py |  23 ---
 .../unit/test_solvers/test_scikits_solvers.py |  24 ---
 tests/unit/test_solvers/test_scipy_solver.py  |  26 ---
 21 files changed, 61 insertions(+), 758 deletions(-)
 delete mode 100644 tests/integration/test_solvers/test_external_variables.py

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index f2dff27cd7..d2e50d616a 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -63,7 +63,6 @@ def __init__(self, mesh=None, spatial_methods=None):
         self._bcs = {}
         self.y_slices = {}
         self._discretised_symbols = {}
-        self.external_variables = {}
 
     @property
     def mesh(self):
@@ -178,11 +177,6 @@ def process_model(
         pybamm.logger.verbose("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)
-        self.set_external_variables(model)
-
         # set boundary conditions (only need key ids for boundary_conditions)
         pybamm.logger.verbose(
             "Discretise boundary conditions for {}".format(model.name)
@@ -242,11 +236,6 @@ def process_model(
             processed_events.append(processed_event)
         model_disc.events = processed_events
 
-        # Set external variables
-        model_disc.external_variables = [
-            self.process_symbol(var) for var in model.external_variables
-        ]
-
         # Create mass matrix
         pybamm.logger.verbose("Create mass matrix for {}".format(model.name))
         model_disc.mass_matrix, model_disc.mass_matrix_inv = self.create_mass_matrix(
@@ -349,70 +338,6 @@ def _get_variable_size(self, variable):
                 size += spatial_method.mesh[dom].npts_for_broadcast_to_nodes * repeats
             return size
 
-    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 != []:
-                new_bcs = self.spatial_methods[
-                    var.domain[0]
-                ].preprocess_external_variables(var)
-
-                model.boundary_conditions.update(new_bcs)
-
-    def set_external_variables(self, model):
-        """
-        Add external variables to the list of variables to account for, being careful
-        about concatenations
-        """
-        for var in model.external_variables:
-            # Find the name in the model variables
-            # Look up dictionary key based on value
-            try:
-                idx = list(model.variables.values()).index(var)
-            except ValueError:
-                raise ValueError(
-                    """
-                    Variable {} must be in the model.variables dictionary to be set
-                    as an external variable
-                    """.format(
-                        var
-                    )
-                )
-            name = list(model.variables.keys())[idx]
-            if isinstance(var, pybamm.Variable):
-                # No need to keep track of the parent
-                self.external_variables[(name, None)] = var
-            elif isinstance(var, pybamm.ConcatenationVariable):
-                start = 0
-                end = 0
-                for child in var.children:
-                    dom = child.domain[0]
-                    if (
-                        self.spatial_methods[dom]._get_auxiliary_domain_repeats(
-                            child.domains
-                        )
-                        > 1
-                    ):
-                        raise NotImplementedError(
-                            "Cannot create 2D external variable with concatenations"
-                        )
-                    end += self._get_variable_size(child)
-                    # Keep a record of the parent
-                    self.external_variables[(name, (var, start, end))] = child
-                    start = end
-
     def set_internal_boundary_conditions(self, model):
         """
         A method to set the internal boundary conditions for the submodel.
@@ -526,16 +451,15 @@ def process_boundary_conditions(self, model):
 
             # check if the boundary condition at the origin for sphere domains is other
             # than no flux
-            if key not in model.external_variables:
-                for subdomain in key.domain:
-                    if self.mesh[subdomain].coord_sys == "spherical polar":
-                        if bcs["left"][0].value != 0 or bcs["left"][1] != "Neumann":
-                            raise pybamm.ModelError(
-                                "Boundary condition at r = 0 must be a homogeneous "
-                                "Neumann condition for {} coordinates".format(
-                                    self.mesh[subdomain].coord_sys
-                                )
+            for subdomain in key.domain:
+                if self.mesh[subdomain].coord_sys == "spherical polar":
+                    if bcs["left"][0].value != 0 or bcs["left"][1] != "Neumann":
+                        raise pybamm.ModelError(
+                            "Boundary condition at r = 0 must be a homogeneous "
+                            "Neumann condition for {} coordinates".format(
+                                self.mesh[subdomain].coord_sys
                             )
+                        )
 
             # Handle any boundary conditions applied on the tabs
             if any("tab" in side for side in list(bcs.keys())):
@@ -956,52 +880,26 @@ def _process_symbol(self, symbol):
             )
 
         elif isinstance(symbol, pybamm.Variable):
-            # Check if variable is a standard variable or an external variable
-            if any(symbol == var for var in self.external_variables.values()):
-                # Look up dictionary key based on value
-                idx = list(self.external_variables.values()).index(symbol)
-                name, parent_and_slice = list(self.external_variables.keys())[idx]
-                if parent_and_slice is None:
-                    # Variable didn't come from a concatenation so we can just create a
-                    # normal external variable using the symbol's name
-                    return pybamm.ExternalVariable(
-                        symbol.name,
-                        size=self._get_variable_size(symbol),
-                        domains=symbol.domains,
-                    )
-                else:
-                    # We have to use a special name since the concatenation doesn't have
-                    # a very informative name. Needs improving
-                    parent, start, end = parent_and_slice
-                    ext = pybamm.ExternalVariable(
-                        name,
-                        size=self._get_variable_size(parent),
-                        domains=parent.domains,
-                    )
-                    out = pybamm.Index(ext, slice(start, end))
-                    out.copy_domains(symbol)
-                    return out
-            else:
-                # add a try except block for a more informative error if a variable
-                # can't be found. This should usually be caught earlier by
-                # model.check_well_posedness, but won't be if debug_mode is False
-                try:
-                    y_slices = self.y_slices[symbol]
-                except KeyError:
-                    raise pybamm.ModelError(
-                        """
-                        No key set for variable '{}'. Make sure it is included in either
-                        model.rhs, model.algebraic, or model.external_variables in an
-                        unmodified form (e.g. not Broadcasted)
-                        """.format(
-                            symbol.name
-                        )
+            # add a try except block for a more informative error if a variable
+            # can't be found. This should usually be caught earlier by
+            # model.check_well_posedness, but won't be if debug_mode is False
+            try:
+                y_slices = self.y_slices[symbol]
+            except KeyError:
+                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)
+                    """.format(
+                        symbol.name
                     )
-                # Add symbol's reference and multiply by the symbol's scale
-                # so that the state vector is of order 1
-                return symbol.reference + symbol.scale * pybamm.StateVector(
-                    *y_slices, domains=symbol.domains
                 )
+            # Add symbol's reference and multiply by the symbol's scale
+            # so that the state vector is of order 1
+            return symbol.reference + symbol.scale * pybamm.StateVector(
+                *y_slices, domains=symbol.domains
+            )
 
         elif isinstance(symbol, pybamm.SpatialVariable):
             return spatial_method.spatial_variable(symbol)
@@ -1083,14 +981,7 @@ 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
             unpacked_variables_set = set(unpacked_variables)
-            external_vars = set(self.external_variables.values())
-            for var in self.external_variables.values():
-                child_vars = set(var.children)
-                external_vars = external_vars.union(child_vars)
-            y_slices_with_external_removed = set(self.y_slices.keys()).difference(
-                external_vars
-            )
-            if unpacked_variables_set != y_slices_with_external_removed:
+            if unpacked_variables_set != set(self.y_slices.keys()):
                 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 45e32641bf..2663a28f0b 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -109,7 +109,6 @@ def __init__(self, name="Unnamed model"):
         self._mass_matrix_inv = None
         self._jacobian = None
         self._jacobian_algebraic = None
-        self.external_variables = []
         self._parameters = None
         self._input_parameters = None
         self._parameter_info = None
@@ -419,7 +418,6 @@ def new_copy(self):
         new_model._boundary_conditions = self.boundary_conditions.copy()
         new_model._variables = self.variables.copy()
         new_model._events = self.events.copy()
-        new_model.external_variables = self.external_variables.copy()
         new_model._variables_casadi = self._variables_casadi.copy()
         return new_model
 
@@ -459,18 +457,6 @@ def build_fundamental_and_external(self):
         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(
-                "Getting external variables for {} submodel ({})".format(
-                    submodel_name, self.name
-                )
-            )
-            external_variables = submodel.get_external_variables()
-
-            self.external_variables += external_variables
-
         self._built_fundamental_and_external = True
 
     def build_coupled_variables(self):
@@ -813,16 +799,7 @@ def check_well_determined(self, post_discretisation):
         all_vars_in_keys = all_vars_in_rhs_keys.union(all_vars_in_algebraic_keys)
         extra_variables_in_equations = all_vars_in_eqns.difference(all_vars_in_keys)
 
-        # get external variables
-        external_vars = set(self.external_variables)
-        for var in self.external_variables:
-            if isinstance(var, pybamm.Concatenation):
-                child_vars = set(var.children)
-                external_vars = external_vars.union(child_vars)
-
-        extra_variables = extra_variables_in_equations.difference(external_vars)
-
-        if extra_variables:
+        if extra_variables_in_equations:
             raise pybamm.ModelError("model is underdetermined (too many variables)")
 
     def check_algebraic_equations(self, post_discretisation):
@@ -858,13 +835,12 @@ def check_variables(self):
 
         vars_in_keys = set()
 
-        model_and_external_variables = (
+        model_keys = (
             list(self.rhs.keys())
             + list(self.algebraic.keys())
-            + self.external_variables
         )
 
-        for var in model_and_external_variables:
+        for var in model_keys:
             if isinstance(var, pybamm.Variable):
                 vars_in_keys.add(var)
             # Key can be a concatenation
@@ -876,8 +852,9 @@ def check_variables(self):
                 raise pybamm.ModelError(
                     """
                     No key set for variable '{}'. Make sure it is included in either
-                    model.rhs, model.algebraic, or model.external_variables in an
-                    unmodified form (e.g. not Broadcasted)
+                    model.rhs or model.algebraic, in an unmodified form
+                    (e.g. not Broadcasted)
+                    
                     """.format(
                         var
                     )
@@ -970,21 +947,13 @@ def export_casadi_objects(self, variable_names, input_parameter_order=None):
         for input_param in self.input_parameters:
             name = input_param.name
             inputs_wrong_order[name] = casadi.MX.sym(name, input_param._expected_size)
-        # Read external variables
-        external_casadi = {}
-        for external_varaiable in self.external_variables:
-            name = external_varaiable.name
-            ev_size = external_varaiable._evaluate_for_shape().shape[0]
-            external_casadi[name] = casadi.MX.sym(name, ev_size)
         # Sort according to input_parameter_order
         if input_parameter_order is None:
             inputs = inputs_wrong_order
         else:
             inputs = {name: inputs_wrong_order[name] for name in input_parameter_order}
-        # Set up external variables and inputs
-        # Put external variables first like the integrator expects
-        ext_and_in = {**external_casadi, **inputs}
-        inputs_stacked = casadi.vertcat(*[p for p in ext_and_in.values()])
+        # Set up inputs
+        inputs_stacked = casadi.vertcat(*[p for p in inputs.values()])
 
         # Convert initial conditions to casadi form
         y0 = self.concatenated_initial_conditions.to_casadi(
@@ -994,7 +963,7 @@ def export_casadi_objects(self, variable_names, input_parameter_order=None):
         z0 = y0[self.concatenated_rhs.size :]
 
         # Convert rhs and algebraic to casadi form and calculate jacobians
-        rhs = self.concatenated_rhs.to_casadi(t_casadi, y_casadi, inputs=ext_and_in)
+        rhs = self.concatenated_rhs.to_casadi(t_casadi, y_casadi, inputs=inputs)
         jac_rhs = casadi.jacobian(rhs, y_casadi)
         algebraic = self.concatenated_algebraic.to_casadi(
             t_casadi, y_casadi, inputs=inputs
@@ -1005,7 +974,7 @@ def export_casadi_objects(self, variable_names, input_parameter_order=None):
         variables = OrderedDict()
         for name in variable_names:
             var = self.variables[name]
-            variables[name] = var.to_casadi(t_casadi, y_casadi, inputs=ext_and_in)
+            variables[name] = var.to_casadi(t_casadi, y_casadi, inputs=inputs)
 
         casadi_dict = {
             "t": t_casadi,
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index ac254ce892..78e55f6016 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -904,18 +904,6 @@ def build_fundamental_and_external(self):
         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(
-                "Getting external variables for {} submodel ({})".format(
-                    submodel_name, self.name
-                )
-            )
-            external_variables = submodel.get_external_variables()
-
-            self.external_variables += external_variables
-
         self._built_fundamental_and_external = True
 
     def build_coupled_variables(self):
diff --git a/pybamm/models/submodels/base_submodel.py b/pybamm/models/submodels/base_submodel.py
index bbabc0d4da..a2f450e04b 100644
--- a/pybamm/models/submodels/base_submodel.py
+++ b/pybamm/models/submodels/base_submodel.py
@@ -150,52 +150,6 @@ def get_fundamental_variables(self):
         """
         return {}
 
-    def get_external_variables(self):
-        """
-        A public method that returns the variables in a submodel which are
-        supplied by an external source.
-
-        Returns
-        -------
-        list :
-            A list of the external variables in the model.
-        """
-
-        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):
-                    if all(
-                        isinstance(child, pybamm.Variable) for child in var.children
-                    ):
-                        list_of_vars += [var]
-
-            # 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]
-                    # also add the ids of the children to unique ids
-                    for child in var.children:
-                        unique_ids += [child]
-
-            # 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]
-                    unique_ids += [var]
-
-        return external_variables
-
     def get_coupled_variables(self, variables):
         """
         A public method that creates and returns the variables in a submodel which
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 0694d79103..506de123e8 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -441,11 +441,6 @@ def process_model(self, unprocessed_model, inplace=True):
 
         model.events = new_events
 
-        # Set external variables
-        model.external_variables = [
-            self.process_symbol(var) for var in unprocessed_model.external_variables
-        ]
-
         # Process timescale
         new_timescale = self.process_symbol(unprocessed_model.timescale)
         if isinstance(new_timescale, pybamm.Scalar):
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 0d7ba94978..d784eddaaf 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -654,7 +654,6 @@ def solve(
         self,
         model,
         t_eval=None,
-        external_variables=None,
         inputs=None,
         initial_conditions=None,
         nproc=None,
@@ -671,9 +670,6 @@ def solve(
             initial_conditions
         t_eval : numeric type
             The times (in seconds) at which to compute the solution
-        external_variables : dict
-            A dictionary of external variables and their corresponding
-            values at the current time
         inputs : dict or list, optional
             A dictionary or list of dictionaries describing any input parameters to
             pass to the model when solving
@@ -753,8 +749,8 @@ def solve(
         # with variable "input_list", which is a list of dictionaries.
         # If "inputs" is a single dict, "inputs_list" is a list of only one dict.
         inputs_list = inputs if isinstance(inputs, list) else [inputs]
-        ext_and_inputs_list = [
-            self._set_up_ext_and_inputs(model, external_variables, inputs)
+        model_inputs_list = [
+            self._set_up_model_inputs(model, inputs)
             for inputs in inputs_list
         ]
 
@@ -788,7 +784,7 @@ def solve(
             # not depend on input parameters. Thefore only `ext_and_inputs[0]`
             # is passed to `set_up`.
             # See https://github.com/pybamm-team/PyBaMM/pull/1261
-            self.set_up(model, ext_and_inputs_list[0], t_eval)
+            self.set_up(model, model_inputs_list[0], t_eval)
             self.models_set_up.update(
                 {model: {"initial conditions": model.concatenated_initial_conditions}}
             )
@@ -804,7 +800,7 @@ def solve(
                 else:
                     # If the new initial conditions are different
                     # and cannot be evaluated directly, set up again
-                    self.set_up(model, ext_and_inputs_list[0], t_eval, ics_only=True)
+                    self.set_up(model, model_inputs_list[0], t_eval, ics_only=True)
                 self.models_set_up[model][
                     "initial conditions"
                 ] = model.concatenated_initial_conditions
@@ -820,7 +816,7 @@ def solve(
         if len(inputs_list) > 1:
             all_inputs_names = set(
                 itertools.chain.from_iterable(
-                    [ext_and_inputs.keys() for ext_and_inputs in ext_and_inputs_list]
+                    [model_inputs.keys() for model_inputs in model_inputs_list]
                 )
             )
             initial_conditions_node_names = set(
@@ -836,12 +832,12 @@ def solve(
         t_eval_dimensionless = t_eval / model.timescale_eval
 
         self._set_initial_conditions(
-            model, t_eval_dimensionless[0], ext_and_inputs_list[0], update_rhs=True
+            model, t_eval_dimensionless[0], model_inputs_list[0], update_rhs=True
         )
 
         # Check initial conditions don't violate events
         self._check_events_with_initial_conditions(
-            t_eval_dimensionless, model, ext_and_inputs_list[0]
+            t_eval_dimensionless, model, model_inputs_list[0]
         )
 
         # Process discontinuities
@@ -865,12 +861,12 @@ def solve(
                     t_eval_dimensionless[end_index - 1] * model.timescale_eval,
                 )
             )
-            ninputs = len(ext_and_inputs_list)
+            ninputs = len(model_inputs_list)
             if ninputs == 1:
                 new_solution = self._integrate(
                     model,
                     t_eval_dimensionless[start_index:end_index],
-                    ext_and_inputs_list[0],
+                    model_inputs_list[0],
                 )
                 new_solutions = [new_solution]
             else:
@@ -880,7 +876,7 @@ def solve(
                         zip(
                             [model] * ninputs,
                             [t_eval_dimensionless[start_index:end_index]] * ninputs,
-                            ext_and_inputs_list,
+                            model_inputs_list,
                         ),
                     )
                     p.close()
@@ -907,7 +903,7 @@ def solve(
                 model.y0 = last_state
                 if len(model.algebraic) > 0:
                     model.y0 = self.calculate_consistent_state(
-                        model, t_eval_dimensionless[end_index], ext_and_inputs_list[0]
+                        model, t_eval_dimensionless[end_index], model_inputs_list[0]
                     )
         solve_time = timer.time()
 
@@ -1067,7 +1063,6 @@ def step(
         model,
         dt,
         npts=2,
-        external_variables=None,
         inputs=None,
         save=True,
     ):
@@ -1088,9 +1083,6 @@ def step(
         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
         inputs : dict, optional
             Any input parameters to pass to the model when solving
         save : bool
@@ -1133,7 +1125,7 @@ def step(
         timer = pybamm.Timer()
 
         # Set up external variables and inputs
-        ext_and_inputs = self._set_up_ext_and_inputs(model, external_variables, inputs)
+        model_inputs = self._set_up_model_inputs(model, inputs)
 
         if (
             isinstance(old_solution, pybamm.EmptySolution)
@@ -1144,14 +1136,14 @@ def step(
                 "Start stepping {} with {}".format(model.name, self.name)
             )
 
-            self.set_up(model, ext_and_inputs)
+            self.set_up(model, model_inputs)
             self.models_set_up.update(
                 {model: {"initial conditions": model.concatenated_initial_conditions}}
             )
         elif model not in self.models_set_up:
             # Run set up if the model has changed
 
-            self.set_up(model, ext_and_inputs)
+            self.set_up(model, model_inputs)
             self.models_set_up.update(
                 {model: {"initial conditions": model.concatenated_initial_conditions}}
             )
@@ -1166,19 +1158,19 @@ def step(
                     old_solution, return_type="ics"
                 )
                 model.y0 = concatenated_initial_conditions.evaluate(
-                    0, inputs=ext_and_inputs
+                    0, inputs=model_inputs
                 )
         set_up_time = timer.time()
 
         # (Re-)calculate consistent initial conditions
-        self._set_initial_conditions(model, t, ext_and_inputs, update_rhs=False)
+        self._set_initial_conditions(model, t, model_inputs, update_rhs=False)
 
         # Non-dimensionalise dt
         dt_dimensionless = dt / model.timescale_eval
         t_eval = np.linspace(t, t + dt_dimensionless, npts)
 
         # Check initial conditions don't violate events
-        self._check_events_with_initial_conditions(t_eval, model, ext_and_inputs)
+        self._check_events_with_initial_conditions(t_eval, model, model_inputs)
 
         # Step
         pybamm.logger.verbose(
@@ -1188,7 +1180,7 @@ def step(
             )
         )
         timer.reset()
-        solution = self._integrate(model, t_eval, ext_and_inputs)
+        solution = self._integrate(model, t_eval, model_inputs)
         solution.solve_time = timer.time()
 
         # Check if extrapolation occurred
@@ -1346,7 +1338,7 @@ def check_extrapolation(self, solution, events):
                         "outside these bounds."
                     )
 
-    def _set_up_ext_and_inputs(self, model, external_variables, inputs):
+    def _set_up_model_inputs(self, model, inputs):
         """Set up external variables and input parameters"""
         inputs = inputs or {}
 
@@ -1362,14 +1354,11 @@ def _set_up_ext_and_inputs(self, model, external_variables, inputs):
                 raise pybamm.SolverError(f"No value provided for input '{name}'")
         inputs = inputs_in_model
 
-        external_variables = external_variables or {}
-
         ordered_inputs_names = list(inputs.keys())
         ordered_inputs_names.sort()
         ordered_inputs = {name: inputs[name] for name in ordered_inputs_names}
 
-        ext_and_inputs = {**external_variables, **ordered_inputs}
-        return ext_and_inputs
+        return ordered_inputs
 
 
 def process(symbol, name, vars_for_processing, use_jacobian=None):
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 8b069ff626..423cdbe356 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -65,32 +65,6 @@ def spatial_variable(self, symbol):
             entries = np.tile(symbol_mesh.nodes, repeats)
         return pybamm.Vector(entries, domains=symbol.domains)
 
-    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
-        """
-
-        new_bcs = {
-            var: {
-                "left": (pybamm.BoundaryGradient(var, "left"), "Neumann"),
-                "right": (pybamm.BoundaryGradient(var, "right"), "Neumann"),
-            }
-        }
-
-        return new_bcs
-
     def gradient(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the gradient operator.
         See :meth:`pybamm.SpatialMethod.gradient`
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 6826414bab..844fb9559e 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -330,9 +330,6 @@ def internal_neumann_condition(
 
         raise NotImplementedError
 
-    def preprocess_external_variables(self, var):
-        return {}
-
     def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
         """
         Returns the boundary value or flux using the approriate expression for the
diff --git a/pybamm/spatial_methods/spectral_volume.py b/pybamm/spatial_methods/spectral_volume.py
index 17fe70e040..af62570cc8 100644
--- a/pybamm/spatial_methods/spectral_volume.py
+++ b/pybamm/spatial_methods/spectral_volume.py
@@ -14,7 +14,7 @@ class SpectralVolume(pybamm.FiniteVolume):
     subdivision of any 1D mesh, so it shouldn't be a problem).
 
     For broadcast and mass_matrix, we follow the default behaviour from
-    SpatialMethod. For spatial_variable, preprocess_external_variables,
+    SpatialMethod. For spatial_variable,
     divergence, divergence_matrix, laplacian, integral,
     definite_integral_matrix, indefinite_integral,
     indefinite_integral_matrix, indefinite_integral_matrix_nodes,
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 977c46a618..d45d7c1560 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
@@ -8,6 +8,7 @@
 
 class TestExternalCC(unittest.TestCase):
     def test_2p1d(self):
+        # To Do: replace external variable with input
         model_options = {
             "current collector": "potential pair",
             "dimensionality": 2,
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 f6b00d0564..de1b803e11 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,6 +7,7 @@
 
 
 class TestExternalThermalModels(unittest.TestCase):
+    # To Do: replace external variable with input
     def test_external_lumped_temperature(self):
         model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
         model = pybamm.lithium_ion.SPMe(model_options)
diff --git a/tests/integration/test_solvers/test_external_variables.py b/tests/integration/test_solvers/test_external_variables.py
deleted file mode 100644
index 130023f3b4..0000000000
--- a/tests/integration/test_solvers/test_external_variables.py
+++ /dev/null
@@ -1,59 +0,0 @@
-#
-# Test solvers with external variables
-#
-import pybamm
-import numpy as np
-import os
-import unittest
-
-
-class TestExternalVariables(unittest.TestCase):
-    def test_on_dfn(self):
-        e_height = 0.25
-
-        model = pybamm.lithium_ion.DFN()
-        geometry = model.default_geometry
-        param = model.default_parameter_values
-        param.update({"Negative electrode conductivity [S.m-1]": "[input]"})
-        inputs = {"Negative electrode conductivity [S.m-1]": e_height}
-        var_pts = {"x_n": 5, "x_s": 5, "x_p": 5, "r_n": 10, "r_p": 10}
-        spatial_methods = model.default_spatial_methods
-
-        solver = pybamm.CasadiSolver()
-        sim = pybamm.Simulation(
-            model=model,
-            geometry=geometry,
-            parameter_values=param,
-            var_pts=var_pts,
-            spatial_methods=spatial_methods,
-            solver=solver,
-        )
-        sim.solve(t_eval=np.linspace(0, 3600, 100), inputs=inputs)
-
-    def test_external_variables_SPMe_thermal(self):
-        model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
-        model = pybamm.lithium_ion.SPMe(model_options)
-        sim = pybamm.Simulation(model)
-        t_eval = np.linspace(0, 100, 3)
-        T_av = 0
-        for i in np.arange(1, len(t_eval) - 1):
-            dt = t_eval[i + 1] - t_eval[i]
-            external_variables = {"Volume-averaged cell temperature": T_av}
-            T_av += 1
-            sim.step(dt, external_variables=external_variables)
-        V = sim.solution["Terminal voltage [V]"].data
-        np.testing.assert_array_less(np.diff(V), 0)
-        # test generate with external variable
-        sim.built_model.generate("test.c", ["Volume-averaged cell temperature"])
-        os.remove("test.c")
-
-
-if __name__ == "__main__":
-    import sys
-
-    print("Add -v for more debug output")
-
-    if "-v" in sys.argv:
-        debug = True
-    pybamm.settings.debug_mode = True
-    unittest.main()
diff --git a/tests/integration/test_solvers/test_solution.py b/tests/integration/test_solvers/test_solution.py
index 3493f66e46..f2837e1a7b 100644
--- a/tests/integration/test_solvers/test_solution.py
+++ b/tests/integration/test_solvers/test_solution.py
@@ -50,64 +50,6 @@ def test_append(self):
             decimal=4,
         )
 
-    def test_append_external_variables(self):
-        model = pybamm.lithium_ion.SPM(
-            {
-                "thermal": "lumped",
-                "external submodels": ["thermal", "negative primary particle"],
-            }
-        )
-        # 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
-        mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
-
-        # discretise model
-        disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-        disc.process_model(model)
-
-        # solve model
-        solver = model.default_solver
-
-        Nr = model.default_var_pts["r_n"]
-
-        T_av = 0
-        c_s_n_av = np.ones((Nr, 1)) * 0.6
-        external_variables = {
-            "Volume-averaged cell temperature": T_av,
-            "X-averaged negative particle concentration": c_s_n_av,
-        }
-
-        # Step
-        dt = 0.1
-        sol_step = None
-        for _ in range(5):
-            sol_step = solver.step(
-                sol_step, model, dt, external_variables=external_variables
-            )
-        np.testing.assert_array_equal(
-            sol_step.all_inputs[0]["Volume-averaged cell temperature"], 0
-        )
-        np.testing.assert_array_equal(
-            sol_step.all_inputs[0]["X-averaged negative particle concentration"], 0.6
-        )
-
-        # Solve
-        t_eval = np.linspace(0, 3600)
-        sol = solver.solve(model, t_eval, external_variables=external_variables)
-        np.testing.assert_array_equal(
-            sol.all_inputs[0]["Volume-averaged cell temperature"], 0
-        )
-        np.testing.assert_array_equal(
-            sol.all_inputs[0]["X-averaged negative particle concentration"], 0.6
-        )
-
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 01dd02a2b4..295d8a2780 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -67,184 +67,6 @@ def test_add_internal_boundary_conditions(self):
         for child in c_e.children:
             self.assertTrue(child in disc.bcs.keys())
 
-    def test_adding_0D_external_variable(self):
-        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)
-
-        self.assertIsInstance(model.variables["b"], pybamm.ExternalVariable)
-        self.assertEqual(model.variables["b"].evaluate(inputs={"b": np.array([1])}), 1)
-
-    def test_adding_0D_external_variable_fail(self):
-        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]
-
-        disc = pybamm.Discretisation()
-        with self.assertRaisesRegex(ValueError, "Variable b must be in the model"):
-            disc.process_model(model)
-
-    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": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}}
-
-        submesh_types = {"test": 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(disc.y_slices[a][0], slice(0, 10, None))
-
-        self.assertEqual(model.y_slices[a][0], slice(0, 10, None))
-        self.assertEqual(model.bounds, disc.bounds)
-
-        b_test = np.ones((10, 1))
-        np.testing.assert_array_equal(
-            model.variables["b"].evaluate(inputs={"b": b_test}), b_test
-        )
-
-        # check that b is added to the boundary conditions
-        model.bcs[b]["left"]
-        model.bcs[b]["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_concatenation_external_variables(self):
-        model = pybamm.BaseModel()
-
-        a = pybamm.Variable("a", domain=["test", "test1"])
-        b1 = pybamm.Variable("b1", domain=["test"])
-        b2 = pybamm.Variable("b2", domain=["test1"])
-        b = pybamm.concatenation(b1, b2)
-
-        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,
-            "b1": b1,
-            "b2": b2,
-            "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")
-        y = pybamm.SpatialVariable("y", domain="test1", coord_sys="cartesian")
-        geometry = {
-            "test": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}},
-            "test1": {y: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}},
-        }
-
-        submesh_types = {
-            "test": pybamm.Uniform1DSubMesh,
-            "test1": pybamm.Uniform1DSubMesh,
-        }
-        var_pts = {x: 10, y: 5}
-        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
-
-        spatial_methods = {
-            "test": pybamm.FiniteVolume(),
-            "test1": pybamm.FiniteVolume(),
-        }
-        disc = pybamm.Discretisation(mesh, spatial_methods)
-        disc.process_model(model)
-
-        self.assertEqual(disc.y_slices[a][0], slice(0, 15, None))
-
-        b_test = np.linspace(0, 1, 15)[:, np.newaxis]
-        np.testing.assert_array_equal(
-            model.variables["b"].evaluate(inputs={"b": b_test}), b_test
-        )
-        np.testing.assert_array_equal(
-            model.variables["b1"].evaluate(inputs={"b": b_test}), b_test[:10]
-        )
-        np.testing.assert_array_equal(
-            model.variables["b2"].evaluate(inputs={"b": b_test}), b_test[10:]
-        )
-
-        # check that b is added to the boundary conditions
-        model.bcs[b]["left"]
-        model.bcs[b]["right"]
-
-        # check that grad and div(grad ) produce the correct shapes
-        self.assertEqual(model.variables["b"].shape_for_testing, (15, 1))
-        self.assertEqual(model.variables["grad b"].shape_for_testing, (16, 1))
-        self.assertEqual(model.variables["div grad b"].shape_for_testing, (15, 1))
-        self.assertEqual(model.variables["b1"].shape_for_testing, (10, 1))
-        self.assertEqual(model.variables["b2"].shape_for_testing, (5, 1))
-
-    def test_adding_2D_external_variable_fail(self):
-        model = pybamm.BaseModel()
-        a = pybamm.Variable(
-            "a",
-            domain=["negative electrode", "separator"],
-            auxiliary_domains={"secondary": "current collector"},
-        )
-        b1 = pybamm.Variable(
-            "b",
-            domain="negative electrode",
-            auxiliary_domains={"secondary": "current collector"},
-        )
-        b2 = pybamm.Variable(
-            "b",
-            domain="separator",
-            auxiliary_domains={"secondary": "current collector"},
-        )
-        b = pybamm.concatenation(b1, b2)
-
-        model.rhs = {a: a * b}
-        model.initial_conditions = {a: 0}
-        model.external_variables = [b]
-        model.variables = {"b": b}
-
-        disc = get_1p1d_discretisation_for_testing()
-        with self.assertRaisesRegex(
-            NotImplementedError, "Cannot create 2D external variable"
-        ):
-            disc.process_model(model)
-
     def test_discretise_slicing(self):
         # create discretisation
         mesh = get_mesh_for_testing()
diff --git a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
index ba91ae8f08..7e3874eb07 100644
--- a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
+++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
@@ -360,31 +360,6 @@ def test_convert_input_parameter(self):
             casadi_inputs["Input 2"] * casadi_y,
         )
 
-    def test_convert_external_variable(self):
-        casadi_t = casadi.MX.sym("t")
-        casadi_y = casadi.MX.sym("y", 10)
-        casadi_inputs = {
-            "External 1": casadi.MX.sym("External 1", 3),
-            "External 2": casadi.MX.sym("External 2", 10),
-        }
-
-        pybamm_y = pybamm.StateVector(slice(0, 10))
-        pybamm_u1 = pybamm.ExternalVariable("External 1", 3)
-        pybamm_u2 = pybamm.ExternalVariable("External 2", 10)
-
-        # External only
-        self.assert_casadi_equal(
-            pybamm_u1.to_casadi(casadi_t, casadi_y, inputs=casadi_inputs),
-            casadi_inputs["External 1"],
-        )
-
-        # More complex
-        expr = pybamm_u2 + pybamm_y
-        self.assert_casadi_equal(
-            expr.to_casadi(casadi_t, casadi_y, inputs=casadi_inputs),
-            casadi_inputs["External 2"] + casadi_y,
-        )
-
     def test_errors(self):
         y = pybamm.StateVector(slice(0, 10))
         with self.assertRaisesRegex(
diff --git a/tests/unit/test_expression_tree/test_variable.py b/tests/unit/test_expression_tree/test_variable.py
index 1dbe3b67db..ee2b97e8d1 100644
--- a/tests/unit/test_expression_tree/test_variable.py
+++ b/tests/unit/test_expression_tree/test_variable.py
@@ -90,42 +90,6 @@ def test_variable_diff(self):
         self.assertEqual(a.diff(b).evaluate(), 0)
 
 
-class TestExternalVariable(unittest.TestCase):
-    def test_external_variable_scalar(self):
-        a = pybamm.ExternalVariable("a", 1)
-        self.assertEqual(a.size, 1)
-
-        self.assertEqual(a.evaluate(inputs={"a": 3}), 3)
-
-        with self.assertRaisesRegex(KeyError, "External variable"):
-            a.evaluate()
-        with self.assertRaisesRegex(TypeError, "inputs should be a dictionary"):
-            a.evaluate(inputs="not a dictionary")
-
-    def test_external_variable_vector(self):
-        a = pybamm.ExternalVariable("a", 10)
-        self.assertEqual(a.size, 10)
-
-        a_test = 2 * np.ones((10, 1))
-        np.testing.assert_array_equal(a.evaluate(inputs={"a": a_test}), a_test)
-        np.testing.assert_array_equal(
-            a.evaluate(inputs={"a": a_test.flatten()}), a_test
-        )
-
-        np.testing.assert_array_equal(a.evaluate(inputs={"a": 2}), a_test)
-
-        with self.assertRaisesRegex(ValueError, "External variable"):
-            a.evaluate(inputs={"a": np.ones((5, 1))})
-
-    def test_external_variable_diff(self):
-        a = pybamm.ExternalVariable("a", 10)
-        b = pybamm.Variable("b")
-        self.assertIsInstance(a.diff(a), pybamm.Scalar)
-        self.assertEqual(a.diff(a).evaluate(), 1)
-        self.assertIsInstance(a.diff(b), pybamm.Scalar)
-        self.assertEqual(a.diff(b).evaluate(), 0)
-
-
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
index af208bd432..b13bae504f 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
@@ -11,6 +11,7 @@ def setUp(self):
         self.model = pybamm.lithium_ion.SPMe
 
     def test_external_variables(self):
+        # To Do: replace external variable with input
         # a concatenation
         model_options = {"external submodels": ["electrolyte diffusion"]}
         model = pybamm.lithium_ion.SPMe(model_options)
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index af9e93499b..8272679e6f 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -149,34 +149,6 @@ def test_set_crate(self):
         self.assertEqual(sim.parameter_values["Current function [A]"], 2 * current_1C)
         self.assertEqual(sim.C_rate, 2)
 
-    def test_set_external_variable(self):
-        model_options = {
-            "thermal": "lumped",
-            "external submodels": ["thermal", "negative primary particle"],
-        }
-        model = pybamm.lithium_ion.SPMe(model_options)
-        sim = pybamm.Simulation(model)
-
-        Nr = model.default_var_pts["r_n"]
-
-        T_av = 0
-        c_s_n_av = np.ones((Nr, 1)) * 0.5
-        external_variables = {
-            "Volume-averaged cell temperature": T_av,
-            "X-averaged negative particle concentration": c_s_n_av,
-        }
-
-        # Step
-        dt = 0.1
-        for _ in range(5):
-            sim.step(dt, external_variables=external_variables)
-        sim.plot(testing=True)
-
-        # Solve
-        t_eval = np.linspace(0, 3600)
-        sim.solve(t_eval, external_variables=external_variables)
-        sim.plot(testing=True)
-
     def test_step(self):
 
         dt = 0.001
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 6b9d073ab6..f4062c0df6 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -440,29 +440,6 @@ def test_model_solver_dae_inputs_in_initial_conditions(self):
             solution.y.full()[-1], 1 * np.exp(-0.1 * solution.t), decimal=5
         )
 
-    def test_model_solver_with_external(self):
-        # Create model
-        model = pybamm.BaseModel()
-        domain = ["negative electrode", "separator", "positive electrode"]
-        var1 = pybamm.Variable("var1", domain=domain)
-        var2 = pybamm.Variable("var2", domain=domain)
-        model.rhs = {var1: -var2}
-        model.initial_conditions = {var1: 1}
-        model.external_variables = [var2]
-        model.variables = {"var1": var1, "var2": var2}
-        # create discretisation
-        mesh = get_mesh_for_testing()
-        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
-        disc = pybamm.Discretisation(mesh, spatial_methods)
-        disc.process_model(model)
-        # Solve
-        solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8)
-        t_eval = np.linspace(0, 10, 100)
-        solution = solver.solve(model, t_eval, external_variables={"var2": 0.5})
-        np.testing.assert_allclose(
-            solution.y.full()[0], 1 - 0.5 * solution.t, rtol=1e-06
-        )
-
     def test_model_solver_with_non_identity_mass(self):
         model = pybamm.BaseModel()
         var1 = pybamm.Variable("var1", domain="negative electrode")
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 9d186a9072..d483d4bcea 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -702,30 +702,6 @@ def test_model_solver_dae_inputs_in_initial_conditions(self):
             solution.y[-1], 1 * np.exp(-0.1 * solution.t), decimal=5
         )
 
-    def test_model_solver_dae_with_external(self):
-        # Create model
-        model = pybamm.BaseModel()
-        domain = ["negative electrode", "separator", "positive electrode"]
-        var1 = pybamm.Variable("var1", domain=domain)
-        var2 = pybamm.Variable("var2", domain=domain)
-        model.rhs = {var1: -var2}
-        model.initial_conditions = {var1: 1}
-        model.external_variables = [var2]
-        model.variables = {"var1": var1, "var2": var2}
-        # No need to set parameters; can use base discretisation (no spatial
-        # operators)
-
-        # create discretisation
-        mesh = get_mesh_for_testing()
-        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
-        disc = pybamm.Discretisation(mesh, spatial_methods)
-        disc.process_model(model)
-        # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
-        t_eval = np.linspace(0, 10, 100)
-        solution = solver.solve(model, t_eval, external_variables={"var2": 0.5})
-        np.testing.assert_allclose(solution.y[0], 1 - 0.5 * solution.t, rtol=1e-06)
-
     def test_solve_ode_model_with_dae_solver_casadi(self):
         model = pybamm.BaseModel()
         model.convert_to_format = "casadi"
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 6363fea6ce..bb1dba09d5 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -371,32 +371,6 @@ def test_model_solver_multiple_inputs_jax_format_error(self):
         ):
             solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
 
-    def test_model_solver_with_external(self):
-        # Create model
-        model = pybamm.BaseModel()
-        domain = ["negative electrode", "separator", "positive electrode"]
-        var1 = pybamm.Variable("var1", domain=domain)
-        var2 = pybamm.Variable("var2", domain=domain)
-        model.rhs = {var1: -var2}
-        model.initial_conditions = {var1: 1}
-        model.external_variables = [var2]
-        model.variables = {"var2": var2}
-        # No need to set parameters; can use base discretisation (no spatial
-        # operators)
-
-        # create discretisation
-        mesh = get_mesh_for_testing()
-        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
-        disc = pybamm.Discretisation(mesh, spatial_methods)
-        disc.process_model(model)
-        # Solve
-        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8)
-        t_eval = np.linspace(0, 10, 100)
-        solution = solver.solve(
-            model, t_eval, external_variables={"var2": 0.5 * np.ones(100)}
-        )
-        np.testing.assert_allclose(solution.y[0], 1 - 0.5 * solution.t, rtol=1e-06)
-
     def test_model_solver_with_event_with_casadi(self):
         # Create model
         model = pybamm.BaseModel()

From 665b4a2f059981c31468f01eaf3234af3b9337c4 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Tue, 29 Nov 2022 14:07:21 +0000
Subject: [PATCH 02/15] style: pre-commit fixes

---
 pybamm/models/base_model.py   | 5 +----
 pybamm/solvers/base_solver.py | 3 +--
 2 files changed, 2 insertions(+), 6 deletions(-)

diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 2663a28f0b..638200de8a 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -835,10 +835,7 @@ def check_variables(self):
 
         vars_in_keys = set()
 
-        model_keys = (
-            list(self.rhs.keys())
-            + list(self.algebraic.keys())
-        )
+        model_keys = list(self.rhs.keys()) + list(self.algebraic.keys())
 
         for var in model_keys:
             if isinstance(var, pybamm.Variable):
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index d784eddaaf..c11a8e3e37 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -750,8 +750,7 @@ def solve(
         # If "inputs" is a single dict, "inputs_list" is a list of only one dict.
         inputs_list = inputs if isinstance(inputs, list) else [inputs]
         model_inputs_list = [
-            self._set_up_model_inputs(model, inputs)
-            for inputs in inputs_list
+            self._set_up_model_inputs(model, inputs) for inputs in inputs_list
         ]
 
         # Cannot use multiprocessing with model in "jax" format

From ccba82c94f40b6ecf659fad103d062e806bfe643 Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Tue, 29 Nov 2022 17:17:10 +0000
Subject: [PATCH 03/15] style

---
 pybamm/models/base_model.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 638200de8a..ede0f69ff1 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -851,7 +851,6 @@ def check_variables(self):
                     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)
-                    
                     """.format(
                         var
                     )

From b725674cb36fb586a4eced45c31f87e9e70d0798 Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Thu, 1 Dec 2022 16:44:36 +0000
Subject: [PATCH 04/15] remove external submodel functionality and test on spme

---
 pybamm/models/base_model.py                   | 55 +++++++++----------
 .../full_battery_models/base_battery_model.py | 55 +++++++++----------
 .../test_lithium_ion/test_spme.py             | 34 ++++++------
 3 files changed, 67 insertions(+), 77 deletions(-)

diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index ede0f69ff1..11d0a5144d 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -453,10 +453,6 @@ def build_fundamental_and_external(self):
             )
             self.variables.update(submodel.get_fundamental_variables())
 
-        # Set the submodels that are external
-        for sub in self.options["external submodels"]:
-            self.submodels[sub].external = True
-
         self._built_fundamental_and_external = True
 
     def build_coupled_variables(self):
@@ -506,38 +502,37 @@ def build_coupled_variables(self):
     def build_model_equations(self):
         # Set model equations
         for submodel_name, submodel in self.submodels.items():
-            if submodel.external is False:
-                pybamm.logger.verbose(
-                    "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
-                )
+            pybamm.logger.verbose(
+                "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
+            )
 
-                submodel.set_rhs(self.variables)
-                pybamm.logger.verbose(
-                    "Setting algebraic for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
+            submodel.set_rhs(self.variables)
+            pybamm.logger.verbose(
+                "Setting algebraic for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
+            )
 
-                submodel.set_algebraic(self.variables)
-                pybamm.logger.verbose(
-                    "Setting boundary conditions for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
+            submodel.set_algebraic(self.variables)
+            pybamm.logger.verbose(
+                "Setting boundary conditions for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
+            )
 
-                submodel.set_boundary_conditions(self.variables)
-                pybamm.logger.verbose(
-                    "Setting initial conditions for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
-                )
-                submodel.set_initial_conditions(self.variables)
-                submodel.set_events(self.variables)
-                pybamm.logger.verbose(
-                    "Updating {} submodel ({})".format(submodel_name, self.name)
+            submodel.set_boundary_conditions(self.variables)
+            pybamm.logger.verbose(
+                "Setting initial conditions for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
-                self.update(submodel)
-                self.check_no_repeated_keys()
+            )
+            submodel.set_initial_conditions(self.variables)
+            submodel.set_events(self.variables)
+            pybamm.logger.verbose(
+                "Updating {} submodel ({})".format(submodel_name, self.name)
+            )
+            self.update(submodel)
+            self.check_no_repeated_keys()
 
     def build_model(self):
         self._build_model()
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 78e55f6016..2ca62e1c05 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -900,10 +900,6 @@ def build_fundamental_and_external(self):
             )
             self.variables.update(submodel.get_fundamental_variables())
 
-        # Set the submodels that are external
-        for sub in self.options["external submodels"]:
-            self.submodels[sub].external = True
-
         self._built_fundamental_and_external = True
 
     def build_coupled_variables(self):
@@ -953,38 +949,37 @@ def build_coupled_variables(self):
     def build_model_equations(self):
         # Set model equations
         for submodel_name, submodel in self.submodels.items():
-            if submodel.external is False:
-                pybamm.logger.verbose(
-                    "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
-                )
+            pybamm.logger.verbose(
+                "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
+            )
 
-                submodel.set_rhs(self.variables)
-                pybamm.logger.verbose(
-                    "Setting algebraic for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
+            submodel.set_rhs(self.variables)
+            pybamm.logger.verbose(
+                "Setting algebraic for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
+            )
 
-                submodel.set_algebraic(self.variables)
-                pybamm.logger.verbose(
-                    "Setting boundary conditions for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
+            submodel.set_algebraic(self.variables)
+            pybamm.logger.verbose(
+                "Setting boundary conditions for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
+            )
 
-                submodel.set_boundary_conditions(self.variables)
-                pybamm.logger.verbose(
-                    "Setting initial conditions for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
-                )
-                submodel.set_initial_conditions(self.variables)
-                submodel.set_events(self.variables)
-                pybamm.logger.verbose(
-                    "Updating {} submodel ({})".format(submodel_name, self.name)
+            submodel.set_boundary_conditions(self.variables)
+            pybamm.logger.verbose(
+                "Setting initial conditions for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
-                self.update(submodel)
-                self.check_no_repeated_keys()
+            )
+            submodel.set_initial_conditions(self.variables)
+            submodel.set_events(self.variables)
+            pybamm.logger.verbose(
+                "Updating {} submodel ({})".format(submodel_name, self.name)
+            )
+            self.update(submodel)
+            self.check_no_repeated_keys()
 
     def build_model(self):
 
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
index b13bae504f..1d3a7dd954 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
@@ -10,23 +10,23 @@ class TestSPMe(BaseUnitTestLithiumIon, unittest.TestCase):
     def setUp(self):
         self.model = pybamm.lithium_ion.SPMe
 
-    def test_external_variables(self):
-        # To Do: replace external variable with input
-        # a concatenation
-        model_options = {"external submodels": ["electrolyte diffusion"]}
-        model = pybamm.lithium_ion.SPMe(model_options)
-        self.assertEqual(
-            model.external_variables[0],
-            model.variables["Porosity times concentration"],
-        )
-
-        # a variable
-        model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
-        model = pybamm.lithium_ion.SPMe(model_options)
-        self.assertEqual(
-            model.external_variables[0],
-            model.variables["Volume-averaged cell temperature"],
-        )
+    # def test_external_variables(self):
+    #     # To Do: replace external variable with input
+    #     # a concatenation
+    #     model_options = {"external submodels": ["electrolyte diffusion"]}
+    #     model = pybamm.lithium_ion.SPMe(model_options)
+    #     self.assertEqual(
+    #         model.external_variables[0],
+    #         model.variables["Porosity times concentration"],
+    #     )
+
+    #     # a variable
+    #     model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
+    #     model = pybamm.lithium_ion.SPMe(model_options)
+    #     self.assertEqual(
+    #         model.external_variables[0],
+    #         model.variables["Volume-averaged cell temperature"],
+    #     )
 
     def test_electrolyte_options(self):
         options = {"electrolyte conductivity": "full"}

From 428c7cbd3a496cd6cc34166e1ebc9a8508fb1cfb Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Mon, 5 Dec 2022 16:55:21 +0000
Subject: [PATCH 05/15] Update thermal test

---
 .../test_external_temperature.py              | 157 ++++++++++--------
 1 file changed, 84 insertions(+), 73 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 de1b803e11..8da30aadff 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,87 +7,98 @@
 
 
 class TestExternalThermalModels(unittest.TestCase):
-    # To Do: replace external variable with input
-    def test_external_lumped_temperature(self):
-        model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
-        model = pybamm.lithium_ion.SPMe(model_options)
+    def test_input_lumped_temperature(self):
+        model = pybamm.lithium_ion.SPMe()
+        parameter_values = model.default_parameter_values
+        # in the default isothermal model, the temperature is everywhere equal
+        # to the ambient temperature
+        parameter_values["Ambient temperature [K]"] = pybamm.InputParameter(
+            "Volume-averaged cell temperature [K]"
+        )
         sim = pybamm.Simulation(model)
 
         t_eval = np.linspace(0, 100, 3)
 
-        T_av = 0
+        T_av = 298
 
         for i in np.arange(1, len(t_eval) - 1):
             dt = t_eval[i + 1] - t_eval[i]
-            external_variables = {"Volume-averaged cell temperature": T_av}
+            inputs = {"Volume-averaged cell temperature [K]": 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)
-
-        t_eval = np.linspace(0, 100, 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):
-
-        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.CasadiSolver()
-        sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver)
-        sim.build()
-
-        t_eval = np.linspace(0, 100, 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)
+            sim.step(dt, inputs=inputs)  # works
+
+    # def test_external_temperature(self):
+
+    #     model = pybamm.lithium_ion.SPM()
+
+    #     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,
+    #     }
+
+    #     parameter_values = model.default_parameter_values
+    #     # in the default isothermal model, the temperature is everywhere equal
+    #     # to the ambient temperature
+    #     parameter_values["Ambient temperature [K]"] = pybamm.InputParameter(
+    #         "Cell temperature [K]",
+    #         domain=["negative electrode", "separator", "positive electrode"],
+    #     )
+    #     sim = pybamm.Simulation(
+    #         model, var_pts=var_pts, parameter_values=parameter_values
+    #     )
+
+    #     t_eval = np.linspace(0, 100, 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
+    #         ]
+    #         inputs = {"Cell temperature [K]": T}
+    #         sim.step(dt, inputs=inputs)  # fails because of broadcasting error
+
+    # 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.CasadiSolver()
+    #     sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver)
+    #     sim.build()
+
+    #     t_eval = np.linspace(0, 100, 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)
 
 
 if __name__ == "__main__":

From 2e1d559c985eafb2dfdb4d4a25d3ca18b915e67d Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Mon, 5 Dec 2022 16:57:18 +0000
Subject: [PATCH 06/15] remove external current collector test for now

---
 .../test_external_current_collector.py        | 81 ++++++++++---------
 1 file changed, 41 insertions(+), 40 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 d45d7c1560..5e12b249cf 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,46 +7,47 @@
 
 
 class TestExternalCC(unittest.TestCase):
-    def test_2p1d(self):
-        # To Do: replace external variable with input
-        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)
-
-        # Simulate 100 seconds
-        t_eval = np.linspace(0, 100, 3)
-
-        for i in np.arange(1, len(t_eval) - 1):
-            dt = t_eval[i + 1] - 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.solution["Positive current collector potential"].data[:, -1]
-
-        self.assertTrue(phi_s_p.shape, (yz_pts**2, 1))
+    pass
+    # def test_2p1d(self):
+    #     # To Do: replace external variable with input
+    #     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)
+
+    #     # Simulate 100 seconds
+    #     t_eval = np.linspace(0, 100, 3)
+
+    #     for i in np.arange(1, len(t_eval) - 1):
+    #         dt = t_eval[i + 1] - 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.solution["Positive current collector potential"].data[:, -1]
+
+    #     self.assertTrue(phi_s_p.shape, (yz_pts**2, 1))
 
 
 if __name__ == "__main__":

From 26ed9952e787cd7e93c99b734d2ad85c613a3b07 Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Mon, 5 Dec 2022 17:07:26 +0000
Subject: [PATCH 07/15] style

---
 .../test_external/test_external_current_collector.py          | 4 ++--
 1 file changed, 2 insertions(+), 2 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 5e12b249cf..df86f9f4dd 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
@@ -1,9 +1,9 @@
 #
 # Tests for external submodels
 #
-import pybamm
+# import pybamm
 import unittest
-import numpy as np
+# import numpy as np
 
 
 class TestExternalCC(unittest.TestCase):

From ec68fe1c23123bee8976e83afabd36303e143c8c Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 5 Dec 2022 17:08:09 +0000
Subject: [PATCH 08/15] style: pre-commit fixes

---
 .../test_external/test_external_current_collector.py             | 1 +
 1 file changed, 1 insertion(+)

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 df86f9f4dd..4b652aeff8 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
@@ -3,6 +3,7 @@
 #
 # import pybamm
 import unittest
+
 # import numpy as np
 
 

From d6d12cd3e8d156fce797b4aab698628de76da447 Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Tue, 6 Dec 2022 16:30:57 +0000
Subject: [PATCH 09/15] remove tests no longer needed

---
 .../test_external_current_collector.py        | 60 ---------------
 .../test_external_temperature.py              | 74 -------------------
 2 files changed, 134 deletions(-)
 delete mode 100644 tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py

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
deleted file mode 100644
index 4b652aeff8..0000000000
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
+++ /dev/null
@@ -1,60 +0,0 @@
-#
-# Tests for external submodels
-#
-# import pybamm
-import unittest
-
-# import numpy as np
-
-
-class TestExternalCC(unittest.TestCase):
-    pass
-    # def test_2p1d(self):
-    #     # To Do: replace external variable with input
-    #     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)
-
-    #     # Simulate 100 seconds
-    #     t_eval = np.linspace(0, 100, 3)
-
-    #     for i in np.arange(1, len(t_eval) - 1):
-    #         dt = t_eval[i + 1] - 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.solution["Positive current collector potential"].data[:, -1]
-
-    #     self.assertTrue(phi_s_p.shape, (yz_pts**2, 1))
-
-
-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
index 8da30aadff..5b4b99ff5e 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
@@ -27,80 +27,6 @@ def test_input_lumped_temperature(self):
             T_av += 1
             sim.step(dt, inputs=inputs)  # works
 
-    # def test_external_temperature(self):
-
-    #     model = pybamm.lithium_ion.SPM()
-
-    #     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,
-    #     }
-
-    #     parameter_values = model.default_parameter_values
-    #     # in the default isothermal model, the temperature is everywhere equal
-    #     # to the ambient temperature
-    #     parameter_values["Ambient temperature [K]"] = pybamm.InputParameter(
-    #         "Cell temperature [K]",
-    #         domain=["negative electrode", "separator", "positive electrode"],
-    #     )
-    #     sim = pybamm.Simulation(
-    #         model, var_pts=var_pts, parameter_values=parameter_values
-    #     )
-
-    #     t_eval = np.linspace(0, 100, 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
-    #         ]
-    #         inputs = {"Cell temperature [K]": T}
-    #         sim.step(dt, inputs=inputs)  # fails because of broadcasting error
-
-    # 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.CasadiSolver()
-    #     sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver)
-    #     sim.build()
-
-    #     t_eval = np.linspace(0, 100, 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)
-
-
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys

From 917edca8dd0f59b34941aca699726ffefdc89c7a Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Tue, 6 Dec 2022 16:31:45 +0000
Subject: [PATCH 10/15] style: pre-commit fixes

---
 .../test_lithium_ion/test_external/test_external_temperature.py  | 1 +
 1 file changed, 1 insertion(+)

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 5b4b99ff5e..ca619f9c14 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
@@ -27,6 +27,7 @@ def test_input_lumped_temperature(self):
             T_av += 1
             sim.step(dt, inputs=inputs)  # works
 
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys

From d5d6bcd47f0aae289a91126e196aa63de2de7bff Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Tue, 6 Dec 2022 16:32:46 +0000
Subject: [PATCH 11/15] Remove ExternalVariable Class

---
 .../operations/convert_to_casadi.py           |   1 -
 pybamm/expression_tree/variable.py            | 104 ------------------
 2 files changed, 105 deletions(-)

diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py
index 661771a123..b3a048b1f1 100644
--- a/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -55,7 +55,6 @@ def _convert(self, symbol, t, y, y_dot, inputs):
                 pybamm.Array,
                 pybamm.Time,
                 pybamm.InputParameter,
-                pybamm.ExternalVariable,
             ),
         ):
             return casadi.MX(symbol.evaluate(t, y, y_dot, inputs))
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index 61f6ed8f55..278e535460 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -223,107 +223,3 @@ def diff(self, variable):
             raise pybamm.ModelError("cannot take second time derivative of a Variable")
         else:
             return pybamm.Scalar(0)
-
-
-class ExternalVariable(Variable):
-    """
-    A node in the expression tree representing an external variable variable.
-
-    This node will be discretised by :class:`.Discretisation` and converted
-    to a :class:`.Vector` node.
-
-    Parameters
-    ----------
-
-    name : str
-        name of the node
-    domain : iterable of str
-        list of domains that this variable is valid over
-    auxiliary_domains : dict
-        dictionary of auxiliary domains ({'secondary': ..., 'tertiary': ...,
-        'quaternary': ...}). For example, for the single particle model, the particle
-        concentration would be a Variable with domain 'negative particle' and secondary
-        auxiliary domain 'current collector'. For the DFN, the particle concentration
-        would be a Variable with domain 'negative particle', secondary domain
-        'negative electrode' and tertiary domain 'current collector'
-    domains : dict
-        A dictionary equivalent to {'primary': domain, auxiliary_domains}. Either
-        'domain' and 'auxiliary_domains', or just 'domains', should be provided
-        (not both). In future, the 'domain' and 'auxiliary_domains' arguments may be
-        deprecated.
-    scale : float or :class:`pybamm.Symbol`, optional
-        The scale of the variable, used for scaling the model when solving. The state
-        vector representing this variable will be multiplied by this scale.
-        Default is 1.
-    reference : float or :class:`pybamm.Symbol`, optional
-        The reference value of the variable, used for scaling the model when solving.
-        This value will be added to the state vector representing this variable.
-        Default is 0.
-
-    *Extends:* :class:`pybamm.Variable`
-    """
-
-    def __init__(
-        self,
-        name,
-        size,
-        domain=None,
-        auxiliary_domains=None,
-        domains=None,
-        scale=1,
-        reference=0,
-    ):
-        self._size = size
-        super().__init__(
-            name, domain, auxiliary_domains, domains, scale=scale, reference=reference
-        )
-
-    @property
-    def size(self):
-        return self._size
-
-    def create_copy(self):
-        """See :meth:`pybamm.Symbol.new_copy()`."""
-        return ExternalVariable(self.name, self.size, domains=self.domains)
-
-    def _evaluate_for_shape(self):
-        """See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()`"""
-        return np.nan * np.ones((self.size, 1))
-
-    def _base_evaluate(self, t=None, y=None, y_dot=None, inputs=None):
-        # inputs should be a dictionary
-        # convert 'None' to empty dictionary for more informative error
-        if inputs is None:
-            inputs = {}
-        if not isinstance(inputs, dict):
-            # if the special input "shape test" is passed, just return 1
-            if inputs == "shape test":
-                return self.evaluate_for_shape()
-            raise TypeError("inputs should be a dictionary")
-        try:
-            out = inputs[self.name]
-            if isinstance(out, numbers.Number) or out.shape[0] == 1:
-                return out * np.ones((self.size, 1))
-            elif out.shape[0] != self.size:
-                raise ValueError(
-                    "External variable input has size {} but should be {}".format(
-                        out.shape[0], self.size
-                    )
-                )
-            else:
-                if isinstance(out, np.ndarray) and out.ndim == 1:
-                    out = out[:, np.newaxis]
-                return out
-        # raise more informative error if can't find name in dict
-        except KeyError:
-            raise KeyError("External variable '{}' not found".format(self.name))
-
-    def diff(self, variable):
-        if variable == self:
-            return pybamm.Scalar(1)
-        elif variable == pybamm.t:
-            raise pybamm.ModelError(
-                "cannot take time derivative of an external variable"
-            )
-        else:
-            return pybamm.Scalar(0)

From a10ced0690f8b0e88159374a72533b2b55c80b48 Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Tue, 6 Dec 2022 16:40:11 +0000
Subject: [PATCH 12/15] remove numbers

---
 pybamm/expression_tree/variable.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index 278e535460..836bd81bca 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -1,7 +1,6 @@
 #
 # Variable class
 #
-import numbers
 
 import numpy as np
 import sympy

From 466a875495441150d86dab3086455a2d545a29a9 Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Tue, 13 Dec 2022 11:00:59 +0000
Subject: [PATCH 13/15] remove more tests

---
 docs/source/expression_tree/variable.rst                    | 3 ---
 tests/unit/test_expression_tree/test_d_dt.py                | 3 ---
 .../unit/test_expression_tree/test_operations/test_copy.py  | 6 ------
 3 files changed, 12 deletions(-)

diff --git a/docs/source/expression_tree/variable.rst b/docs/source/expression_tree/variable.rst
index 510e366cf1..ab42f76bcf 100644
--- a/docs/source/expression_tree/variable.rst
+++ b/docs/source/expression_tree/variable.rst
@@ -7,6 +7,3 @@ Variable
 .. autoclass:: pybamm.VariableDot
   :members:
 
-.. autoclass:: pybamm.ExternalVariable
-  :members:
-
diff --git a/tests/unit/test_expression_tree/test_d_dt.py b/tests/unit/test_expression_tree/test_d_dt.py
index ce4b1d920b..9171452406 100644
--- a/tests/unit/test_expression_tree/test_d_dt.py
+++ b/tests/unit/test_expression_tree/test_d_dt.py
@@ -40,9 +40,6 @@ def test_time_derivative_of_variable(self):
         with self.assertRaises(pybamm.ModelError):
             a = (pybamm.Variable("a")).diff(pybamm.t).diff(pybamm.t)
 
-        with self.assertRaises(pybamm.ModelError):
-            a = pybamm.ExternalVariable("a", 1).diff(pybamm.t)
-
     def test_time_derivative_of_state_vector(self):
 
         sv = pybamm.StateVector(slice(0, 10))
diff --git a/tests/unit/test_expression_tree/test_operations/test_copy.py b/tests/unit/test_expression_tree/test_operations/test_copy.py
index 684037b28a..d9558de3ae 100644
--- a/tests/unit/test_expression_tree/test_operations/test_copy.py
+++ b/tests/unit/test_expression_tree/test_operations/test_copy.py
@@ -55,12 +55,6 @@ def test_symbol_new_copy(self):
             pybamm.t,
             pybamm.Index(vec, 1),
             pybamm.NotConstant(a),
-            pybamm.ExternalVariable(
-                "external variable",
-                20,
-                domain="test",
-                auxiliary_domains={"secondary": "test2"},
-            ),
             pybamm.minimum(a, b),
             pybamm.maximum(a, b),
             pybamm.SparseStack(mat, mat),

From 3ce54390f25ea3bc900f67984eeaefa7b4899300 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Tue, 13 Dec 2022 10:55:45 -0500
Subject: [PATCH 14/15] flake8 and remove other mentions of external
 variables/submodels

---
 CHANGELOG.md                                         |  1 +
 pybamm/expression_tree/exceptions.py                 |  6 ------
 pybamm/models/base_model.py                          | 12 ++++++------
 .../models/full_battery_models/base_battery_model.py | 12 +++---------
 .../equivalent_circuit/thevenin.py                   |  6 ------
 pybamm/solvers/base_solver.py                        | 12 ++++++------
 pybamm/solvers/casadi_solver.py                      |  4 ++--
 pybamm/solvers/idaklu_solver.py                      |  2 +-
 .../test_external/test_external_temperature.py       |  4 ++--
 tests/unit/test_models/test_base_model.py            |  9 +--------
 .../test_base_battery_model.py                       |  1 -
 11 files changed, 22 insertions(+), 47 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index a6ea119383..8ef61adc8e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,7 @@
 
 ## Breaking changes
 
+- Removed external variables and submodels. InputParameter should now be used in all cases ([#2502](https://github.com/pybamm-team/PyBaMM/pull/2502))
 - Trying to use a solver to solve multiple models results in a RuntimeError exception ([#2481](https://github.com/pybamm-team/PyBaMM/pull/2481))
 - Inputs for the `ElectrodeSOH` solver are now (i) "Q_Li", the total cyclable capacity of lithium in the electrodes (previously "n_Li", the total number of moles, n_Li = 3600/F \* Q_Li) (ii) "Q_n", the capacity of the negative electrode (previously "C_n"), and "Q_p", the capacity of the positive electrode (previously "C_p") ([#2508](https://github.com/pybamm-team/PyBaMM/pull/2508))
 
diff --git a/pybamm/expression_tree/exceptions.py b/pybamm/expression_tree/exceptions.py
index 1a66d230c3..3d064df520 100644
--- a/pybamm/expression_tree/exceptions.py
+++ b/pybamm/expression_tree/exceptions.py
@@ -56,12 +56,6 @@ class ModelWarning(UserWarning):
     pass
 
 
-class InputError(Exception):
-    """An external variable has been input incorrectly into PyBaMM."""
-
-    pass
-
-
 class DiscretisationError(Exception):
     """A model could not be discretised."""
 
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index b0e5b174d1..937383a512 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -90,9 +90,9 @@ class BaseModel:
 
     def __init__(self, name="Unnamed model"):
         self.name = name
-        self._options = {"external submodels": []}
+        self._options = {}
         self._built = False
-        self._built_fundamental_and_external = False
+        self._built_fundamental = False
 
         # Initialise empty model
         self.submodels = {}
@@ -443,7 +443,7 @@ def update(self, *submodels):
             self.variables.update(submodel.variables)  # keys are strings so no check
             self._events += submodel.events
 
-    def build_fundamental_and_external(self):
+    def build_fundamental(self):
         # Get the fundamental variables
         for submodel_name, submodel in self.submodels.items():
             pybamm.logger.debug(
@@ -453,7 +453,7 @@ def build_fundamental_and_external(self):
             )
             self.variables.update(submodel.get_fundamental_variables())
 
-        self._built_fundamental_and_external = True
+        self._built_fundamental = True
 
     def build_coupled_variables(self):
         # Note: pybamm will try to get the coupled variables for the submodels in the
@@ -549,8 +549,8 @@ def _build_model(self):
 
         pybamm.logger.info("Start building {}".format(self.name))
 
-        if self._built_fundamental_and_external is False:
-            self.build_fundamental_and_external()
+        if self._built_fundamental is False:
+            self.build_fundamental()
 
         self.build_coupled_variables()
 
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 2ca62e1c05..b9d177c8ee 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -50,11 +50,6 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             * "electrolyte conductivity" : str
                 Can be "default" (default), "full", "leading order", "composite" or
                 "integrated".
-            * "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`.
             * "hydrolysis" : str
                 Whether to include hydrolysis in the model. Only implemented for
                 lead-acid models. Can be "false" (default) or "true". If "true", then
@@ -278,7 +273,6 @@ def __init__(self, extra_options):
         default_options = {
             name: options[0] for name, options in self.possible_options.items()
         }
-        default_options["external submodels"] = []
         default_options["timescale"] = "default"
 
         # Change the default for cell geometry based on which thermal option is provided
@@ -530,7 +524,7 @@ def __init__(self, extra_options):
 
         # Check options are valid
         for option, value in options.items():
-            if option in ["external submodels", "working electrode"]:
+            if option in ["working electrode"]:
                 pass
             else:
                 if isinstance(value, str) or option in [
@@ -890,7 +884,7 @@ def set_standard_output_variables(self):
                 {"y": var.y, "y [m]": var.y * L_z, "z": var.z, "z [m]": var.z * L_z}
             )
 
-    def build_fundamental_and_external(self):
+    def build_fundamental(self):
         # Get the fundamental variables
         for submodel_name, submodel in self.submodels.items():
             pybamm.logger.debug(
@@ -900,7 +894,7 @@ def build_fundamental_and_external(self):
             )
             self.variables.update(submodel.get_fundamental_variables())
 
-        self._built_fundamental_and_external = True
+        self._built_fundamental = True
 
     def build_coupled_variables(self):
         # Note: pybamm will try to get the coupled variables for the submodels in the
diff --git a/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py b/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py
index 078bb5de29..0414f3fde3 100644
--- a/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py
+++ b/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py
@@ -44,11 +44,6 @@ class Thevenin(pybamm.BaseModel):
                 - callable : if a callable is given as this option, the function \
                     defines the residual of an algebraic equation. The applied current \
                     will be solved for such that the algebraic constraint is satisfied.
-            * "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`.
     build :  bool, optional
         Whether to build the model on instantiation. Default is True. Setting this
         option to False allows users to change any number of the submodels before
@@ -89,7 +84,6 @@ def set_options(self, extra_options=None):
             "calculate discharge energy": ["false", "true"],
             "operating mode": OperatingModes("current"),
             "number of rc elements": NaturalNumberOption(1),
-            "external submodels": [[]],
         }
 
         default_options = {
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 676d5cdf53..7c4b20cd75 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -745,7 +745,7 @@ def solve(
         if (np.diff(t_eval) < 0).any():
             raise pybamm.SolverError("t_eval must increase monotonically")
 
-        # Set up external variables and inputs
+        # Set up inputs
         #
         # Argument "inputs" can be either a list of input dicts or
         # a single dict. The remaining of this function is only working
@@ -790,7 +790,7 @@ def solve(
                 )
             # It is assumed that when len(inputs_list) > 1, model set
             # up (initial condition, time-scale and length-scale) does
-            # not depend on input parameters. Thefore only `ext_and_inputs[0]`
+            # not depend on input parameters. Thefore only `model_inputs[0]`
             # is passed to `set_up`.
             # See https://github.com/pybamm-team/PyBaMM/pull/1261
             self.set_up(model, model_inputs_list[0], t_eval)
@@ -819,7 +819,7 @@ def solve(
 
         # (Re-)calculate consistent initial conditions
         # Assuming initial conditions do not depend on input parameters
-        # when len(inputs_list) > 1, only `ext_and_inputs_list[0]`
+        # when len(inputs_list) > 1, only `model_inputs_list[0]`
         # is passed to `_set_initial_conditions`.
         # See https://github.com/pybamm-team/PyBaMM/pull/1261
         if len(inputs_list) > 1:
@@ -1133,7 +1133,7 @@ def step(
         # Set timer
         timer = pybamm.Timer()
 
-        # Set up external variables and inputs
+        # Set up inputs
         model_inputs = self._set_up_model_inputs(model, inputs)
 
         t = old_solution.t[-1]
@@ -1148,7 +1148,7 @@ def step(
                     f'"{existing_model.name}". Please create a separate '
                     "solver for this model"
                 )
-            self.set_up(model, ext_and_inputs)
+            self.set_up(model, model_inputs)
             self._model_set_up.update(
                 {model: {"initial conditions": model.concatenated_initial_conditions}}
             )
@@ -1356,7 +1356,7 @@ def check_extrapolation(self, solution, events):
                     )
 
     def _set_up_model_inputs(self, model, inputs):
-        """Set up external variables and input parameters"""
+        """Set up input parameters"""
         inputs = inputs or {}
 
         # Go through all input parameters that can be found in the model
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index eea853eeda..12627b7120 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -143,7 +143,7 @@ def _integrate(self, model, t_eval, inputs_dict=None):
         t_eval : numeric type
             The times at which to compute the solution
         inputs_dict : dict, optional
-            Any external variables or input parameters to pass to the model when solving
+            Any input parameters to pass to the model when solving
         """
 
         # Record whether there are any symbolic inputs
@@ -605,7 +605,7 @@ def _run_integrator(
         y0:
             casadi vector of initial conditions
         inputs_dict : dict, optional
-            Any external variables or input parameters to pass to the model when solving
+            Any input parameters to pass to the model when solving
         inputs:
             Casadi vector of inputs
         t_eval : numeric type
diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py
index 8a7671f084..1c2822d968 100644
--- a/pybamm/solvers/idaklu_solver.py
+++ b/pybamm/solvers/idaklu_solver.py
@@ -510,7 +510,7 @@ def _integrate(self, model, t_eval, inputs_dict=None):
         t_eval : numeric type
             The times at which to compute the solution
         inputs_dict : dict, optional
-            Any external variables or input parameters to pass to the model when solving
+            Any input parameters to pass to the model when solving
         """
         inputs_dict = inputs_dict or {}
         # stack inputs
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 ca619f9c14..9a1c8ff732 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
@@ -1,12 +1,12 @@
 #
-# Tests for external submodels
+# Tests for inputting a temperature profile
 #
 import pybamm
 import unittest
 import numpy as np
 
 
-class TestExternalThermalModels(unittest.TestCase):
+class TestInputLumpedTemperature(unittest.TestCase):
     def test_input_lumped_temperature(self):
         model = pybamm.lithium_ion.SPMe()
         parameter_values = model.default_parameter_values
diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py
index 161f347a26..ea593642bd 100644
--- a/tests/unit/test_models/test_base_model.py
+++ b/tests/unit/test_models/test_base_model.py
@@ -531,15 +531,8 @@ def test_export_casadi(self):
         np.testing.assert_array_equal(np.array(jac_alg_fn(5, 6, 7, [9, 8])), [[1, -1]])
         self.assertEqual(var_fn(6, 3, 2, [2, 7]), -1)
 
-        # Test model with external variable runs
-        model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
-        model = pybamm.lithium_ion.SPMe(model_options)
-        sim = pybamm.Simulation(model)
-        sim.build()
-        variable_names = ["Volume-averaged cell temperature"]
-        out = sim.built_model.export_casadi_objects(variable_names)
-
         # Test fails if not discretised
+        model = pybamm.lithium_ion.SPMe()
         with self.assertRaisesRegex(
             pybamm.DiscretisationError, "Cannot automatically discretise model"
         ):
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index b2fa082f7e..6ef7a02760 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -44,7 +44,6 @@
 'total interfacial current density as a state': 'false' (possible: ['false', 'true'])
 'working electrode': 'both' (possible: ['both', 'negative', 'positive'])
 'x-average side reactions': 'false' (possible: ['false', 'true'])
-'external submodels': []
 'timescale': 'default'
 """  # noqa: E501
 

From bff1f468f6c125ceb0786ad64d156d7f204544ae Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Tue, 13 Dec 2022 11:59:45 -0500
Subject: [PATCH 15/15] coverage and error message for wrong order

---
 pybamm/models/base_model.py               | 17 +++++++++++----
 pybamm/simulation.py                      |  2 +-
 tests/unit/test_models/test_base_model.py | 25 +++++++++++++++++++++++
 3 files changed, 39 insertions(+), 5 deletions(-)

diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 937383a512..38b957f40f 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -916,8 +916,10 @@ def export_casadi_objects(self, variable_names, input_parameter_order=None):
         variable_names : list
             Variables to be exported alongside the model structure
         input_parameter_order : list, optional
-            Order in which the input parameters should be stacked. If None, the order
-            returned by :meth:`BaseModel.input_parameters` is used
+            Order in which the input parameters should be stacked.
+            If input_parameter_order=None and len(self.input_parameters) > 1, a
+            ValueError is raised (this helps to avoid accidentally using the wrong
+            order)
 
         Returns
         -------
@@ -940,6 +942,11 @@ def export_casadi_objects(self, variable_names, input_parameter_order=None):
             inputs_wrong_order[name] = casadi.MX.sym(name, input_param._expected_size)
         # Sort according to input_parameter_order
         if input_parameter_order is None:
+            if len(inputs_wrong_order) > 1:
+                raise ValueError(
+                    "input_parameter_order must be specified if there is more than one "
+                    "input parameter"
+                )
             inputs = inputs_wrong_order
         else:
             inputs = {name: inputs_wrong_order[name] for name in input_parameter_order}
@@ -996,8 +1003,10 @@ def generate(
         variable_names : list
             Variables to be exported alongside the model structure
         input_parameter_order : list, optional
-            Order in which the input parameters should be stacked. If None, the order
-            returned by :meth:`BaseModel.input_parameters` is used
+            Order in which the input parameters should be stacked.
+            If input_parameter_order=None and len(self.input_parameters) > 1, a
+            ValueError is raised (this helps to avoid accidentally using the wrong
+            order)
         cg_options : dict
             Options to pass to the code generator.
             See https://web.casadi.org/docs/#generating-c-code
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 972debf03f..b41b48581b 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -108,7 +108,7 @@ def __init__(
         self.submesh_types = submesh_types or self.model.default_submesh_types
         self.var_pts = var_pts or self.model.default_var_pts
         self.spatial_methods = spatial_methods or self.model.default_spatial_methods
-        self._solver = solver or self.model.default_solver
+        self.solver = solver or self.model.default_solver
         self.output_variables = output_variables
 
         # Initialize empty built states
diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py
index ea593642bd..1d00cca1b4 100644
--- a/tests/unit/test_models/test_base_model.py
+++ b/tests/unit/test_models/test_base_model.py
@@ -531,6 +531,31 @@ def test_export_casadi(self):
         np.testing.assert_array_equal(np.array(jac_alg_fn(5, 6, 7, [9, 8])), [[1, -1]])
         self.assertEqual(var_fn(6, 3, 2, [2, 7]), -1)
 
+        # Test fails if order not specified
+        with self.assertRaisesRegex(
+            ValueError, "input_parameter_order must be specified"
+        ):
+            model.export_casadi_objects(["a+b"])
+
+        # Fine if order is not specified if there is only one input parameter
+        model = pybamm.BaseModel()
+        p = pybamm.InputParameter("p")
+        model.rhs = {a: -a}
+        model.algebraic = {b: a - b}
+        model.initial_conditions = {a: 1, b: 1}
+        model.variables = {"a+b": a + b - p}
+
+        out = model.export_casadi_objects(["a+b"])
+
+        # Try making a function from the outputs
+        t, x, z, p = out["t"], out["x"], out["z"], out["inputs"]
+        var = out["variables"]["a+b"]
+        var_fn = casadi.Function("var", [t, x, z, p], [var])
+
+        # Test that function values are as expected
+        # a + b - p = 3 + 2 - 7 = -2
+        self.assertEqual(var_fn(6, 3, 2, [7]), -2)
+
         # Test fails if not discretised
         model = pybamm.lithium_ion.SPMe()
         with self.assertRaisesRegex(