Repository: pybamm-team/PyBaMM

Issue:
Changable parameters 
**Summary**
Add ability to have parameter values that change without recreating the tree. 

- Add an extra subclass of symbol `ChangableParameter`. 
- Do not insert any values for `ChangableParameters` when we `params.process(model)` 
- During `params.process(model)` replace each `ChangableParameter` by `changable_parameters[ChangableParameters]` (a dict with keys that are the changable paramerters).
- add an extra optional input to `model.concatenated_rhs.evaluate(t, y, changable_params)` which is a dict of the changed parameters with their values which will update the `changable_parameters` dict. 
- create a model.set_param() which takes in a dict of parameter values and puts them into the correct location in the changable_params vector

 **Motivation**
 Allows parameters to be changed easily without rebuilding the tree. Efficiently changing parameter values is a requirement of the control community. 

Time dependent functions could also be implemented within this approach. 

**Questions** 
It would be possible to have this as the default approach for setting all parameters. We should consider this as it means someone doesn't need to rebuild the discretisation if they realise afterwards that they want to change a parameter value. I also don't think that this will slow thing down to any great extent. We would also only pass in the parameters values that are changed an update the parameters dict. 
What are peoples thought on this? 


Discussion:

Patch:
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 60b5a7a99..0b4902ccf 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -37,8 +37,6 @@ class Discretisation(object):
         self._bcs = {}
         self._y_slices = {}
 
-        mesh.add_ghost_meshes()
-
     @property
     def mesh(self):
         return self._mesh
@@ -305,10 +303,6 @@ class Discretisation(object):
         elif isinstance(symbol, pybamm.Concatenation):
             new_children = [self.process_symbol(child) for child in symbol.children]
             new_symbol = pybamm.DomainConcatenation(new_children, self.mesh)
-
-            if new_symbol.is_constant():
-                value = new_symbol.evaluate()
-                return pybamm.Vector(value)
             return new_symbol
 
         else:
diff --git a/pybamm/expression_tree/broadcasts.py b/pybamm/expression_tree/broadcasts.py
index 5397f6df1..9c1771e27 100644
--- a/pybamm/expression_tree/broadcasts.py
+++ b/pybamm/expression_tree/broadcasts.py
@@ -86,6 +86,13 @@ class NumpyBroadcast(Broadcast):
         # domain)
         self.broadcasting_vector = np.ones(self.broadcasting_vector_size)
 
+        # store mesh
+        self._mesh = mesh
+
+    @property
+    def mesh(self):
+        return self._mesh
+
     def evaluate(self, t=None, y=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         child = self.children[0]
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index a6dcf5ab0..068973c6b 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -111,6 +111,9 @@ class DomainConcatenation(Concatenation):
         # Allow the base class to sort the domains into the correct order
         super().__init__(*children, name="domain concatenation")
 
+        # store mesh
+        self._mesh = mesh
+
         # Check that there is a domain, otherwise the functionality won't work and we
         # should raise a DomainError
         if self.domain == []:
@@ -122,7 +125,7 @@ class DomainConcatenation(Concatenation):
             )
 
         # create dict of domain => slice of final vector
-        self._slices = self.create_slices(self, mesh)
+        self._slices = self.create_slices(self)
 
         # store size of final vector
         self._size = self._slices[self.domain[-1]].stop
@@ -130,15 +133,27 @@ class DomainConcatenation(Concatenation):
         # create disc of domain => slice for each child
         self._children_slices = []
         for child in self.children:
-            self._children_slices.append(self.create_slices(child, mesh))
+            self._children_slices.append(self.create_slices(child))
+
+    @property
+    def mesh(self):
+        return self._mesh
+
+    @property
+    def size(self):
+        return self._size
+
+    @property
+    def shape(self):
+        return (self.size,)
 
-    def create_slices(self, node, mesh):
+    def create_slices(self, node):
         slices = {}
         start = 0
         end = 0
         for dom in node.domain:
-            prim_pts = mesh[dom][0].npts
-            second_pts = len(mesh[dom])
+            prim_pts = self.mesh[dom][0].npts
+            second_pts = len(self.mesh[dom])
             end += prim_pts * second_pts
             slices[dom] = slice(start, end)
             start = end
