Repository: pybamm-team/PyBaMM (python)

Issue:
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`

Pull request description:
Issue 963 misc bugs
# Description

- Add some extra checks when creating model
- Make grad (broadcast) and div(broadcast) automatically return broadcast(0)
- don't allow rediscretising a model once discretised (if allowed, it usually fails anyway)

Fixes #963 

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/master/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [ ] New feature (non-breaking change which adds functionality)
- [ ] Optimization (back-end change that speeds up the code)
- [x] Bug fix (non-breaking change which fixes an issue)


# Key checklist:

- [x] No style issues: `$ flake8`
- [ ] All tests pass: `$ python run-tests.py --unit`
- [x] The documentation builds: `$ cd docs` and then `$ make clean; make html`

You can run all three at once, using `$ python run-tests.py --quick`.

## Further checks:

- [x] Code is commented, particularly in hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works


Patch:
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(


Test patch:
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))