diff --git a/.github/workflows/test_on_push.yml b/.github/workflows/test_on_push.yml
index 962c4b6ad6..385914bc3a 100644
--- a/.github/workflows/test_on_push.yml
+++ b/.github/workflows/test_on_push.yml
@@ -40,7 +40,7 @@ jobs:
       uses: actions/setup-python@v1
       with:
         python-version: ${{ matrix.python-version }}
-    
+  
     - name: Install Linux system dependencies
       if: matrix.os == 'ubuntu-latest'
       run: |
diff --git a/.gitignore b/.gitignore
index ea13d52adc..adf3743680 100644
--- a/.gitignore
+++ b/.gitignore
@@ -110,8 +110,5 @@ test.c
 # tox
 .tox/
 
-# julia
-Manifest.toml
-
 # vcpkg
-vcpkg_installed/
\ No newline at end of file
+vcpkg_installed/
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 76f036cef7..60fb615ad1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added `print_name` to some symbols ([#1495](https://github.com/pybamm-team/PyBaMM/pull/1497), [#1495](https://github.com/pybamm-team/PyBaMM/pull/1497))
 -   Added Base Parameters class and SymPy in dependencies ([#1495](https://github.com/pybamm-team/PyBaMM/pull/1495))
 -   Added a new "reaction-driven" model for LAM from Reniers et al (2019) ([#1490](https://github.com/pybamm-team/PyBaMM/pull/1490))
 -   Some features ("loss of active material" and "particle mechanics") can now be specified separately for the negative electrode and positive electrode by passing a 2-tuple ([#1490](https://github.com/pybamm-team/PyBaMM/pull/1490))
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index cdf440fced..e5f476195a 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -95,7 +95,7 @@ def version(formatted=False):
 from .expression_tree.exceptions import *
 
 # Operations
-from .expression_tree.operations.evaluate import (
+from .expression_tree.operations.evaluate_python import (
     find_symbols,
     id_to_python_variable,
     to_python,
@@ -103,8 +103,8 @@ def version(formatted=False):
 )
 
 if system() != "Windows":
-    from .expression_tree.operations.evaluate import EvaluatorJax
-    from .expression_tree.operations.evaluate import JaxCooMatrix
+    from .expression_tree.operations.evaluate_python import EvaluatorJax
+    from .expression_tree.operations.evaluate_python import JaxCooMatrix
 
 from .expression_tree.operations.jacobian import Jacobian
 from .expression_tree.operations.convert_to_casadi import CasadiConverter
@@ -115,7 +115,7 @@ def version(formatted=False):
 # Model classes
 #
 from .models.base_model import BaseModel
-from .models import standard_variables
+from .models.standard_variables import standard_variables
 from .models.event import Event
 from .models.event import EventType
 
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index e9ac8152ca..2d72b0efcc 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -349,6 +349,16 @@ def __init__(self, *children):
             np.min([child.bounds[1] for child in children]),
         )
 
+        if not any(c.print_name is None for c in children):
+            print_name = intersect(children[0].print_name, children[1].print_name)
+            for child in children[2:]:
+                print_name = intersect(print_name, child.print_name)
+            if print_name.endswith("_"):
+                print_name = print_name[:-1]
+        else:
+            print_name = None
+        self.print_name = print_name
+
 
 def substrings(s):
     for i in range(len(s)):
diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index 4c4c5269f4..e1da19b902 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -358,12 +358,26 @@ def log10(child):
     return log(child, base=10)
 
 
+class Max(SpecificFunction):
+    """ Max function """
+
+    def __init__(self, child):
+        super().__init__(np.max, child)
+
+
 def max(child):
     """
     Returns max function of child. Not to be confused with :meth:`pybamm.maximum`, which
     returns the larger of two objects.
     """
-    return pybamm.simplify_if_constant(Function(np.max, child))
+    return pybamm.simplify_if_constant(Max(child))
+
+
+class Min(SpecificFunction):
+    """ Min function """
+
+    def __init__(self, child):
+        super().__init__(np.min, child)
 
 
 def min(child):
@@ -371,7 +385,7 @@ def min(child):
     Returns min function of child. Not to be confused with :meth:`pybamm.minimum`, which
     returns the smaller of two objects.
     """
-    return pybamm.simplify_if_constant(Function(np.min, child))
+    return pybamm.simplify_if_constant(Min(child))
 
 
 def sech(child):
diff --git a/pybamm/expression_tree/operations/evaluate.py b/pybamm/expression_tree/operations/evaluate_python.py
similarity index 98%
rename from pybamm/expression_tree/operations/evaluate.py
rename to pybamm/expression_tree/operations/evaluate_python.py
index 08b40eb605..638dc33418 100644
--- a/pybamm/expression_tree/operations/evaluate.py
+++ b/pybamm/expression_tree/operations/evaluate_python.py
@@ -313,23 +313,23 @@ def find_symbols(symbol, constant_symbols, variable_symbols, output_jax=False):
 
     elif isinstance(symbol, pybamm.Concatenation):
 
-        # don't bother to concatenate if there is only a single child
+        # no need to concatenate if there is only a single child
         if isinstance(symbol, pybamm.NumpyConcatenation):
-            if len(children_vars) > 1:
-                symbol_str = "np.concatenate(({}))".format(",".join(children_vars))
+            if len(children_vars) == 1:
+                symbol_str = children_vars[0]
             else:
-                symbol_str = "{}".format(",".join(children_vars))
+                symbol_str = "np.concatenate(({}))".format(",".join(children_vars))
 
         elif isinstance(symbol, pybamm.SparseStack):
-            if len(children_vars) > 1:
+            if len(children_vars) == 1:
+                symbol_str = children_vars[0]
+            else:
                 if output_jax:
                     raise NotImplementedError
                 else:
                     symbol_str = "scipy.sparse.vstack(({}))".format(
                         ",".join(children_vars)
                     )
-            else:
-                symbol_str = "{}".format(",".join(children_vars))
 
         # DomainConcatenation specifies a particular ordering for the concatenation,
         # which we must follow
@@ -376,7 +376,9 @@ def find_symbols(symbol, constant_symbols, variable_symbols, output_jax=False):
 
     else:
         raise NotImplementedError(
-            "Not implemented for a symbol of type '{}'".format(type(symbol))
+            "Conversion to python not implemented for a symbol of type '{}'".format(
+                type(symbol)
+            )
         )
 
     variable_symbols[symbol.id] = symbol_str
diff --git a/pybamm/expression_tree/parameter.py b/pybamm/expression_tree/parameter.py
index f38b03e7fc..d7b885e9ae 100644
--- a/pybamm/expression_tree/parameter.py
+++ b/pybamm/expression_tree/parameter.py
@@ -4,6 +4,7 @@
 import numbers
 import numpy as np
 import pybamm
+import sys
 
 
 class Parameter(pybamm.Symbol):
@@ -62,6 +63,9 @@ class FunctionParameter(pybamm.Symbol):
         if diff_variable is specified, the FunctionParameter node will be replaced by a
         :class:`pybamm.Function` and then differentiated with respect to diff_variable.
         Default is None.
+    print_name : str, optional
+        The name to show when printing. Default is 'calculate', in which case the name
+        is calculated using sys._getframe().
     """
 
     def __init__(
@@ -69,6 +73,7 @@ def __init__(
         name,
         inputs,
         diff_variable=None,
+        print_name="calculate",
     ):
         # assign diff variable
         self.diff_variable = diff_variable
@@ -90,6 +95,23 @@ def __init__(
 
         self.input_names = list(inputs.keys())
 
+        # Use the inspect module to find the function's "short name" from the
+        # Parameters module that called it
+        if print_name != "calculate":
+            self.print_name = print_name
+        else:
+            frame = sys._getframe().f_back
+            print_name = frame.f_code.co_name
+            if print_name.startswith("_"):
+                self.print_name = None
+            else:
+                if print_name.endswith("_dimensional"):
+                    self.print_name = print_name[: -len("_dimensional")]
+                elif print_name.endswith("_dim"):
+                    self.print_name = print_name[: -len("_dim")]
+                else:
+                    self.print_name = print_name
+
     @property
     def input_names(self):
         return self._input_names
@@ -159,9 +181,14 @@ def diff(self, variable):
 
     def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
-        return self._function_parameter_new_copy(self._input_names, self.orphans)
+        out = self._function_parameter_new_copy(
+            self._input_names, self.orphans, print_name=self.print_name
+        )
+        return out
 
-    def _function_parameter_new_copy(self, input_names, children):
+    def _function_parameter_new_copy(
+        self, input_names, children, print_name="calculate"
+    ):
         """Returns a new copy of the function parameter.
 
         Inputs
@@ -173,14 +200,17 @@ def _function_parameter_new_copy(self, input_names, children):
 
         Returns
         -------
-            : :pybamm.FunctionParameter
+        :class:`pybamm.FunctionParameter`
             A new copy of the function parameter
         """
 
         input_dict = {input_names[i]: children[i] for i in range(len(input_names))}
 
         return FunctionParameter(
-            self.name, input_dict, diff_variable=self.diff_variable
+            self.name,
+            input_dict,
+            diff_variable=self.diff_variable,
+            print_name=print_name,
         )
 
     def _evaluate_for_shape(self):
diff --git a/pybamm/expression_tree/scalar.py b/pybamm/expression_tree/scalar.py
index 5bde1f147e..333844fcdb 100644
--- a/pybamm/expression_tree/scalar.py
+++ b/pybamm/expression_tree/scalar.py
@@ -31,6 +31,9 @@ def __init__(self, value, name=None, domain=[]):
 
         super().__init__(name, domain=domain)
 
+    def __str__(self):
+        return str(self.value)
+
     @property
     def value(self):
         """the value returned by the node when evaluated"""
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index ee54d7eda8..c38dd50c63 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -47,13 +47,16 @@ def __init__(self, name, domain=None, auxiliary_domains=None, bounds=None):
                     + "Lower bound should be strictly less than upper bound."
                 )
         self.bounds = bounds
+        self.print_name = None
 
     def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
 
-        return self.__class__(
+        out = self.__class__(
             self.name, self.domain, self.auxiliary_domains, self.bounds
         )
+        out.print_name = self.print_name
+        return out
 
     def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 2a0b140f0c..36bd5b862e 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -774,6 +774,22 @@ def info(self, symbol_name):
 
         print(div)
 
+    def check_discretised_or_discretise_inplace_if_0D(self):
+        """
+        Discretise model if it isn't already discretised
+        This only works with purely 0D models, as otherwise the mesh and spatial
+        method should be specified by the user
+        """
+        if self.is_discretised is False:
+            try:
+                disc = pybamm.Discretisation()
+                disc.process_model(self)
+            except pybamm.DiscretisationError as e:
+                raise pybamm.DiscretisationError(
+                    "Cannot automatically discretise model, model should be "
+                    "discretised before exporting casadi functions ({})".format(e)
+                )
+
     def export_casadi_objects(self, variable_names, input_parameter_order=None):
         """
         Export the constituent parts of the model (rhs, algebraic, initial conditions,
@@ -793,18 +809,7 @@ def export_casadi_objects(self, variable_names, input_parameter_order=None):
             Dictionary of {str: casadi object} pairs representing the model in casadi
             format
         """
-        # Discretise model if it isn't already discretised
-        # This only works with purely 0D models, as otherwise the mesh and spatial
-        # method should be specified by the user
-        if self.is_discretised is False:
-            try:
-                disc = pybamm.Discretisation()
-                disc.process_model(self)
-            except pybamm.DiscretisationError as e:
-                raise pybamm.DiscretisationError(
-                    "Cannot automatically discretise model, model should be "
-                    "discretised before exporting casadi functions ({})".format(e)
-                )
+        self.check_discretised_or_discretise_inplace_if_0D()
 
         # Create casadi functions for the model
         t_casadi = casadi.MX.sym("t")
diff --git a/pybamm/models/standard_variables.py b/pybamm/models/standard_variables.py
index cb05482b4b..ea1e4e4391 100644
--- a/pybamm/models/standard_variables.py
+++ b/pybamm/models/standard_variables.py
@@ -4,306 +4,316 @@
 import pybamm
 import numpy as np
 
-# Electrolyte concentration
-c_e_n = pybamm.Variable(
-    "Negative electrolyte concentration",
-    domain="negative electrode",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, np.inf),
-)
-c_e_s = pybamm.Variable(
-    "Separator electrolyte concentration",
-    domain="separator",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, np.inf),
-)
-c_e_p = pybamm.Variable(
-    "Positive electrolyte concentration",
-    domain="positive electrode",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, np.inf),
-)
-c_e = pybamm.concatenation(c_e_n, c_e_s, c_e_p)
 
-c_e_av = pybamm.Variable(
-    "X-averaged electrolyte concentration",
-    domain="current collector",
-    bounds=(0, np.inf),
-)
+class StandardVariables:
+    def __init__(self):
+        # Discharge capacity
+        self.Q = pybamm.Variable("Discharge capacity [A.h]")
 
-# Electrolyte porosity times concentration
-eps_c_e_n = pybamm.Variable(
-    "Negative electrode porosity times concentration",
-    domain="negative electrode",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, np.inf),
-)
-eps_c_e_s = pybamm.Variable(
-    "Separator porosity times concentration",
-    domain="separator",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, np.inf),
-)
-eps_c_e_p = pybamm.Variable(
-    "Positive electrode porosity times concentration",
-    domain="positive electrode",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, np.inf),
-)
-eps_c_e = pybamm.concatenation(eps_c_e_n, eps_c_e_s, eps_c_e_p)
+        # Electrolyte concentration
+        self.c_e_n = pybamm.Variable(
+            "Negative electrolyte concentration",
+            domain="negative electrode",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, np.inf),
+        )
+        self.c_e_s = pybamm.Variable(
+            "Separator electrolyte concentration",
+            domain="separator",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, np.inf),
+        )
+        self.c_e_p = pybamm.Variable(
+            "Positive electrolyte concentration",
+            domain="positive electrode",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, np.inf),
+        )
 
-eps_c_e_av = pybamm.Variable(
-    "X-averaged porosity times concentration",
-    domain="current collector",
-    bounds=(0, np.inf),
-)
+        self.c_e_av = pybamm.Variable(
+            "X-averaged electrolyte concentration",
+            domain="current collector",
+            bounds=(0, np.inf),
+        )
 
-# Electrolyte potential
-phi_e_n = pybamm.Variable(
-    "Negative electrolyte potential",
-    domain="negative electrode",
-    auxiliary_domains={"secondary": "current collector"},
-)
-phi_e_s = pybamm.Variable(
-    "Separator electrolyte potential",
-    domain="separator",
-    auxiliary_domains={"secondary": "current collector"},
-)
-phi_e_p = pybamm.Variable(
-    "Positive electrolyte potential",
-    domain="positive electrode",
-    auxiliary_domains={"secondary": "current collector"},
-)
-phi_e = pybamm.concatenation(phi_e_n, phi_e_s, phi_e_p)
+        # Electrolyte porosity times concentration
+        self.eps_c_e_n = pybamm.Variable(
+            "Negative electrode porosity times concentration",
+            domain="negative electrode",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, np.inf),
+        )
+        self.eps_c_e_s = pybamm.Variable(
+            "Separator porosity times concentration",
+            domain="separator",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, np.inf),
+        )
+        self.eps_c_e_p = pybamm.Variable(
+            "Positive electrode porosity times concentration",
+            domain="positive electrode",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, np.inf),
+        )
 