diff --git a/pybamm/expression_tree/scalar.py b/pybamm/expression_tree/scalar.py
index c03ea53fc..1a3cbc70a 100644
--- a/pybamm/expression_tree/scalar.py
+++ b/pybamm/expression_tree/scalar.py
@@ -29,11 +29,11 @@ class Scalar(pybamm.Symbol):
 
         """
         # set default name if not provided
+        self.value = value
         if name is None:
-            name = str(value)
+            name = str(self.value)
 
         super().__init__(name, domain=domain)
-        self.value = value
 
     @property
     def id(self):
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 7b22fafd2..24c586278 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -354,13 +354,12 @@ class Symbol(anytree.NodeMixin):
         evaluate : evaluate the expression
 
         """
-
         # if any of the nodes are instances of any of these types, then the whole
         # expression depends on either t or y
         search_types = (pybamm.Variable, pybamm.StateVector, pybamm.IndependentVariable)
 
         # do the search, return true if no relevent nodes are found
-        return all([not isinstance(n, search_types) for n in self.pre_order()])
+        return all([not (isinstance(n, search_types)) for n in self.pre_order()])
 
     def evaluates_to_value(self, value):
         """
@@ -379,7 +378,11 @@ class Symbol(anytree.NodeMixin):
         evaluate : evaluate the expression
 
         """
-        return self.evaluates_to_number() and self.evaluate() == value
+        return (
+            self.is_constant()
+            and self.evaluates_to_number()
+            and self.evaluate() == value
+        )
 
     def evaluates_to_number(self):
         """
diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py
index 64909a1d0..f5b82ff9d 100644
--- a/pybamm/meshes/meshes.py
+++ b/pybamm/meshes/meshes.py
@@ -46,6 +46,7 @@ class Mesh(dict):
             self[domain] = [
                 submesh_types[domain](geometry[domain]["primary"], submesh_pts[domain])
             ] * repeats
+        self.add_ghost_meshes()
 
     def combine_submeshes(self, *submeshnames):
         """Combine submeshes into a new submesh, using self.submeshclass
diff --git a/pybamm/models/reaction_diffusion.py b/pybamm/models/reaction_diffusion.py
index 4fe39ff24..2cae8a8e4 100644
--- a/pybamm/models/reaction_diffusion.py
+++ b/pybamm/models/reaction_diffusion.py
@@ -38,3 +38,6 @@ class ReactionDiffusionModel(pybamm.BaseModel):
 
         # Create own model from diffusion model
         self.update(diffusion_model)
+
+        # Add j to variables dict
+        self.variables.update({"Interfacial current density": j})
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index b3d66c27d..3f3eb09dd 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -114,6 +114,28 @@ class ParameterValues(dict):
         for idx, equation in enumerate(model.events):
             model.events[idx] = self.process_symbol(equation)
 
+    def process_discretised_model(self, model, disc):
+        """Process a discretised model.
+        Currently inplace, could be changed to return a new model.
+
+        Parameters
+        ----------
+        model : :class:`pybamm.models.core.BaseModel` (or subclass) instance
+            Model to assign parameter values for
+        disc : :class:`pybamm.Discretisation`
+            The class that was used to discretise
+
+        """
+        # process parameter values for the model
+        self.process_model(model)
+
+        # update discretised quantities using disc
+        model.concatenated_rhs = disc.concatenate(*model.rhs.values())
+        model.concatenated_algebraic = disc.concatenate(*model.algebraic.values())
+        model.concatenated_initial_conditions = disc._concatenate_init(
+            model.initial_conditions
+        ).evaluate(0, None)
+
     def process_geometry(self, geometry):
         """Assign parameter values to a geometry.
             Currently inplace, could be changed to return a new model.
@@ -134,6 +156,7 @@ class ParameterValues(dict):
 
     def process_symbol(self, symbol):
         """Walk through the symbol and replace any Parameter with a Value.
+        Can process a model either before or after discretisation
 
         Parameters
         ----------
