Repository: pybamm-team/PyBaMM

Issue:
Evaluate all variables using known_evals
Quick test suggests this speeds up the evaluation a fair bit

Discussion:

Patch:
diff --git a/pybamm/processed_variable.py b/pybamm/processed_variable.py
index a2787d336..f00e0fb3d 100644
--- a/pybamm/processed_variable.py
+++ b/pybamm/processed_variable.py
@@ -3,6 +3,7 @@
 #
 import numbers
 import numpy as np
+import pybamm
 import scipy.interpolate as interp
 
 
@@ -31,10 +32,14 @@ def post_process_variables(variables, t_sol, y_sol, mesh=None, interp_kind="line
         Dictionary of processed variables
     """
     processed_variables = {}
+    known_evals = {t: {} for t in t_sol}
     for var, eqn in variables.items():
+        pybamm.logger.debug("Post-processing {}".format(var))
         processed_variables[var] = ProcessedVariable(
-            eqn, t_sol, y_sol, mesh, interp_kind
+            eqn, t_sol, y_sol, mesh, interp_kind, known_evals
         )
+        for t in known_evals:
+            known_evals[t].update(processed_variables[var].known_evals[t])
     return processed_variables
 
 
@@ -62,15 +67,29 @@ class ProcessedVariable(object):
         The method to use for interpolation
     """
 
-    def __init__(self, base_variable, t_sol, y_sol, mesh=None, interp_kind="linear"):
+    def __init__(
+        self,
+        base_variable,
+        t_sol,
+        y_sol,
+        mesh=None,
+        interp_kind="linear",
+        known_evals=None,
+    ):
         self.base_variable = base_variable
         self.t_sol = t_sol
         self.y_sol = y_sol
         self.mesh = mesh
         self.interp_kind = interp_kind
         self.domain = base_variable.domain
+        self.known_evals = known_evals
 
-        self.base_eval = base_variable.evaluate(t_sol[0], y_sol[:, 0])
+        if self.known_evals:
+            self.base_eval, self.known_evals[t_sol[0]] = base_variable.evaluate(
+                t_sol[0], y_sol[:, 0], self.known_evals[t_sol[0]]
+            )
+        else:
+            self.base_eval = base_variable.evaluate(t_sol[0], y_sol[:, 0])
 
         if (
             isinstance(self.base_eval, numbers.Number)
@@ -79,7 +98,8 @@ class ProcessedVariable(object):
         ):
             self.initialise_1D()
         else:
-            if len(self.mesh.combine_submeshes(*self.domain)) == 1:
+            n = self.mesh.combine_submeshes(*self.domain)[0].npts
+            if self.base_eval.shape[0] in [n, n + 1]:
                 self.initialise_2D()
             else:
                 self.initialise_3D()
@@ -92,9 +112,13 @@ class ProcessedVariable(object):
         entries = np.empty(len(self.t_sol))
         # Evaluate the base_variable index-by-index
         for idx in range(len(self.t_sol)):
-            entries[idx] = self.base_variable.evaluate(
-                self.t_sol[idx], self.y_sol[:, idx]
-            )
+            t = self.t_sol[idx]
+            if self.known_evals:
+                entries[idx], self.known_evals[t] = self.base_variable.evaluate(
+                    t, self.y_sol[:, idx], self.known_evals[t]
+                )
+            else:
+                entries[idx] = self.base_variable.evaluate(t, self.y_sol[:, idx])
 
         # No discretisation provided, or variable has no domain (function of t only)
         self._interpolation_function = interp.interp1d(
@@ -115,9 +139,16 @@ class ProcessedVariable(object):
 
         # Evaluate the base_variable index-by-index
         for idx in range(len(self.t_sol)):
-            entries[:, idx] = self.base_variable.evaluate(
-                self.t_sol[idx], self.y_sol[:, idx]
-            )[:, 0]
+            t = self.t_sol[idx]
+            y = self.y_sol[:, idx]
+            if self.known_evals:
+                eval_and_known_evals = self.base_variable.evaluate(
+                    t, y, self.known_evals[t]
+                )
+                entries[:, idx] = eval_and_known_evals[0][:, 0]
+                self.known_evals[t] = eval_and_known_evals[1]
+            else:
+                entries[:, idx] = self.base_variable.evaluate(t, y)[:, 0]
 
         # Process the discretisation to get x values
         nodes = self.mesh.combine_submeshes(*self.domain)[0].nodes
@@ -126,8 +157,6 @@ class ProcessedVariable(object):
             space = nodes
         elif entries.shape[0] == len(edges):
             space = edges
-        else:
-            raise ValueError("variable shape does not match domain shape")
 
         # add points outside domain for extrapolation to boundaries
         extrap_space_left = np.array([2 * space[0] - space[1]])
@@ -163,10 +192,18 @@ class ProcessedVariable(object):
 
         # Evaluate the base_variable index-by-index
         for idx in range(len(self.t_sol)):
-            entries[:, :, idx] = np.reshape(
-                self.base_variable.evaluate(self.t_sol[idx], self.y_sol[:, idx]),
-                [len_x, len_r],
-            )
+            t = self.t_sol[idx]
+            y = self.y_sol[:, idx]
+            if self.known_evals:
+                eval_and_known_evals = self.base_variable.evaluate(
+                    t, y, self.known_evals[t]
+                )
+                entries[:, :, idx] = np.reshape(eval_and_known_evals[0], [len_x, len_r])
+                self.known_evals[t] = eval_and_known_evals[1]
+            else:
+                entries[:, :, idx] = np.reshape(
+                    self.base_variable.evaluate(t, y), [len_x, len_r]
+                )
         # Process the discretisation to get x values
         nodes = self.mesh.combine_submeshes(*self.domain)[0].nodes
         edges = self.mesh.combine_submeshes(*self.domain)[0].edges
@@ -175,7 +212,7 @@ class ProcessedVariable(object):
         elif entries.shape[1] == len(edges):
             r_sol = edges
         else:
-            raise ValueError("variable shape does not match domain shape")
+            raise ValueError("3D variable shape does not match domain shape")
 
         # Get x values
         if self.domain == ["negative particle"]:


Test patch:
diff --git a/tests/integration/test_models/standard_output_comparison.py b/tests/integration/test_models/standard_output_comparison.py
index 6fa9bb604..fd237d295 100644
--- a/tests/integration/test_models/standard_output_comparison.py
+++ b/tests/integration/test_models/standard_output_comparison.py
@@ -9,6 +9,14 @@ class StandardOutputComparison(object):
     "Calls all the tests comparing standard output variables."
 
     def __init__(self, models, discs, solvers):
+        # Process variables
+        for model in models:
+            disc = discs[model]
+            solver = solvers[model]
+            model.variables = pybamm.post_process_variables(
+                model.variables, solver.t, solver.y, disc.mesh
+            )
+
         self.models = models
         self.discs = discs
         self.solvers = solvers
@@ -79,22 +87,10 @@ class BaseOutputComparison(object):
         self.mesh = mesh
         self.solvers = solvers
 
-    def get_vars(self, var):
-        "Helper function to reduce repeated code."
-        return [
-            pybamm.ProcessedVariable(
-                model.variables[var],
-                self.solvers[model].t,
-                self.solvers[model].y,
-                mesh=self.mesh,
-            )
-            for model in self.models
-        ]
-
     def compare(self, var, tol=1e-2):
         "Compare variables from different models"
         # Get variable for each model
-        model_variables = self.get_vars(var)
+        model_variables = [model.variables[var] for model in self.models]
         var0 = model_variables[0]
 
         if var0.domain == []:
diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py
index 330b70efa..00228a18a 100644
--- a/tests/integration/test_models/standard_output_tests.py
+++ b/tests/integration/test_models/standard_output_tests.py
@@ -9,6 +9,12 @@ class StandardOutputTests(object):
     "Calls all the tests on the standard output variables."
 
     def __init__(self, model, parameter_values, disc, solver):
+        # Process variables
+        model.variables = pybamm.post_process_variables(
+            model.variables, solver.t, solver.y, disc.mesh
+        )
+
+        # Assign attributes
         self.model = model
         self.parameter_values = parameter_values
         self.disc = disc
@@ -27,6 +33,9 @@ class StandardOutputTests(object):
         else:
             self.operating_condition = "off"
 
+    def process_variables(self):
+        return
+
     def run_test_class(self, ClassName):
         "Run all tests from a class 'ClassName'"
         tests = ClassName(
@@ -83,35 +92,29 @@ class BaseOutputTest(object):
             pybamm.electrical_parameters.current_with_time
         ).evaluate(self.t)
 
-    def get_var(self, var):
-        "Helper function to reduce repeated code."
-        pybamm.logger.debug("Processing {} for {}".format(var, self.model.name))
-        return pybamm.ProcessedVariable(
-            self.model.variables[var], self.solver.t, self.solver.y, mesh=self.disc.mesh
-        )
-
 
 class VoltageTests(BaseOutputTest):
     def __init__(self, model, param, disc, solver, operating_condition):
         super().__init__(model, param, disc, solver, operating_condition)
+        variables = self.model.variables
 
-        self.eta_n = self.get_var("Negative reaction overpotential [V]")
-        self.eta_p = self.get_var("Positive reaction overpotential [V]")
-        self.eta_r_n_av = self.get_var("Average negative reaction overpotential [V]")
-        self.eta_r_p_av = self.get_var("Average positive reaction overpotential [V]")
-        self.eta_r_av = self.get_var("Average reaction overpotential [V]")
+        self.eta_n = variables["Negative reaction overpotential [V]"]
+        self.eta_p = variables["Positive reaction overpotential [V]"]
+        self.eta_r_n_av = variables["Average negative reaction overpotential [V]"]
+        self.eta_r_p_av = variables["Average positive reaction overpotential [V]"]
+        self.eta_r_av = variables["Average reaction overpotential [V]"]
 
-        self.eta_e_av = self.get_var("Average electrolyte overpotential [V]")
-        self.delta_phi_s_av = self.get_var("Average solid phase ohmic losses [V]")
+        self.eta_e_av = variables["Average electrolyte overpotential [V]"]
+        self.delta_phi_s_av = variables["Average solid phase ohmic losses [V]"]
 
-        self.ocp_n_av = self.get_var(
+        self.ocp_n_av = variables[
             "Average negative electrode open circuit potential [V]"
-        )
-        self.ocp_p_av = self.get_var(
+        ]
+        self.ocp_p_av = variables[
             "Average positive electrode open circuit potential [V]"
-        )
-        self.ocv_av = self.get_var("Average open circuit voltage [V]")
-        self.voltage = self.get_var("Terminal voltage [V]")
+        ]
+        self.ocv_av = variables["Average open circuit voltage [V]"]
+        self.voltage = variables["Terminal voltage [V]"]
 
     def test_each_reaction_overpotential(self):
         """Testing that:
@@ -233,15 +236,16 @@ class VoltageTests(BaseOutputTest):
 class ParticleConcentrationTests(BaseOutputTest):
     def __init__(self, model, param, disc, solver, operating_condition):
         super().__init__(model, param, disc, solver, operating_condition)
+        variables = self.model.variables
 
-        self.c_s_n = self.get_var("Negative particle concentration")
-        self.c_s_p = self.get_var("Positive particle concentration")
+        self.c_s_n = variables["Negative particle concentration"]
+        self.c_s_p = variables["Positive particle concentration"]
 
-        self.c_s_n_surf = self.get_var("Negative particle surface concentration")
-        self.c_s_p_surf = self.get_var("Positive particle surface concentration")
+        self.c_s_n_surf = variables["Negative particle surface concentration"]
+        self.c_s_p_surf = variables["Positive particle surface concentration"]
 
-        self.N_s_n = self.get_var("Negative particle flux")
-        self.N_s_p = self.get_var("Positive particle flux")
+        self.N_s_n = variables["Negative particle flux"]
+        self.N_s_p = variables["Positive particle flux"]
 
     def test_concentration_increase_decrease(self):
         """Test all concentrations in negative particles decrease and all
@@ -314,21 +318,22 @@ class ParticleConcentrationTests(BaseOutputTest):
 class ElectrolyteConcentrationTests(BaseOutputTest):
     def __init__(self, model, param, disc, solver, operating_condition):
         super().__init__(model, param, disc, solver, operating_condition)
+        variables = self.model.variables
 
-        self.c_e = self.get_var("Electrolyte concentration")
+        self.c_e = variables["Electrolyte concentration"]
 
-        self.c_e_n = self.get_var("Negative electrolyte concentration")
-        self.c_e_s = self.get_var("Separator electrolyte concentration")
-        self.c_e_p = self.get_var("Positive electrolyte concentration")
+        self.c_e_n = variables["Negative electrolyte concentration"]
+        self.c_e_s = variables["Separator electrolyte concentration"]
+        self.c_e_p = variables["Positive electrolyte concentration"]
 
         # TODO: output average electrolyte concentration
-        # self.c_e_av = self.get_var("Average electrolyte concentration")
-        # self.c_e_n_av = self.get_var("Average negative electrolyte concentration")
-        # self.c_e_s_av = self.get_var("Average separator electrolyte concentration")
-        # self.c_e_p_av = self.get_var("Average positive electrolyte concentration")
+        # self.c_e_av = variables["Average electrolyte concentration"]
+        # self.c_e_n_av = variables["Average negative electrolyte concentration"]
+        # self.c_e_s_av = variables["Average separator electrolyte concentration"]
+        # self.c_e_p_av = variables["Average positive electrolyte concentration"]
 
-        # self.N_e = self.get_var("Electrolyte flux")
-        self.N_e_hat = self.get_var("Reduced cation flux")
+        # self.N_e = variables["Electrolyte flux"]
+        self.N_e_hat = variables["Reduced cation flux"]
 
     def test_concentration_limit(self):
         "Test that the electrolyte concentration is always greater than zero."
@@ -399,15 +404,16 @@ class ElectrolyteConcentrationTests(BaseOutputTest):
 class PotentialTests(BaseOutputTest):
     def __init__(self, model, param, disc, solver, operating_condition):
         super().__init__(model, param, disc, solver, operating_condition)
+        variables = self.model.variables
 
-        self.phi_s_n = self.get_var("Negative electrode potential [V]")
-        self.phi_s_p = self.get_var("Positive electrode potential [V]")
+        self.phi_s_n = variables["Negative electrode potential [V]"]
+        self.phi_s_p = variables["Positive electrode potential [V]"]
 
-        self.phi_e = self.get_var("Electrolyte potential [V]")
+        self.phi_e = variables["Electrolyte potential [V]"]
 
-        self.phi_e_n = self.get_var("Negative electrolyte potential [V]")
-        self.phi_e_s = self.get_var("Separator electrolyte potential [V]")
-        self.phi_e_p = self.get_var("Positive electrolyte potential [V]")
+        self.phi_e_n = variables["Negative electrolyte potential [V]"]
+        self.phi_e_s = variables["Separator electrolyte potential [V]"]
+        self.phi_e_p = variables["Positive electrolyte potential [V]"]
 
     def test_negative_electrode_potential_profile(self):
         """Test that negative electrode potential is zero on left boundary. Test
@@ -440,26 +446,27 @@ class PotentialTests(BaseOutputTest):
 class CurrentTests(BaseOutputTest):
     def __init__(self, model, param, disc, solver, operating_condition):
         super().__init__(model, param, disc, solver, operating_condition)
+        variables = self.model.variables
 
-        self.j = self.get_var("Interfacial current density")
-        self.j0 = self.get_var("Exchange-current density")
+        self.j = variables["Interfacial current density"]
+        self.j0 = variables["Exchange-current density"]
 
-        self.j_n = self.get_var("Negative electrode interfacial current density")
-        self.j_p = self.get_var("Positive electrode interfacial current density")
-        self.j_n_av = self.get_var(
+        self.j_n = variables["Negative electrode interfacial current density"]
+        self.j_p = variables["Positive electrode interfacial current density"]
+        self.j_n_av = variables[
             "Average negative electrode interfacial current density"
-        )
-        self.j_p_av = self.get_var(
+        ]
+        self.j_p_av = variables[
             "Average positive electrode interfacial current density"
-        )
+        ]
 
-        self.j0_n = self.get_var("Negative electrode exchange-current density")
-        self.j0_p = self.get_var("Positive electrode exchange-current density")
+        self.j0_n = variables["Negative electrode exchange-current density"]
+        self.j0_p = variables["Positive electrode exchange-current density"]
 
-        self.i_s_n = self.get_var("Negative electrode current density")
-        self.i_s_p = self.get_var("Positive electrode current density")
-        self.i_s = self.get_var("Electrode current density")
-        self.i_e = self.get_var("Electrolyte current density")
+        self.i_s_n = variables["Negative electrode current density"]
+        self.i_s_p = variables["Positive electrode current density"]
+        self.i_s = variables["Electrode current density"]
+        self.i_e = variables["Electrolyte current density"]
 
     def test_interfacial_current_average(self):
         """Test that average of the interfacial current density is equal to the true
@@ -512,10 +519,11 @@ class CurrentTests(BaseOutputTest):
 class VelocityTests(BaseOutputTest):
     def __init__(self, model, param, disc, solver, operating_condition):
         super().__init__(model, param, disc, solver, operating_condition)
+        variables = self.model.variables
 
-        self.v_box = self.get_var("Volume-averaged velocity")
-        self.i_e = self.get_var("Electrolyte current density")
-        self.dVbox_dz = self.get_var("Vertical volume-averaged acceleration")
+        self.v_box = variables["Volume-averaged velocity"]
+        self.i_e = variables["Electrolyte current density"]
+        self.dVbox_dz = variables["Vertical volume-averaged acceleration"]
 
     def test_velocity_boundaries(self):
         """Test the boundary values of the current densities"""
diff --git a/tests/integration/test_models/test_lead_acid/test_composite.py b/tests/integration/test_models/test_lead_acid/test_composite.py
index 05c58022f..f48da169b 100644
--- a/tests/integration/test_models/test_lead_acid/test_composite.py
+++ b/tests/integration/test_models/test_lead_acid/test_composite.py
@@ -47,22 +47,8 @@ class TestLeadAcidCompositeCapacitance(unittest.TestCase):
         options = {"capacitance": "algebraic"}
         model = pybamm.lead_acid.Composite(options)
         modeltest = tests.StandardModelTest(model)
-
         modeltest.test_all()
 
-    def test_optimisations(self):
-        options = {"capacitance": "differential"}
-        model = pybamm.lead_acid.Composite(options)
-        optimtest = tests.OptimisationsTest(model)
-
-        original = optimtest.evaluate_model()
-        simplified = optimtest.evaluate_model(simplify=True)
-        using_known_evals = optimtest.evaluate_model(use_known_evals=True)
-        simp_and_known = optimtest.evaluate_model(simplify=True, use_known_evals=True)
-        np.testing.assert_array_almost_equal(original, simplified)
-        np.testing.assert_array_almost_equal(original, using_known_evals)
-        np.testing.assert_array_almost_equal(original, simp_and_known)
-
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/integration/test_models/test_lithium_ion/test_lithium_ion_dfn.py b/tests/integration/test_models/test_lithium_ion/test_lithium_ion_dfn.py
index 9138ffb69..33dfd1978 100644
--- a/tests/integration/test_models/test_lithium_ion/test_lithium_ion_dfn.py
+++ b/tests/integration/test_models/test_lithium_ion/test_lithium_ion_dfn.py
@@ -14,7 +14,7 @@ class TestDFN(unittest.TestCase):
     def test_basic_processing(self):
         model = pybamm.lithium_ion.DFN()
         var = pybamm.standard_spatial_vars
-        var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.r_n: 10, var.r_p: 10}
+        var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.r_n: 5, var.r_p: 5}
         modeltest = tests.StandardModelTest(model, var_pts=var_pts)
         modeltest.test_all()
 