-# Electrode potential
-phi_s_n = pybamm.Variable(
-    "Negative electrode potential",
-    domain="negative electrode",
-    auxiliary_domains={"secondary": "current collector"},
-)
-phi_s_p = pybamm.Variable(
-    "Positive electrode potential",
-    domain="positive electrode",
-    auxiliary_domains={"secondary": "current collector"},
-)
+        self.eps_c_e_av = pybamm.Variable(
+            "X-averaged porosity times concentration",
+            domain="current collector",
+            bounds=(0, np.inf),
+        )
 
-# Potential difference
-delta_phi_n = pybamm.Variable(
-    "Negative electrode surface potential difference",
-    domain="negative electrode",
-    auxiliary_domains={"secondary": "current collector"},
-)
-delta_phi_p = pybamm.Variable(
-    "Positive electrode surface potential difference",
-    domain="positive electrode",
-    auxiliary_domains={"secondary": "current collector"},
-)
+        # Electrolyte potential
+        self.phi_e_n = pybamm.Variable(
+            "Negative electrolyte potential",
+            domain="negative electrode",
+            auxiliary_domains={"secondary": "current collector"},
+        )
+        self.phi_e_s = pybamm.Variable(
+            "Separator electrolyte potential",
+            domain="separator",
+            auxiliary_domains={"secondary": "current collector"},
+        )
+        self.phi_e_p = pybamm.Variable(
+            "Positive electrolyte potential",
+            domain="positive electrode",
+            auxiliary_domains={"secondary": "current collector"},
+        )
 
