diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4d7f16265..adb52c99f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -16,20 +16,21 @@
 
 ## Bug fixes
 
--   Fixed `Interpolant` ids to allow processing ([#962](https://github.com/pybamm-team/PyBaMM/pull/962)
+-   Added extra checks when creating a model, for clearer errors ([#971](https://github.com/pybamm-team/PyBaMM/pull/971))
+-   Fixed `Interpolant` ids to allow processing ([#962](https://github.com/pybamm-team/PyBaMM/pull/962))
 -   Fixed a bug in the initial conditions of the potential pair model ([#954](https://github.com/pybamm-team/PyBaMM/pull/954))
--   Changed simulation attributes to assign copies rather than the objects themselves ([#952](https://github.com/pybamm-team/PyBaMM/pull/952)
--   Added default values to base model so that it works with the `Simulation` class ([#952](https://github.com/pybamm-team/PyBaMM/pull/952)
--   Fixed solver to recompute initial conditions when inputs are changed ([#951](https://github.com/pybamm-team/PyBaMM/pull/951)
--   Reformatted thermal submodels ([#938](https://github.com/pybamm-team/PyBaMM/pull/938)
+-   Changed simulation attributes to assign copies rather than the objects themselves ([#952](https://github.com/pybamm-team/PyBaMM/pull/952))
+-   Added default values to base model so that it works with the `Simulation` class ([#952](https://github.com/pybamm-team/PyBaMM/pull/952))
+-   Fixed solver to recompute initial conditions when inputs are changed ([#951](https://github.com/pybamm-team/PyBaMM/pull/951))
+-   Reformatted thermal submodels ([#938](https://github.com/pybamm-team/PyBaMM/pull/938))
 -   Reformatted electrolyte submodels ([#927](https://github.com/pybamm-team/PyBaMM/pull/927))
 -   Reformatted convection submodels ([#635](https://github.com/pybamm-team/PyBaMM/pull/635))
 
 ## Breaking changes
 
 -   Removed some inputs like `T_inf`, `R_g` and activation energies to some of the standard function parameters. This is because each of those inputs is specific to a particular function (e.g. the reference temperature at which the function was measured). To change a property such as the activation energy, users should create a new function, specifying the relevant property as a `Parameter` or `InputParameter` ([#942](https://github.com/pybamm-team/PyBaMM/pull/942))
--   The thermal option 'xyz-lumped' has been removed. The option 'thermal current collector' has also been removed ([#938](https://github.com/pybamm-team/PyBaMM/pull/938)
--   The 'C-rate' parameter has been deprecated. Use 'Current function [A]' instead. The cell capacity can be accessed as 'Cell capacity [A.h]', and used to calculate current from C-rate ([#952](https://github.com/pybamm-team/PyBaMM/pull/952)
+-   The thermal option 'xyz-lumped' has been removed. The option 'thermal current collector' has also been removed ([#938](https://github.com/pybamm-team/PyBaMM/pull/938))
+-   The 'C-rate' parameter has been deprecated. Use 'Current function [A]' instead. The cell capacity can be accessed as 'Cell capacity [A.h]', and used to calculate current from C-rate ([#952](https://github.com/pybamm-team/PyBaMM/pull/952))
 
 # [v0.2.1](https://github.com/pybamm-team/PyBaMM/tree/v0.2.1) - 2020-03-31
 
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 5549b121d..e9cae42dd 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -117,6 +117,12 @@ class Discretisation(object):
             `model.variables = {}`)
 
         """
+        if model.is_discretised is True:
+            raise pybamm.ModelError(
+                "Cannot re-discretise a model. "
+                "Set 'inplace=False' when first discretising a model to then be able "
+                "to discretise it more times (e.g. for convergence studies)."
+            )
 
         pybamm.logger.info("Start discretising {}".format(model.name))
 
diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index 1d5b3580a..32bc0ca4c 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -169,6 +169,10 @@ class Function(pybamm.Symbol):
             ]
             return self._function_evaluate(evaluated_children)
 
+    def evaluates_on_edges(self):
+        """ See :meth:`pybamm.Symbol.evaluates_on_edges()`. """
+        return any(child.evaluates_on_edges() for child in self.children)
+
     def _evaluate_for_shape(self):
         """
         Default behaviour: has same shape as all child
diff --git a/pybamm/expression_tree/operations/simplify.py b/pybamm/expression_tree/operations/simplify.py
index df033cea4..60282c893 100644
--- a/pybamm/expression_tree/operations/simplify.py
+++ b/pybamm/expression_tree/operations/simplify.py
@@ -567,7 +567,7 @@ class Simplification(object):
     def __init__(self, simplified_symbols=None):
         self._simplified_symbols = simplified_symbols or {}
 
-    def simplify(self, symbol):
+    def simplify(self, symbol, clear_domains=True):
         """
         This function recurses down the tree, applying any simplifications defined in
         classes derived from pybamm.Symbol. E.g. any expression multiplied by a
@@ -577,7 +577,9 @@ class Simplification(object):
         Parameters
         ----------
         symbol : :class:`pybamm.Symbol`
-        The symbol to simplify
+            The symbol to simplify
+        clear_domains : bool
+            Whether to remove a symbol's domain when simplifying. Default is True.
 
         Returns
         -------
@@ -588,15 +590,16 @@ class Simplification(object):
         try:
             return self._simplified_symbols[symbol.id]
         except KeyError:
-            simplified_symbol = self._simplify(symbol)
+            simplified_symbol = self._simplify(symbol, clear_domains)
 
             self._simplified_symbols[symbol.id] = simplified_symbol
 
             return simplified_symbol
 
-    def _simplify(self, symbol):
+    def _simplify(self, symbol, clear_domains=True):
         """ See :meth:`Simplification.simplify()`. """
-        symbol.clear_domains()
+        if clear_domains:
+            symbol.clear_domains()
 
         if isinstance(symbol, pybamm.BinaryOperator):
             left, right = symbol.children
@@ -607,7 +610,11 @@ class Simplification(object):
             new_symbol = symbol._binary_simplify(new_left, new_right)
 
         elif isinstance(symbol, pybamm.UnaryOperator):
-            new_child = self.simplify(symbol.child)
+            # Reassign domain for gradient and divergence
+            if isinstance(symbol, (pybamm.Gradient, pybamm.Divergence)):
+                new_child = self.simplify(symbol.child, clear_domains=False)
+            else:
+                new_child = self.simplify(symbol.child)
             # _unary_simplify defined in derived classes for specific rules
             new_symbol = symbol._unary_simplify(new_child)
 
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index bc30e48cf..63f10e618 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -498,10 +498,16 @@ class Symbol(anytree.NodeMixin):
 
     def jac(self, variable, known_jacs=None, clear_domain=True):
         """
-        Differentiate a symbol with respect to a (slice of) a State Vector.
+        Differentiate a symbol with respect to a (slice of) a StateVector
+        or StateVectorDot.
         See :class:`pybamm.Jacobian`.
         """
         jac = pybamm.Jacobian(known_jacs, clear_domain=clear_domain)
+        if not isinstance(variable, (pybamm.StateVector, pybamm.StateVectorDot)):
+            raise TypeError(
+                "Jacobian can only be taken with respect to a 'StateVector' "
+                "or 'StateVectorDot', but {} is a {}".format(variable, type(variable))
+            )
         return jac.jac(self, variable)
 
     def _jac(self, variable):
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index de123d836..b1d785c39 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -310,6 +310,16 @@ class Gradient(SpatialOperator):
     """
 
     def __init__(self, child):
+        if child.domain == []:
+            raise pybamm.DomainError(
+                "Cannot take gradient of '{}' since its domain is empty. ".format(child)
+                + "Try broadcasting the object first, e.g.\n\n"
+                "\tpybamm.grad(pybamm.PrimaryBroadcast(symbol, 'domain'))"
+            )
+        if child.evaluates_on_edges() is True:
+            raise TypeError(
+                "Cannot take gradient of '{}' since it evaluates on edges".format(child)
+            )
         super().__init__("grad", child)
 
     def evaluates_on_edges(self):
@@ -324,6 +334,20 @@ class Divergence(SpatialOperator):
     """
 
     def __init__(self, child):
+        if child.domain == []:
+            raise pybamm.DomainError(
+                "Cannot take divergence of '{}' since its domain is empty. ".format(
+                    child
+                )
+                + "Try broadcasting the object first, e.g.\n\n"
+                "\tpybamm.div(pybamm.PrimaryBroadcast(symbol, 'domain'))"
+            )
+        if child.evaluates_on_edges() is False:
+            raise TypeError(
+                "Cannot take divergence of '{}' since it does not ".format(child)
+                + "evaluates on nodes. Usually, a gradient should be taken before the "
+                "divergence."
+            )
         super().__init__("div", child)
 
     def evaluates_on_edges(self):
@@ -814,8 +838,12 @@ def grad(expression):
     :class:`Gradient`
         the gradient of ``expression``
     """
-
-    return Gradient(expression)
+    # Gradient of a broadcast is zero
+    if isinstance(expression, pybamm.PrimaryBroadcast):
+        new_child = pybamm.PrimaryBroadcast(0, expression.child.domain)
+        return pybamm.PrimaryBroadcastToEdges(new_child, expression.domain)
+    else:
+        return Gradient(expression)
 
 
 def div(expression):
@@ -833,8 +861,12 @@ def div(expression):
     :class:`Divergence`
         the divergence of ``expression``
     """
-
-    return Divergence(expression)
+    # Divergence of a broadcast is zero
+    if isinstance(expression, pybamm.PrimaryBroadcastToEdges):
+        new_child = pybamm.PrimaryBroadcast(0, expression.child.domain)
+        return pybamm.PrimaryBroadcast(new_child, expression.domain)
+    else:
+        return Divergence(expression)
 
 
 def laplacian(expression):
diff --git a/pybamm/input/parameters/lithium-ion/experiments/1C_charge_from_empty_Mohtat2020/parameters.csv b/pybamm/input/parameters/lithium-ion/experiments/1C_charge_from_empty_Mohtat2020/parameters.csv
index ae32f40bb..05800bb53 100644
--- a/pybamm/input/parameters/lithium-ion/experiments/1C_charge_from_empty_Mohtat2020/parameters.csv
+++ b/pybamm/input/parameters/lithium-ion/experiments/1C_charge_from_empty_Mohtat2020/parameters.csv
@@ -10,7 +10,7 @@ Ambient temperature [K], 298.15,,
 Number of electrodes connected in parallel to make a cell,1,,
 Number of cells connected in series to make a battery,1,,
 Lower voltage cut-off [V],2.5,,
-Upper voltage cut-off [V],4.2,,
+Upper voltage cut-off [V],4.25,,just above 4.2
 ,,,
 # Initial conditions
 Initial concentration in negative electrode [mol.m-3],48.8682,Peyman MPM, x0 (0.0017) * Csmax_n
diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py
index 7cfa06ad2..5c4a6d9cd 100644
--- a/pybamm/meshes/meshes.py
+++ b/pybamm/meshes/meshes.py
@@ -142,6 +142,8 @@ class Mesh(dict):
         submesh: :class:`self.submeshclass`
             A new submesh with the class defined by self.submeshclass
         """
+        if submeshnames == ():
+            raise ValueError("Submesh domains being combined cannot be empty")
         # Check that the final edge of each submesh is the same as the first edge of the
         # next submesh
         for i in range(len(submeshnames) - 1):
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 131c5acd3..2a16c024a 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -101,12 +101,12 @@ class BaseModel(object):
         self.options = {}
 
         # Initialise empty model
-        self._rhs = {}
-        self._algebraic = {}
-        self._initial_conditions = {}
-        self._boundary_conditions = {}
-        self._variables = pybamm.FuzzyDict()
-        self._events = []
+        self.rhs = {}
+        self.algebraic = {}
+        self.initial_conditions = {}
+        self.boundary_conditions = {}
+        self.variables = {}
+        self.events = []
         self._concatenated_rhs = None
         self._concatenated_algebraic = None
         self._concatenated_initial_conditions = None
@@ -128,30 +128,6 @@ class BaseModel(object):
         # Default timescale is 1 second
         self.timescale = pybamm.Scalar(1)
 
-    def _set_dictionary(self, dict, name):
-        """
-        Convert any scalar equations in dict to 'pybamm.Scalar'
-        and check that domains are consistent
-        """
-        # Convert any numbers to a pybamm.Scalar
-        for var, eqn in dict.items():
-            if isinstance(eqn, numbers.Number):
-                dict[var] = pybamm.Scalar(eqn)
-
-        if not all(
-            [
-                variable.domain == equation.domain
-                or variable.domain == []
-                or equation.domain == []
-                for variable, equation in dict.items()
-            ]
-        ):
-            raise pybamm.DomainError(
-                "variable and equation in '{}' must have the same domain".format(name)
-            )
-
-        return dict
-
     @property
     def name(self):
         return self._name
@@ -166,7 +142,7 @@ class BaseModel(object):
 
     @rhs.setter
     def rhs(self, rhs):
-        self._rhs = self._set_dictionary(rhs, "rhs")
+        self._rhs = EquationDict("rhs", rhs)
 
     @property
     def algebraic(self):
@@ -174,7 +150,7 @@ class BaseModel(object):
 
     @algebraic.setter
     def algebraic(self, algebraic):
-        self._algebraic = self._set_dictionary(algebraic, "algebraic")
+        self._algebraic = EquationDict("algebraic", algebraic)
 
     @property
     def initial_conditions(self):
@@ -182,8 +158,8 @@ class BaseModel(object):
 
     @initial_conditions.setter
     def initial_conditions(self, initial_conditions):
-        self._initial_conditions = self._set_dictionary(
-            initial_conditions, "initial_conditions"
+        self._initial_conditions = EquationDict(
+            "initial_conditions", initial_conditions
         )
 
     @property
@@ -192,23 +168,7 @@ class BaseModel(object):
 
     @boundary_conditions.setter
     def boundary_conditions(self, boundary_conditions):
-        # Convert any numbers to a pybamm.Scalar
-        for var, bcs in boundary_conditions.items():
-            for side, bc in bcs.items():
-                if isinstance(bc[0], numbers.Number):
-                    # typ is the type of the bc, e.g. "Dirichlet" or "Neumann"
-                    eqn, typ = boundary_conditions[var][side]
-                    boundary_conditions[var][side] = (pybamm.Scalar(eqn), typ)
-                # Check types
-                if bc[1] not in ["Dirichlet", "Neumann"]:
-                    raise pybamm.ModelError(
-                        """
-                        boundary condition types must be Dirichlet or Neumann, not '{}'
-                        """.format(
-                            bc[1]
-                        )
-                    )
-        self._boundary_conditions = boundary_conditions
+        self._boundary_conditions = BoundaryConditionsDict(boundary_conditions)
 
     @property
     def variables(self):
@@ -299,9 +259,12 @@ class BaseModel(object):
 
     @param.setter
     def param(self, values):
-        # convert module into a class
-        # (StackOverflow: https://tinyurl.com/yk3euon3)
-        self._param = ParamClass(values)
+        if values is None:
+            self._param = None
+        else:
+            # convert module into a class
+            # (StackOverflow: https://tinyurl.com/yk3euon3)
+            self._param = ParamClass(values)
 
     @property
     def options(self):
@@ -723,3 +686,95 @@ def find_symbol_in_model(model, name):
         dic_return = find_symbol_in_dict(dic, name)
         if dic_return:
             return dic_return
+
+
+class EquationDict(dict):
+    def __init__(self, name, equations):
+        self.name = name
+        equations = self.check_and_convert_equations(equations)
+        super().__init__(equations)
+
+    def __setitem__(self, key, value):
+        "Call the update functionality when doing a setitem"
+        self.update({key: value})
+
+    def update(self, equations):
+        equations = self.check_and_convert_equations(equations)
+        super().update(equations)
+
+    def check_and_convert_equations(self, equations):
+        """
+        Convert any scalar equations in dict to 'pybamm.Scalar'
+        and check that domains are consistent
+        """
+        # Convert any numbers to a pybamm.Scalar
+        for var, eqn in equations.items():
+            if isinstance(eqn, numbers.Number):
+                equations[var] = pybamm.Scalar(eqn)
+
+        if not all(
+            [
+                variable.domain == equation.domain
+                or variable.domain == []
+                or equation.domain == []
+                for variable, equation in equations.items()
+            ]
+        ):
+            raise pybamm.DomainError(
+                "variable and equation in '{}' must have the same domain".format(
+                    self.name
+                )
+            )
+
+        # For initial conditions, check that the equation doesn't contain any
+        # Variable objects
+        # skip this if the dictionary has no "name" attribute (which will be the case
+        # after pickling)
+        if hasattr(self, "name") and self.name == "initial_conditions":
+            for var, eqn in equations.items():
+                if eqn.has_symbol_of_classes(pybamm.Variable):
+                    unpacker = pybamm.SymbolUnpacker(pybamm.Variable)
+                    variable_in_equation = list(unpacker.unpack_symbol(eqn).values())[0]
+                    raise TypeError(
+                        "Initial conditions cannot contain 'Variable' objects, "
+                        "but '{!r}' found in initial conditions for '{}'".format(
+                            variable_in_equation, var
+                        )
+                    )
+
+        return equations
+
+
+class BoundaryConditionsDict(dict):
+    def __init__(self, bcs):
+        bcs = self.check_and_convert_bcs(bcs)
+        super().__init__(bcs)
+
+    def __setitem__(self, key, value):
+        "Call the update functionality when doing a setitem"
+        self.update({key: value})
+
+    def update(self, bcs):
+        bcs = self.check_and_convert_bcs(bcs)
+        super().update(bcs)
+
+    def check_and_convert_bcs(self, boundary_conditions):
+        """ Convert any scalar bcs in dict to 'pybamm.Scalar', and check types """
+        # Convert any numbers to a pybamm.Scalar
+        for var, bcs in boundary_conditions.items():
+            for side, bc in bcs.items():
+                if isinstance(bc[0], numbers.Number):
+                    # typ is the type of the bc, e.g. "Dirichlet" or "Neumann"
+                    eqn, typ = boundary_conditions[var][side]
+                    boundary_conditions[var][side] = (pybamm.Scalar(eqn), typ)
+                # Check types
+                if bc[1] not in ["Dirichlet", "Neumann"]:
+                    raise pybamm.ModelError(
+                        """
+                        boundary condition types must be Dirichlet or Neumann, not '{}'
+                        """.format(
+                            bc[1]
+                        )
+                    )
+
+        return boundary_conditions
diff --git a/pybamm/models/submodels/base_submodel.py b/pybamm/models/submodels/base_submodel.py
index 335fdb8e9..73d2c2965 100644
--- a/pybamm/models/submodels/base_submodel.py
+++ b/pybamm/models/submodels/base_submodel.py
@@ -4,7 +4,7 @@
 import pybamm
 
 
-class BaseSubModel:
+class BaseSubModel(pybamm.BaseModel):
     """
     The base class for all submodels. All submodels inherit from this class and must
     only provide public methods which overwrite those in this base class. Any methods
@@ -54,16 +54,8 @@ class BaseSubModel:
         name="Unnamed submodel",
         external=False,
     ):
-        super().__init__()
+        super().__init__(name)
         self.param = param
-        # Initialise empty variables (to avoid overwriting with 'None')
-
-        self.rhs = {}
-        self.algebraic = {}
-        self.boundary_conditions = {}
-        self.initial_conditions = {}
-        self.variables = {}
-        self.events = []
 
         self.domain = domain
         self.set_domain_for_broadcast()
diff --git a/pybamm/models/submodels/convection/base_convection.py b/pybamm/models/submodels/convection/base_convection.py
index 65f123093..95828ab11 100644
--- a/pybamm/models/submodels/convection/base_convection.py
+++ b/pybamm/models/submodels/convection/base_convection.py
@@ -110,49 +110,3 @@ class BaseModel(pybamm.BaseSubModel):
         variables = {"Pressure": p}
 
         return variables
-
-    def _separator_velocity(self, variables):
-        """
-        A private method to calculate x- and z-components of velocity in the separator
-        Parameters
-        ----------
-        variables : dict
-            Dictionary of variables in the whole model.
-
-        Returns
-        -------
-        v_box_s : :class:`pybamm.Symbol`
-            The x-component of velocity in the separator
-        dVbox_dz : :class:`pybamm.Symbol`
-            The z-component of velocity in the separator
-        """
-        # Set up
-        param = self.param
-        l_n = pybamm.geometric_parameters.l_n
-        l_s = pybamm.geometric_parameters.l_s
-        x_s = pybamm.standard_spatial_vars.x_s
-
-        # Difference in negative and positive electrode velocities determines the
-        # velocity in the separator
-        i_boundary_cc = variables["Current collector current density"]
-        v_box_n_right = param.beta_n * i_boundary_cc
-        v_box_p_left = param.beta_p * i_boundary_cc
-        d_vbox_s__dx = (v_box_p_left - v_box_n_right) / l_s
-
-        # Simple formula for velocity in the separator
-        dVbox_dz = pybamm.Concatenation(
-            pybamm.FullBroadcast(
-                0,
-                "negative electrode",
-                auxiliary_domains={"secondary": "current collector"},
-            ),
-            pybamm.PrimaryBroadcast(-d_vbox_s__dx, "separator"),
-            pybamm.FullBroadcast(
-                0,
-                "positive electrode",
-                auxiliary_domains={"secondary": "current collector"},
-            ),
-        )
-        v_box_s = d_vbox_s__dx * (x_s - l_n) + v_box_n_right
-
-        return v_box_s, dVbox_dz
diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py
index 98ed2a050..e9b2e2e73 100644
--- a/pybamm/models/submodels/particle/base_particle.py
+++ b/pybamm/models/submodels/particle/base_particle.py
@@ -72,7 +72,7 @@ class BaseParticle(pybamm.BaseSubModel):
 
     def set_events(self, variables):
         c_s_surf = variables[self.domain + " particle surface concentration"]
-        tol = 0.01
+        tol = 1e-4
 
         self.events.append(
             pybamm.Event(
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 7e7668d65..e09d99fc9 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -493,9 +493,21 @@ class ParameterValues:
             elif isinstance(function_name, pybamm.InputParameter):
                 # Replace the function with an input parameter
                 function = function_name
-            else:
+            elif (
+                isinstance(function_name, pybamm.Symbol)
+                and function_name.evaluates_to_number()
+            ):
+                # If the "function" provided is a pybamm scalar-like, use ones_like to
+                # get the right shape
+                function = function_name * pybamm.ones_like(*new_children)
+            elif callable(function_name):
                 # otherwise evaluate the function to create a new PyBaMM object
                 function = function_name(*new_children)
+            else:
+                raise TypeError(
+                    "Parameter provided for '{}' ".format(symbol.name)
+                    + "is of the wrong type (should either be scalar-like or callable)"
+                )
             # Differentiate if necessary
             if symbol.diff_variable is None:
                 function_out = function
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 0b88a125c..9b9f58a08 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -446,14 +446,14 @@ class Simulation:
                     or "[experiment]" in self._solution.termination
                 ):
                     pybamm.logger.warning(
-                        """
-                        Experiment is infeasible: '{}' was triggered during '{}'. Try
-                        reducing current, shortening the time interval, or reducing
-                        the period.
-                        """.format(
+                        "\n\n\tExperiment is infeasible: '{}' ".format(
                             self._solution.termination,
+                        )
+                        + "was triggered during '{}'. ".format(
                             self.experiment.operating_conditions_strings[idx],
                         )
+                        + "Try reducing current, shortening the time interval, "
+                        "or reducing the period.\n\n"
                     )
                     break
             pybamm.logger.info(

Better error messages
Some small changes / better error messages needed for increased usability:
- Raise error when taking gradient of an object with no domain (but gradient of broadcast should work and return zero)
- Raise error when trying to do boundary value of a gradient object
- ~~Allow re-processing discretised model~~ Raise error when trying to re-process discretised model, and recommend using `inplace=False`
diff --git a/tests/integration/test_spatial_methods/test_finite_volume.py b/tests/integration/test_spatial_methods/test_finite_volume.py
index 5f46b4d59..75505c5af 100644
--- a/tests/integration/test_spatial_methods/test_finite_volume.py
+++ b/tests/integration/test_spatial_methods/test_finite_volume.py
@@ -9,6 +9,23 @@ import unittest
 
 
 class TestFiniteVolumeConvergence(unittest.TestCase):
+    def test_grad_div_broadcast(self):
+        # create mesh and discretisation
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        mesh = get_mesh_for_testing()
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        a = pybamm.PrimaryBroadcast(1, "negative electrode")
+        grad_a = disc.process_symbol(pybamm.grad(a))
+        np.testing.assert_array_equal(grad_a.evaluate(), 0)
+
+        a_edge = pybamm.PrimaryBroadcastToEdges(1, "negative electrode")
+        div_a = disc.process_symbol(pybamm.div(a_edge))
+        np.testing.assert_array_equal(div_a.evaluate(), 0)
+
+        div_grad_a = disc.process_symbol(pybamm.div(pybamm.grad(a)))
+        np.testing.assert_array_equal(div_grad_a.evaluate(), 0)
+
     def test_cartesian_spherical_grad_convergence(self):
         # note that grad function is the same for cartesian and spherical
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
@@ -71,10 +88,10 @@ class TestFiniteVolumeConvergence(unittest.TestCase):
             disc = pybamm.Discretisation(mesh, spatial_methods)
             combined_submesh = mesh.combine_submeshes(*whole_cell)
             x = combined_submesh[0].nodes
-            x_edge = combined_submesh[0].edges
+            x_edge = pybamm.standard_spatial_vars.x_edge
 
             # Define flux and bcs
-            N = pybamm.Vector(x_edge ** 2 * np.cos(x_edge), domain=whole_cell)
+            N = x_edge ** 2 * pybamm.cos(x_edge)
             div_eqn = pybamm.div(N)
             # Define exact solutions
             # N = x**2 * cos(x) --> dNdx = x*(2cos(x) - xsin(x))
@@ -106,12 +123,10 @@ class TestFiniteVolumeConvergence(unittest.TestCase):
             disc = pybamm.Discretisation(mesh, spatial_methods)
             submesh = mesh["negative particle"]
             r = submesh[0].nodes
-            r_edge = submesh[0].edges
+            r_edge = pybamm.standard_spatial_vars.r_n_edge
 
             # Define flux and bcs
-            N = pybamm.Vector(
-                r_edge ** 2 * np.sin(r_edge), domain=["negative particle"]
-            )
+            N = r_edge ** 2 * pybamm.sin(r_edge)
             div_eqn = pybamm.div(N)
             # Define exact solutions
             # N = r**3 --> div(N) = 5 * r**2
@@ -143,10 +158,10 @@ class TestFiniteVolumeConvergence(unittest.TestCase):
             disc = pybamm.Discretisation(mesh, spatial_methods)
             submesh = mesh["negative particle"]
             r = submesh[0].nodes
-            r_edge = submesh[0].edges
+            r_edge = pybamm.standard_spatial_vars.r_n_edge
 
             # Define flux and bcs
-            N = pybamm.Vector(r_edge * np.sin(r_edge), domain=["negative particle"])
+            N = r_edge * pybamm.sin(r_edge)
             div_eqn = pybamm.div(N)
             # Define exact solutions
             # N = r*sin(r) --> div(N) = 3*sin(r) + r*cos(r)
@@ -178,12 +193,9 @@ class TestFiniteVolumeConvergence(unittest.TestCase):
             disc = pybamm.Discretisation(mesh, spatial_methods)
             submesh = mesh["negative particle"]
             r = submesh[0].nodes
-            r_edge = submesh[0].edges
+            r_edge = pybamm.standard_spatial_vars.r_n_edge
 
-            N = pybamm.Matrix(
-                np.kron(np.ones(len(submesh)), r_edge ** 2 * np.sin(r_edge)),
-                domain=["negative particle"],
-            )
+            N = r_edge ** 2 * pybamm.sin(r_edge)
             div_eqn = pybamm.div(N)
             # Define exact solutions
             # N = r**2*sin(r) --> div(N) = 4*r*sin(r) - r**2*cos(r)
@@ -206,7 +218,10 @@ class TestFiniteVolumeConvergence(unittest.TestCase):
 
     def test_p2d_with_x_dep_bcs_spherical_convergence(self):
         # test div_r( (r**2 * sin(r)) * x ) == (4*r*sin(r) - r**2*cos(r)) * x
-        spatial_methods = {"negative particle": pybamm.FiniteVolume()}
+        spatial_methods = {
+            "negative particle": pybamm.FiniteVolume(),
+            "negative electrode": pybamm.FiniteVolume(),
+        }
 
         # Function for convergence testing
         def get_error(m):
@@ -215,18 +230,17 @@ class TestFiniteVolumeConvergence(unittest.TestCase):
             disc = pybamm.Discretisation(mesh, spatial_methods)
             submesh_r = mesh["negative particle"]
             r = submesh_r[0].nodes
-            r_edge = submesh_r[0].edges
-            x = pybamm.Vector(mesh["negative electrode"][0].nodes)
+            r_edge = pybamm.standard_spatial_vars.r_n_edge
+            x = pybamm.standard_spatial_vars.x_n
 
-            N = pybamm.Matrix(
-                np.kron(x.entries[:, 0], r_edge ** 2 * np.sin(r_edge)),
-                domain=["negative particle"],
+            N = pybamm.PrimaryBroadcast(x, "negative particle") * (
+                r_edge ** 2 * pybamm.sin(r_edge)
             )
             div_eqn = pybamm.div(N)
             # Define exact solutions
             # N = r**2*sin(r) --> div(N) = 4*r*sin(r) - r**2*cos(r)
             div_exact = 4 * r * np.sin(r) + r ** 2 * np.cos(r)
-            div_exact = np.kron(x.entries[:, 0], div_exact)
+            div_exact = np.kron(mesh["negative electrode"][0].nodes, div_exact)
 
             # Discretise and evaluate
             div_eqn_disc = disc.process_symbol(div_eqn)
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 537d0d9a6..06fb95811 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -475,12 +475,11 @@ class TestDiscretise(unittest.TestCase):
         disc.set_variable_slices(variables)
 
         # Simple expressions
-        for eqn in [pybamm.grad(var), pybamm.div(var)]:
+        for eqn in [pybamm.grad(var), pybamm.div(pybamm.grad(var))]:
             eqn_disc = disc.process_symbol(eqn)
 
             self.assertIsInstance(eqn_disc, pybamm.MatrixMultiplication)
             self.assertIsInstance(eqn_disc.children[0], pybamm.Matrix)
-            self.assertIsInstance(eqn_disc.children[1], pybamm.StateVector)
 
             combined_submesh = mesh.combine_submeshes(*whole_cell)
             y = combined_submesh[0].nodes ** 2
@@ -491,14 +490,13 @@ class TestDiscretise(unittest.TestCase):
             )
 
         # More complex expressions
-        for eqn in [var * pybamm.grad(var), var * pybamm.div(var)]:
+        for eqn in [var * pybamm.grad(var), var * pybamm.div(pybamm.grad(var))]:
             eqn_disc = disc.process_symbol(eqn)
 
             self.assertIsInstance(eqn_disc, pybamm.Multiplication)
             self.assertIsInstance(eqn_disc.children[0], pybamm.StateVector)
             self.assertIsInstance(eqn_disc.children[1], pybamm.MatrixMultiplication)
             self.assertIsInstance(eqn_disc.children[1].children[0], pybamm.Matrix)
-            self.assertIsInstance(eqn_disc.children[1].children[1], pybamm.StateVector)
 
             y = combined_submesh[0].nodes ** 2
             var_disc = disc.process_symbol(var)
@@ -693,6 +691,10 @@ class TestDiscretise(unittest.TestCase):
         model_jacobian = model.jacobian.evaluate(0, y0)
         np.testing.assert_array_equal(model_jacobian.toarray(), jacobian.toarray())
 
+        # test that discretising again gives an error
+        with self.assertRaisesRegex(pybamm.ModelError, "Cannot re-discretise a model"):
+            disc.process_model(model)
+
         # test that not enough initial conditions raises an error
         model = pybamm.BaseModel()
         model.rhs = {c: pybamm.div(N), T: pybamm.div(q), S: pybamm.div(p)}
@@ -1116,7 +1118,7 @@ class TestDiscretise(unittest.TestCase):
 
         # check doesn't raise if concatenation
         model.variables = {c_n.name: pybamm.Concatenation(c_n, c_s)}
-        disc.process_model(model)
+        disc.process_model(model, inplace=False)
 
         # check doesn't raise if broadcast
         model.variables = {
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 e3fdb43c7..a72faf649 100644
--- a/tests/unit/test_expression_tree/test_operations/test_copy.py
+++ b/tests/unit/test_expression_tree/test_operations/test_copy.py
@@ -26,8 +26,8 @@ class TestCopy(unittest.TestCase):
             abs(a),
             pybamm.Function(np.sin, a),
             pybamm.FunctionParameter("function", {"a": a}),
-            pybamm.grad(a),
-            pybamm.div(a),
+            pybamm.grad(v_n),
+            pybamm.div(pybamm.grad(v_n)),
             pybamm.Integral(a, pybamm.t),
             pybamm.BoundaryValue(v_n, "right"),
             pybamm.BoundaryGradient(v_n, "right"),
diff --git a/tests/unit/test_expression_tree/test_operations/test_jac.py b/tests/unit/test_expression_tree/test_operations/test_jac.py
index 56895d26b..f7374068d 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac.py
@@ -14,6 +14,13 @@ def test_multi_var_function(arg1, arg2):
 
 
 class TestJacobian(unittest.TestCase):
+    def test_variable_is_statevector(self):
+        a = pybamm.Symbol("a")
+        with self.assertRaisesRegex(
+            TypeError, "Jacobian can only be taken with respect to a 'StateVector'"
+        ):
+            a.jac(a)
+
     def test_linear(self):
         y = pybamm.StateVector(slice(0, 4))
         u = pybamm.StateVector(slice(0, 2))
@@ -233,7 +240,7 @@ class TestJacobian(unittest.TestCase):
         a = pybamm.Scalar(1)
         b = pybamm.Scalar(2)
 
-        y = pybamm.Variable("y")
+        y = pybamm.StateVector(slice(0, 1))
 
         self.assertEqual(a.jac(y).evaluate(), 0)
 
@@ -261,14 +268,16 @@ class TestJacobian(unittest.TestCase):
     def test_spatial_operator(self):
         a = pybamm.Variable("a")
         b = pybamm.SpatialOperator("Operator", a)
+        y = pybamm.StateVector(slice(0, 1))
         with self.assertRaises(NotImplementedError):
-            b.jac(None)
+            b.jac(y)
 
     def test_jac_of_unary_operator(self):
         a = pybamm.Scalar(1)
         b = pybamm.UnaryOperator("Operator", a)
+        y = pybamm.StateVector(slice(0, 1))
         with self.assertRaises(NotImplementedError):
-            b.jac(None)
+            b.jac(y)
 
     def test_jac_of_independent_variable(self):
         a = pybamm.IndependentVariable("Variable")
diff --git a/tests/unit/test_expression_tree/test_operations/test_simplify.py b/tests/unit/test_expression_tree/test_operations/test_simplify.py
index d898709a7..6ce32eb3c 100644
--- a/tests/unit/test_expression_tree/test_operations/test_simplify.py
+++ b/tests/unit/test_expression_tree/test_operations/test_simplify.py
@@ -10,7 +10,7 @@ from tests import get_discretisation_for_testing
 
 class TestSimplify(unittest.TestCase):
     def test_symbol_simplify(self):
-        a = pybamm.Scalar(0)
+        a = pybamm.Scalar(0, domain="domain")
         b = pybamm.Scalar(1)
         c = pybamm.Parameter("c")
         d = pybamm.Scalar(-1)
@@ -58,13 +58,17 @@ class TestSimplify(unittest.TestCase):
         # Gradient
         self.assertIsInstance((pybamm.grad(a)).simplify(), pybamm.Scalar)
         self.assertEqual((pybamm.grad(a)).simplify().evaluate(), 0)
-        v = pybamm.Variable("v")
-        self.assertIsInstance((pybamm.grad(v)).simplify(), pybamm.Gradient)
+        v = pybamm.Variable("v", domain="domain")
+        grad_v = pybamm.grad(v)
+        self.assertIsInstance(grad_v.simplify(), pybamm.Gradient)
 
         # Divergence
-        self.assertIsInstance((pybamm.div(a)).simplify(), pybamm.Scalar)
-        self.assertEqual((pybamm.div(a)).simplify().evaluate(), 0)
-        self.assertIsInstance((pybamm.div(v)).simplify(), pybamm.Divergence)
+        div_b = pybamm.div(pybamm.PrimaryBroadcastToEdges(b, "domain"))
+        self.assertIsInstance(div_b.simplify(), pybamm.PrimaryBroadcast)
+        self.assertEqual(div_b.simplify().child.child.evaluate(), 0)
+        self.assertIsInstance(
+            (pybamm.div(pybamm.grad(v))).simplify(), pybamm.Divergence
+        )
 
         # Integral
         self.assertIsInstance(
diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py
index 52638aebb..6aa58f033 100644
--- a/tests/unit/test_expression_tree/test_symbol.py
+++ b/tests/unit/test_expression_tree/test_symbol.py
@@ -271,12 +271,6 @@ class TestSymbol(unittest.TestCase):
             + r", \*, children\=\['c', 'd'\], domain=\['test'\]"
             + r", auxiliary_domains\=\{'sec': \"\['other test'\]\"\}\)",
         )
-        self.assertRegex(
-            pybamm.grad(a).__repr__(),
-            r"Gradient\("
-            + hex_regex
-            + r", grad, children\=\['a'\], domain=\[\], auxiliary_domains\=\{\}\)",
-        )
         self.assertRegex(
             pybamm.grad(c).__repr__(),
             r"Gradient\("
@@ -332,9 +326,9 @@ class TestSymbol(unittest.TestCase):
             rhs.visualise("StefanMaxwell_test")
 
     def test_has_spatial_derivatives(self):
-        var = pybamm.Variable("var")
+        var = pybamm.Variable("var", domain="test")
         grad_eqn = pybamm.grad(var)
-        div_eqn = pybamm.div(var)
+        div_eqn = pybamm.div(pybamm.standard_spatial_vars.x_edge)
         grad_div_eqn = pybamm.div(grad_eqn)
         algebraic_eqn = 2 * var + 3
         self.assertTrue(grad_eqn.has_symbol_of_classes(pybamm.Gradient))
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 3fad2f3c0..28d206b92 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -52,9 +52,58 @@ class TestUnaryOperators(unittest.TestCase):
         )
 
     def test_gradient(self):
+        # gradient of scalar symbol should fail
         a = pybamm.Symbol("a")
+        with self.assertRaisesRegex(
+            pybamm.DomainError, "Cannot take gradient of 'a' since its domain is empty"
+        ):
+            pybamm.Gradient(a)
+
+        # gradient of variable evaluating on edges should fail
+        a = pybamm.PrimaryBroadcastToEdges(pybamm.Scalar(1), "test")
+        with self.assertRaisesRegex(TypeError, "evaluates on edges"):
+            pybamm.Gradient(a)
+
+        # gradient of broadcast should return broadcasted zero
+        a = pybamm.PrimaryBroadcast(pybamm.Variable("a"), "test domain")
+        grad = pybamm.grad(a)
+        self.assertIsInstance(grad, pybamm.PrimaryBroadcastToEdges)
+        self.assertIsInstance(grad.child, pybamm.PrimaryBroadcast)
+        self.assertIsInstance(grad.child.child, pybamm.Scalar)
+        self.assertEqual(grad.child.child.value, 0)
+
+        # otherwise gradient should work
+        a = pybamm.Symbol("a", domain="test domain")
         grad = pybamm.Gradient(a)
         self.assertEqual(grad.children[0].name, a.name)
+        self.assertEqual(grad.domain, a.domain)
+
+    def test_div(self):
+        # divergence of scalar symbol should fail
+        a = pybamm.Symbol("a")
+        with self.assertRaisesRegex(
+            pybamm.DomainError,
+            "Cannot take divergence of 'a' since its domain is empty",
+        ):
+            pybamm.Divergence(a)
+
+        # divergence of variable evaluating on edges should fail
+        a = pybamm.PrimaryBroadcast(pybamm.Scalar(1), "test")
+        with self.assertRaisesRegex(TypeError, "evaluates on nodes"):
+            pybamm.Divergence(a)
+
+        # divergence of broadcast should return broadcasted zero
+        a = pybamm.PrimaryBroadcastToEdges(pybamm.Variable("a"), "test domain")
+        div = pybamm.div(a)
+        self.assertIsInstance(div, pybamm.PrimaryBroadcast)
+        self.assertIsInstance(div.child, pybamm.PrimaryBroadcast)
+        self.assertIsInstance(div.child.child, pybamm.Scalar)
+        self.assertEqual(div.child.child.value, 0)
+
+        # otherwise divergence should work
+        a = pybamm.Symbol("a", domain="test domain")
+        div = pybamm.Divergence(pybamm.Gradient(a))
+        self.assertEqual(div.domain, a.domain)
 
     def test_integral(self):
         # time integral
@@ -186,7 +235,7 @@ class TestUnaryOperators(unittest.TestCase):
             spatial_a.diff(a)
 
     def test_printing(self):
-        a = pybamm.Symbol("a")
+        a = pybamm.Symbol("a", domain="test")
         self.assertEqual(str(-a), "-a")
         grad = pybamm.Gradient(a)
         self.assertEqual(grad.name, "grad")
diff --git a/tests/unit/test_meshes/test_meshes.py b/tests/unit/test_meshes/test_meshes.py
index 068001023..c8ce228ca 100644
--- a/tests/unit/test_meshes/test_meshes.py
+++ b/tests/unit/test_meshes/test_meshes.py
@@ -213,6 +213,11 @@ class TestMesh(unittest.TestCase):
         with self.assertRaisesRegex(pybamm.DomainError, "trying"):
             mesh.combine_submeshes("negative electrode", "negative particle")
 
+        with self.assertRaisesRegex(
+            ValueError, "Submesh domains being combined cannot be empty"
+        ):
+            mesh.combine_submeshes()
+
     def test_ghost_cells(self):
         param = pybamm.ParameterValues(
             values={
diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py
index 1e19c0dfa..6e1d642e9 100644
--- a/tests/unit/test_models/test_base_model.py
+++ b/tests/unit/test_models/test_base_model.py
@@ -50,10 +50,16 @@ class TestBaseModel(unittest.TestCase):
 
         # Test number input
         c0 = pybamm.Symbol("c0")
-        model.initial_conditions = {c0: 34}
+        model.initial_conditions[c0] = 34
         self.assertIsInstance(model.initial_conditions[c0], pybamm.Scalar)
         self.assertEqual(model.initial_conditions[c0].value, 34)
 
+        # Variable in initial conditions should fail
+        with self.assertRaisesRegex(
+            TypeError, "Initial conditions cannot contain 'Variable' objects"
+        ):
+            model.initial_conditions = {c0: pybamm.Variable("v")}
+
         # non-matching domains should fail
         with self.assertRaises(pybamm.DomainError):
             model.initial_conditions = {
@@ -72,8 +78,9 @@ class TestBaseModel(unittest.TestCase):
 
         # Test number input
         c0 = pybamm.Symbol("c0")
-        model.boundary_conditions = {
-            c0: {"left": (-2, "Dirichlet"), "right": (4, "Dirichlet")}
+        model.boundary_conditions[c0] = {
+            "left": (-2, "Dirichlet"),
+            "right": (4, "Dirichlet"),
         }
         self.assertIsInstance(model.boundary_conditions[c0]["left"][0], pybamm.Scalar)
         self.assertIsInstance(model.boundary_conditions[c0]["right"][0], pybamm.Scalar)
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte_conductivity/test_full_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte_conductivity/test_full_conductivity.py
index 4f16b15b1..7056c1bb3 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte_conductivity/test_full_conductivity.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte_conductivity/test_full_conductivity.py
@@ -13,7 +13,11 @@ class TestFull(unittest.TestCase):
         a = pybamm.Scalar(0)
         variables = {
             "Electrolyte tortuosity": a,
-            "Electrolyte concentration": a,
+            "Electrolyte concentration": pybamm.FullBroadcast(
+                a,
+                ["negative electrode", "separator", "positive electrode"],
+                "current collector",
+            ),
             "Negative electrode interfacial current density": pybamm.FullBroadcast(
                 a, "negative electrode", "current collector"
             ),
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index 351f7b23e..2c04971e3 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -52,6 +52,11 @@ class TestParameterValues(unittest.TestCase):
         ):
             pybamm.ParameterValues(values=1, chemistry={})
 
+    def test_repr(self):
+        param = pybamm.ParameterValues({"a": 1})
+        self.assertEqual(repr(param), "{'a': 1}")
+        self.assertEqual(param._ipython_key_completions_(), ["a"])
+
     def test_update_from_chemistry(self):
         # incomplete chemistry
         with self.assertRaisesRegex(KeyError, "must provide 'cell' parameters"):
@@ -127,11 +132,11 @@ class TestParameterValues(unittest.TestCase):
         self.assertEqual(processed_integ.integration_variable[0].id, x.id)
 
         # process unary operation
-        grad = pybamm.Gradient(a)
+        v = pybamm.Variable("v", domain="test")
+        grad = pybamm.Gradient(v)
         processed_grad = parameter_values.process_symbol(grad)
         self.assertIsInstance(processed_grad, pybamm.Gradient)
-        self.assertIsInstance(processed_grad.children[0], pybamm.Scalar)
-        self.assertEqual(processed_grad.children[0].value, 1)
+        self.assertIsInstance(processed_grad.children[0], pybamm.Variable)
 
         # process delta function
         aa = pybamm.Parameter("a")
@@ -252,6 +257,8 @@ class TestParameterValues(unittest.TestCase):
                 "func": pybamm.load_function("process_symbol_test_function.py"),
                 "const": 254,
                 "float_func": lambda x: 42,
+                "mult": pybamm.InputParameter("b") * 5,
+                "bad type": np.array([1, 2, 3]),
             }
         )
         a = pybamm.InputParameter("a")
@@ -267,6 +274,12 @@ class TestParameterValues(unittest.TestCase):
         self.assertIsInstance(processed_const, pybamm.Scalar)
         self.assertEqual(processed_const.evaluate(), 254)
 
+        # process case where parameter provided is a pybamm symbol
+        # (e.g. a multiplication)
+        mult = pybamm.FunctionParameter("mult", {"a": a})
+        processed_mult = parameter_values.process_symbol(mult)
+        self.assertEqual(processed_mult.evaluate(inputs={"a": 14, "b": 63}), 63 * 5)
+
         # process differentiated function parameter
         diff_func = func.diff(a)
         processed_diff_func = parameter_values.process_symbol(diff_func)
@@ -277,13 +290,24 @@ class TestParameterValues(unittest.TestCase):
         processed_func = parameter_values.process_symbol(func)
         self.assertEqual(processed_func.evaluate(), 42)
 
+        # weird type raises error
+        func = pybamm.FunctionParameter("bad type", {"a": a})
+        with self.assertRaisesRegex(TypeError, "Parameter provided for"):
+            parameter_values.process_symbol(func)
+
         # function itself as input (different to the variable being an input)
-        parameter_values = pybamm.ParameterValues({"func": "[input]"})
+        parameter_values = pybamm.ParameterValues(
+            {"func": "[input]", "vector func": pybamm.InputParameter("vec", "test")}
+        )
         a = pybamm.Scalar(3)
         func = pybamm.FunctionParameter("func", {"a": a})
         processed_func = parameter_values.process_symbol(func)
         self.assertEqual(processed_func.evaluate(inputs={"func": 13}), 13)
 
+        func = pybamm.FunctionParameter("vector func", {"a": a})
+        processed_func = parameter_values.process_symbol(func)
+        self.assertEqual(processed_func.evaluate(inputs={"vec": 13}), 13)
+
     def test_process_inline_function_parameters(self):
         def D(c):
             return c ** 2
@@ -364,7 +388,13 @@ class TestParameterValues(unittest.TestCase):
                 "interpolation": "[data]lico2_data_example",
             },
             path=os.path.join(
-                "input", "parameters", "lithium-ion", "cathodes", "lico2_Marquis2019"
+                pybamm.root_dir(),
+                "pybamm",
+                "input",
+                "parameters",
+                "lithium-ion",
+                "cathodes",
+                "lico2_Marquis2019",
             ),
             check_already_exists=False,
         )
@@ -435,8 +465,8 @@ class TestParameterValues(unittest.TestCase):
         b = pybamm.Parameter("b")
         c = pybamm.Parameter("c")
         d = pybamm.Parameter("d")
-        var1 = pybamm.Variable("var1")
-        var2 = pybamm.Variable("var2")
+        var1 = pybamm.Variable("var1", domain="test")
+        var2 = pybamm.Variable("var2", domain="test")
         model.rhs = {var1: a * pybamm.grad(var1)}
         model.algebraic = {var2: c * var2}
         model.initial_conditions = {var1: b, var2: d}
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 47dab28e8..1d7e0f264 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -23,11 +23,11 @@ class TestCasadiSolver(unittest.TestCase):
 
         # create discretisation
         disc = pybamm.Discretisation()
-        disc.process_model(model)
+        model_disc = disc.process_model(model, inplace=False)
         # Solve
         solver = pybamm.CasadiSolver(mode="fast", rtol=1e-8, atol=1e-8)
         t_eval = np.linspace(0, 1, 100)
-        solution = solver.solve(model, t_eval)
+        solution = solver.solve(model_disc, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
         np.testing.assert_array_almost_equal(
             solution.y[0], np.exp(0.1 * solution.t), decimal=5
@@ -79,7 +79,7 @@ class TestCasadiSolver(unittest.TestCase):
 
         # create discretisation
         disc = pybamm.Discretisation()
-        disc.process_model(model)
+        model_disc = disc.process_model(model, inplace=False)
 
         solver = pybamm.CasadiSolver(extra_options_call={"regularity_check": False})
         solver_old = pybamm.CasadiSolver(
@@ -88,15 +88,15 @@ class TestCasadiSolver(unittest.TestCase):
         # Solve with failure at t=2
         t_eval = np.linspace(0, 20, 100)
         with self.assertRaises(pybamm.SolverError):
-            solver.solve(model, t_eval)
+            solver.solve(model_disc, t_eval)
         with self.assertRaises(pybamm.SolverError):
-            solver_old.solve(model, t_eval)
+            solver_old.solve(model_disc, t_eval)
         # Solve with failure at t=0
         model.initial_conditions = {var: 0}
-        disc.process_model(model)
+        model_disc = disc.process_model(model, inplace=False)
         t_eval = np.linspace(0, 20, 100)
         with self.assertRaises(pybamm.SolverError):
-            solver.solve(model, t_eval)
+            solver.solve(model_disc, t_eval)
 
     def test_model_solver_events(self):
         # Create model
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 1bab2b4d4..985ebf0f6 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -529,12 +529,12 @@ class TestScikitsSolvers(unittest.TestCase):
                 pybamm.Event("var2 = 2.5", pybamm.min(var2 - 2.5)),
             ]
             disc = get_discretisation_for_testing()
-            disc.process_model(model)
+            model_disc = disc.process_model(model, inplace=False)
 
             # Solve
             solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
             t_eval = np.linspace(0, 5, 100)
-            solution = solver.solve(model, t_eval)
+            solution = solver.solve(model_disc, t_eval)
             np.testing.assert_array_less(solution.y[0], 1.5)
             np.testing.assert_array_less(solution.y[-1], 2.5)
             np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index fd21be080..eadb6fcbd 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -257,11 +257,11 @@ class TestScipySolver(unittest.TestCase):
             mesh = get_mesh_for_testing()
             spatial_methods = {"macroscale": pybamm.FiniteVolume()}
             disc = pybamm.Discretisation(mesh, spatial_methods)
-            disc.process_model(model)
+            model_disc = disc.process_model(model, inplace=False)
             # Solve
             solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
             t_eval = np.linspace(0, 10, 100)
-            solution = solver.solve(model, t_eval)
+            solution = solver.solve(model_disc, t_eval)
             self.assertLess(len(solution.t), len(t_eval))
             np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
             np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))

pip install -e .[dev,docs]
parse_log_pytest
requirements.txt
pytest
apt-get update apt-get install -y gcc gfortran
3.7
requirements.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
aiofiles==22.1.0
aiosqlite==0.19.0
alabaster==0.7.13
anyio==3.7.1
anytree==2.12.1
argon2-cffi==23.1.0
argon2-cffi-bindings==21.2.0
arrow==1.2.3
attrs==24.2.0
autograd==1.6.2
Babel==2.14.0
backcall==0.2.0
beautifulsoup4==4.13.3
black==23.3.0
bleach==6.0.0
cached-property==1.5.2
casadi==3.7.0
certifi @ file:///croot/certifi_1671487769961/work/certifi
cffi==1.15.1
charset-normalizer==3.4.1
click==8.1.8
comm==0.1.4
cycler==0.11.0
debugpy==1.7.0
decorator==5.1.1
defusedxml==0.7.1
docutils==0.17.1
entrypoints==0.4
exceptiongroup==1.2.2
fastjsonschema==2.21.1
flake8==5.0.4
fonttools==4.38.0
fqdn==1.5.1
future==1.0.0
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
importlib-metadata==4.2.0
importlib-resources==5.12.0
iniconfig==2.0.0
ipykernel==6.16.2
ipython==7.34.0
ipython-genutils==0.2.0
ipywidgets==8.1.5
isoduration==20.11.0
jedi==0.19.2
Jinja2==3.1.6
json5==0.9.16
jsonpointer==3.0.0
jsonschema==4.17.3
jupyter==1.1.1
jupyter-console==6.6.3
jupyter-events==0.6.3
jupyter-server==1.24.0
jupyter-ydoc==0.2.5
jupyter_client==7.4.9
jupyter_core==4.12.0
jupyter_server_fileid==0.9.3
jupyter_server_ydoc==0.8.0
jupyterlab==3.5.3
jupyterlab-pygments==0.2.2
jupyterlab_server==2.15.2
jupyterlab_widgets==3.0.13
kiwisolver==1.4.5
MarkupSafe==2.1.5
matplotlib==3.5.3
matplotlib-inline==0.1.6
mccabe==0.7.0
mistune==3.0.2
mypy-extensions==1.0.0
nbclassic==1.2.0
nbclient==0.7.4
nbconvert==7.6.0
nbformat==5.8.0
nest-asyncio==1.6.0
notebook==6.5.7
notebook_shim==0.2.4
numpy==1.21.6
packaging==24.0
pandas==1.3.5
pandocfilters==1.5.1
parso==0.8.4
pathspec==0.11.2
pexpect==4.9.0
pickleshare==0.7.5
Pillow==9.5.0
pkgutil_resolve_name==1.3.10
platformdirs==4.0.0
pluggy==1.2.0
prometheus-client==0.17.1
prompt_toolkit==3.0.48
psutil==7.0.0
ptyprocess==0.7.0
-e git+https://github.com/pybamm-team/PyBaMM.git@6852e0e518157e6802ce83a2549562e7d0ed4b9f#egg=pybamm
pycodestyle==2.9.1
pycparser==2.21
pyflakes==2.5.0
Pygments==2.17.2
pyparsing==3.1.4
pyrsistent==0.19.3
pytest==7.4.4
python-dateutil==2.9.0.post0
python-json-logger==3.0.1
pytz==2025.2
PyYAML==6.0.1
pyzmq==26.2.1
requests==2.31.0
rfc3339-validator==0.1.4
rfc3986-validator==0.1.1
scikit-fem==8.1.0
scipy==1.7.3
Send2Trash==1.8.3
six==1.17.0
sniffio==1.3.1
snowballstemmer==2.2.0
soupsieve==2.4.1
Sphinx==4.3.2
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
terminado==0.17.1
tinycss2==1.2.1
tomli==2.0.1
tornado==6.2
traitlets==5.9.0
typed-ast==1.5.5
typing_extensions==4.7.1
uri-template==1.3.0
urllib3==2.0.7
wcwidth==0.2.13
webcolors==1.13
webencodings==0.5.1
websocket-client==1.6.1
widgetsnbextension==4.0.13
y-py==0.6.2
ypy-websocket==0.8.4
zipp==3.15.0

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2022.12.7=py37h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.4.4=h6a678d5_1
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - pip=22.3.1=py37h06a4308_0
  - python=3.7.16=h7a1cb2a_0
  - readline=8.2=h5eee18b_0
  - setuptools=65.6.3=py37h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.38.4=py37h06a4308_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - aiofiles==22.1.0
      - aiosqlite==0.19.0
      - alabaster==0.7.13
      - anyio==3.7.1
      - anytree==2.12.1
      - argon2-cffi==23.1.0
      - argon2-cffi-bindings==21.2.0
      - arrow==1.2.3
      - attrs==24.2.0
      - autograd==1.6.2
      - babel==2.14.0
      - backcall==0.2.0
      - beautifulsoup4==4.13.3
      - black==23.3.0
      - bleach==6.0.0
      - cached-property==1.5.2
      - casadi==3.7.0
      - cffi==1.15.1
      - charset-normalizer==3.4.1
      - click==8.1.8
      - comm==0.1.4
      - cycler==0.11.0
      - debugpy==1.7.0
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.17.1
      - entrypoints==0.4
      - exceptiongroup==1.2.2
      - fastjsonschema==2.21.1
      - flake8==5.0.4
      - fonttools==4.38.0
      - fqdn==1.5.1
      - future==1.0.0
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - importlib-metadata==4.2.0
      - importlib-resources==5.12.0
      - iniconfig==2.0.0
      - ipykernel==6.16.2
      - ipython==7.34.0
      - ipython-genutils==0.2.0
      - ipywidgets==8.1.5
      - isoduration==20.11.0
      - jedi==0.19.2
      - jinja2==3.1.6
      - json5==0.9.16
      - jsonpointer==3.0.0
      - jsonschema==4.17.3
      - jupyter==1.1.1
      - jupyter-client==7.4.9
      - jupyter-console==6.6.3
      - jupyter-core==4.12.0
      - jupyter-events==0.6.3
      - jupyter-server==1.24.0
      - jupyter-server-fileid==0.9.3
      - jupyter-server-ydoc==0.8.0
      - jupyter-ydoc==0.2.5
      - jupyterlab==3.5.3
      - jupyterlab-pygments==0.2.2
      - jupyterlab-server==2.15.2
      - jupyterlab-widgets==3.0.13
      - kiwisolver==1.4.5
      - markupsafe==2.1.5
      - matplotlib==3.5.3
      - matplotlib-inline==0.1.6
      - mccabe==0.7.0
      - mistune==3.0.2
      - mypy-extensions==1.0.0
      - nbclassic==1.2.0
      - nbclient==0.7.4
      - nbconvert==7.6.0
      - nbformat==5.8.0
      - nest-asyncio==1.6.0
      - notebook==6.5.7
      - notebook-shim==0.2.4
      - numpy==1.21.6
      - packaging==24.0
      - pandas==1.3.5
      - pandocfilters==1.5.1
      - parso==0.8.4
      - pathspec==0.11.2
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==9.5.0
      - pkgutil-resolve-name==1.3.10
      - platformdirs==4.0.0
      - pluggy==1.2.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.48
      - psutil==7.0.0
      - ptyprocess==0.7.0
      - pycodestyle==2.9.1
      - pycparser==2.21
      - pyflakes==2.5.0
      - pygments==2.17.2
      - pyparsing==3.1.4
      - pyrsistent==0.19.3
      - pytest==7.4.4
      - python-dateutil==2.9.0.post0
      - python-json-logger==3.0.1
      - pytz==2025.2
      - pyyaml==6.0.1
      - pyzmq==26.2.1
      - requests==2.31.0
      - rfc3339-validator==0.1.4
      - rfc3986-validator==0.1.1
      - scikit-fem==8.1.0
      - scipy==1.7.3
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.3.1
      - snowballstemmer==2.2.0
      - soupsieve==2.4.1
      - sphinx==4.3.2
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - terminado==0.17.1
      - tinycss2==1.2.1
      - tomli==2.0.1
      - tornado==6.2
      - traitlets==5.9.0
      - typed-ast==1.5.5
      - typing-extensions==4.7.1
      - uri-template==1.3.0
      - urllib3==2.0.7
      - wcwidth==0.2.13
      - webcolors==1.13
      - webencodings==0.5.1
      - websocket-client==1.6.1
      - widgetsnbextension==4.0.13
      - y-py==0.6.2
      - ypy-websocket==0.8.4
      - zipp==3.15.0
prefix: /opt/conda/envs/PyBaMM

tests/integration/test_spatial_methods/test_finite_volume.py::TestFiniteVolumeConvergence::test_cartesian_div_convergence tests/integration/test_spatial_methods/test_finite_volume.py::TestFiniteVolumeConvergence::test_grad_div_broadcast tests/integration/test_spatial_methods/test_finite_volume.py::TestFiniteVolumeConvergence::test_p2d_spherical_convergence_quadratic tests/integration/test_spatial_methods/test_finite_volume.py::TestFiniteVolumeConvergence::test_p2d_with_x_dep_bcs_spherical_convergence tests/integration/test_spatial_methods/test_finite_volume.py::TestFiniteVolumeConvergence::test_spherical_div_convergence_linear tests/integration/test_spatial_methods/test_finite_volume.py::TestFiniteVolumeConvergence::test_spherical_div_convergence_quadratic tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_ode tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_variable_is_statevector tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_symbol_simplify tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_div tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_gradient tests/unit/test_meshes/test_meshes.py::TestMesh::test_combine_submeshes tests/unit/test_models/test_base_model.py::TestBaseModel::test_boundary_conditions_set_get tests/unit/test_models/test_base_model.py::TestBaseModel::test_initial_conditions_set_get tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_function_parameter
tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_mass_matirx_inverse tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_visualise
tests/integration/test_spatial_methods/test_finite_volume.py::TestFiniteVolumeConvergence::test_cartesian_spherical_grad_convergence tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_add_internal_boundary_conditions tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_adding_0D_external_variable tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_adding_0D_external_variable_fail tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_adding_1D_external_variable tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_adding_2D_external_variable_fail tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_broadcast tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_broadcast_2D tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_check_tab_bcs_error tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenate_in_order tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation_2D tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation_external_variables tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation_of_scalars tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_discretise_slicing tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_discretise_spatial_operator tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_exceptions tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_no_mesh tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_complex_expression tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_dict tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_empty_model tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_input_variable tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_concatenation tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_dae tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_fail tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_not_inplace tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_symbol_base tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_variables_dict tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_with_no_check tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_secondary_broadcast_2D tests/unit/test_expression_tree/test_operations/test_copy.py::TestCopy::test_symbol_new_copy tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_functions tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_index tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_jac_of_abs tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_jac_of_domain_concatenation tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_jac_of_heaviside tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_jac_of_independent_variable tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_jac_of_inner tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_jac_of_minimum_maximum tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_jac_of_number tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_jac_of_sign tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_jac_of_symbol tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_jac_of_unary_operator tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_linear tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_linear_ydot tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_multislice_raises tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_nonlinear tests/unit/test_expression_tree/test_operations/test_jac.py::TestJacobian::test_spatial_operator tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_domain_concatenation_simplify tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_matrix_divide_simplify tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_matrix_simplifications tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_simplify_broadcast tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_simplify_concatenation_state_vectors tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_simplify_heaviside tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_simplify_inner tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_vector_zero_simplify tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_cached_children tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_evaluate_ignoring_errors tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_has_spatial_derivatives tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_multiple_symbols tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_orphans tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_shape tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_shape_and_size_for_testing tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_auxiliary_domains tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_diff tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_domains tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_evaluates_to_number tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_evaluation tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_init tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_is_constant tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_methods tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_repr tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_test_shape tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_absolute tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_boundary_operators tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_boundary_value tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_delta_function tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_diff tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_evaluates_on_edges tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_id tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_index tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_integral tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_negation tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_printing tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_r_average tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_sign tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_unary_operator tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_x_average tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_yz_average tests/unit/test_meshes/test_meshes.py::TestMesh::test_1plus1D_tabs_left_right tests/unit/test_meshes/test_meshes.py::TestMesh::test_1plus1D_tabs_right_left tests/unit/test_meshes/test_meshes.py::TestMesh::test_ghost_cells tests/unit/test_meshes/test_meshes.py::TestMesh::test_init_failure tests/unit/test_meshes/test_meshes.py::TestMesh::test_mesh_coord_sys tests/unit/test_meshes/test_meshes.py::TestMesh::test_mesh_creation tests/unit/test_meshes/test_meshes.py::TestMesh::test_mesh_creation_no_parameters tests/unit/test_meshes/test_meshes.py::TestMesh::test_mesh_sizes tests/unit/test_meshes/test_meshes.py::TestMesh::test_multiple_meshes tests/unit/test_meshes/test_meshes.py::TestMesh::test_multiple_meshes_macro tests/unit/test_meshes/test_meshes.py::TestMesh::test_submesh_creation tests/unit/test_meshes/test_meshes.py::TestMesh::test_unimplemented_meshes tests/unit/test_meshes/test_meshes.py::TestMeshGenerator::test_init_name tests/unit/test_models/test_base_model.py::TestBaseModel::test_algebraic_set_get tests/unit/test_models/test_base_model.py::TestBaseModel::test_check_well_posedness_initial_boundary_conditions tests/unit/test_models/test_base_model.py::TestBaseModel::test_check_well_posedness_output_variables tests/unit/test_models/test_base_model.py::TestBaseModel::test_check_well_posedness_variables tests/unit/test_models/test_base_model.py::TestBaseModel::test_jac_set_get tests/unit/test_models/test_base_model.py::TestBaseModel::test_model_dict_behaviour tests/unit/test_models/test_base_model.py::TestBaseModel::test_read_input_parameters tests/unit/test_models/test_base_model.py::TestBaseModel::test_rhs_set_get tests/unit/test_models/test_base_model.py::TestBaseModel::test_update tests/unit/test_models/test_base_model.py::TestBaseModel::test_variables_set_get tests/unit/test_models/test_base_model.py::TestStandardBatteryBaseModel::test_default_parameters tests/unit/test_models/test_base_model.py::TestStandardBatteryBaseModel::test_default_solver tests/unit/test_models/test_base_model.py::TestStandardBatteryBaseModel::test_timescale tests/unit/test_models/test_submodels/test_electrolyte_conductivity/test_full_conductivity.py::TestFull::test_public_functions tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_check_parameter_values tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_evaluate tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_init tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_inplace tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_interpolant_against_function tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_multi_var_function_parameter tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_multi_var_function_with_parameters tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_complex_expression tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_empty_model tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_inline_function_parameters tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_input_parameter tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_interpolant tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_model tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_symbol tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_read_parameters_csv tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_repr tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_update tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_update_from_chemistry tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_update_model tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_bad_mode tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_dae_solver_algebraic_model tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_dae_inputs_in_initial_conditions tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_events tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_failure tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_python tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_with_external tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_with_inputs tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_with_non_identity_mass tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_step tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_step_events tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_step_with_input tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_failure tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_inputs_in_initial_conditions tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_ode_with_jacobian_python tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_python tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_event_python tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_event_with_casadi tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_external tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_inputs tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_inputs_with_casadi tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_step_python
pybamm/parameters/parameter_values.py pybamm/expression_tree/functions.py pybamm/expression_tree/unary_operators.py pybamm/simulation.py pybamm/discretisations/discretisation.py pybamm/meshes/meshes.py CHANGELOG.md pybamm/expression_tree/operations/simplify.py pybamm/models/submodels/base_submodel.py pybamm/input/parameters/lithium-ion/experiments/1C_charge_from_empty_Mohtat2020/parameters.csv pybamm/expression_tree/symbol.py pybamm/models/base_model.py pybamm/models/submodels/particle/base_particle.py pybamm/models/submodels/convection/base_convection.py
pybamm/parameters/parameter_values.py pybamm/expression_tree/functions.py pybamm/expression_tree/unary_operators.py pybamm/simulation.py pybamm/discretisations/discretisation.py pybamm/meshes/meshes.py CHANGELOG.md pybamm/expression_tree/operations/simplify.py pybamm/models/submodels/base_submodel.py pybamm/input/parameters/lithium-ion/experiments/1C_charge_from_empty_Mohtat2020/parameters.csv pybamm/expression_tree/symbol.py pybamm/models/base_model.py pybamm/models/submodels/particle/base_particle.py pybamm/models/submodels/convection/base_convection.py