diff --git a/tests/unit/test_processed_variable.py b/tests/unit/test_processed_variable.py
index 013962580..26922ab60 100644
--- a/tests/unit/test_processed_variable.py
+++ b/tests/unit/test_processed_variable.py
@@ -25,6 +25,7 @@ class TestProcessedVariable(unittest.TestCase):
         x = pybamm.SpatialVariable("x", domain=["negative electrode", "separator"])
         eqn = t * var + x
 
+        # On nodes
         disc = tests.get_discretisation_for_testing()
         disc.set_variable_slices([var])
         x_sol = disc.process_symbol(x).entries[:, 0]
@@ -47,6 +48,13 @@ class TestProcessedVariable(unittest.TestCase):
             processed_var.entries[1], 2 * y_sol[-1] - y_sol[-2]
         )
 
+        # On edges
+        x_s_edge = pybamm.Matrix(disc.mesh["separator"][0].edges, domain="separator")
+        processed_x_s_edge = pybamm.ProcessedVariable(x_s_edge, t_sol, y_sol, disc.mesh)
+        np.testing.assert_array_equal(
+            x_s_edge.entries[:, 0], processed_x_s_edge.entries[1:-1, 0]
+        )
+
     def test_processed_variable_3D(self):
         var = pybamm.Variable("var", domain=["negative particle"])
         x = pybamm.SpatialVariable("x", domain=["negative electrode"])