-delta_phi_n_av = pybamm.Variable(
-    "X-averaged negative electrode surface potential difference",
-    domain="current collector",
-)
-delta_phi_p_av = pybamm.Variable(
-    "X-averaged positive electrode surface potential difference",
-    domain="current collector",
-)
+        # Electrode potential
+        self.phi_s_n = pybamm.Variable(
+            "Negative electrode potential",
+            domain="negative electrode",
+            auxiliary_domains={"secondary": "current collector"},
+        )
+        self.phi_s_p = pybamm.Variable(
+            "Positive electrode potential",
+            domain="positive electrode",
+            auxiliary_domains={"secondary": "current collector"},
+        )
 
-# current collector variables
-phi_s_cn = pybamm.Variable(
-    "Negative current collector potential", domain="current collector"
-)
-phi_s_cp = pybamm.Variable(
-    "Positive current collector potential", domain="current collector"
-)
-i_boundary_cc = pybamm.Variable(
-    "Current collector current density", domain="current collector"
-)
-phi_s_cn_composite = pybamm.Variable(
-    "Composite negative current collector potential", domain="current collector"
-)
-phi_s_cp_composite = pybamm.Variable(
-    "Composite positive current collector potential", domain="current collector"
-)
-i_boundary_cc_composite = pybamm.Variable(
-    "Composite current collector current density", domain="current collector"
-)
+        # Potential difference
+        self.delta_phi_n = pybamm.Variable(
+            "Negative electrode surface potential difference",
+            domain="negative electrode",
+            auxiliary_domains={"secondary": "current collector"},
+        )
+        self.delta_phi_p = pybamm.Variable(
+            "Positive electrode surface potential difference",
+            domain="positive electrode",
+            auxiliary_domains={"secondary": "current collector"},
+        )
 