@@ -148,7 +171,8 @@ class ParameterValues(dict):
         """
         if isinstance(symbol, pybamm.Parameter):
             value = self.get_parameter_value(symbol)
-            return pybamm.Scalar(value, domain=symbol.domain)
+            # Scalar inherits name (for updating parameters) and domain (for Broadcast)
+            return pybamm.Scalar(value, name=symbol.name, domain=symbol.domain)
 
         elif isinstance(symbol, pybamm.FunctionParameter):
             new_child = self.process_symbol(symbol.children[0])
@@ -161,38 +185,49 @@ class ParameterValues(dict):
                 new_diff_variable = self.process_symbol(symbol.children[0])
                 return function.diff(new_diff_variable)
 
+        elif isinstance(symbol, pybamm.Scalar):
+            # update any Scalar nodes if their name is in the parameter dict (no error)
+            try:
+                value = self.get_parameter_value(symbol)
+            except KeyError:
+                # KeyError -> name not in parameter dict, don't update, return old value
+                value = symbol.value
+            return pybamm.Scalar(value, name=symbol.name, domain=symbol.domain)
+
         elif isinstance(symbol, pybamm.BinaryOperator):
             left, right = symbol.children
             new_left = self.process_symbol(left)
             new_right = self.process_symbol(right)
             return symbol.__class__(new_left, new_right)
 
-        elif isinstance(symbol, pybamm.Broadcast):
-            new_child = self.process_symbol(symbol.children[0])
-            return pybamm.Broadcast(new_child, symbol.domain)
-
-        elif isinstance(symbol, pybamm.Function):
-            new_child = self.process_symbol(symbol.children[0])
-            return pybamm.Function(symbol.func, new_child)
-
-        elif isinstance(symbol, pybamm.Integral):
-            new_child = self.process_symbol(symbol.children[0])
-            return pybamm.Integral(new_child, symbol.integration_variable)
-
-        elif isinstance(symbol, pybamm.BoundaryValue):
-            new_child = self.process_symbol(symbol.children[0])
-            return pybamm.BoundaryValue(new_child, symbol.side)
-
         elif isinstance(symbol, pybamm.UnaryOperator):
             new_child = self.process_symbol(symbol.children[0])
-            return symbol.__class__(new_child)
+            if isinstance(symbol, pybamm.NumpyBroadcast):
+                return pybamm.NumpyBroadcast(new_child, symbol.domain, symbol.mesh)
+            if isinstance(symbol, pybamm.Broadcast):
+                return pybamm.Broadcast(new_child, symbol.domain)
+            elif isinstance(symbol, pybamm.Function):
+                return pybamm.Function(symbol.func, new_child)
+            elif isinstance(symbol, pybamm.Integral):
+                return pybamm.Integral(new_child, symbol.integration_variable)
+            elif isinstance(symbol, pybamm.NodeToEdge):
+                return pybamm.NodeToEdge(new_child, symbol.node_to_edge_function)
+            elif isinstance(symbol, pybamm.BoundaryValue):
+                return pybamm.BoundaryValue(new_child, symbol.side)
+            else:
+                return symbol.__class__(new_child)
 
+        # Concatenations
         elif isinstance(symbol, pybamm.Concatenation):
             new_children = []
             for child in symbol.children:
                 new_child = self.process_symbol(child)
                 new_children.append(new_child)
-            return pybamm.Concatenation(*new_children)
+            if isinstance(symbol, pybamm.DomainConcatenation):
+                return pybamm.DomainConcatenation(new_children, symbol.mesh)
+            else:
+                # Concatenation or NumpyConcatenation
+                return symbol.__class__(*new_children)
 
         else:
             new_symbol = copy.deepcopy(symbol)
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index ec0e7f10d..1354795e6 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -60,17 +60,7 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         See :meth: `pybamm.SpatialMethod.broadcast`
         """
-
-        broadcasted_symbol = pybamm.NumpyBroadcast(symbol, domain, self.mesh)
-
-        # if the broadcasted symbol evaluates to a constant value, replace the
-        # symbol-Vector multiplication with a single array
-        if broadcasted_symbol.is_constant():
-            broadcasted_symbol = pybamm.Array(
-                broadcasted_symbol.evaluate(), domain=broadcasted_symbol.domain
-            )
-
-        return broadcasted_symbol
+        return pybamm.NumpyBroadcast(symbol, domain, self.mesh)
 
     def gradient(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the gradient operator.
@@ -538,6 +528,10 @@ class NodeToEdge(pybamm.SpatialOperator):
         )
         self._node_to_edge_function = node_to_edge_function
 
