Repository: pybamm-team/PyBaMM (python)

Issue:
cache symbol shapes for testing
For some bigger problems discretising all of the variables can be slow. It would be good to have the option to only discretise the variables actually used by the model and then be able to discretise a user specified list of variables after solving for post-processing etc.

Pull request:
Issue 749 cache symbol shape
# Description

Save symbol shapes when they have already been evaluated

Fixes #749 

## 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)
- [x] Optimization (back-end change that speeds up the code)
- [ ] Bug fix (non-breaking change which fixes an issue)


# Key checklist:

- [ ] No style issues: `$ flake8`
- [ ] All tests pass: `$ python run-tests.py --unit`
- [ ] 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:

- [ ] Code is commented, particularly in hard-to-understand areas
- [ ] Tests added that prove fix is effective or that feature works


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9a15d3259..ee7535c19 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -31,6 +31,7 @@
 
 ## Optimizations
 
+-   Added caching for shape evaluation, used during discretisation ([#780](https://github.com/pybamm-team/PyBaMM/pull/780))
 -   Added an option to skip model checks during discretisation, which could be slow for large models ([#739](https://github.com/pybamm-team/PyBaMM/pull/739))
 -   Use CasADi's automatic differentation algorithms by default when solving a model ([#714](https://github.com/pybamm-team/PyBaMM/pull/714))
 -   Avoid re-checking size when making a copy of an `Index` object ([#656](https://github.com/pybamm-team/PyBaMM/pull/656))
@@ -53,6 +54,7 @@
 
 ## Breaking changes
 
+-   Removed `Outer` and `Kron` nodes as no longer used ([#777](https://github.com/pybamm-team/PyBaMM/pull/777))
 -   Moved `results` to separate repositories ([#761](https://github.com/pybamm-team/PyBaMM/pull/761))
 -   The parameters "Bruggeman coefficient" must now be specified separately as "Bruggeman coefficient (electrolyte)" and "Bruggeman coefficient (electrode)"
 -   The current classes (`GetConstantCurrent`, `GetUserCurrent` and `GetUserData`) have now been removed. Please refer to the [`change-input-current` notebook](https://github.com/pybamm-team/PyBaMM/blob/master/examples/notebooks/change-input-current.ipynb) for information on how to specify an input current
diff --git a/docs/source/expression_tree/binary_operator.rst b/docs/source/expression_tree/binary_operator.rst
index 822864bbc..d8082d8c4 100644
--- a/docs/source/expression_tree/binary_operator.rst
+++ b/docs/source/expression_tree/binary_operator.rst
@@ -25,15 +25,7 @@ Binary Operators
 .. autoclass:: pybamm.Inner
   :members:
 
-.. autoclass:: pybamm.Outer
-  :members:
-
-.. autoclass:: pybamm.Kron
-  :members:
-
 .. autoclass:: pybamm.Heaviside
   :members:
 
-.. autofunction:: pybamm.outer
-
 .. autofunction:: pybamm.source
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index b47efcfa2..663540e3c 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -81,10 +81,7 @@ from .expression_tree.binary_operators import (
     Division,
     Inner,
     inner,
-    Outer,
-    Kron,
     Heaviside,
-    outer,
     source,
 )
 from .expression_tree.concatenations import (
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 1503e3c94..f99234378 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -989,7 +989,7 @@ class Discretisation(object):
         """
         Check variables in variable list against rhs
         Be lenient with size check if the variable in model.variables is broadcasted, or
-        a concatenation, or an outer product
+        a concatenation
         (if broadcasted, variable is a multiplication with a vector of ones)
         """
         for rhs_var in model.rhs.keys():
@@ -1001,7 +1001,6 @@ class Discretisation(object):
                 )
 
                 not_concatenation = not isinstance(var, pybamm.Concatenation)
-                not_outer = not isinstance(var, pybamm.Outer)
 
                 not_mult_by_one_vec = not (
                     isinstance(var, pybamm.Multiplication)
@@ -1012,7 +1011,6 @@ class Discretisation(object):
                 if (
                     different_shapes
                     and not_concatenation
-                    and not_outer
                     and not_mult_by_one_vec
                 ):
                     raise pybamm.ModelError(
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index f5f2cb024..6ac3cbc72 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -5,7 +5,7 @@ import pybamm
 
 import numpy as np
 import numbers
-from scipy.sparse import issparse, csr_matrix, kron
+from scipy.sparse import issparse, csr_matrix
 
 
 def is_scalar_zero(expr):
@@ -79,17 +79,8 @@ class BinaryOperator(pybamm.Symbol):
     def __init__(self, name, left, right):
         left, right = self.format(left, right)
 
-        # Check and process domains, except for Outer symbol which takes the outer
-        # product of two smbols in different domains, and gives it the domain of the
-        # right child.
-        if isinstance(self, (pybamm.Outer, pybamm.Kron)):
-            domain = right.domain
-            auxiliary_domains = {}
-            if domain != []:
-                auxiliary_domains["secondary"] = left.domain
-        else:
-            domain = self.get_children_domains(left.domain, right.domain)
-            auxiliary_domains = self.get_children_auxiliary_domains([left, right])
+        domain = self.get_children_domains(left.domain, right.domain)
+        auxiliary_domains = self.get_children_auxiliary_domains([left, right])
         super().__init__(
             name,
             children=[left, right],
@@ -116,11 +107,7 @@ class BinaryOperator(pybamm.Symbol):
             )
 
         # Do some broadcasting in special cases, to avoid having to do this manually
-        if (
-            not isinstance(self, (Outer, Kron))
-            and left.domain != []
-            and right.domain != []
-        ):
+        if left.domain != [] and right.domain != []:
             if (
                 left.domain != right.domain
                 and "secondary" in right.auxiliary_domains
@@ -192,7 +179,7 @@ class BinaryOperator(pybamm.Symbol):
             right = self.right.evaluate(t, y, u)
             return self._binary_evaluate(left, right)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape()`. """
         left = self.children[0].evaluate_for_shape()
         right = self.children[1].evaluate_for_shape()
@@ -654,86 +641,6 @@ def inner(left, right):
     return pybamm.Inner(left, right)
 
 
-class Outer(BinaryOperator):
-    """A node in the expression tree representing an outer product.
-    This takes a 1D vector in the current collector domain of size (n,1) and a 1D
-    variable of size (m,1), takes their outer product, and reshapes this into a vector
-    of size (nm,1). It can also take in a vector in a single particle and a vector
-    of the electrolyte domain to repeat that particle.
-    Note: this class might be a bit dangerous, so at the moment it is very restrictive
-    in what symbols can be passed to it
-
-    **Extends:** :class:`BinaryOperator`
-    """
-
-    def __init__(self, left, right):
-        """ See :meth:`pybamm.BinaryOperator.__init__()`. """
-        # cannot have certain types of objects in the right symbol, as these
-        # can already be 2D objects (so we can't take an outer product with them)
-        if right.has_symbol_of_classes(
-            (pybamm.Variable, pybamm.StateVector, pybamm.Matrix, pybamm.SpatialVariable)
-        ):
-            raise TypeError("right child must only contain Vectors and Scalars" "")
-
-        super().__init__("outer product", left, right)
-
-    def __str__(self):
-        """ See :meth:`pybamm.Symbol.__str__()`. """
-        return "outer({!s}, {!s})".format(self.left, self.right)
-
-    def diff(self, variable):
-        """ See :meth:`pybamm.Symbol.diff()`. """
-        raise NotImplementedError("diff not implemented for symbol of type 'Outer'")
-
-    def _outer_jac(self, left_jac, right_jac, variable):
-        """
-        Calculate jacobian of outer product.
-        See :meth:`pybamm.Jacobian._jac()`.
-        """
-        # right cannot be a StateVector, so no need for product rule
-        left, right = self.orphans
-        if left.evaluates_to_number():
-            # Return zeros of correct size
-            return pybamm.Matrix(
-                csr_matrix((self.size, variable.evaluation_array.count(True)))
-            )
-        else:
-            return pybamm.Kron(left_jac, right)
-
-    def _binary_evaluate(self, left, right):
-        """ See :meth:`pybamm.BinaryOperator._binary_evaluate()`. """
-
-        return np.outer(left, right).reshape(-1, 1)
-
-
-class Kron(BinaryOperator):
-    """A node in the expression tree representing a (sparse) kronecker product operator
-
-    **Extends:** :class:`BinaryOperator`
-    """
-
-    def __init__(self, left, right):
-        """ See :meth:`pybamm.BinaryOperator.__init__()`. """
-
-        super().__init__("kronecker product", left, right)
-
-    def __str__(self):
-        """ See :meth:`pybamm.Symbol.__str__()`. """
-        return "kron({!s}, {!s})".format(self.left, self.right)
-
-    def diff(self, variable):
-        """ See :meth:`pybamm.Symbol.diff()`. """
-        raise NotImplementedError("diff not implemented for symbol of type 'Kron'")
-
-    def _binary_jac(self, left_jac, right_jac):
-        """ See :meth:`pybamm.BinaryOperator._binary_jac()`. """
-        raise NotImplementedError("jac not implemented for symbol of type 'Kron'")
-
-    def _binary_evaluate(self, left, right):
-        """ See :meth:`pybamm.BinaryOperator._binary_evaluate()`. """
-        return csr_matrix(kron(left, right))
-
-
 class Heaviside(BinaryOperator):
     """A node in the expression tree representing a heaviside step function.
 
@@ -783,18 +690,6 @@ class Heaviside(BinaryOperator):
         return Heaviside(left, right, self.equal)
 
 
-def outer(left, right):
-    """
-    Return outer product of two symbols. If the symbols have the same domain, the outer
-    product is just a multiplication. If they have different domains, make a copy of the
-    left child with same domain as right child, and then take outer product.
-    """
-    try:
-        return left * right
-    except pybamm.DomainError:
-        return pybamm.Outer(left, right)
-
-
 def source(left, right, boundary=False):
     """A convinience function for creating (part of) an expression tree representing
     a source term. This is necessary for spatial methods where the mass matrix
diff --git a/pybamm/expression_tree/broadcasts.py b/pybamm/expression_tree/broadcasts.py
index 4011c80bd..31228d808 100644
--- a/pybamm/expression_tree/broadcasts.py
+++ b/pybamm/expression_tree/broadcasts.py
@@ -130,7 +130,7 @@ class PrimaryBroadcast(Broadcast):
         """ See :meth:`pybamm.UnaryOperator.simplify()`. """
         return PrimaryBroadcast(child, self.broadcast_domain)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """
         Returns a vector of NaNs to represent the shape of a Broadcast.
         See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()`
@@ -210,7 +210,7 @@ class SecondaryBroadcast(Broadcast):
         """ See :meth:`pybamm.UnaryOperator.simplify()`. """
         return SecondaryBroadcast(child, self.broadcast_domain)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """
         Returns a vector of NaNs to represent the shape of a Broadcast.
         See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()`
@@ -253,7 +253,7 @@ class FullBroadcast(Broadcast):
         """ See :meth:`pybamm.UnaryOperator.simplify()`. """
         return FullBroadcast(child, self.broadcast_domain, self.auxiliary_domains)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """
         Returns a vector of NaNs to represent the shape of a Broadcast.
         See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()`
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index af4e0a72d..8b0cfacc2 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -90,7 +90,7 @@ class Concatenation(pybamm.Symbol):
         new_symbol.clear_domains()
         return new_symbol
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape` """
         if len(self.children) == 0:
             return np.array([])
diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index 242d2756d..be98461a7 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -167,7 +167,7 @@ class Function(pybamm.Symbol):
             evaluated_children = [child.evaluate(t, y, u) for child in self.children]
             return self._function_evaluate(evaluated_children)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """
         Default behaviour: has same shape as all child
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
diff --git a/pybamm/expression_tree/independent_variable.py b/pybamm/expression_tree/independent_variable.py
index dceb84ea5..985b15fd1 100644
--- a/pybamm/expression_tree/independent_variable.py
+++ b/pybamm/expression_tree/independent_variable.py
@@ -27,7 +27,7 @@ class IndependentVariable(pybamm.Symbol):
     def __init__(self, name, domain=None, auxiliary_domains=None):
         super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
         return pybamm.evaluate_for_shape_using_domain(
             self.domain, self.auxiliary_domains
@@ -57,7 +57,7 @@ class Time(IndependentVariable):
             raise ValueError("t must be provided")
         return t
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """
         Return the scalar '0' to represent the shape of the independent variable `Time`.
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
diff --git a/pybamm/expression_tree/input_parameter.py b/pybamm/expression_tree/input_parameter.py
index 67b8a8755..8d0cbf1e6 100644
--- a/pybamm/expression_tree/input_parameter.py
+++ b/pybamm/expression_tree/input_parameter.py
@@ -25,7 +25,7 @@ class InputParameter(pybamm.Symbol):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
         return InputParameter(self.name)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """
         Returns the scalar 'NaN' to represent the shape of a parameter.
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py
index 983008d70..706960d3c 100644
--- a/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -59,11 +59,8 @@ class CasadiConverter(object):
             # process children
             converted_left = self.convert(left, t, y, u)
             converted_right = self.convert(right, t, y, u)
-            if isinstance(symbol, pybamm.Outer):
-                return casadi.kron(converted_left, converted_right)
-            else:
-                # _binary_evaluate defined in derived classes for specific rules
-                return symbol._binary_evaluate(converted_left, converted_right)
+            # _binary_evaluate defined in derived classes for specific rules
+            return symbol._binary_evaluate(converted_left, converted_right)
 
         elif isinstance(symbol, pybamm.UnaryOperator):
             converted_child = self.convert(symbol.child, t, y, u)
diff --git a/pybamm/expression_tree/operations/evaluate.py b/pybamm/expression_tree/operations/evaluate.py
index 6558357c3..a11769b06 100644
--- a/pybamm/expression_tree/operations/evaluate.py
+++ b/pybamm/expression_tree/operations/evaluate.py
@@ -92,14 +92,6 @@ def find_symbols(symbol, constant_symbols, variable_symbols):
                 "if scipy.sparse.issparse({1}) else "
                 "{0} * {1}".format(children_vars[0], children_vars[1])
             )
-        elif isinstance(symbol, pybamm.Outer):
-            symbol_str = "np.outer({}, {}).reshape(-1, 1)".format(
-                children_vars[0], children_vars[1]
-            )
-        elif isinstance(symbol, pybamm.Kron):
-            symbol_str = "scipy.sparse.csr_matrix(scipy.sparse.kron({}, {}))".format(
-                children_vars[0], children_vars[1]
-            )
         else:
             symbol_str = children_vars[0] + " " + symbol.name + " " + children_vars[1]
 
diff --git a/pybamm/expression_tree/operations/jacobian.py b/pybamm/expression_tree/operations/jacobian.py
index 629d498b9..42bd3683d 100644
--- a/pybamm/expression_tree/operations/jacobian.py
+++ b/pybamm/expression_tree/operations/jacobian.py
@@ -46,15 +46,8 @@ class Jacobian(object):
             # process children
             left_jac = self.jac(left, variable)
             right_jac = self.jac(right, variable)
-            # Need to treat outer differently. If the left child of an Outer
-            # evaluates to number then we need to return a matrix of zeros
-            # of the correct size, which requires variable.evaluation_array
-            if isinstance(symbol, pybamm.Outer):
-                # _outer_jac defined in pybamm.Outer
-                jac = symbol._outer_jac(left_jac, right_jac, variable)
-            else:
-                # _binary_jac defined in derived classes for specific rules
-                jac = symbol._binary_jac(left_jac, right_jac)
+            # _binary_jac defined in derived classes for specific rules
+            jac = symbol._binary_jac(left_jac, right_jac)
 
         elif isinstance(symbol, pybamm.UnaryOperator):
             child_jac = self.jac(symbol.child, variable)
diff --git a/pybamm/expression_tree/parameter.py b/pybamm/expression_tree/parameter.py
index f99bb0991..7000bca24 100644
--- a/pybamm/expression_tree/parameter.py
+++ b/pybamm/expression_tree/parameter.py
@@ -27,7 +27,7 @@ class Parameter(pybamm.Symbol):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
         return Parameter(self.name, self.domain)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """
         Returns the scalar 'NaN' to represent the shape of a parameter.
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
@@ -118,7 +118,7 @@ class FunctionParameter(pybamm.Symbol):
         """
         return FunctionParameter(self.name, *children, diff_variable=self.diff_variable)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """
         Returns the sum of the evaluated children
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index 688f07173..dbb1d32d3 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -173,7 +173,7 @@ class StateVector(pybamm.Symbol):
             evaluation_array=self.evaluation_array,
         )
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """
         Returns a vector of NaNs to represent the shape of a StateVector.
         The size of a StateVector is the number of True elements in its evaluation_array
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 89b1a8f89..37f063243 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -531,6 +531,14 @@ class Symbol(anytree.NodeMixin):
         shape is returned instead, using the symbol's domain.
         See :meth:`pybamm.Symbol.evaluate()`
         """
+        try:
+            return self._saved_evaluate_for_shape
+        except AttributeError:
+            self._saved_evaluate_for_shape = self._evaluate_for_shape()
+            return self._saved_evaluate_for_shape
+
+    def _evaluate_for_shape(self):
+        "See :meth:`Symbol.evaluate_for_shape`"
         return self.evaluate()
 
     def is_constant(self):
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index d8d13ed57..0b4f68e22 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -74,7 +74,7 @@ class UnaryOperator(pybamm.Symbol):
             child = self.child.evaluate(t, y, u)
             return self._unary_evaluate(child)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """
         Default behaviour: unary operator has same shape as child
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
@@ -228,7 +228,7 @@ class Index(UnaryOperator):
 
         return self.__class__(child, self.index, check_size=False)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         return self._unary_evaluate(self.children[0].evaluate_for_shape())
 
     def evaluates_on_edges(self):
@@ -347,7 +347,7 @@ class Mass(SpatialOperator):
     def __init__(self, child):
         super().__init__("mass", child)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         return pybamm.evaluate_for_shape_using_domain(self.domain, typ="matrix")
 
 
@@ -361,7 +361,7 @@ class BoundaryMass(SpatialOperator):
     def __init__(self, child):
         super().__init__("boundary mass", child)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         return pybamm.evaluate_for_shape_using_domain(self.domain, typ="matrix")
 
 
@@ -455,7 +455,7 @@ class Integral(SpatialOperator):
 
         return self.__class__(child, self.integration_variable)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
         return pybamm.evaluate_for_shape_using_domain(self.domain)
 
@@ -501,7 +501,7 @@ class IndefiniteIntegral(Integral):
         if isinstance(integration_variable, pybamm.SpatialVariable):
             self.name += " on {}".format(integration_variable.domain)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         return self.children[0].evaluate_for_shape()
 
 
@@ -550,7 +550,7 @@ class DefiniteIntegralVector(SpatialOperator):
 
         return self.__class__(child, vector_type=self.vector_type)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
         return pybamm.evaluate_for_shape_using_domain(self.domain)
 
@@ -611,7 +611,7 @@ class BoundaryIntegral(SpatialOperator):
 
         return self.__class__(child, region=self.region)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
         return pybamm.evaluate_for_shape_using_domain(self.domain)
 
@@ -722,7 +722,7 @@ class BoundaryOperator(SpatialOperator):
         """ See :meth:`UnaryOperator._unary_new_copy()`. """
         return self.__class__(child, self.side)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
         return pybamm.evaluate_for_shape_using_domain(
             self.domain, self.auxiliary_domains
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index c0a350717..6abd0f7ad 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -39,7 +39,7 @@ class Variable(pybamm.Symbol):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
         return Variable(self.name, self.domain, self.auxiliary_domains)
 
-    def evaluate_for_shape(self):
+    def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
         return pybamm.evaluate_for_shape_using_domain(
             self.domain, self.auxiliary_domains
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index ec9ee9433..1976b806c 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -98,9 +98,15 @@ class SpatialMethod:
         )
 
         if broadcast_type == "primary":
-            out = pybamm.Outer(
-                symbol, pybamm.Vector(np.ones(primary_domain_size), domain=domain)
-            )
+            # Make copies of the child stacked on top of each other
+            sub_vector = np.ones((primary_domain_size, 1))
+            if symbol.shape_for_testing == ():
+                out = symbol * pybamm.Vector(sub_vector)
+            else:
+                # Repeat for secondary points
+                matrix = csr_matrix(kron(eye(symbol.shape_for_testing[0]), sub_vector))
+                out = pybamm.Matrix(matrix) @ symbol
+            out.domain = domain
         elif broadcast_type == "secondary":
             secondary_domain_size = sum(
                 self.mesh[dom][0].npts_for_broadcast


Test patch:
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 8ff690c17..5e823c2a0 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -781,13 +781,18 @@ class TestDiscretise(unittest.TestCase):
 
         disc.set_variable_slices([var])
         broad_disc = disc.process_symbol(broad)
-        self.assertIsInstance(broad_disc, pybamm.Outer)
-        self.assertIsInstance(broad_disc.children[0], pybamm.StateVector)
-        self.assertIsInstance(broad_disc.children[1], pybamm.Vector)
+        self.assertIsInstance(broad_disc, pybamm.MatrixMultiplication)
+        self.assertIsInstance(broad_disc.children[0], pybamm.Matrix)
+        self.assertIsInstance(broad_disc.children[1], pybamm.StateVector)
         self.assertEqual(
             broad_disc.shape,
             (mesh["separator"][0].npts * mesh["current collector"][0].npts, 1),
         )
+        y_test = np.linspace(0, 1, mesh["current collector"][0].npts)
+        np.testing.assert_array_equal(
+            broad_disc.evaluate(y=y_test),
+            np.outer(y_test, np.ones(mesh["separator"][0].npts)).reshape(-1, 1),
+        )
 
     def test_secondary_broadcast_2D(self):
         # secondary broadcast in 2D --> Matrix multiplication
@@ -808,29 +813,6 @@ class TestDiscretise(unittest.TestCase):
             (mesh["negative particle"][0].npts * mesh["negative electrode"][0].npts, 1),
         )
 
-    def test_outer(self):
-
-        # create discretisation
-        disc = get_1p1d_discretisation_for_testing()
-        mesh = disc.mesh
-
-        var_z = pybamm.Variable("var_z", ["current collector"])
-        var_x = pybamm.Vector(
-            np.linspace(0, 1, mesh["separator"][0].npts), domain="separator"
-        )
-
-        # process Outer variable
-        disc.set_variable_slices([var_z, var_x])
-        outer = pybamm.outer(var_z, var_x)
-        outer_disc = disc.process_symbol(outer)
-        self.assertIsInstance(outer_disc, pybamm.Outer)
-        self.assertIsInstance(outer_disc.children[0], pybamm.StateVector)
-        self.assertIsInstance(outer_disc.children[1], pybamm.Vector)
-        self.assertEqual(
-            outer_disc.shape,
-            (mesh["separator"][0].npts * mesh["current collector"][0].npts, 1),
-        )
-
     def test_concatenation(self):
         a = pybamm.Symbol("a")
         b = pybamm.Symbol("b")
diff --git a/tests/unit/test_expression_tree/test_binary_operators.py b/tests/unit/test_expression_tree/test_binary_operators.py
index 2381bcab6..8127b41b2 100644
--- a/tests/unit/test_expression_tree/test_binary_operators.py
+++ b/tests/unit/test_expression_tree/test_binary_operators.py
@@ -58,54 +58,6 @@ class TestBinaryOperators(unittest.TestCase):
         pow2 = pybamm.Power(a, b)
         self.assertEqual(pow2.evaluate(), 16)
 
-    def test_outer(self):
-        # Outer class
-        v = pybamm.Vector(np.ones(5), domain="current collector")
-        w = pybamm.Vector(2 * np.ones(3), domain="test")
-        outer = pybamm.Outer(v, w)
-        np.testing.assert_array_equal(outer.evaluate(), 2 * np.ones((15, 1)))
-        self.assertEqual(outer.domain, w.domain)
-        self.assertEqual(
-            str(outer), "outer(Column vector of length 5, Column vector of length 3)"
-        )
-
-        # outer function
-        # if there is no domain clash, normal multiplication is retured
-        u = pybamm.Vector(np.linspace(0, 1, 5))
-        outer = pybamm.outer(u, v)
-        self.assertIsInstance(outer, pybamm.Multiplication)
-        np.testing.assert_array_equal(outer.evaluate(), u.evaluate())
-        # otherwise, Outer class is returned
-        outer_fun = pybamm.outer(v, w)
-        outer_class = pybamm.Outer(v, w)
-        self.assertEqual(outer_fun.id, outer_class.id)
-
-        # failures
-        y = pybamm.StateVector(slice(10))
-        with self.assertRaisesRegex(
-            TypeError, "right child must only contain Vectors and Scalars"
-        ):
-            pybamm.Outer(v, y)
-        with self.assertRaises(NotImplementedError):
-            outer_fun.diff(None)
-
-    def test_kron(self):
-        # Kron class
-        A = pybamm.Matrix(np.eye(2))
-        b = pybamm.Vector(np.array([[4], [5]]))
-        kron = pybamm.Kron(A, b)
-        np.testing.assert_array_equal(
-            kron.evaluate().toarray(), np.kron(A.entries, b.entries)
-        )
-
-        # failures
-        with self.assertRaises(NotImplementedError):
-            kron.diff(None)
-
-        y = pybamm.StateVector(slice(0, 2))
-        with self.assertRaises(NotImplementedError):
-            kron.jac(y)
-
     def test_known_eval(self):
         # Scalars
         a = pybamm.Scalar(4)
diff --git a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
index d0ee4fcf1..a5fc04748 100644
--- a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
+++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
@@ -71,12 +71,6 @@ class TestCasadiConverter(unittest.TestCase):
         # State Vector
         self.assert_casadi_equal(pybamm_y.to_casadi(casadi_t, casadi_y), casadi_y)
 
-        # outer product
-        outer = pybamm.Outer(pybamm_a, pybamm_a)
-        self.assert_casadi_equal(
-            outer.to_casadi(), casadi.MX(outer.evaluate()), evalf=True
-        )
-
     def test_special_functions(self):
         a = pybamm.Array(np.array([1, 2, 3, 4, 5]))
         self.assert_casadi_equal(pybamm.max(a).to_casadi(), casadi.MX(5), evalf=True)
diff --git a/tests/unit/test_expression_tree/test_operations/test_evaluate.py b/tests/unit/test_expression_tree/test_operations/test_evaluate.py
index 52addcb63..cd6aae7dd 100644
--- a/tests/unit/test_expression_tree/test_operations/test_evaluate.py
+++ b/tests/unit/test_expression_tree/test_operations/test_evaluate.py
@@ -391,14 +391,6 @@ class TestEvaluate(unittest.TestCase):
             result = evaluator.evaluate(t=t, y=y).toarray()
             np.testing.assert_allclose(result, expr.evaluate(t=t, y=y).toarray())
 
-        # test Outer
-        v = pybamm.Vector(np.ones(5), domain="current collector")
-        w = pybamm.Vector(2 * np.ones(3), domain="test")
-        expr = pybamm.Outer(v, w)
-        evaluator = pybamm.EvaluatorPython(expr)
-        result = evaluator.evaluate()
-        np.testing.assert_allclose(result, expr.evaluate())
-
         # test Inner
         v = pybamm.Vector(np.ones(5), domain="test")
         w = pybamm.Vector(2 * np.ones(5), domain="test")
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 526f31390..270b030c4 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac.py
@@ -56,32 +56,6 @@ class TestJacobian(unittest.TestCase):
         du_dv = u.jac(v).evaluate().toarray()
         np.testing.assert_array_equal(du_dv, jacobian)
 
-        # test Jacobian of Outer (must set domain to be 'current collector')
-        u.domain = ["current collector"]
-        func = pybamm.Outer(u, pybamm.Scalar(4))
-        jacobian = np.array([[4, 0, 0, 0], [0, 4, 0, 0]])
-        dfunc_dy = func.jac(y).evaluate(y=y0)
-        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
-
-        func = pybamm.Outer(u, pybamm.Vector(np.array([1, 2, 3])))
-        jacobian = np.array(
-            [
-                [1, 0, 0, 0],
-                [2, 0, 0, 0],
-                [3, 0, 0, 0],
-                [0, 1, 0, 0],
-                [0, 2, 0, 0],
-                [0, 3, 0, 0],
-            ]
-        )
-        dfunc_dy = func.jac(y).evaluate(y=y0)
-        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
-
-        # test jac of outer if left evaluates to number
-        func = pybamm.Outer(pybamm.Scalar(1), pybamm.Scalar(4))
-        dfunc_dy = func.jac(y).evaluate(y=y0)
-        np.testing.assert_array_equal(0, dfunc_dy.toarray())
-
     def test_nonlinear(self):
         y = pybamm.StateVector(slice(0, 4))
         u = pybamm.StateVector(slice(0, 2))
diff --git a/tests/unit/test_expression_tree/test_operations/test_jac_2D.py b/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
index 0d3679715..d08e2b313 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
@@ -86,40 +86,6 @@ class TestJacobian(unittest.TestCase):
         with self.assertRaises(NotImplementedError):
             u.jac(v)
 
-        # test Jacobian of Outer (must set domain to be 'current collector')
-        u.domain = ["current collector"]
-        func = pybamm.Outer(u, pybamm.Scalar(4))
-        jacobian = np.array(
-            [
-                [4, 0, 0, 0, 0, 0, 0, 0],
-                [0, 4, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, 4, 0, 0, 0],
-                [0, 0, 0, 0, 0, 4, 0, 0],
-            ]
-        )
-        dfunc_dy = func.jac(y).evaluate(y=y0)
-        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
-
-        func = pybamm.Outer(u, pybamm.Vector(np.array([1, 2, 3])))
-        jacobian = np.array(
-            [
-                [1, 0, 0, 0, 0, 0, 0, 0],
-                [2, 0, 0, 0, 0, 0, 0, 0],
-                [3, 0, 0, 0, 0, 0, 0, 0],
-                [0, 1, 0, 0, 0, 0, 0, 0],
-                [0, 2, 0, 0, 0, 0, 0, 0],
-                [0, 3, 0, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, 1, 0, 0, 0],
-                [0, 0, 0, 0, 2, 0, 0, 0],
-                [0, 0, 0, 0, 3, 0, 0, 0],
-                [0, 0, 0, 0, 0, 1, 0, 0],
-                [0, 0, 0, 0, 0, 2, 0, 0],
-                [0, 0, 0, 0, 0, 3, 0, 0],
-            ]
-        )
-        dfunc_dy = func.jac(y).evaluate(y=y0)
-        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
-
     def test_nonlinear(self):
         y = pybamm.StateVector(slice(0, 8))
         u = pybamm.StateVector(slice(0, 2), slice(4, 6))
@@ -251,21 +217,16 @@ class TestJacobian(unittest.TestCase):
         a_dom = ["negative electrode"]
         b_dom = ["separator"]
         c_dom = ["positive electrode"]
-        a_npts = mesh[a_dom[0]][0].npts
-        b_npts = mesh[b_dom[0]][0].npts
-        c_npts = mesh[c_dom[0]][0].npts
         cc_npts = mesh["current collector"][0].npts
         curr_coll_vector = pybamm.Vector(np.ones(cc_npts), domain="current collector")
-        a = 2 * pybamm.Outer(
-            curr_coll_vector, pybamm.Vector(np.ones(a_npts), domain=a_dom)
-        )
-        b = pybamm.Outer(curr_coll_vector, pybamm.Vector(np.ones(b_npts), domain=b_dom))
-        c = 3 * pybamm.Outer(
-            curr_coll_vector, pybamm.Vector(np.ones(c_npts), domain=c_dom)
-        )
+        a = 2 * pybamm.PrimaryBroadcast(curr_coll_vector, a_dom)
+        b = pybamm.PrimaryBroadcast(curr_coll_vector, b_dom)
+        c = 3 * pybamm.PrimaryBroadcast(curr_coll_vector, c_dom)
 
-        conc = pybamm.DomainConcatenation([a, b, c], mesh)
-        jac = conc.jac(y).evaluate().toarray()
+        conc = pybamm.Concatenation(a, b, c)
+        disc.set_variable_slices([conc])
+        conc_disc = disc.process_symbol(conc)
+        jac = conc_disc.jac(y).evaluate().toarray()
         np.testing.assert_array_equal(jac, np.zeros((1500, 1500)))
 
         # Jacobian of a DomainConcatenation of StateVectors
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 9625a090d..9c47da440 100644
--- a/tests/unit/test_expression_tree/test_operations/test_simplify.py
+++ b/tests/unit/test_expression_tree/test_operations/test_simplify.py
@@ -570,40 +570,6 @@ class TestSimplify(unittest.TestCase):
         broad_simp = broad.simplify()
         self.assertEqual(broad_simp.id, broad.id)
 
-    def test_simplify_outer(self):
-        v = pybamm.Vector(np.ones(5), domain="current collector")
-        w = pybamm.Vector(2 * np.ones(3), domain="test")
-        outer_simp = pybamm.Outer(v, w).simplify()
-        self.assertIsInstance(outer_simp, pybamm.Vector)
-        np.testing.assert_array_equal(outer_simp.evaluate(), 2 * np.ones((15, 1)))
-
-    def test_simplify_divide_outer(self):
-        u = pybamm.Scalar(1)
-        v = pybamm.StateVector(slice(0, 5), domain="current collector")
-        outer = pybamm.Outer(v, u)
-
-        exp1 = pybamm.Division(pybamm.Division(outer, u), u)
-        self.assertIsInstance(exp1.simplify(), pybamm.Outer)
-
-        exp2 = pybamm.Division(pybamm.Division(outer, 2 * u), u)
-        self.assertIsInstance(exp2.simplify(), pybamm.Multiplication)
-
-        exp3 = pybamm.Division(pybamm.Division(outer, u), 2 * u)
-        self.assertIsInstance(exp3.simplify(), pybamm.Multiplication)
-
-        exp4 = pybamm.Division(pybamm.Division(outer, 2 * u), 2 * u)
-        self.assertIsInstance(exp4.simplify(), pybamm.Multiplication)
-
-    def test_simplify_kron(self):
-        A = pybamm.Matrix(np.eye(2))
-        b = pybamm.Vector(np.array([[4], [5]]))
-        kron = pybamm.Kron(A, b)
-        kron_simp = kron.simplify()
-        self.assertIsInstance(kron_simp, pybamm.Matrix)
-        np.testing.assert_array_equal(
-            kron_simp.evaluate().toarray(), np.kron(A.entries, b.entries)
-        )
-
     def test_simplify_heaviside(self):
         a = pybamm.Scalar(1)
         b = pybamm.Scalar(2)
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index 75bc3a4a7..a41b3e5ea 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -235,12 +235,6 @@ class TestParameterValues(unittest.TestCase):
             processed_g.evaluate(y=np.ones(10)), np.ones((10, 1))
         )
 
-        # process outer
-        c = pybamm.Parameter("c", domain="current collector")
-        outer = pybamm.Outer(c, b)
-        processed_outer = parameter_values.process_symbol(outer)
-        self.assertIsInstance(processed_outer, pybamm.Outer)
-
         # not implemented
         sym = pybamm.Symbol("sym")
         with self.assertRaises(NotImplementedError):