+        self.delta_phi_n_av = pybamm.Variable(
+            "X-averaged negative electrode surface potential difference",
+            domain="current collector",
+        )
+        self.delta_phi_p_av = pybamm.Variable(
+            "X-averaged positive electrode surface potential difference",
+            domain="current collector",
+        )
 
-# Particle concentration
-c_s_n = pybamm.Variable(
-    "Negative particle concentration",
-    domain="negative particle",
-    auxiliary_domains={
-        "secondary": "negative electrode",
-        "tertiary": "current collector",
-    },
-    bounds=(0, 1),
-)
-c_s_p = pybamm.Variable(
-    "Positive particle concentration",
-    domain="positive particle",
-    auxiliary_domains={
-        "secondary": "positive electrode",
-        "tertiary": "current collector",
-    },
-    bounds=(0, 1),
-)
-c_s_n_xav = pybamm.Variable(
-    "X-averaged negative particle concentration",
-    domain="negative particle",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, 1),
-)
-c_s_p_xav = pybamm.Variable(
-    "X-averaged positive particle concentration",
-    domain="positive particle",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, 1),
-)
-c_s_n_rav = pybamm.Variable(
-    "R-averaged negative particle concentration",
-    domain="negative electrode",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, 1),
-)
-c_s_p_rav = pybamm.Variable(
-    "R-averaged positive particle concentration",
-    domain="positive electrode",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, 1),
-)
-c_s_n_rxav = pybamm.Variable(
-    "R-X-averaged negative particle concentration",
-    domain="current collector",
-    bounds=(0, 1),
-)
-c_s_p_rxav = pybamm.Variable(
-    "R-X-averaged positive particle concentration",
-    domain="current collector",
-    bounds=(0, 1),
-)
-c_s_n_surf = pybamm.Variable(
-    "Negative particle surface concentration",
-    domain="negative electrode",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, 1),
-)
-c_s_p_surf = pybamm.Variable(
-    "Positive particle surface concentration",
-    domain="positive electrode",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, 1),
-)
-c_s_n_surf_xav = pybamm.Variable(
-    "X-averaged negative particle surface concentration",
-    domain="current collector",
-    bounds=(0, 1),
-)
-c_s_p_surf_xav = pybamm.Variable(
-    "X-averaged positive particle surface concentration",
-    domain="current collector",
-    bounds=(0, 1),
-)
-# Average particle concentration gradient (for polynomial particle concentration
-# models). Note: we make the distinction here between the flux defined as
-# N = -D*dc/dr and the concentration gradient q = dc/dr
-q_s_n_rav = pybamm.Variable(
-    "R-averaged negative particle concentration gradient",
-    domain="negative electrode",
-    auxiliary_domains={"secondary": "current collector"},
-)
-q_s_p_rav = pybamm.Variable(
-    "R-averaged positive particle concentration gradient",
-    domain="positive electrode",
-    auxiliary_domains={"secondary": "current collector"},
-)
-q_s_n_rxav = pybamm.Variable(
-    "R-X-averaged negative particle concentration gradient", domain="current collector"
-)
-q_s_p_rxav = pybamm.Variable(
-    "R-X-averaged positive particle concentration gradient", domain="current collector"
-)
+        # current collector variables
+        self.phi_s_cn = pybamm.Variable(
+            "Negative current collector potential", domain="current collector"
+        )
+        self.phi_s_cp = pybamm.Variable(
+            "Positive current collector potential", domain="current collector"
+        )
+        self.i_boundary_cc = pybamm.Variable(
+            "Current collector current density", domain="current collector"
+        )
+        self.phi_s_cn_composite = pybamm.Variable(
+            "Composite negative current collector potential", domain="current collector"
+        )
+        self.phi_s_cp_composite = pybamm.Variable(
+            "Composite positive current collector potential", domain="current collector"
+        )
+        self.i_boundary_cc_composite = pybamm.Variable(
+            "Composite current collector current density", domain="current collector"
+        )
 
