Repository: pybamm-team/PyBaMM

Issue:
Quick Plot functionality
Add more functionality to `QuickPlot`:
- [x] User can define labels (otherwise defaults to model names)
- [x] Can provide list of list of variables, which plots several variables on the same graph
- [ ] Adjust x ticks to fit for lithium ion plots
- [x] Automatically adjust label sizes to avoid overlap

Discussion:

Patch:
diff --git a/.gitignore b/.gitignore
index 4c691e285..e567fff73 100644
--- a/.gitignore
+++ b/.gitignore
@@ -45,6 +45,7 @@ dist/
 # coverage
 .coverage
 coverage.xml
+htmlcov/
 
 # black setup file seems to make Travis CI fail
 pyproject.toml
@@ -61,6 +62,3 @@ pyvenv.cfg
 
 # sundials
 sundials
-
-# code coverage
-htmlcov/
diff --git a/examples/scripts/compare_comsol/compare_comsol_DFN.py b/examples/scripts/compare_comsol/compare_comsol_DFN.py
index ba34e461e..0a6b982c9 100644
--- a/examples/scripts/compare_comsol/compare_comsol_DFN.py
+++ b/examples/scripts/compare_comsol/compare_comsol_DFN.py
@@ -2,8 +2,7 @@ import pybamm
 import numpy as np
 import os
 import pickle
-import matplotlib.pyplot as plt
-from matplotlib.widgets import Slider
+import scipy.interpolate as interp
 
 # change working directory to the root of pybamm
 os.chdir(pybamm.__path__[0] + "/..")
@@ -26,234 +25,93 @@ comsol_variables = pickle.load(
 "Create and solve pybamm model"
 
 # load model and geometry
-model = pybamm.lithium_ion.DFN()
-geometry = model.default_geometry
+pybamm.set_logging_level("INFO")
+pybamm_model = pybamm.lithium_ion.DFN()
+geometry = pybamm_model.default_geometry
 
 # load parameters and process model and geometry
-param = model.default_parameter_values
+param = pybamm_model.default_parameter_values
 param["Electrode depth [m]"] = 1
 param["Electrode height [m]"] = 1
 param["Typical current [A]"] = 24 * C_rates[C_rate]
-param.process_model(model)
+param.process_model(pybamm_model)
 param.process_geometry(geometry)
 
 # create mesh
 var = pybamm.standard_spatial_vars
 var_pts = {var.x_n: 31, var.x_s: 11, var.x_p: 31, var.r_n: 11, var.r_p: 11}
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
+mesh = pybamm.Mesh(geometry, pybamm_model.default_submesh_types, var_pts)
 
 # discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
+disc = pybamm.Discretisation(mesh, pybamm_model.default_spatial_methods)
+disc.process_model(pybamm_model)
 
 # discharge timescale
-tau = param.process_symbol(
-    pybamm.standard_parameters_lithium_ion.tau_discharge
-).evaluate(0, 0)
+tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge)
 
 # solve model at comsol times
