Repository: tinosulzer/PyBaMM

Issue:
Move all test statments of the discretisation in to tests.shared
Reformulate tests in test_discretisation and test_finite_volume to call `tests.shared.testing_disc` and tests.shared.finite_volume_1D_macro_disc` etc. instead of defining the disc from the defaults

Discussion:

Patch:
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 8e1d0694..0d6637d2 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -24,12 +24,17 @@ class Discretisation(object):
 
     def __init__(self, mesh, spatial_methods):
         self._mesh = mesh
+        # Unpack macroscale to the constituent subdomains
+        if "macroscale" in spatial_methods.keys():
+            method = spatial_methods["macroscale"]
+            spatial_methods["negative electrode"] = method
+            spatial_methods["separator"] = method
+            spatial_methods["positive electrode"] = method
         self._spatial_methods = {
             dom: method(mesh) for dom, method in spatial_methods.items()
         }
         self._bcs = {}
         self._y_slices = {}
-        self._variables = []
 
     @property
     def mesh(self):
@@ -49,10 +54,10 @@ class Discretisation(object):
         # set boundary conditions (only need key ids for boundary_conditions)
         self._bcs = {key.id: value for key, value in model.boundary_conditions.items()}
         # set variables (we require the full variable not just id)
-        self._variables = list(model.rhs.keys()) + list(model.algebraic.keys())
+        variables = list(model.rhs.keys()) + list(model.algebraic.keys())
 
         # Set the y split for variables
-        self.set_variable_slices()
+        self.set_variable_slices(variables)
 
         # Process initial condtions
         self.process_initial_conditions(model)
@@ -68,13 +73,16 @@ class Discretisation(object):
         # Check that resulting model makes sense
         self.check_model(model)
 
-    def set_variable_slices(self):
+    def set_variable_slices(self, variables):
         """Sets the slicing for variables.
+
+        variables : iterable of :class:`pybamm.Variables`
+            The variables for which to set slices
         """
-        y_slices = {variable.id: None for variable in self._variables}
+        y_slices = {variable.id: None for variable in variables}
         start = 0
         end = 0
-        for variable in self._variables:
+        for variable in variables:
             # If domain is empty then variable has size 1
             if variable.domain == []:
                 end += 1
diff --git a/pybamm/geometry/geometry.py b/pybamm/geometry/geometry.py
index 8b47ee8c..019dbc4f 100644
--- a/pybamm/geometry/geometry.py
+++ b/pybamm/geometry/geometry.py
@@ -18,8 +18,20 @@ class Geometry(dict):
      custom_geometry : dict containing any extra user defined geometry
      """
 
-    # left empty for now but I think we should have methods to
-    # add Geometries together (i.e. choose 3D macro with 1D micro etc)
+    def __init__(self, *geometries, custom_geometry={}):
+        for geometry in geometries:
+            if geometry == "1D macro":
+                geometry = Geometry1DMacro()
+            elif geometry == "3D macro":
+                geometry = Geometry3DMacro()
+            elif geometry == "1D micro":
+                geometry = Geometry1DMicro()
+            # avoid combining geometries that clash
+            if any([k in self.keys() for k in geometry.keys()]):
+                raise ValueError("trying to overwrite existing geometry")
+            self.update(geometry)
+        # Allow overwriting with a custom geometry
+        self.update(custom_geometry)
 
 
 class Geometry1DMacro(Geometry):
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index b81d42ac..e928f582 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -48,35 +48,30 @@ class BaseModel(object):
         self._concatenated_rhs = None
         self._concatenated_initial_conditions = None
 
-        # Default parameter values, discretisation and solver
+        # Default parameter values, geometry, submesh, spatial methods and solver
         self.default_parameter_values = pybamm.ParameterValues(
             "input/parameters/lithium-ion/parameters/LCO.csv"
         )
-
-        self.default_geometry = pybamm.Geometry1DMacro()
-
-        self.default_parameter_values.process_geometry(self.default_geometry)
-        # provide mesh properties
-        submesh_pts = {
+        self.default_geometry = pybamm.Geometry("1D macro", "1D micro")
+        self.default_submesh_pts = {
             "negative electrode": {"x": 40},
             "separator": {"x": 25},
             "positive electrode": {"x": 35},
+            "negative particle": {"r": 10},
+            "positive particle": {"r": 10},
         }
-        submesh_types = {
+        self.default_submesh_types = {
             "negative electrode": pybamm.Uniform1DSubMesh,
             "separator": pybamm.Uniform1DSubMesh,
             "positive electrode": pybamm.Uniform1DSubMesh,
+            "negative particle": pybamm.Uniform1DSubMesh,
+            "positive particle": pybamm.Uniform1DSubMesh,
         }
-
         self.default_spatial_methods = {
-            "negative electrode": pybamm.FiniteVolume,
-            "separator": pybamm.FiniteVolume,
-            "positive electrode": pybamm.FiniteVolume,
+            "macroscale": pybamm.FiniteVolume,
+            "negative particle": pybamm.FiniteVolume,
+            "positive particle": pybamm.FiniteVolume,
         }
-        self.mesh = pybamm.Mesh(self.default_geometry, submesh_types, submesh_pts)
-        self.default_discretisation = pybamm.Discretisation(
-            self.mesh, self.default_spatial_methods
-        )
         self.default_solver = pybamm.ScipySolver(method="RK45")
 
     def _set_dict(self, dict, name):
diff --git a/pybamm/models/li_ion/single_particle_model.py b/pybamm/models/li_ion/single_particle_model.py
index fcf2e535..830b0c12 100644
--- a/pybamm/models/li_ion/single_particle_model.py
+++ b/pybamm/models/li_ion/single_particle_model.py
@@ -30,27 +30,6 @@ class SPM(pybamm.BaseModel):
     def __init__(self):
         super().__init__()
 
-        # Overwrite default geometry
-        # NOTE: Is this the best way/place to do this?
-        self.default_geometry = pybamm.Geometry1DMicro()
-        self.default_parameter_values.process_geometry(self.default_geometry)
-        submesh_pts = {
-            "negative particle": {"r": 10},
-            "positive particle": {"r": 10},
-        }
-        submesh_types = {
-            "negative particle": pybamm.Uniform1DSubMesh,
-            "positive particle": pybamm.Uniform1DSubMesh,
-        }
-        self.default_spatial_methods = {
-            "negative particle": pybamm.FiniteVolume,
-            "positive particle": pybamm.FiniteVolume,
-        }
-        self.mesh = pybamm.Mesh(self.default_geometry, submesh_types, submesh_pts)
-        self.default_discretisation = pybamm.Discretisation(
-            self.mesh, self.default_spatial_methods
-        )
-
         # Variables
         cn = pybamm.Variable("cn", domain="negative particle")
         cp = pybamm.Variable("cp", domain="positive particle")
@@ -76,26 +55,24 @@ class SPM(pybamm.BaseModel):
         cp_init = pybamm.standard_parameters.cp0
 
         # PDE RHS
-        Nn = - gamma_n * D_n(cn) * pybamm.grad(cn)
-        dcndt = - pybamm.div(Nn)
-        Np = - gamma_p * D_p(cp) * pybamm.grad(cp)
-        dcpdt = - pybamm.div(Np)
+        Nn = -gamma_n * D_n(cn) * pybamm.grad(cn)
+        dcndt = -pybamm.div(Nn)
+        Np = -gamma_p * D_p(cp) * pybamm.grad(cp)
+        dcpdt = -pybamm.div(Np)
         self.rhs = {cn: dcndt, cp: dcpdt}
 
         # Boundary conditions
         # Note: this is for constant current discharge only
         self.boundary_conditions = {
-            Nn: {"left": pybamm.Scalar(0),
-                 "right": pybamm.Scalar(1) / ln / beta_n},
-            Np: {"left": pybamm.Scalar(0),
-                 "right": - pybamm.Scalar(1) / lp / beta_p / C_hat_p},
+            Nn: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(1) / ln / beta_n},
+            Np: {
+                "left": pybamm.Scalar(0),
+                "right": -pybamm.Scalar(1) / lp / beta_p / C_hat_p,
+            },
         }
 
         # Initial conditions
-        self.initial_conditions = {
-            cn: cn_init,
-            cp: cp_init,
-        }
+        self.initial_conditions = {cn: cn_init, cp: cp_init}
 
         # Variables
         cn_surf = pybamm.surf(cn)
@@ -103,9 +80,12 @@ class SPM(pybamm.BaseModel):
         gn = m_n * cn_surf ** 0.5 * (1 - cn_surf) ** 0.5
         gp = m_p * C_hat_p * cp_surf ** 0.5 * (1 - cp_surf) ** 0.5
         # linearise BV for now
-        V = (U_p(cp_surf) - U_n(cn_surf)
-             - (2 / Lambda) * (1 / (gp * lp))
-             - (2 / Lambda) * (1 / (gn * ln)))
+        V = (
+            U_p(cp_surf)
+            - U_n(cn_surf)
+            - (2 / Lambda) * (1 / (gp * lp))
+            - (2 / Lambda) * (1 / (gn * ln))
+        )
 
         self.variables = {
             "cn": cn,
@@ -114,3 +94,6 @@ class SPM(pybamm.BaseModel):
             "cp_surf": cp_surf,
             "V": V,
         }
+
+        # Overwrite default solver
+        self.default_solver = pybamm.ScipySolver(method="BDF")


Test patch:
diff --git a/tests/__init__.py b/tests/__init__.py
index b5d1a622..37653849 100644
--- a/tests/__init__.py
+++ b/tests/__init__.py
@@ -6,3 +6,4 @@ from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 
 from .test_models.standard_model_tests import StandardModelTest
+from .shared import get_mesh_for_testing, get_discretisation_for_testing
diff --git a/tests/shared.py b/tests/shared.py
index 6fffc90d..3c7fb6cd 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -6,61 +6,6 @@ import pybamm
 import numpy as np
 
 
-class TestDefaults1DMacro:
-    def __init__(self, npts=0):
-        self.param = pybamm.ParameterValues(
-            base_parameters={"Ln": 0.3, "Ls": 0.3, "Lp": 0.3}
-        )
-
-        self.geometry = pybamm.Geometry1DMacro()
-        self.param.process_geometry(self.geometry)
-
-        self.submesh_pts = {
-            "negative electrode": {"x": 40},
-            "separator": {"x": 25},
-            "positive electrode": {"x": 35},
-        }
-
-        self.submesh_types = {
-            "negative electrode": pybamm.Uniform1DSubMesh,
-            "separator": pybamm.Uniform1DSubMesh,
-            "positive electrode": pybamm.Uniform1DSubMesh,
-        }
-
-        if npts != 0:
-            n = 3 * round(npts / 3)
-            self.submesh_pts = {
-                "negative electrode": {"x": n},
-                "separator": {"x": n},
-                "positive electrode": {"x": n},
-            }
-
-        self.mesh = pybamm.Mesh(self.geometry, self.submesh_types, self.submesh_pts)
-
-        self.spatial_methods = {
-            "negative electrode": SpatialMethodForTesting,
-            "separator": SpatialMethodForTesting,
-            "positive electrode": SpatialMethodForTesting,
-        }
-
-
-class TestDefaults1DParticle:
-    def __init__(self, n):
-        self.geometry = {
-            "negative particle": {
-                "r": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}
-            }
-        }
-        self.param = pybamm.ParameterValues(base_parameters={})
-        self.param.process_geometry(self.geometry)
-        self.submesh_pts = {"negative particle": {"r": n}}
-        self.submesh_types = {"negative particle": pybamm.Uniform1DSubMesh}
-
-        self.mesh = pybamm.Mesh(self.geometry, self.submesh_types, self.submesh_pts)
-
-        self.spatial_methods = {"negative particle": pybamm.FiniteVolume}
-
-
 class SpatialMethodForTesting(pybamm.SpatialMethod):
     """Identity operators, no boundary conditions."""
 
@@ -103,3 +48,48 @@ class SpatialMethodForTesting(pybamm.SpatialMethod):
 
     def compute_diffusivity(self, symbol):
         return symbol
+
+
+def get_mesh_for_testing(npts=None):
+    param = pybamm.ParameterValues(base_parameters={"Ln": 0.3, "Ls": 0.3, "Lp": 0.3})
+
+    geometry = pybamm.Geometry("1D macro", "1D micro")
+    param.process_geometry(geometry)
+
+    submesh_types = {
+        "negative electrode": pybamm.Uniform1DSubMesh,
+        "separator": pybamm.Uniform1DSubMesh,
+        "positive electrode": pybamm.Uniform1DSubMesh,
+        "negative particle": pybamm.Uniform1DSubMesh,
+        "positive particle": pybamm.Uniform1DSubMesh,
+    }
+
+    if npts is None:
+        submesh_pts = {
+            "negative electrode": {"x": 40},
+            "separator": {"x": 25},
+            "positive electrode": {"x": 35},
+            "negative particle": {"r": 10},
+            "positive particle": {"r": 10},
+        }
+    else:
+        n = 3 * round(npts / 3)
+        submesh_pts = {
+            "negative electrode": {"x": n},
+            "separator": {"x": n},
+            "positive electrode": {"x": n},
+            "negative particle": {"r": npts},
+            "positive particle": {"r": npts},
+        }
+    return pybamm.Mesh(geometry, submesh_types, submesh_pts)
+
+
+def get_discretisation_for_testing(npts=None):
+    mesh = get_mesh_for_testing(npts)
+    spatial_methods = {
+        "macroscale": SpatialMethodForTesting,
+        "negative particle": SpatialMethodForTesting,
+        "positive particle": SpatialMethodForTesting,
+    }
+
+    return pybamm.Discretisation(mesh, spatial_methods)
diff --git a/tests/test_discretisations/test_base_discretisations.py b/tests/test_discretisations/test_base_discretisations.py
index 3eda96e3..9bd1e0ff 100644
--- a/tests/test_discretisations/test_base_discretisations.py
+++ b/tests/test_discretisations/test_base_discretisations.py
@@ -5,7 +5,7 @@ import pybamm
 
 import numpy as np
 import unittest
-import tests.shared as shared
+from tests import get_mesh_for_testing, get_discretisation_for_testing
 
 
 class TestDiscretise(unittest.TestCase):
@@ -21,8 +21,7 @@ class TestDiscretise(unittest.TestCase):
         }
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
 
         disc._y_slices = {c.id: slice(0, 1), a.id: slice(2, 3), b.id: slice(3, 4)}
         result = disc._concatenate_init(initial_conditions)
@@ -38,15 +37,14 @@ class TestDiscretise(unittest.TestCase):
 
     def test_discretise_slicing(self):
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
+        mesh = get_mesh_for_testing()
         spatial_methods = {}
-        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-        mesh = disc.mesh
+        disc = pybamm.Discretisation(mesh, spatial_methods)
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         c = pybamm.Variable("c", domain=whole_cell)
-        disc._variables = [c]
-        disc.set_variable_slices()
+        variables = [c]
+        disc.set_variable_slices(variables)
 
         self.assertEqual(disc._y_slices, {c.id: slice(0, 100)})
 
@@ -59,8 +57,8 @@ class TestDiscretise(unittest.TestCase):
         # Several variables
         d = pybamm.Variable("d", domain=whole_cell)
         jn = pybamm.Variable("jn", domain=["negative electrode"])
-        disc._variables = [c, d, jn]
-        disc.set_variable_slices()
+        variables = [c, d, jn]
+        disc.set_variable_slices(variables)
 
         self.assertEqual(
             disc._y_slices,
@@ -75,9 +73,9 @@ class TestDiscretise(unittest.TestCase):
 
     def test_process_symbol_base(self):
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
+        mesh = get_mesh_for_testing()
         spatial_methods = {}
-        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
+        disc = pybamm.Discretisation(mesh, spatial_methods)
 
         # variable
         var = pybamm.Variable("var")
@@ -135,8 +133,7 @@ class TestDiscretise(unittest.TestCase):
         expression = (scal1 * (scal3 + var2)) / ((var1 - scal4) + scal2)
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
 
         disc._y_slices = {var1.id: slice(53), var2.id: slice(53, 59)}
         exp_disc = disc.process_symbol(expression)
@@ -176,11 +173,10 @@ class TestDiscretise(unittest.TestCase):
         # create discretisation
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=whole_cell)
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         mesh = disc.mesh
-        disc._variables = [var]
-        disc.set_variable_slices()
+        variables = [var]
+        disc.set_variable_slices(variables)
 
         # Simple expressions
         for eqn in [pybamm.grad(var), pybamm.div(var)]:
@@ -238,17 +234,15 @@ class TestDiscretise(unittest.TestCase):
         }
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
 
         y = combined_submesh.nodes ** 2
-        disc._variables = list(rhs.keys())
         disc._bcs = boundary_conditions
 
-        disc.set_variable_slices()
+        disc.set_variable_slices(list(rhs.keys()))
         # rhs - grad and div are identity operators here
         processed_rhs = disc.process_dict(rhs)
         np.testing.assert_array_equal(y, processed_rhs[c].evaluate(None, y))
@@ -273,8 +267,8 @@ class TestDiscretise(unittest.TestCase):
             [combined_submesh.nodes ** 2, mesh["negative electrode"].nodes ** 4]
         )
 
-        disc._variables = list(rhs.keys())
-        disc.set_variable_slices()
+        variables = list(rhs.keys())
+        disc.set_variable_slices(variables)
         # rhs
         processed_rhs = disc.process_dict(rhs)
         np.testing.assert_array_equal(
@@ -306,8 +300,7 @@ class TestDiscretise(unittest.TestCase):
         model.variables = {"c": c, "N": N}
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -385,8 +378,7 @@ class TestDiscretise(unittest.TestCase):
         model.variables = {"c": c, "N": N, "d": d}
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
         disc.process_model(model)
@@ -444,9 +436,7 @@ class TestDiscretise(unittest.TestCase):
         var = pybamm.Variable("var")
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -487,8 +477,7 @@ class TestDiscretise(unittest.TestCase):
         c = pybamm.Symbol("c")
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
 
         conc = disc.concatenate(a, b, c)
         self.assertIsInstance(conc, pybamm.NumpyModelConcatenation)
@@ -499,12 +488,11 @@ class TestDiscretise(unittest.TestCase):
         b = pybamm.Scalar(4, domain=["positive electrode"])
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
-        disc._variables = [pybamm.Variable("var", domain=whole_cell)]
-        disc.set_variable_slices()
+        variables = [pybamm.Variable("var", domain=whole_cell)]
+        disc.set_variable_slices(variables)
 
         eqn = pybamm.Concatenation(a, b)
         eqn_disc = disc.process_symbol(eqn)
@@ -524,8 +512,7 @@ class TestDiscretise(unittest.TestCase):
         x3 = 3 * pybamm.Space(["negative electrode"])
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
 
         # space
         x1_disc = disc.process_symbol(x1)
diff --git a/tests/test_expression_tree/test_concatenations.py b/tests/test_expression_tree/test_concatenations.py
index 21cca07d..a6c89734 100644
--- a/tests/test_expression_tree/test_concatenations.py
+++ b/tests/test_expression_tree/test_concatenations.py
@@ -2,7 +2,7 @@
 # Tests for the Concatenation class and subclasses
 #
 import pybamm
-import tests.shared as shared
+from tests import get_discretisation_for_testing
 import numpy as np
 import unittest
 
@@ -79,8 +79,7 @@ class TestConcatenations(unittest.TestCase):
 
     def test_numpy_domain_concatenation(self):
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
         a_dom = ["negative electrode"]
diff --git a/tests/test_expression_tree/test_unary_operators.py b/tests/test_expression_tree/test_unary_operators.py
index 304ba568..452833d6 100644
--- a/tests/test_expression_tree/test_unary_operators.py
+++ b/tests/test_expression_tree/test_unary_operators.py
@@ -4,7 +4,7 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 import pybamm
-import tests.shared as shared
+from tests import get_discretisation_for_testing
 import unittest
 import numpy as np
 
@@ -99,8 +99,7 @@ class TestUnaryOperators(unittest.TestCase):
 
     def test_numpy_broadcast(self):
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
diff --git a/tests/test_geometry/test_geometry.py b/tests/test_geometry/test_geometry.py
index cf7236ea..3cea6591 100644
--- a/tests/test_geometry/test_geometry.py
+++ b/tests/test_geometry/test_geometry.py
@@ -77,3 +77,87 @@ class TestGeometry3DMacro(unittest.TestCase):
                 self.assertIsInstance(spatial_var, pybamm.IndependentVariable)
                 for spatial_var in spatial_vars.keys()
             )
+
+
+class TestGeometry(unittest.TestCase):
+    def test_combine_geometries(self):
+        geometry1Dmacro = pybamm.Geometry1DMacro()
+        geometry1Dmicro = pybamm.Geometry1DMicro()
+        geometry = pybamm.Geometry(geometry1Dmacro, geometry1Dmicro)
+        self.assertEqual(
+            list(geometry.keys()),
+            [
+                "negative electrode",
+                "separator",
+                "positive electrode",
+                "negative particle",
+                "positive particle",
+            ],
+        )
+
+        # update with custom geometry
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        x = pybamm.IndependentVariable("x", whole_cell)
+        custom_geometry = {
+            "negative electrode": {
+                x: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}
+            }
+        }
+        geometry = pybamm.Geometry(
+            geometry1Dmacro, geometry1Dmicro, custom_geometry=custom_geometry
+        )
+        self.assertEqual(
+            geometry["negative electrode"], custom_geometry["negative electrode"]
+        )
+
+    def test_combine_geometries_3D(self):
+        geometry3Dmacro = pybamm.Geometry3DMacro()
+        geometry1Dmicro = pybamm.Geometry1DMicro()
+        geometry = pybamm.Geometry(geometry3Dmacro, geometry1Dmicro)
+        self.assertEqual(
+            list(geometry.keys()),
+            [
+                "negative electrode",
+                "separator",
+                "positive electrode",
+                "negative particle",
+                "positive particle",
+            ],
+        )
+
+        with self.assertRaises(ValueError):
+            geometry1Dmacro = pybamm.Geometry1DMacro()
+            geometry = pybamm.Geometry(geometry3Dmacro, geometry1Dmacro)
+
+    def test_combine_geometries_strings(self):
+        geometry = pybamm.Geometry("1D macro", "1D micro")
+        self.assertEqual(
+            list(geometry.keys()),
+            [
+                "negative electrode",
+                "separator",
+                "positive electrode",
+                "negative particle",
+                "positive particle",
+            ],
+        )
+        geometry = pybamm.Geometry("3D macro", "1D micro")
+        self.assertEqual(
+            list(geometry.keys()),
+            [
+                "negative electrode",
+                "separator",
+                "positive electrode",
+                "negative particle",
+                "positive particle",
+            ],
+        )
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/test_models/standard_model_tests.py b/tests/test_models/standard_model_tests.py
index 2c303aa2..3917eba8 100644
--- a/tests/test_models/standard_model_tests.py
+++ b/tests/test_models/standard_model_tests.py
@@ -3,16 +3,25 @@
 #
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
+import pybamm
+
 import numpy as np
 
 
 class StandardModelTest(object):
     def __init__(self, model):
         self.model = model
-        # Set defaults
+        # Set default parameters
         self.param = model.default_parameter_values
-        self.geometry = model.default_geometry
-        self.disc = model.default_discretisation
+        # Process geometry
+        self.param.process_geometry(model.default_geometry)
+        geometry = model.default_geometry
+        # Set default discretisation
+        mesh = pybamm.Mesh(
+            geometry, model.default_submesh_types, model.default_submesh_pts
+        )
+        self.disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+        # Set default solver
         self.solver = model.default_solver
 
     def test_processing_parameters(self, param=None):
diff --git a/tests/test_models/test_interface.py b/tests/test_models/test_interface.py
index 2b36824d..ad2f61b9 100644
--- a/tests/test_models/test_interface.py
+++ b/tests/test_models/test_interface.py
@@ -4,7 +4,7 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 import pybamm
-from tests import shared
+from tests import get_discretisation_for_testing
 
 import unittest
 import numpy as np
@@ -53,8 +53,7 @@ class TestHomogeneousReaction(unittest.TestCase):
         param = pybamm.ParameterValues(
             "input/parameters/lithium-ion/parameters/LCO.csv"
         )
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
 
         rxn = pybamm.interface.homogeneous_reaction()
 
@@ -162,11 +161,10 @@ class TestButlerVolmerLeadAcid(unittest.TestCase):
         param_bv_p = param.process_symbol(bv_p)
 
         # discretise
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         mesh = disc.mesh
-        disc._variables = [self.cn, self.cp, self.phin, self.phip]
-        disc.set_variable_slices()
+        variables = [self.cn, self.cp, self.phin, self.phip]
+        disc.set_variable_slices(variables)
         processed_bv_n = disc.process_symbol(param_bv_n)
         processed_bv_p = disc.process_symbol(param_bv_p)
 
@@ -196,12 +194,11 @@ class TestButlerVolmerLeadAcid(unittest.TestCase):
         param_bv_whole = param.process_symbol(bv_whole)
 
         # discretise
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
-        disc._variables = [self.cn, self.cp, self.phin, self.phip]
-        disc.set_variable_slices()
+        variables = [self.cn, self.cp, self.phin, self.phip]
+        disc.set_variable_slices(variables)
         processed_bv_whole = disc.process_symbol(param_bv_whole)
 
         # test
@@ -278,12 +275,11 @@ class TestExchangeCurrentDensity(unittest.TestCase):
         param_j0p = param.process_symbol(j0p)
 
         # discretise
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
-        disc._variables = [self.cn, self.cp]
-        disc.set_variable_slices()
+        variables = [self.cn, self.cp]
+        disc.set_variable_slices(variables)
         processed_j0n = disc.process_symbol(param_j0n)
         processed_j0p = disc.process_symbol(param_j0p)
 
diff --git a/tests/test_models/test_lead_acid.py b/tests/test_models/test_lead_acid.py
index 9677f20d..b58460d9 100644
--- a/tests/test_models/test_lead_acid.py
+++ b/tests/test_models/test_lead_acid.py
@@ -19,15 +19,9 @@ class TestLeadAcidLOQS(unittest.TestCase):
 
     def test_solution(self):
         model = pybamm.lead_acid.LOQS()
-
-        # process parameter values, discretise and solve
-        model.default_parameter_values.process_model(model)
-        disc = model.default_discretisation
-        disc.process_model(model)
-        t_eval = np.linspace(0, 1, 100)
-        solver = model.default_solver
-        solver.solve(model, t_eval)
-        T, Y = solver.t, solver.y
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+        T, Y = modeltest.solver.t, modeltest.solver.y
 
         # check output
         # make sure concentration and voltage are monotonically decreasing
diff --git a/tests/test_models/test_li_ion.py b/tests/test_models/test_li_ion.py
index e0a1e59c..c289c5f8 100644
--- a/tests/test_models/test_li_ion.py
+++ b/tests/test_models/test_li_ion.py
@@ -19,14 +19,9 @@ class TestLiIonSPM(unittest.TestCase):
 
     def test_surface_concentrartion(self):
         model = pybamm.li_ion.SPM()
-        params = model.default_parameter_values
-        params.process_model(model)
-        disc = model.default_discretisation
-        disc.process_model(model)
-        t_eval = np.linspace(0, 1, 100)
-        solver = model.default_solver
-        solver.solve(model, t_eval)
-        T, Y = solver.t, solver.y
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+        T, Y = modeltest.solver.t, modeltest.solver.y
 
         # check surface concentration decreases in negative particle and
         # increases in positive particle for discharge
@@ -38,3 +33,12 @@ class TestLiIonSPM(unittest.TestCase):
             model.variables["cp_surf"].evaluate(T, Y)[:, :-1],
             model.variables["cp_surf"].evaluate(T, Y)[:, 1:],
         )
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/test_models/test_simple_ode_model.py b/tests/test_models/test_simple_ode_model.py
index b6523c80..20d41292 100644
--- a/tests/test_models/test_simple_ode_model.py
+++ b/tests/test_models/test_simple_ode_model.py
@@ -18,17 +18,13 @@ class TestSimpleODEModel(unittest.TestCase):
         modeltest.test_all()
 
     def test_solution(self):
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
         model = pybamm.SimpleODEModel()
-
-        # discretise and solve
-        disc = model.default_discretisation
-        combined_submesh = disc.mesh.combine_submeshes(*whole_cell)
-        disc.process_model(model)
-        t_eval = np.linspace(0, 1, 100)
-        solver = model.default_solver
-        solver.solve(model, t_eval)
-        T, Y = solver.t, solver.y
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+        T, Y = modeltest.solver.t, modeltest.solver.y
+        mesh = modeltest.disc.mesh
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        combined_submesh = mesh.combine_submeshes(*whole_cell)
 
         # check output
         np.testing.assert_array_almost_equal(
@@ -40,9 +36,9 @@ class TestSimpleODEModel(unittest.TestCase):
         )
         np.testing.assert_array_almost_equal(
             model.variables["c broadcasted"].evaluate(T, Y),
-            np.ones(
-                sum([disc.mesh[d].npts for d in ["negative electrode", "separator"]])
-            )[:, np.newaxis]
+            np.ones(sum([mesh[d].npts for d in ["negative electrode", "separator"]]))[
+                :, np.newaxis
+            ]
             * np.exp(-T),
         )
 
diff --git a/tests/test_parameters/test_standard_parameters_lead_acid.py b/tests/test_parameters/test_standard_parameters_lead_acid.py
index 39363eb1..10530990 100644
--- a/tests/test_parameters/test_standard_parameters_lead_acid.py
+++ b/tests/test_parameters/test_standard_parameters_lead_acid.py
@@ -4,7 +4,7 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 import pybamm
-from tests import shared
+from tests import get_discretisation_for_testing
 
 import unittest
 
@@ -60,8 +60,7 @@ class TestStandardParametersLeadAcid(unittest.TestCase):
         parameter_values = pybamm.ParameterValues(
             "input/parameters/lead-acid/default.csv", {"current scale": 1}
         )
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        disc = get_discretisation_for_testing()
         processed_s = disc.process_symbol(parameter_values.process_symbol(s))
 
         # test output
diff --git a/tests/test_solvers/test_scikits_solvers.py b/tests/test_solvers/test_scikits_solvers.py
index 6d0e37b8..c2bbcbd6 100644
--- a/tests/test_solvers/test_scikits_solvers.py
+++ b/tests/test_solvers/test_scikits_solvers.py
@@ -5,7 +5,7 @@ from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 import pybamm
 from pybamm.solvers.scikits_ode_solver import scikits_odes_spec
-
+from tests import StandardModelTest
 
 import unittest
 import numpy as np
@@ -73,7 +73,7 @@ class TestScikitsSolver(unittest.TestCase):
         model.rhs = {var: 0.1 * var}
         model.initial_conditions = {var: 1}
         model.initial_conditions_ydot = {var: 0.1}
-        disc = model.default_discretisation
+        disc = StandardModelTest(model).disc
         disc.process_model(model)
 
         # Solve
@@ -93,7 +93,7 @@ class TestScikitsSolver(unittest.TestCase):
         model.algebraic = {var2: 2 * var1 - var2}
         model.initial_conditions = {var1: 1, var2: 2}
         model.initial_conditions_ydot = {var1: 0.1, var2: 0.2}
-        disc = model.default_discretisation
+        disc = StandardModelTest(model).disc
         disc.process_model(model)
 
         # Solve
diff --git a/tests/test_solvers/test_scipy_solver.py b/tests/test_solvers/test_scipy_solver.py
index 469fe512..b7901174 100644
--- a/tests/test_solvers/test_scipy_solver.py
+++ b/tests/test_solvers/test_scipy_solver.py
@@ -7,7 +7,7 @@ import pybamm
 
 import unittest
 import numpy as np
-import tests.shared as shared
+from tests import get_mesh_for_testing
 
 
 class TestScipySolver(unittest.TestCase):
@@ -45,13 +45,9 @@ class TestScipySolver(unittest.TestCase):
         # No need to set parameters; can use base discretisation (no spatial operators)
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        spatial_methods = {
-            "negative electrode": pybamm.FiniteVolume,
-            "separator": pybamm.FiniteVolume,
-            "positive electrode": pybamm.FiniteVolume,
-        }
-        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
         disc.process_model(model)
         # Solve
         solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
diff --git a/tests/test_spatial_methods/test_finite_volume.py b/tests/test_spatial_methods/test_finite_volume.py
index ce84a225..76590990 100644
--- a/tests/test_spatial_methods/test_finite_volume.py
+++ b/tests/test_spatial_methods/test_finite_volume.py
@@ -2,7 +2,7 @@
 # Test for the operator class
 #
 import pybamm
-import tests.shared as shared
+from tests import get_mesh_for_testing
 
 import numpy as np
 import unittest
@@ -34,18 +34,16 @@ class TestFiniteVolume(unittest.TestCase):
 
     def test_surface_value(self):
         # create discretisation
-        defaults = shared.TestDefaults1DParticle(10)
+        mesh = get_mesh_for_testing()
         spatial_methods = {"negative particle": pybamm.FiniteVolume}
-        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-        mesh = disc.mesh
+        disc = pybamm.Discretisation(mesh, spatial_methods)
 
         combined_submesh = mesh.combine_submeshes("negative particle")
 
         # create variable
         var = pybamm.Variable("var", domain="negative particle")
         surf_eqn = pybamm.surf(var)
-        disc._variables = [var]
-        disc.set_variable_slices()
+        disc.set_variable_slices([var])
         surf_eqn_disc = disc.process_symbol(surf_eqn)
 
         # check constant extrapolates to constant
@@ -62,21 +60,15 @@ class TestFiniteVolume(unittest.TestCase):
         whole_cell = ["negative electrode", "separator", "positive electrode"]
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        spatial_methods = {
-            "negative electrode": pybamm.FiniteVolume,
-            "separator": pybamm.FiniteVolume,
-            "positive electrode": pybamm.FiniteVolume,
-        }
-        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-        mesh = disc.mesh
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
 
         # Discretise some equations where averaging is needed
         var = pybamm.Variable("var", domain=whole_cell)
-        disc._variables = [var]
-        disc.set_variable_slices()
+        disc.set_variable_slices([var])
         y_test = np.ones_like(combined_submesh.nodes)
         for eqn in [
             var * pybamm.grad(var),
@@ -121,24 +113,18 @@ class TestFiniteVolume(unittest.TestCase):
         # Set up
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        spatial_methods = {
-            "negative electrode": pybamm.FiniteVolume,
-            "separator": pybamm.FiniteVolume,
-            "positive electrode": pybamm.FiniteVolume,
-        }
-        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-        mesh = disc.mesh
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
 
         # Add ghost nodes
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=whole_cell)
-        disc._variables = [var]
-        disc.set_variable_slices()
+        disc.set_variable_slices([var])
         discretised_symbol = pybamm.StateVector(disc._y_slices[var.id])
         lbc = pybamm.Scalar(0)
         rbc = pybamm.Scalar(3)
-        symbol_plus_ghost = pybamm.FiniteVolume(defaults.mesh).add_ghost_nodes(
+        symbol_plus_ghost = pybamm.FiniteVolume(mesh).add_ghost_nodes(
             discretised_symbol, lbc, rbc
         )
 
@@ -172,14 +158,9 @@ class TestFiniteVolume(unittest.TestCase):
         """
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        spatial_methods = {
-            "negative electrode": pybamm.FiniteVolume,
-            "separator": pybamm.FiniteVolume,
-            "positive electrode": pybamm.FiniteVolume,
-        }
-        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-        mesh = disc.mesh
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
 
@@ -194,8 +175,7 @@ class TestFiniteVolume(unittest.TestCase):
         }
         disc._bcs = boundary_conditions
 
-        disc._variables = [var]
-        disc.set_variable_slices()
+        disc.set_variable_slices([var])
         grad_eqn_disc = disc.process_symbol(grad_eqn)
 
         constant_y = np.ones_like(combined_submesh.nodes)
@@ -229,10 +209,9 @@ class TestFiniteVolume(unittest.TestCase):
         Test grad and div with Dirichlet boundary conditions (applied by grad on var)
         """
         # create discretisation
-        defaults = shared.TestDefaults1DParticle(10)
+        mesh = get_mesh_for_testing()
         spatial_methods = {"negative particle": pybamm.FiniteVolume}
-        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-        mesh = disc.mesh
+        disc = pybamm.Discretisation(mesh, spatial_methods)
 
         combined_submesh = mesh.combine_submeshes("negative particle")
 
@@ -249,8 +228,7 @@ class TestFiniteVolume(unittest.TestCase):
 
         disc._bcs = boundary_conditions
 
-        disc._variables = [var]
-        disc.set_variable_slices()
+        disc.set_variable_slices([var])
         grad_eqn_disc = disc.process_symbol(grad_eqn)
 
         constant_y = np.ones_like(combined_submesh.nodes)
@@ -290,22 +268,16 @@ class TestFiniteVolume(unittest.TestCase):
         whole_cell = ["negative electrode", "separator", "positive electrode"]
 
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        spatial_methods = {
-            "negative electrode": pybamm.FiniteVolume,
-            "separator": pybamm.FiniteVolume,
-            "positive electrode": pybamm.FiniteVolume,
-        }
-        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-        mesh = disc.mesh
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
 
         # grad
         var = pybamm.Variable("var", domain=whole_cell)
         grad_eqn = pybamm.grad(var)
-        disc._variables = [var]
-        disc.set_variable_slices()
+        disc.set_variable_slices([var])
         grad_eqn_disc = disc.process_symbol(grad_eqn)
 
         constant_y = np.ones_like(combined_submesh.nodes)
@@ -337,18 +309,16 @@ class TestFiniteVolume(unittest.TestCase):
         """Test grad and div with Neumann boundary conditions (applied by div on N)"""
 
         # create discretisation
-        defaults = shared.TestDefaults1DParticle(10)
+        mesh = get_mesh_for_testing()
         spatial_methods = {"negative particle": pybamm.FiniteVolume}
-        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-        mesh = disc.mesh
+        disc = pybamm.Discretisation(mesh, spatial_methods)
 
         combined_submesh = mesh.combine_submeshes("negative particle")
 
         # grad
         var = pybamm.Variable("var", domain="negative particle")
         grad_eqn = pybamm.grad(var)
-        disc._variables = [var]
-        disc.set_variable_slices()
+        disc.set_variable_slices([var])
         grad_eqn_disc = disc.process_symbol(grad_eqn)
 
         constant_y = np.ones_like(combined_submesh.nodes)
@@ -384,14 +354,9 @@ class TestFiniteVolume(unittest.TestCase):
         Test grad and div with Dirichlet boundary conditions (applied by grad on var)
         """
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        spatial_methods = {
-            "negative electrode": pybamm.FiniteVolume,
-            "separator": pybamm.FiniteVolume,
-            "positive electrode": pybamm.FiniteVolume,
-        }
-        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-        mesh = disc.mesh
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
 
         mesh.add_ghost_meshes()
         disc.mesh.add_ghost_meshes()
@@ -404,8 +369,7 @@ class TestFiniteVolume(unittest.TestCase):
         }
         disc._bcs = boundary_conditions
 
-        disc._variables = [var]
-        disc.set_variable_slices()
+        disc.set_variable_slices([var])
 
         grad_eqn_disc = disc.process_symbol(grad_eqn)
 
@@ -449,14 +413,9 @@ class TestFiniteVolume(unittest.TestCase):
             # Set up discretisation
             n = 3 * round(n / 3)
             # create discretisation
-            defaults = shared.TestDefaults1DMacro(n)
-            spatial_methods = {
-                "negative electrode": pybamm.FiniteVolume,
-                "separator": pybamm.FiniteVolume,
-                "positive electrode": pybamm.FiniteVolume,
-            }
-            disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-            mesh = disc.mesh
+            mesh = get_mesh_for_testing(n)
+            spatial_methods = {"macroscale": pybamm.FiniteVolume}
+            disc = pybamm.Discretisation(mesh, spatial_methods)
 
             combined_submesh = mesh.combine_submeshes(*whole_cell)
             # Define exact solutions
@@ -464,8 +423,7 @@ class TestFiniteVolume(unittest.TestCase):
             grad_exact = np.cos(combined_submesh.edges[1:-1])
 
             # Discretise and evaluate
-            disc._variables = [var]
-            disc.set_variable_slices()
+            disc.set_variable_slices([var])
             grad_eqn_disc = disc.process_symbol(grad_eqn)
             grad_approx = grad_eqn_disc.evaluate(None, y)
 
@@ -495,14 +453,9 @@ class TestFiniteVolume(unittest.TestCase):
         # Function for convergence testing
         def get_l2_error(n):
             # create discretisation
-            defaults = shared.TestDefaults1DMacro(n)
-            spatial_methods = {
-                "negative electrode": pybamm.FiniteVolume,
-                "separator": pybamm.FiniteVolume,
-                "positive electrode": pybamm.FiniteVolume,
-            }
-            disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-            mesh = disc.mesh
+            mesh = get_mesh_for_testing(n)
+            spatial_methods = {"macroscale": pybamm.FiniteVolume}
+            disc = pybamm.Discretisation(mesh, spatial_methods)
             disc._bcs = boundary_conditions
 
             # Set up discretisation
@@ -517,8 +470,7 @@ class TestFiniteVolume(unittest.TestCase):
             grad_exact = np.cos(combined_submesh.edges)
 
             # Discretise and evaluate
-            disc._variables = [var]
-            disc.set_variable_slices()
+            disc.set_variable_slices([var])
             grad_eqn_disc = disc.process_symbol(grad_eqn)
             grad_approx = grad_eqn_disc.evaluate(None, y)
 
@@ -548,15 +500,10 @@ class TestFiniteVolume(unittest.TestCase):
         def get_l2_error(n):
 
             # create discretisation
-            defaults = shared.TestDefaults1DMacro(n)
-            spatial_methods = {
-                "negative electrode": pybamm.FiniteVolume,
-                "separator": pybamm.FiniteVolume,
-                "positive electrode": pybamm.FiniteVolume,
-            }
-            disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
+            mesh = get_mesh_for_testing(n)
+            spatial_methods = {"macroscale": pybamm.FiniteVolume}
+            disc = pybamm.Discretisation(mesh, spatial_methods)
             disc._bcs = boundary_conditions
-            mesh = disc.mesh
 
             whole_cell = ["negative electrode", "separator", "positive electrode"]
             combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -566,8 +513,7 @@ class TestFiniteVolume(unittest.TestCase):
             div_exact_internal = -np.sin(combined_submesh.nodes[1:-1])
 
             # Discretise and evaluate
-            disc._variables = [var]
-            disc.set_variable_slices()
+            disc.set_variable_slices([var])
             div_eqn_disc = disc.process_symbol(div_eqn)
             div_approx_internal = div_eqn_disc.evaluate(None, y)[1:-1]
 
@@ -599,14 +545,9 @@ class TestFiniteVolume(unittest.TestCase):
         def get_l2_error(n):
             whole_cell = ["negative electrode", "separator", "positive electrode"]
             # create discretisation
-            defaults = shared.TestDefaults1DMacro(n)
-            spatial_methods = {
-                "negative electrode": pybamm.FiniteVolume,
-                "separator": pybamm.FiniteVolume,
-                "positive electrode": pybamm.FiniteVolume,
-            }
-            disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
-            mesh = disc.mesh
+            mesh = get_mesh_for_testing(n)
+            spatial_methods = {"macroscale": pybamm.FiniteVolume}
+            disc = pybamm.Discretisation(mesh, spatial_methods)
             disc._bcs = boundary_conditions
 
             combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -616,8 +557,7 @@ class TestFiniteVolume(unittest.TestCase):
             div_exact_internal = -np.sin(combined_submesh.nodes[1:-1])
 
             # Discretise and evaluate
-            disc._variables = [var]
-            disc.set_variable_slices()
+            disc.set_variable_slices([var])
             div_eqn_disc = disc.process_symbol(div_eqn)
             div_approx_internal = div_eqn_disc.evaluate(None, y)[1:-1]
 
@@ -647,10 +587,9 @@ class TestFiniteVolume(unittest.TestCase):
         }
 
         def get_l2_error(n):
-            defaults = shared.TestDefaults1DParticle(n)
-
+            mesh = get_mesh_for_testing(n)
             spatial_methods = {"negative particle": pybamm.FiniteVolume}
-            disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
+            disc = pybamm.Discretisation(mesh, spatial_methods)
             disc._bcs = boundary_conditions
             mesh = disc.mesh["negative particle"]
             r = mesh.nodes
@@ -661,8 +600,8 @@ class TestFiniteVolume(unittest.TestCase):
             exact_internal = exact[1:-1]
 
             # discretise and evaluate
-            disc._variables = [c]
-            disc.set_variable_slices()
+            variables = [c]
+            disc.set_variable_slices(variables)
             eqn_disc = disc.process_symbol(eqn)
             approx_internal = eqn_disc.evaluate(None, y)[1:-1]