-# Porosity
-eps_n = pybamm.Variable(
-    "Negative electrode porosity",
-    domain="negative electrode",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, 1),
-)
-eps_s = pybamm.Variable(
-    "Separator porosity",
-    domain="separator",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, 1),
-)
-eps_p = pybamm.Variable(
-    "Positive electrode porosity",
-    domain="positive electrode",
-    auxiliary_domains={"secondary": "current collector"},
-    bounds=(0, 1),
-)
-eps = pybamm.concatenation(eps_n, eps_s, eps_p)
+        # Particle concentration
+        self.c_s_n = pybamm.Variable(
+            "Negative particle concentration",
+            domain="negative particle",
+            auxiliary_domains={
+                "secondary": "negative electrode",
+                "tertiary": "current collector",
+            },
+            bounds=(0, 1),
+        )
+        self.c_s_p = pybamm.Variable(
+            "Positive particle concentration",
+            domain="positive particle",
+            auxiliary_domains={
+                "secondary": "positive electrode",
+                "tertiary": "current collector",
+            },
+            bounds=(0, 1),
+        )
+        self.c_s_n_xav = pybamm.Variable(
+            "X-averaged negative particle concentration",
+            domain="negative particle",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, 1),
+        )
+        self.c_s_p_xav = pybamm.Variable(
+            "X-averaged positive particle concentration",
+            domain="positive particle",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, 1),
+        )
+        self.c_s_n_rav = pybamm.Variable(
+            "R-averaged negative particle concentration",
+            domain="negative electrode",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, 1),
+        )
+        self.c_s_p_rav = pybamm.Variable(
+            "R-averaged positive particle concentration",
+            domain="positive electrode",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, 1),
+        )
+        self.c_s_n_rxav = pybamm.Variable(
+            "R-X-averaged negative particle concentration",
+            domain="current collector",
+            bounds=(0, 1),
+        )
+        self.c_s_p_rxav = pybamm.Variable(
+            "R-X-averaged positive particle concentration",
+            domain="current collector",
+            bounds=(0, 1),
+        )
+        self.c_s_n_surf = pybamm.Variable(
+            "Negative particle surface concentration",
+            domain="negative electrode",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, 1),
+        )
+        self.c_s_p_surf = pybamm.Variable(
+            "Positive particle surface concentration",
+            domain="positive electrode",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, 1),
+        )
+        self.c_s_n_surf_xav = pybamm.Variable(
+            "X-averaged negative particle surface concentration",
+            domain="current collector",
+            bounds=(0, 1),
+        )
+        self.c_s_p_surf_xav = pybamm.Variable(
+            "X-averaged positive particle surface concentration",
+            domain="current collector",
+            bounds=(0, 1),
+        )
+        # Average particle concentration gradient (for polynomial particle concentration
+        # models). Note: we make the distinction here between the flux defined as
+        # N = -D*dc/dr and the concentration gradient q = dc/dr
+        self.q_s_n_rav = pybamm.Variable(
+            "R-averaged negative particle concentration gradient",
+            domain="negative electrode",
+            auxiliary_domains={"secondary": "current collector"},
+        )
+        self.q_s_p_rav = pybamm.Variable(
+            "R-averaged positive particle concentration gradient",
+            domain="positive electrode",
+            auxiliary_domains={"secondary": "current collector"},
+        )
+        self.q_s_n_rxav = pybamm.Variable(
+            "R-X-averaged negative particle concentration gradient",
+            domain="current collector",
+        )
+        self.q_s_p_rxav = pybamm.Variable(
+            "R-X-averaged positive particle concentration gradient",
+            domain="current collector",
+        )
 