-time = comsol_variables["time"] / tau
-solution = model.default_solver.solve(model, time)
-
-"-----------------------------------------------------------------------------"
-"Get variables for comparison"
-
-# discharge capacity
-discharge_capacity = pybamm.ProcessedVariable(
-    model.variables["Discharge capacity [A.h]"], solution.t, solution.y, mesh=mesh
-)
-discharge_capacity_sol = discharge_capacity(solution.t)
-comsol_discharge_capacity = (
-    comsol_variables["time"] * param["Typical current [A]"] / 3600
-)
-
-# spatial points
-l_n = param.process_symbol(pybamm.geometric_parameters.l_n).evaluate(0, 0)
-l_s = param.process_symbol(pybamm.geometric_parameters.l_s).evaluate(0, 0)
-l_p = param.process_symbol(pybamm.geometric_parameters.l_p).evaluate(0, 0)
-x_n = np.linspace(0, l_n, 40)
-x_s = np.linspace(l_n, l_n + l_s, 20)
-x_p = np.linspace(l_n + l_s, 1, 40)
-x = np.linspace(0, 1, 100)
-
-# output variables
-output_variables = {
-    "c_n_surf": model.variables["Negative particle surface concentration [mol.m-3]"],
-    "c_e": model.variables["Electrolyte concentration [mol.m-3]"],
-    "c_p_surf": model.variables["Positive particle surface concentration [mol.m-3]"],
-    "phi_n": model.variables["Negative electrode potential [V]"],
-    "phi_e": model.variables["Electrolyte potential [V]"],
-    "phi_p": model.variables["Positive electrode potential [V]"],
-    "voltage": model.variables["Terminal voltage [V]"],
+time = comsol_variables["time"] / tau.evaluate(0)
+solution = pybamm_model.default_solver.solve(pybamm_model, time)
+
+
+# Make Comsol 'model' for comparison
+whole_cell = ["negative electrode", "separator", "positive electrode"]
+comsol_t = comsol_variables["time"]
+L_x = param.process_symbol(pybamm.standard_parameters_lithium_ion.L_x).evaluate()
+
+
+def get_interp_fun(variable, domain):
+    """
+    Create a :class:`pybamm.Function` object using the variable, to allow plotting with
+    :class:`'pybamm.QuickPlot'` (interpolate in space to match edges, and then create
+    function to interpolate in time)
+    """
+    if domain == ["negative electrode"]:
+        comsol_x = comsol_variables["x_n"]
+    elif domain == ["positive electrode"]:
+        comsol_x = comsol_variables["x_p"]
+    elif domain == whole_cell:
+        comsol_x = comsol_variables["x"]
+    # Make sure to use dimensional space
+    pybamm_x = mesh.combine_submeshes(*domain)[0].nodes * L_x
+    variable = interp.interp1d(comsol_x, variable, axis=0)(pybamm_x)
+
+    def myinterp(t):
+        return interp.interp1d(comsol_t, variable)(t)[:, np.newaxis]
+
+    # Make sure to use dimensional time
+    fun = pybamm.Function(myinterp, pybamm.t * tau)
+    fun.domain = domain
+    return fun
+
+
+comsol_c_n_surf = get_interp_fun(comsol_variables["c_n_surf"], ["negative electrode"])
+comsol_c_e = get_interp_fun(comsol_variables["c_e"], whole_cell)
+comsol_c_p_surf = get_interp_fun(comsol_variables["c_p_surf"], ["positive electrode"])
+comsol_phi_n = get_interp_fun(comsol_variables["phi_n"], ["negative electrode"])
+comsol_phi_e = get_interp_fun(comsol_variables["phi_e"], whole_cell)
+comsol_phi_p = get_interp_fun(comsol_variables["phi_p"], ["positive electrode"])
+comsol_voltage = interp.interp1d(comsol_t, comsol_variables["voltage"])
+
+# Create comsol model with dictionary of Matrix variables
+comsol_model = pybamm.BaseModel()
+comsol_model.variables = {
+    "Negative particle surface concentration [mol.m-3]": comsol_c_n_surf,
+    "Electrolyte concentration [mol.m-3]": comsol_c_e,
+    "Positive particle surface concentration [mol.m-3]": comsol_c_p_surf,
+    "Current [A]": pybamm_model.variables["Current [A]"],
+    "Negative electrode potential [V]": comsol_phi_n,
+    "Electrolyte potential [V]": comsol_phi_e,
+    "Positive electrode potential [V]": comsol_phi_p,
+    "Terminal voltage [V]": pybamm.Function(comsol_voltage, pybamm.t * tau),
 }
 
-processed_variables = pybamm.post_process_variables(
-    output_variables, solution.t, solution.y, mesh=mesh
-)
-
-"-----------------------------------------------------------------------------"
-"Make plots"
-
-fig, ax = plt.subplots(figsize=(15, 8))
-plt.tight_layout()
-plt.subplots_adjust(left=-0.1)
-
-
-# negative particle surface concentration
-c_n_surf_min = 0.9 * np.min(comsol_variables["c_n_surf"])
-c_n_surf_max = 1.1 * np.max(comsol_variables["c_n_surf"])
-plt.subplot(241)
-c_n_surf_plot, = plt.plot(x_n, processed_variables["c_n_surf"](time[0], x_n), "b-")
-comsol_c_n_surf_plot, = plt.plot(
-    comsol_variables["x_n"] / comsol_variables["x"][-1],
-    comsol_variables["c_n_surf"][:, 0],
-    "r:",
-)
-plt.axis([0, l_n, c_n_surf_min, c_n_surf_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"Surface $c_n$ (mol/m$^3$)")
-
-# electrolyte concentration
-c_e_min = 0.9 * np.min(comsol_variables["c_e"])
-c_e_max = 1.1 * np.max(comsol_variables["c_e"])
-plt.subplot(242)
-c_e_plot, = plt.plot(x, processed_variables["c_e"](time[0], x), "b-")
-comsol_c_e_plot, = plt.plot(
-    comsol_variables["x"] / comsol_variables["x"][-1],
-    comsol_variables["c_e"][:, 0],
-    "r:",
-)
-plt.axis([0, 1, c_e_min, c_e_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"$c_e$ (mol/m$^3$)")
-
-# negative particle surface concentration
-c_p_surf_min = 0.9 * np.min(comsol_variables["c_p_surf"])
-c_p_surf_max = 1.1 * np.max(comsol_variables["c_p_surf"])
-plt.subplot(243)
-c_p_surf_plot, = plt.plot(x_p, processed_variables["c_n_surf"](time[0], x_p), "b-")
-comsol_c_p_surf_plot, = plt.plot(
-    comsol_variables["x_p"] / comsol_variables["x"][-1],
-    comsol_variables["c_p_surf"][:, 0],
-    "r:",
-)
-plt.axis([l_n + l_s, 1, c_p_surf_min, c_p_surf_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"Surface $c_p$ (mol/m$^3$)")
-
-# current
-I_min = 0.9 * param["Typical current [A]"]
-I_max = 1.1 * param["Typical current [A]"]
-plt.subplot(244)
-time_tracer, = plt.plot(
-    [comsol_variables["time"][0], comsol_variables["time"][0]], [I_min, I_max], "k--"
-)
-plt.plot(
-    [comsol_variables["time"][0], comsol_variables["time"][-1]],
-    [param["Typical current [A]"], param["Typical current [A]"]],
-    "b-",
-)
-plt.axis(
-    [comsol_variables["time"][0], comsol_variables["time"][-1] * 1.1, I_min, I_max]
-)
-plt.xlabel("Time (s)")
-plt.ylabel("Applied current (A)")
-
-# negative electrode potential
-phi_n_min = 1.1 * np.min(comsol_variables["phi_n"])
-phi_n_max = 0.9 * np.max(comsol_variables["phi_n"])
-plt.subplot(245)
-phi_n_plot, = plt.plot(x_n, processed_variables["phi_n"](time[0], x_n), "b-")
-comsol_phi_n_plot, = plt.plot(
-    comsol_variables["x_n"] / comsol_variables["x"][-1],
-    comsol_variables["phi_n"][:, 0],
-    "r:",
-)
-plt.axis([0, l_n, phi_n_min, phi_n_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"$\phi_n$ (V)")
-
-# electrolyte potential
-phi_e_min = 1.1 * np.min(comsol_variables["phi_e"])
-phi_e_max = 0.9 * np.max(comsol_variables["phi_e"])
-plt.subplot(246)
-phi_e_plot, = plt.plot(x, processed_variables["phi_e"](time[0], x), "b-")
-comsol_phi_e_plot, = plt.plot(
-    comsol_variables["x"] / comsol_variables["x"][-1],
-    comsol_variables["phi_e"][:, 0],
-    "r:",
-)
-plt.axis([0, 1, phi_e_min, phi_e_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"$\phi_e$ (V)")
-
-# positive electrode potential
-phi_p_min = 0.9 * np.min(comsol_variables["phi_p"])
-phi_p_max = 1.1 * np.max(comsol_variables["phi_p"])
-plt.subplot(247)
-phi_p_plot, = plt.plot(x_p, processed_variables["phi_p"](time[0], x_p), "b-")
-comsol_phi_p_plot, = plt.plot(
-    comsol_variables["x_p"] / comsol_variables["x"][-1],
-    comsol_variables["phi_p"][:, 0],
-    "r:",
+# plot
+plot = pybamm.QuickPlot(
+    [pybamm_model, comsol_model],
+    mesh,
+    [solution, solution],
+    output_variables=comsol_model.variables.keys(),
+    labels=["PyBaMM", "Comsol"],
 )
-plt.axis([l_n + l_s, 1, phi_p_min, phi_p_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"$\phi_p$ (V)")
-
-# discharge curve
-v_min = 3.2
-v_max = 3.9
-plt.subplot(248)
-discharge_capacity_tracer, = plt.plot(
-    [comsol_discharge_capacity[0], comsol_discharge_capacity[0]], [v_min, v_max], "k--"
-)
-plt.plot(
-    discharge_capacity_sol,
-    processed_variables["voltage"](solution.t),
-    "b-",
-    label="PyBaMM",
-)
-plt.plot(comsol_discharge_capacity, comsol_variables["voltage"], "r:", label="Comsol")
-plt.axis(
-    [comsol_discharge_capacity[0], comsol_discharge_capacity[-1] * 1.1, v_min, v_max]
-)
-plt.xlabel("Discharge Capacity (Ah)")
-plt.ylabel("Voltage (V)")
-plt.legend(loc="best")
-
-axcolor = "lightgoldenrodyellow"
-axfreq = plt.axes([0.315, 0.02, 0.37, 0.03], facecolor=axcolor)
-sfreq = Slider(axfreq, "Time", 0, comsol_variables["time"][-1], valinit=0)
-
-
-def update_plot(t):
-    # find t index
-    ind = (np.abs(comsol_variables["time"] - t)).argmin()
-    # update time
-    time_tracer.set_xdata(comsol_variables["time"][ind])
-    discharge_capacity_tracer.set_xdata(comsol_discharge_capacity[ind])
-    # update negative particle surface concentration
-    c_n_surf_plot.set_ydata(processed_variables["c_n_surf"](time[ind], x_n))
-    comsol_c_n_surf_plot.set_ydata(comsol_variables["c_n_surf"][:, ind])
-    # update electrolyte concentration
-    c_e_plot.set_ydata(processed_variables["c_e"](time[ind], x))
-    comsol_c_e_plot.set_ydata(comsol_variables["c_e"][:, ind])
-    # update negative particle surface concentration
-    c_p_surf_plot.set_ydata(processed_variables["c_p_surf"](time[ind], x_p))
-    comsol_c_p_surf_plot.set_ydata(comsol_variables["c_p_surf"][:, ind])
-    # update negative electrode potential
-    phi_n_plot.set_ydata(processed_variables["phi_n"](time[ind], x_n))
-    comsol_phi_n_plot.set_ydata(comsol_variables["phi_n"][:, ind])
-    # update electrolyte potential
-    phi_e_plot.set_ydata(processed_variables["phi_e"](time[ind], x))
-    comsol_phi_e_plot.set_ydata(comsol_variables["phi_e"][:, ind])
-    # update positive electrode potential
-    phi_p_plot.set_ydata(processed_variables["phi_p"](time[ind], x_p))
-    comsol_phi_p_plot.set_ydata(comsol_variables["phi_p"][:, ind])
-    fig.canvas.draw_idle()
-
-
-sfreq.on_changed(update_plot)
-plt.subplots_adjust(top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5)
-plt.show()
+plot.dynamic_plot()
diff --git a/examples/scripts/compare_lead_acid.py b/examples/scripts/compare_lead_acid.py
index cdf22b3db..4ecccb998 100644
--- a/examples/scripts/compare_lead_acid.py
+++ b/examples/scripts/compare_lead_acid.py
@@ -47,7 +47,10 @@ for i, model in enumerate(models):
 
 # plot
 output_variables = [
-    "Interfacial current density [A.m-2]",
+    [
+        "Average negative electrode interfacial current density [A.m-2]",
+        "Average positive electrode interfacial current density [A.m-2]",
+    ],
     "Electrolyte concentration [mol.m-3]",
     "Porosity",
     "Electrolyte current density [A.m-2]",
diff --git a/pandas/_libs/parsers.pyx b/pandas/_libs/parsers.pyx
deleted file mode 100644
index e69de29bb..000000000
diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index 8701785c6..efc634fa6 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -25,7 +25,10 @@ class Function(pybamm.Symbol):
 
     def __init__(self, function, *children):
 
-        name = "function ({})".format(function.__name__)
+        try:
+            name = "function ({})".format(function.__name__)
+        except AttributeError:
+            name = "function ({})".format(function.__class__)
         children_list = list(children)
         domain = self.get_children_domains(children_list)
 
diff --git a/pybamm/quick_plot.py b/pybamm/quick_plot.py
index 79a23a55d..5465d6be7 100644
--- a/pybamm/quick_plot.py
+++ b/pybamm/quick_plot.py
@@ -3,6 +3,7 @@
 #
 import numpy as np
 import pybamm
+from collections import defaultdict
 
 
 def ax_min(data):
@@ -23,6 +24,18 @@ def ax_max(data):
         return 1.1 * data_max
 
 
+def split_long_string(title, max_words=4):
+    "Get title in a nice format"
+    words = title.split()
+    # Don't split if fits on one line, don't split just for units
+    if len(words) <= max_words or words[max_words].startswith("["):
+        return title
+    else:
+        first_line = (" ").join(words[:max_words])
+        second_line = (" ").join(words[max_words:])
+        return first_line + "\n" + second_line
+
+
 class QuickPlot(object):
     """
     Generates a quick plot of a subset of key outputs of the model so that the model
@@ -39,11 +52,14 @@ class QuickPlot(object):
     solutions: (iter of) :class:`pybamm.Solver`
         The numerical solution(s) for the model(s) which contained the solution to the
         model(s).
-    output_variables : list of str
+    output_variables : list of str, optional
         List of variables to plot
+    labels : list of str, optional
+        Labels for the different models. Defaults to model names
     """
 
-    def __init__(self, models, mesh, solutions, output_variables=None):
+    def __init__(self, models, mesh, solutions, output_variables=None, labels=None):
+        # Pre-process models and solutions
         if isinstance(models, pybamm.BaseModel):
             models = [models]
         elif not isinstance(models, list):
@@ -57,14 +73,17 @@ class QuickPlot(object):
         else:
             raise ValueError("must provide the same number of models and solutions")
 
+        # Set labels
+        self.labels = labels or [model.name for model in models]
+
         # Scales (default to 1 if information not in model)
-        vars = models[0].variables
+        variables = models[0].variables
         self.x_scale = 1
         self.time_scale = 1
-        if "x [m]" and "x" in vars:
-            self.x_scale = (vars["x [m]"] / vars["x"]).evaluate()[-1]
-        if "Time [m]" and "Time" in vars:
-            self.time_scale = (vars["Time [h]"] / vars["Time"]).evaluate(t=1)
+        if "x [m]" and "x" in variables:
+            self.x_scale = (variables["x [m]"] / variables["x"]).evaluate()[-1]
+        if "Time [m]" and "Time" in variables:
+            self.time_scale = (variables["Time [h]"] / variables["Time"]).evaluate(t=1)
 
         # Time parameters
         self.ts = [solution.t for solution in solutions]
@@ -107,55 +126,83 @@ class QuickPlot(object):
 
         # Calculate subplot positions based on number of variables supplied
         self.subplot_positions = {}
-        n = int(len(output_variables) // np.sqrt(len(output_variables)))
-        m = np.ceil(len(output_variables) / n)
+        self.n_rows = int(len(output_variables) // np.sqrt(len(output_variables)))
+        self.n_cols = int(np.ceil(len(output_variables) / self.n_rows))
 
         # Process output variables into a form that can be plotted
-        for k, var in enumerate(output_variables):
-            self.variables[var] = [
-                pybamm.ProcessedVariable(
-                    models[i].variables[var], solutions[i].t, solutions[i].y, mesh
-                )
-                for i in range(len(models))
-            ]
-            if self.variables[var][0].dimensions == 2:
-                domain = models[0].variables[var].domain
-                self.x_values[var] = mesh.combine_submeshes(*domain)[0].edges
-            self.subplot_positions[var] = (n, m, k + 1)
+        for k, variable_list in enumerate(output_variables):
+            # Make sure we always have a list of lists of variables
+            if isinstance(variable_list, str):
+                variable_list = [variable_list]
 
-        # Set labels
-        self.labels = [model.name for model in models]
+            # Prepare list of variables
+            key = tuple(variable_list)
+            self.variables[key] = [None] * len(models)
+
+            # process each variable in variable_list for each model
+            for i, model in enumerate(models):
+                # self.variables is a dictionary of lists of lists
+                self.variables[key][i] = [
+                    pybamm.ProcessedVariable(
+                        model.variables[var], solutions[i].t, solutions[i].y, mesh
+                    )
+                    for var in variable_list
+                ]
+
+            # Make sure variables have the same dimensions and domain
+            domain = self.variables[key][0][0].domain
+            for variable in self.variables[key][0]:
+                if variable.domain != domain:
+                    raise ValueError("mismatching variable domains")
+
+            # Set the x variable for any two-dimensional variables
+            if self.variables[key][0][0].dimensions == 2:
+                self.x_values[key] = mesh.combine_submeshes(*domain)[0].edges
+
+            # Don't allow 3D variables
+            elif any(var.dimensions == 3 for var in self.variables[key][0]):
+                raise NotImplementedError("cannot plot 3D variables")
 
-        # Don't allow 3D variables
-        if self.variables[var][0].dimensions == 3:
-            raise NotImplementedError("cannot plot 3D variables")
+            # Define subplot position
+            self.subplot_positions[key] = (self.n_rows, self.n_cols, k + 1)
 
     def reset_axis(self):
         """
         Reset the axis limits to the default values.
+        These are calculated to fit around the minimum and maximum values of all the
+        variables in each subplot
         """
         self.axis = {}
-        for name, variable in self.variables.items():
-            if variable[0].dimensions == 1:
-                y_min = np.min([ax_min(v(self.ts[i])) for i, v in enumerate(variable)])
-                y_max = np.max([ax_max(v(self.ts[i])) for i, v in enumerate(variable)])
-                if y_min == y_max:
-                    y_min -= 1
-                    y_max += 1
-                self.axis[name] = [self.min_t, self.max_t, y_min, y_max]
-            elif variable[0].dimensions == 2:
-                x = self.x_values[name]
+        for key, variable_lists in self.variables.items():
+            if variable_lists[0][0].dimensions == 1:
+                x = None
+                x_min = self.min_t
+                x_max = self.max_t
+            elif variable_lists[0][0].dimensions == 2:
+                x = self.x_values[key]
                 x_scaled = x * self.x_scale
-                y_min = np.min(
-                    [ax_min(v(self.ts[i], x)) for i, v in enumerate(variable)]
-                )
-                y_max = np.max(
-                    [ax_max(v(self.ts[i], x)) for i, v in enumerate(variable)]
-                )
-                if y_min == y_max:
-                    y_min -= 1
-                    y_max += 1
-                self.axis[name] = [x_scaled[0], x_scaled[-1], y_min, y_max]
+                x_min = x_scaled[0]
+                x_max = x_scaled[-1]
+
+            # Get min and max y values
+            y_min = np.min(
+                [
+                    ax_min(var(self.ts[i], x))
+                    for i, variable_list in enumerate(variable_lists)
+                    for var in variable_list
+                ]
+            )
+            y_max = np.max(
+                [
+                    ax_max(var(self.ts[i], x))
+                    for i, variable_list in enumerate(variable_lists)
+                    for var in variable_list
+                ]
+            )
+            if y_min == y_max:
+                y_min -= 1
+                y_max += 1
+            self.axis[key] = [x_min, x_max, y_min, y_max]
 
     def plot(self, t):
         """Produces a quick plot with the internal states at time t.
@@ -169,51 +216,81 @@ class QuickPlot(object):
         import matplotlib.pyplot as plt
 
         t /= self.time_scale
-        self.fig, self.ax = plt.subplots(figsize=(15, 8))
+        self.fig, self.ax = plt.subplots(self.n_rows, self.n_cols, figsize=(15, 8))
         plt.tight_layout()
         plt.subplots_adjust(left=-0.1)
         self.plots = {}
         self.time_lines = {}
 
-        for k, name in enumerate(self.variables.keys()):
-            variable = self.variables[name]
-            plt.subplot(*self.subplot_positions[name])
-            plt.ylabel(name, fontsize=14)
-            plt.axis(self.axis[name])
-            self.plots[name] = [None] * self.num_models
-            # Set labels
+        colors = ["r", "b", "k", "g"]
+        linestyles = ["-", ":", "--", "-."]
+        fontsize = 42 // self.n_cols
+
+        for k, (key, variable_lists) in enumerate(self.variables.items()):
+            if len(self.variables) == 1:
+                ax = self.ax
+            else:
+                ax = self.ax.flat[k]
+            ax.set_xlim(self.axis[key][:2])
+            ax.set_ylim(self.axis[key][2:])
+            ax.xaxis.set_major_locator(plt.MaxNLocator(3))
+            self.plots[key] = defaultdict(dict)
+            # Set labels for the first subplot only (avoid repetition)
             if k == 0:
                 labels = self.labels
             else:
                 labels = [None] * len(self.labels)
-            if variable[0].dimensions == 2:
+            if variable_lists[0][0].dimensions == 2:
                 # 2D plot: plot as a function of x at time t
-                plt.xlabel("Position [m]", fontsize=14)
-                x_value = self.x_values[name]
-                for i in range(self.num_models):
-                    self.plots[name][i], = plt.plot(
-                        x_value * self.x_scale,
-                        variable[i](t, x_value),
-                        lw=2,
-                        label=labels[i],
-                    )
+                ax.set_xlabel("Position [m]", fontsize=fontsize)
+                x_value = self.x_values[key]
+                for i, variable_list in enumerate(variable_lists):
+                    for j, variable in enumerate(variable_list):
+                        if j == 0:
+                            label = labels[i]
+                        else:
+                            label = None
+                        self.plots[key][i][j], = ax.plot(
+                            x_value * self.x_scale,
+                            variable(t, x_value),
+                            lw=2,
+                            color=colors[i],
+                            linestyle=linestyles[j],
+                            label=label,
+                        )
             else:
                 # 1D plot: plot as a function of time, indicating time t with a line
-                plt.xlabel("Time [h]", fontsize=14)
-                for i in range(self.num_models):
-                    full_t = self.ts[i]
-                    self.plots[name][i], = plt.plot(
-                        full_t * self.time_scale,
-                        variable[i](full_t),
-                        lw=2,
-                        label=labels[i],
-                    )
-                    y_min, y_max = self.axis[name][2:]
-                    self.time_lines[name], = plt.plot(
-                        [t * self.time_scale, t * self.time_scale],
-                        [y_min, y_max],
-                        "k--",
-                    )
+                ax.set_xlabel("Time [h]", fontsize=fontsize)
+                for i, variable_list in enumerate(variable_lists):
+                    for j, variable in enumerate(variable_list):
+                        full_t = self.ts[i]
+                        if j == 0:
+                            label = labels[i]
+                        else:
+                            label = None
+                        self.plots[key][i][j], = ax.plot(
+                            full_t * self.time_scale,
+                            variable(full_t),
+                            lw=2,
+                            color=colors[i],
+                            linestyle=linestyles[j],
+                            label=label,
+                        )
+                y_min, y_max = self.axis[key][2:]
+                self.time_lines[key], = ax.plot(
+                    [t * self.time_scale, t * self.time_scale], [y_min, y_max], "k--"
+                )
+            # Set either y label or legend entries
+            if len(key) == 1:
+                title = split_long_string(key[0])
+                ax.set_title(title, fontsize=fontsize)
+            else:
+                ax.legend(
+                    [split_long_string(s, 6) for s in key],
+                    bbox_to_anchor=(0.5, 1.2),
+                    fontsize=8,
+                    loc="upper center",
+                )
         self.fig.legend(loc="lower right")
 
     def dynamic_plot(self, testing=False):
@@ -246,12 +323,13 @@ class QuickPlot(object):
         """
         t = self.sfreq.val
         t_dimensionless = t / self.time_scale
-        for var, plot in self.plots.items():
-            if self.variables[var][0].dimensions == 2:
-                x = self.x_values[var]
-                for i in range(self.num_models):
-                    plot[i].set_ydata(self.variables[var][i](t_dimensionless, x))
+        for key, plot in self.plots.items():
+            if self.variables[key][0][0].dimensions == 2:
+                x = self.x_values[key]
+                for i, variable_lists in enumerate(self.variables[key]):
+                    for j, variable in enumerate(variable_lists):
+                        plot[i][j].set_ydata(variable(t_dimensionless, x))
             else:
-                self.time_lines[var].set_xdata([t])
+                self.time_lines[key].set_xdata([t])
 
         self.fig.canvas.draw_idle()


Test patch:
diff --git a/tests/integration/test_quick_plot.py b/tests/integration/test_quick_plot.py
index 05b49d33e..d573bf674 100644
--- a/tests/integration/test_quick_plot.py
+++ b/tests/integration/test_quick_plot.py
@@ -30,12 +30,12 @@ class TestQuickPlot(unittest.TestCase):
 
         # update the axis
         new_axis = [0, 0.5, 0, 1]
-        quick_plot.axis.update({"Electrolyte concentration": new_axis})
-        self.assertEqual(quick_plot.axis["Electrolyte concentration"], new_axis)
+        quick_plot.axis.update({("Electrolyte concentration",): new_axis})
+        self.assertEqual(quick_plot.axis[("Electrolyte concentration",)], new_axis)
 
         # and now reset them
         quick_plot.reset_axis()
-        self.assertNotEqual(quick_plot.axis["Electrolyte concentration"], new_axis)
+        self.assertNotEqual(quick_plot.axis[("Electrolyte concentration",)], new_axis)
 
         # check dynamic plot loads
         quick_plot.dynamic_plot(testing=True)
@@ -61,12 +61,12 @@ class TestQuickPlot(unittest.TestCase):
 
         # update the axis
         new_axis = [0, 0.5, 0, 1]
-        quick_plot.axis.update({"Electrolyte concentration": new_axis})
-        self.assertEqual(quick_plot.axis["Electrolyte concentration"], new_axis)
+        quick_plot.axis.update({("Electrolyte concentration",): new_axis})
+        self.assertEqual(quick_plot.axis[("Electrolyte concentration",)], new_axis)
 
         # and now reset them
         quick_plot.reset_axis()
-        self.assertNotEqual(quick_plot.axis["Electrolyte concentration"], new_axis)
+        self.assertNotEqual(quick_plot.axis[("Electrolyte concentration",)], new_axis)
 
         # check dynamic plot loads
         quick_plot.dynamic_plot(testing=True)
diff --git a/tests/unit/test_expression_tree/test_functions.py b/tests/unit/test_expression_tree/test_functions.py
index 70b1e8a04..0ac27a800 100644
--- a/tests/unit/test_expression_tree/test_functions.py
+++ b/tests/unit/test_expression_tree/test_functions.py
@@ -6,6 +6,7 @@ import pybamm
 import unittest
 import numpy as np
 import autograd.numpy as auto_np
+from scipy.interpolate import interp1d
 
 
 def test_function(arg):
@@ -93,6 +94,15 @@ class TestFunction(unittest.TestCase):
         with self.assertRaises(pybamm.DomainError):
             pybamm.Function(test_multi_var_function, a, b)
 
+    def test_function_unnamed(self):
+        t = np.linspace(0, 1)
+        entries = 2 * t
+        interpfun = interp1d(t, entries)
+        fun = pybamm.Function(interpfun, pybamm.t)
+        self.assertEqual(
+            fun.name, "function (<class 'scipy.interpolate.interpolate.interp1d'>)"
+        )
+
 
 class TestSpecificFunctions(unittest.TestCase):
     def test_cos(self):
diff --git a/tests/unit/test_quick_plot.py b/tests/unit/test_quick_plot.py
index 526a04af4..3913c50f2 100644
--- a/tests/unit/test_quick_plot.py
+++ b/tests/unit/test_quick_plot.py
@@ -21,12 +21,12 @@ class TestQuickPlot(unittest.TestCase):
 
         # update the axis
         new_axis = [0, 0.5, 0, 1]
-        quick_plot.axis.update({"a": new_axis})
-        self.assertEqual(quick_plot.axis["a"], new_axis)
+        quick_plot.axis.update({("a",): new_axis})
+        self.assertEqual(quick_plot.axis[("a",)], new_axis)
 
         # and now reset them
         quick_plot.reset_axis()
-        self.assertNotEqual(quick_plot.axis["a"], new_axis)
+        self.assertNotEqual(quick_plot.axis[("a",)], new_axis)
 
         # check dynamic plot loads
         quick_plot.dynamic_plot(testing=True)
@@ -34,27 +34,66 @@ class TestQuickPlot(unittest.TestCase):
         quick_plot.update(0.01)
 
         # Test with different output variables
+        quick_plot = pybamm.QuickPlot(model, mesh, solution, ["b broadcasted"])
+        self.assertEqual(len(quick_plot.axis), 1)
+        quick_plot.plot(0)
+
         quick_plot = pybamm.QuickPlot(
-            model, mesh, solution, ["b broadcasted", "c broadcasted"]
+            model,
+            mesh,
+            solution,
+            [["a", "a"], ["b broadcasted", "b broadcasted"], "c broadcasted"],
         )
-        self.assertEqual(len(quick_plot.axis), 2)
+        self.assertEqual(len(quick_plot.axis), 3)
         quick_plot.plot(0)
 
         # update the axis
         new_axis = [0, 0.5, 0, 1]
-        var = "c broadcasted"
-        quick_plot.axis.update({var: new_axis})
-        self.assertEqual(quick_plot.axis[var], new_axis)
+        var_key = ("c broadcasted",)
+        quick_plot.axis.update({var_key: new_axis})
+        self.assertEqual(quick_plot.axis[var_key], new_axis)
 
         # and now reset them
         quick_plot.reset_axis()
-        self.assertNotEqual(quick_plot.axis[var], new_axis)
+        self.assertNotEqual(quick_plot.axis[var_key], new_axis)
 
         # check dynamic plot loads
         quick_plot.dynamic_plot(testing=True)
 
         quick_plot.update(0.01)
 
+        # Test longer name
+        model.variables["Variable with a very long name"] = model.variables["a"]
+        quick_plot = pybamm.QuickPlot(model, mesh, solution)
+        quick_plot.plot(0)
+
+        # Test errors
+        with self.assertRaisesRegex(ValueError, "mismatching variable domains"):
+            pybamm.QuickPlot(model, mesh, solution, [["a", "b broadcasted"]])
+        model.variables["3D variable"] = disc.process_symbol(
+            pybamm.Broadcast(1, ["negative particle"])
+        )
+        with self.assertRaisesRegex(NotImplementedError, "cannot plot 3D variables"):
+            pybamm.QuickPlot(model, mesh, solution, ["3D variable"])
+
+    def test_loqs_spm_base(self):
+        t_eval = np.linspace(0, 0.01, 2)
+
+        # SPM
+        for model in [pybamm.lithium_ion.SPM(), pybamm.lead_acid.LOQS()]:
+            geometry = model.default_geometry
+            param = model.default_parameter_values
+            param.process_model(model)
+            param.process_geometry(geometry)
+            mesh = pybamm.Mesh(
+                geometry, model.default_submesh_types, model.default_var_pts
+            )
+            disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+            disc.process_model(model)
+            solver = model.default_solver
+            solution = solver.solve(model, t_eval)
+            pybamm.QuickPlot(model, mesh, solution)
+
     def test_failure(self):
         with self.assertRaisesRegex(TypeError, "'models' must be"):
             pybamm.QuickPlot(1, None, None)