@@ -125,6 +133,16 @@ class TestProcessedVariable(unittest.TestCase):
         # 2 scalars
         self.assertEqual(processed_eqn(0.5, x_sol[-1]).shape, (1,))
 
+        # On microscale
+        r_n = pybamm.Matrix(
+            disc.mesh["negative particle"][0].nodes, domain="negative particle"
+        )
+        processed_r_n = pybamm.ProcessedVariable(r_n, t_sol, y_sol, disc.mesh)
+        np.testing.assert_array_equal(r_n.entries[:, 0], processed_r_n.entries[1:-1, 0])
+        np.testing.assert_array_almost_equal(
+            processed_r_n(0, r=np.linspace(0, 1))[:, 0], np.linspace(0, 1)
+        )
+
     def test_processed_var_3D_interpolation(self):
         var = pybamm.Variable("var", domain=["negative particle"])
         x = pybamm.SpatialVariable("x", domain=["negative electrode"])
@@ -158,6 +176,24 @@ class TestProcessedVariable(unittest.TestCase):
         # 3 scalars
         np.testing.assert_array_equal(processed_var(0.2, 0.2, 0.2).shape, ())
 
+        # positive particle
+        var = pybamm.Variable("var", domain=["positive particle"])
+        x = pybamm.SpatialVariable("x", domain=["positive electrode"])
+        r = pybamm.SpatialVariable("r", domain=["positive particle"])
+
+        disc.set_variable_slices([var])
+        x_sol = disc.process_symbol(x).entries[:, 0]
+        r_sol = disc.process_symbol(r).entries[:, 0]
+        var_sol = disc.process_symbol(var)
+        t_sol = np.linspace(0, 1)
+        y_sol = np.ones(len(x_sol) * len(r_sol))[:, np.newaxis] * np.linspace(0, 5)
+
+        processed_var = pybamm.ProcessedVariable(var_sol, t_sol, y_sol, mesh=disc.mesh)
+        # 3 vectors
+        np.testing.assert_array_equal(
+            processed_var(t_sol, x_sol, r_sol).shape, (35, 10, 50)
+        )
+
     def test_processed_variable_ode_pde_solution(self):
         # without space
         model = pybamm.StandardBatteryBaseModel()