-# Piecewise constant (for asymptotic models)
-eps_n_pc = pybamm.Variable(
-    "X-averaged negative electrode porosity", domain="current collector", bounds=(0, 1)
-)
-eps_s_pc = pybamm.Variable(
-    "X-averaged separator porosity", domain="current collector", bounds=(0, 1)
-)
-eps_p_pc = pybamm.Variable(
-    "X-averaged positive electrode porosity", domain="current collector", bounds=(0, 1)
-)
+        # Porosity
+        self.eps_n = pybamm.Variable(
+            "Negative electrode porosity",
+            domain="negative electrode",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, 1),
+        )
+        self.eps_s = pybamm.Variable(
+            "Separator porosity",
+            domain="separator",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, 1),
+        )
+        self.eps_p = pybamm.Variable(
+            "Positive electrode porosity",
+            domain="positive electrode",
+            auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, 1),
+        )
 
-eps_piecewise_constant = pybamm.concatenation(
-    pybamm.PrimaryBroadcast(eps_n_pc, "negative electrode"),
-    pybamm.PrimaryBroadcast(eps_s_pc, "separator"),
-    pybamm.PrimaryBroadcast(eps_p_pc, "positive electrode"),
-)
+        # Piecewise constant (for asymptotic models)
+        self.eps_n_pc = pybamm.Variable(
+            "X-averaged negative electrode porosity",
+            domain="current collector",
+            bounds=(0, 1),
+        )
+        self.eps_s_pc = pybamm.Variable(
+            "X-averaged separator porosity", domain="current collector", bounds=(0, 1)
+        )
+        self.eps_p_pc = pybamm.Variable(
+            "X-averaged positive electrode porosity",
+            domain="current collector",
+            bounds=(0, 1),
+        )
 