+    @property
+    def node_to_edge_function(self):
+        return self._node_to_edge_function
+
     def evaluate(self, t=None, y=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         evaluated_child = self.children[0].evaluate(t, y)


Test patch:
diff --git a/tests/shared.py b/tests/shared.py
index 6ed712bea..b1106df3b 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -21,17 +21,7 @@ class SpatialMethodForTesting(pybamm.SpatialMethod):
         return pybamm.Vector(symbol_mesh[0].nodes)
 
     def broadcast(self, symbol, domain):
-        # for finite volume we send variables to cells and so use number_of_cells
-        broadcasted_symbol = pybamm.NumpyBroadcast(symbol, domain, self.mesh)
-
-        # if the broadcasted symbol evaluates to a constant value, replace the
-        # symbol-Vector multiplication with a single array
-        if broadcasted_symbol.is_constant():
-            broadcasted_symbol = pybamm.Array(
-                broadcasted_symbol.evaluate(), domain=broadcasted_symbol.domain
-            )
-
-        return broadcasted_symbol
+        return pybamm.NumpyBroadcast(symbol, domain, self.mesh)
 
     def gradient(self, symbol, discretised_symbol, boundary_conditions):
         n = 0
diff --git a/tests/test_discretisations/test_discretisation.py b/tests/test_discretisations/test_discretisation.py
index b09a02430..d1d29ed08 100644
--- a/tests/test_discretisations/test_discretisation.py
+++ b/tests/test_discretisations/test_discretisation.py
@@ -512,15 +512,13 @@ class TestDiscretise(unittest.TestCase):
 
         # scalar
         broad = disc._spatial_methods[whole_cell[0]].broadcast(a, whole_cell)
-        self.assertIsInstance(broad, pybamm.Array)
         np.testing.assert_array_equal(
             broad.evaluate(), 7 * np.ones_like(combined_submesh[0].nodes)
         )
         self.assertEqual(broad.domain, whole_cell)
 
         broad_disc = disc.process_symbol(broad)
-        # type of broad will be array as broad is constant
-        self.assertIsInstance(broad_disc, pybamm.Array)
+        self.assertIsInstance(broad_disc, pybamm.NumpyBroadcast)
 
         # process Broadcast variable
         disc._y_slices = {var.id: slice(53)}
@@ -554,7 +552,6 @@ class TestDiscretise(unittest.TestCase):
 
         eqn = pybamm.Concatenation(a, b)
         eqn_disc = disc.process_symbol(eqn)
-        self.assertIsInstance(eqn_disc, pybamm.Vector)
         expected_vector = np.concatenate(
             [
                 5 * np.ones_like(mesh["negative electrode"][0].nodes),
diff --git a/tests/test_expression_tree/test_broadcasts.py b/tests/test_expression_tree/test_broadcasts.py
index a83047b87..327f3f9bd 100644
--- a/tests/test_expression_tree/test_broadcasts.py
+++ b/tests/test_expression_tree/test_broadcasts.py
@@ -25,7 +25,7 @@ class TestBroadcasts(unittest.TestCase):
         broad_a = pybamm.Broadcast(1, ["negative electrode"])
         self.assertEqual(broad_a.name, "broadcast")
         self.assertIsInstance(broad_a.children[0], pybamm.Symbol)
-        self.assertEqual(broad_a.children[0].name, str(1))
+        self.assertEqual(broad_a.children[0].name, str(1.0))
         self.assertEqual(broad_a.domain, ["negative electrode"])
 
         b = pybamm.Symbol("b", domain=["negative electrode"])
diff --git a/tests/test_expression_tree/test_concatenations.py b/tests/test_expression_tree/test_concatenations.py
index 2bd9b17f7..b5cb7a1a1 100644
--- a/tests/test_expression_tree/test_concatenations.py
+++ b/tests/test_expression_tree/test_concatenations.py
@@ -102,6 +102,9 @@ class TestConcatenations(unittest.TestCase):
                 [np.full(mesh[a_dom[0]][0].npts, 2), np.full(mesh[b_dom[0]][0].npts, 1)]
             ),
         )
+        # test size and shape
+        self.assertEqual(conc.size, mesh[a_dom[0]][0].npts + mesh[b_dom[0]][0].npts)
+        self.assertEqual(conc.shape, (mesh[a_dom[0]][0].npts + mesh[b_dom[0]][0].npts,))
 
         # check the reordering in case a child vector has to be split up
         a_dom = ["separator"]
@@ -125,6 +128,15 @@ class TestConcatenations(unittest.TestCase):
                 ]
             ),
         )
+        # test size and shape
+        self.assertEqual(
+            conc.size,
+            mesh[b_dom[0]][0].npts + mesh[a_dom[0]][0].npts + mesh[b_dom[1]][0].npts,
+        )
+        self.assertEqual(
+            conc.shape,
+            (mesh[b_dom[0]][0].npts + mesh[a_dom[0]][0].npts + mesh[b_dom[1]][0].npts,),
+        )
 
     def test_concatenation_orphans(self):
         a = pybamm.Variable("a")