@@ -176,9 +212,19 @@ class TestProcessedVariable(unittest.TestCase):
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         model = pybamm.StandardBatteryBaseModel()
         c = pybamm.Variable("conc", domain=whole_cell)
-        model.rhs = {c: -c}
-        model.initial_conditions = {c: 1}
-        model.variables = {"c": c}
+        c_s = pybamm.Variable("particle conc", domain="negative particle")
+        model.rhs = {c: -c, c_s: 1 - c_s}
+        model.initial_conditions = {c: 1, c_s: 0.5}
+        model.boundary_conditions = {
+            c: {"left": (0, "Neumann"), "right": (0, "Neumann")},
+            c_s: {"left": (0, "Neumann"), "right": (0, "Neumann")},
+        }
+        model.variables = {
+            "c": c,
+            "N": pybamm.grad(c),
+            "c_s": c_s,
+            "N_s": pybamm.grad(c_s),
+        }
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
         # set up testing
@@ -197,19 +243,11 @@ class TestProcessedVariable(unittest.TestCase):
 
     def test_failure(self):
         t = np.ones(25)
-        y = np.ones((15, 25))
-        mat = pybamm.Vector(np.ones(15), domain=["negative electrode"])
-        disc = tests.get_discretisation_for_testing()
-        with self.assertRaisesRegex(
-            ValueError, "variable shape does not match domain shape"
-        ):
-            pybamm.ProcessedVariable(mat, t, y, disc.mesh)
-
         y = np.ones((120, 25))
         mat = pybamm.Vector(np.ones(120), domain=["negative particle"])
         disc = tests.get_p2d_discretisation_for_testing()
         with self.assertRaisesRegex(
-            ValueError, "variable shape does not match domain shape"
+            ValueError, "3D variable shape does not match domain shape"
         ):
             pybamm.ProcessedVariable(mat, t, y, disc.mesh)