-# Temperature
-T_cn = pybamm.Variable(
-    "Negative currents collector temperature", domain="current collector"
-)
-T_n = pybamm.Variable(
-    "Negative electrode temperature",
-    domain="negative electrode",
-    auxiliary_domains={"secondary": "current collector"},
-)
-T_s = pybamm.Variable(
-    "Separator temperature",
-    domain="separator",
-    auxiliary_domains={"secondary": "current collector"},
-)
-T_p = pybamm.Variable(
-    "Positive electrode temperature",
-    domain="positive electrode",
-    auxiliary_domains={"secondary": "current collector"},
-)
-T_cp = pybamm.Variable(
-    "Positive currents collector temperature", domain="current collector"
-)
-T = pybamm.concatenation(T_n, T_s, T_p)
-T_av = pybamm.Variable("X-averaged cell temperature", domain="current collector")
-T_vol_av = pybamm.Variable("Volume-averaged cell temperature")
+        # Temperature
+        self.T_cn = pybamm.Variable(
+            "Negative currents collector temperature", domain="current collector"
+        )
+        self.T_n = pybamm.Variable(
+            "Negative electrode temperature",
+            domain="negative electrode",
+            auxiliary_domains={"secondary": "current collector"},
+        )
+        self.T_s = pybamm.Variable(
+            "Separator temperature",
+            domain="separator",
+            auxiliary_domains={"secondary": "current collector"},
+        )
+        self.T_p = pybamm.Variable(
+            "Positive electrode temperature",
+            domain="positive electrode",
+            auxiliary_domains={"secondary": "current collector"},
+        )
+        self.T_cp = pybamm.Variable(
+            "Positive currents collector temperature", domain="current collector"
+        )
+        self.T_av = pybamm.Variable(
+            "X-averaged cell temperature", domain="current collector"
+        )
+        self.T_vol_av = pybamm.Variable("Volume-averaged cell temperature")
 
+        # SEI variables
+        self.L_inner_av = pybamm.Variable(
+            "X-averaged inner negative electrode SEI thickness",
+            domain="current collector",
+        )
+        self.L_inner = pybamm.Variable(
+            "Inner negative electrode SEI thickness",
+            domain=["negative electrode"],
+            auxiliary_domains={"secondary": "current collector"},
+        )
+        self.L_outer_av = pybamm.Variable(
+            "X-averaged outer negative electrode SEI thickness",
+            domain="current collector",
+        )
+        self.L_outer = pybamm.Variable(
+            "Outer negative electrode SEI thickness",
+            domain=["negative electrode"],
+            auxiliary_domains={"secondary": "current collector"},
+        )
 
-# SEI variables
-L_inner_av = pybamm.Variable(
-    "X-averaged inner negative electrode SEI thickness", domain="current collector"
-)
-L_inner = pybamm.Variable(
-    "Inner negative electrode SEI thickness",
-    domain=["negative electrode"],
-    auxiliary_domains={"secondary": "current collector"},
-)
-L_outer_av = pybamm.Variable(
-    "X-averaged outer negative electrode SEI thickness", domain="current collector"
-)
-L_outer = pybamm.Variable(
-    "Outer negative electrode SEI thickness",
-    domain=["negative electrode"],
-    auxiliary_domains={"secondary": "current collector"},
-)
+    def __setattr__(self, name, value):
+        value.print_name = name
+        super().__setattr__(name, value)
+
+
+standard_variables = StandardVariables()
diff --git a/pybamm/models/submodels/external_circuit/base_external_circuit.py b/pybamm/models/submodels/external_circuit/base_external_circuit.py
index cd37a82a6f..95df5b7582 100644
--- a/pybamm/models/submodels/external_circuit/base_external_circuit.py
+++ b/pybamm/models/submodels/external_circuit/base_external_circuit.py
@@ -11,7 +11,7 @@ def __init__(self, param):
         super().__init__(param)
 
     def get_fundamental_variables(self):
-        Q = pybamm.Variable("Discharge capacity [A.h]")
+        Q = pybamm.standard_variables.Q
         variables = {"Discharge capacity [A.h]": Q}
         return variables
 
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index d6023e5e8d..c8b1a7ba77 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -571,6 +571,8 @@ def process_geometry(self, geometry):
                             geometry[domain][spatial_variable][
                                 lim
                             ] = self.process_symbol(sym)