diff --git a/tests/test_expression_tree/test_symbol.py b/tests/test_expression_tree/test_symbol.py
index c47dad489..972cb09a4 100644
--- a/tests/test_expression_tree/test_symbol.py
+++ b/tests/test_expression_tree/test_symbol.py
@@ -53,7 +53,7 @@ class TestSymbol(unittest.TestCase):
         self.assertEqual((A @ a).simplify().evaluate(), 0)
 
         # test when other node is a parameter
-        c = pybamm.Parameter('c')
+        c = pybamm.Parameter("c")
         self.assertIsInstance((a + c).simplify(), pybamm.Parameter)
         self.assertIsInstance((c + a).simplify(), pybamm.Parameter)
         self.assertIsInstance((c + b).simplify(), pybamm.Addition)
diff --git a/tests/test_meshes/test_meshes.py b/tests/test_meshes/test_meshes.py
index 349a84377..ed2188983 100644
--- a/tests/test_meshes/test_meshes.py
+++ b/tests/test_meshes/test_meshes.py
@@ -87,8 +87,12 @@ class TestMesh(unittest.TestCase):
         # create mesh
         mesh = mesh_type(geometry, submesh_types, submesh_pts)
         for domain in mesh:
-            self.assertEqual(mesh[domain][0].npts, submesh_pts[domain]["x"])
-            self.assertEqual(len(mesh[domain][0].edges) - 1, submesh_pts[domain]["x"])
+            # ignore ghost cells for this test
+            if "ghost" not in domain:
+                self.assertEqual(mesh[domain][0].npts, submesh_pts[domain]["x"])
+                self.assertEqual(
+                    len(mesh[domain][0].edges) - 1, submesh_pts[domain]["x"]
+                )
 
     def test_combine_submeshes(self):
         param = pybamm.ParameterValues(
diff --git a/tests/test_models/standard_model_tests.py b/tests/test_models/standard_model_tests.py
index 6c5d41005..bb7d00d20 100644
--- a/tests/test_models/standard_model_tests.py
+++ b/tests/test_models/standard_model_tests.py
@@ -53,3 +53,24 @@ class StandardModelTest(object):
         self.test_processing_parameters(param)
         self.test_processing_disc(disc)
         self.test_solving(solver, t_eval)
+
+    def test_update_parameters(self, param):
+        # check if geometry has changed, throw error if so (need to re-discretise)
+        if any(
+            [
+                length in param.keys() and param[length] != self.param[length]
+                for length in [
+                    "Negative electrode width",
+                    "Separator width",
+                    "Positive electrode width",
+                ]
+            ]
+        ):
+            raise ValueError(
+                "geometry has changed, the orginal model needs to be re-discretised"
+            )
+        # otherwise update self.param and change the parameters in the discretised model
+        self.param = param
+        param.process_discretised_model(self.model, self.disc)
+        # Model should still be well-posed after processing
+        self.model.check_well_posedness()
diff --git a/tests/test_parameters/test_parameter_values.py b/tests/test_parameters/test_parameter_values.py
index 4a3320ce8..604278057 100644
--- a/tests/test_parameters/test_parameter_values.py
+++ b/tests/test_parameters/test_parameter_values.py
@@ -234,7 +234,9 @@ class TestParameterValues(unittest.TestCase):
         self.assertTrue(
             isinstance(model.variables["grad_var1"].children[0], pybamm.Variable)
         )
-        self.assertEqual(model.variables["d_var1"].id, (42 * var1).id)
+        self.assertEqual(
+            model.variables["d_var1"].id, (pybamm.Scalar(42, name="d") * var1).id
+        )
         self.assertIsInstance(model.variables["d_var1"].children[0], pybamm.Scalar)
         self.assertTrue(
             isinstance(model.variables["d_var1"].children[1], pybamm.Variable)
diff --git a/tests/test_parameters/test_standard_parameters_lead_acid.py b/tests/test_parameters/test_standard_parameters_lead_acid.py
index 9afaaa151..55802a948 100644
--- a/tests/test_parameters/test_standard_parameters_lead_acid.py
+++ b/tests/test_parameters/test_standard_parameters_lead_acid.py
@@ -69,8 +69,6 @@ class TestStandardParametersLeadAcid(unittest.TestCase):
         processed_s = disc.process_symbol(parameter_values.process_symbol(s))
 
         # test output
-        self.assertIsInstance(processed_s, pybamm.Vector)
-
         combined_submeshes = disc.mesh.combine_submeshes(
             "negative electrode", "separator", "positive electrode"
         )
diff --git a/tests/test_parameters/test_update_parameters.py b/tests/test_parameters/test_update_parameters.py
new file mode 100644
index 000000000..374c40795
--- /dev/null
+++ b/tests/test_parameters/test_update_parameters.py
@@ -0,0 +1,112 @@
+#
+# Tests for updating parameter values
+#
+import pybamm
+
+import unittest
+import numpy as np
+import tests
+
+
+class TestUpdateParameters(unittest.TestCase):
+    def test_update_parameters_eqn(self):
+        a = pybamm.Scalar(1)
+        b = pybamm.Scalar(2, name="test parameter")
+        c = pybamm.Scalar(3)
+        eqn = a + b * c
+        self.assertEqual(eqn.evaluate(), 7)
+
+        parameter_values = pybamm.ParameterValues({"test parameter": 3})
+        eqn_changed = parameter_values.process_symbol(eqn)
+        self.assertEqual(eqn_changed.evaluate(), 10)
+
+    def test_set_and_update_parameters(self):
+        a = pybamm.Scalar(1)
+        b = pybamm.Parameter(name="test parameter")
+        c = pybamm.Scalar(3)
+        eqn = a + b * c
+
+        parameter_values = pybamm.ParameterValues({"test parameter": 2})
+        eqn_changed = parameter_values.process_symbol(eqn)
+        self.assertEqual(eqn_changed.evaluate(), 7)
+
+        parameter_values = pybamm.ParameterValues({"test parameter": 3})
+        eqn_changed_again = parameter_values.process_symbol(eqn_changed)
+        self.assertEqual(eqn_changed_again.evaluate(), 10)
+
+    def test_update_model(self):
+        # standard model
+        model1 = pybamm.ReactionDiffusionModel()
+        modeltest1 = tests.StandardModelTest(model1)
+        modeltest1.test_all()
+        Y1 = modeltest1.solver.y
+
+        # double initial conditions
+        model2 = pybamm.ReactionDiffusionModel()
+        # process and solve the model a first time
+        modeltest2 = tests.StandardModelTest(model2)
+        modeltest2.test_all()
+        # process and solve with updated parameter values
+        parameter_values_update = pybamm.ParameterValues(
+            base_parameters=model2.default_parameter_values,
+            optional_parameters={"Typical current density": 2},
+        )
+        modeltest2.test_update_parameters(parameter_values_update)
+        modeltest2.test_solving()
+        Y2 = modeltest2.solver.y
+
+        # results should be different
+        self.assertNotEqual(np.linalg.norm(Y1 - Y2), 0)
+
+    def test_update_geometry(self):
+        # standard model
+        model1 = pybamm.ReactionDiffusionModel()
+        modeltest1 = tests.StandardModelTest(model1)
+        modeltest1.test_all()
+        T1, Y1 = modeltest1.solver.t, modeltest1.solver.y
+
+        # trying to update the geometry fails
+        parameter_values_update = pybamm.ParameterValues(
+            base_parameters=model1.default_parameter_values,
+            optional_parameters={
+                "Negative electrode width": 0.0002,
+                "Separator width": 0.0003,
+                "Positive electrode width": 0.0004,
+            },
+        )
+        with self.assertRaisesRegex(ValueError, "geometry has changed"):
+            modeltest1.test_update_parameters(parameter_values_update)
+
+        # instead we need to make a new model and re-discretise
+        model2 = pybamm.ReactionDiffusionModel()
+        parameter_values_update = pybamm.ParameterValues(
+            base_parameters=model2.default_parameter_values,
+            optional_parameters={
+                "Negative electrode width": 0.0002,
+                "Separator width": 0.0003,
+                "Positive electrode width": 0.0004,
+            },
+        )
+        modeltest2 = tests.StandardModelTest(model2)
+        modeltest2.test_all(param=parameter_values_update)
+        T2, Y2 = modeltest2.solver.t, modeltest2.solver.y
+        # results should be different
+        for idx in range(len(T1)):
+            j1 = modeltest1.model.variables["Interfacial current density"].evaluate(
+                T1[idx], Y1[:, idx]
+            )
+            j2 = modeltest2.model.variables["Interfacial current density"].evaluate(
+                T2[idx], Y2[:, idx]
+            )
+            self.assertNotEqual(np.linalg.norm(j1 - j2), 0)
+        self.assertNotEqual(np.linalg.norm(Y1 - Y2), 0)
+
+
+if __name__ == "__main__":
+    import sys
+
+    print("Add -v for more debug output")
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()