+                        elif isinstance(sym, numbers.Number):
+                            geometry[domain][spatial_variable][lim] = pybamm.Scalar(sym)
 
     def process_symbol(self, symbol):
         """Walk through the symbol and replace any Parameter with a Value.
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index 69bfa77f79..5117e2930b 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -15,7 +15,7 @@ class CasadiSolver(pybamm.BaseSolver):
     Parameters
     ----------
     mode : str
-            How to solve the model (default is "safe"):
+        How to solve the model (default is "safe"):
 
             - "fast": perform direct integration, without accounting for events. \
             Recommended when simulating a drive cycle or other simulation where \
diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py
index e73f7ea4dd..f6b8362702 100644
--- a/pybamm/solvers/idaklu_solver.py
+++ b/pybamm/solvers/idaklu_solver.py
@@ -10,8 +10,11 @@
 
 idaklu_spec = importlib.util.find_spec("pybamm.solvers.idaklu")
 if idaklu_spec is not None:
-    idaklu = importlib.util.module_from_spec(idaklu_spec)
-    idaklu_spec.loader.exec_module(idaklu)
+    try:
+        idaklu = importlib.util.module_from_spec(idaklu_spec)
+        idaklu_spec.loader.exec_module(idaklu)
+    except ImportError:  # pragma: no cover
+        idaklu_spec = None
 
 
 def have_idaklu():
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 9b1fcd8faa..2015a6e841 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -815,13 +815,12 @@ def add_neumann_values(self, symbol, discretised_gradient, bcs, domain):
         # Need to match the domain. E.g. in the case of the boundary condition
         # on the particle, the gradient has domain particle but the bcs_vector
         # has domain electrode, since it is a function of the macroscopic variables
-        bcs_vector.domain = discretised_gradient.domain
-        bcs_vector.auxiliary_domains = discretised_gradient.auxiliary_domains
+        bcs_vector.copy_domains(discretised_gradient)
 
         # Make matrix which makes "gaps" in the the discretised gradient into
         # which the known Neumann values will be added. E.g. in 1D if the left
         # boundary condition is Dirichlet and the right Neumann, this matrix will
-        # act to append a zero to the end of the discretsied gradient
+        # act to append a zero to the end of the discretised gradient
         if lbc_type == "Neumann":
             left_vector = csr_matrix((1, n))
         else:
diff --git a/pybamm/util.py b/pybamm/util.py
index e407cef78e..27ff4a88cc 100644
--- a/pybamm/util.py
+++ b/pybamm/util.py
@@ -276,8 +276,9 @@ def load_function(filename):
 
 
 def rmse(x, y):
-    """Calculate the root-mean-square-error between
-    two vectors x and y, ignoring NaNs"""
+    """
+    Calculate the root-mean-square-error between two vectors x and y, ignoring NaNs
+    """
     # Check lengths
     if len(x) != len(y):
         raise ValueError("Vectors must have the same length")
diff --git a/tests/unit/test_expression_tree/test_operations/test_evaluate.py b/tests/unit/test_expression_tree/test_operations/test_evaluate_python.py
similarity index 99%
rename from tests/unit/test_expression_tree/test_operations/test_evaluate.py
rename to tests/unit/test_expression_tree/test_operations/test_evaluate_python.py
index 7b30cefb56..3907470234 100644
--- a/tests/unit/test_expression_tree/test_operations/test_evaluate.py
+++ b/tests/unit/test_expression_tree/test_operations/test_evaluate_python.py
@@ -1,5 +1,5 @@
 #
-# Test for the evaluate functions
+# Test for the evaluate-to-python functions
 #
 import pybamm
 
@@ -427,6 +427,10 @@ def test_evaluator_python(self):
         for t, y in zip(t_tests, y_tests):
             result = evaluator.evaluate(t=t, y=y).toarray()
             np.testing.assert_allclose(result, expr.evaluate(t=t, y=y).toarray())
+        expr = pybamm.SparseStack(A)
+        evaluator = pybamm.EvaluatorPython(expr)
+        result = evaluator.evaluate().toarray()
+        np.testing.assert_allclose(result, expr.evaluate().toarray())
 
         # test Inner
         expr = pybamm.Inner(a, b)
diff --git a/tests/unit/test_expression_tree/test_parameter.py b/tests/unit/test_expression_tree/test_parameter.py
index ee3231a0d8..6f51624b55 100644
--- a/tests/unit/test_expression_tree/test_parameter.py
+++ b/tests/unit/test_expression_tree/test_parameter.py
@@ -74,6 +74,25 @@ def test_set_input_names(self):
             new_input_names = [var]
             func.input_names = new_input_names
 
+    def test_print_name(self):
+        def myfun(x):
+            return pybamm.FunctionParameter("my function", {"x": x})
+
+        def myfun_dim(x):
+            return pybamm.FunctionParameter("my function", {"x": x})
+
+        def myfun_dimensional(x):
+            return pybamm.FunctionParameter("my function", {"x": x})
+
+        def _myfun(x):
+            return pybamm.FunctionParameter("my function", {"x": x})
+
+        x = pybamm.Scalar(1)
+        self.assertEqual(myfun(x).print_name, "myfun")
+        self.assertEqual(myfun_dim(x).print_name, "myfun")
+        self.assertEqual(myfun_dimensional(x).print_name, "myfun")
+        self.assertEqual(_myfun(x).print_name, None)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")