Repository: tinosulzer/PyBaMM

Issue:
refactor reaction diffusion model
See #16. 

refactor the `ReactionDiffusionModel` class to inherit from `BaseModel` (see #19), implementing a reaction diffusion model


Discussion:

Patch:
diff --git a/docs/source/models/electrolyte_current.rst b/docs/source/models/electrolyte_current.rst
deleted file mode 100644
index c3bdf616..00000000
--- a/docs/source/models/electrolyte_current.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-Electrolyte Current Model
-=========================
-
-.. autoclass:: pybamm.ElectrolyteCurrentModel
-    :members:
-
diff --git a/docs/source/models/index.rst b/docs/source/models/index.rst
index e8088c1f..a439e9d2 100644
--- a/docs/source/models/index.rst
+++ b/docs/source/models/index.rst
@@ -3,12 +3,11 @@ Models
 
 .. toctree::
 
-  electrolyte_current
   reaction_diffusion
-
+  submodels/index
+  
 Core
 ----
 
 .. autoclass:: pybamm.BaseModel
     :members:
-
diff --git a/docs/source/models/submodels/electrolyte.rst b/docs/source/models/submodels/electrolyte.rst
new file mode 100644
index 00000000..a1519b2a
--- /dev/null
+++ b/docs/source/models/submodels/electrolyte.rst
@@ -0,0 +1,5 @@
+Electrolyte Submodels
+=====================
+
+.. autoclass:: pybamm.electrolyte.StefanMaxwellDiffusion
+    :members:
diff --git a/docs/source/models/submodels/index.rst b/docs/source/models/submodels/index.rst
new file mode 100644
index 00000000..fcde4d5d
--- /dev/null
+++ b/docs/source/models/submodels/index.rst
@@ -0,0 +1,7 @@
+Models
+======
+
+.. toctree::
+
+  electrolyte
+  interface
diff --git a/docs/source/models/submodels/interface.rst b/docs/source/models/submodels/interface.rst
new file mode 100644
index 00000000..fbf3ac53
--- /dev/null
+++ b/docs/source/models/submodels/interface.rst
@@ -0,0 +1,4 @@
+Interface Submodels
+===================
+
+.. automethod:: pybamm.interface.homogeneous_reaction
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 4ed40545..73d18fa6 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -89,7 +89,6 @@ from .expression_tree.exceptions import DomainError, ModelError
 #
 from .models.base_model import BaseModel
 from .models.reaction_diffusion import ReactionDiffusionModel
-from .models.electrolyte_current import ElectrolyteCurrentModel
 
 #
 # Submodel classes
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 7d2a33ff..cca3af14 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -23,6 +23,10 @@ class BaseModel(object):
     boundary_conditions: dict
         A dictionary that maps expressions (variables) to expressions that represent
         the boundary conditions
+    variables: dict
+        A dictionary that maps strings to expressions that represent
+        the useful variables
+
     """
 
     def __init__(self):
@@ -34,6 +38,14 @@ class BaseModel(object):
         self._concatenated_rhs = None
         self._concatenated_initial_conditions = None
 
+        # Default parameter values, discretisation and solver
+        self.default_parameter_values = pybamm.ParameterValues(
+            "input/parameters/lithium-ion/parameters/LCO.csv"
+        )
+        mesh = pybamm.FiniteVolumeMacroMesh(self.default_parameter_values, 2)
+        self.default_discretisation = pybamm.FiniteVolumeDiscretisation(mesh)
+        self.default_solver = pybamm.ScipySolver(method="RK45")
+
     @property
     def rhs(self):
         return self._rhs
@@ -121,6 +133,28 @@ class BaseModel(object):
     def __getitem__(self, key):
         return self.rhs[key]
 
+    def update(self, *submodels):
+        """
+        Update model to add new physics from submodels
+
+        Parameters
+        ----------
+        submodel : iterable of submodels (subclasses of :class:`pybamm.BaseModel`)
+            The submodels from which to create new model
+        """
+        for submodel in submodels:
+            # check for duplicates in keys
+            vars = [var.id for var in submodel.rhs.keys()] + [
+                var.id for var in self.rhs.keys()
+            ]
+            assert len(vars) == len(set(vars)), pybamm.ModelError("duplicate variables")
+
+            # update dicts
+            self._rhs.update(submodel.rhs)
+            self._initial_conditions.update(submodel.initial_conditions)
+            self._boundary_conditions.update(submodel.boundary_conditions)
+            self._variables.update(submodel.variables)
+
     def check_well_posedness(self):
         """
         Check that the model is well-posed by executing the following tests:
diff --git a/pybamm/models/electrolyte_current.py b/pybamm/models/electrolyte_current.py
deleted file mode 100644
index 582133d4..00000000
--- a/pybamm/models/electrolyte_current.py
+++ /dev/null
@@ -1,95 +0,0 @@
-#
-# Electrolyte current model
-#
-from __future__ import absolute_import, division
-from __future__ import print_function, unicode_literals
-import pybamm
-
-import numpy as np
-
-
-class ElectrolyteCurrentModel(pybamm.BaseModel):
-    """Electrolyte current model.
-
-    Parameters
-    ----------
-    tests : dict, optional
-        A dictionary for testing the convergence of the numerical solution:
-            * {} (default): We are not running in test mode, use built-ins.
-            * {'inits': dict of initial conditions,
-               'bcs': dict of boundary conditions,
-               'sources': dict of source terms
-               }: To be used for testing convergence to an exact solution.
-    """
-
-    def __init__(self, tests={}):
-        super()
-        self.name = "Electrolyte Current"
-        # Assign tests as an attribute
-        if tests:
-            assert set(tests.keys()) == {
-                "inits",
-                "bcs",
-                "sources",
-            }, "tests.keys() must include, 'inits', 'bcs' and 'sources'"
-        self.tests = tests
-
-        # Set variables
-        self.variables = [("en", "xcn"), ("ep", "xcp")]
-
-        # Initialise the class(es) that will be called upon for equations
-        self.electrolyte = pybamm.Electrolyte()
-        self.interface = pybamm.Interface()
-
-    def set_simulation(self, param, operators, mesh):
-        """
-        Assign simulation-specific objects as attributes.
-
-        Parameters
-        ----------
-        param : :class:`pybamm.ParameterValues` instance
-            The parameters of the simulation
-        operators : :class:`pybamm.Operators` instance
-            The spatial operators.
-        mesh : :class:`pybamm.Mesh` instance
-            The spatial and temporal discretisation.
-        """
-        self.param = param
-        self.operators = operators
-        self.mesh = mesh
-
-        # Set simulation for the components
-        self.electrolyte.set_simulation(param, operators, mesh)
-        self.interface.set_simulation(param, mesh)
-
-    def initial_conditions(self):
-        """See :meth:`pybamm.BaseModel.initial_conditions`"""
-        if not self.tests:
-            electrolyte_inits = self.electrolyte.initial_conditions()
-            y0 = np.concatenate([electrolyte_inits["en"], electrolyte_inits["ep"]])
-            return y0
-        else:
-            return self.tests["inits"]
-
-    def pdes_rhs(self, vars):
-        """See :meth:`pybamm.BaseModel.pdes_rhs`"""
-        cn = np.ones_like(self.mesh.xcn)
-        cp = np.ones_like(self.mesh.xcp)
-        if not self.tests:
-            current_bcs_n = self.electrolyte.bcs_current("xcn", vars.t)
-            current_bcs_p = self.electrolyte.bcs_current("xcp", vars.t)
-            jn = self.interface.butler_volmer("xcn", cn, vars.en)
-            jp = self.interface.butler_volmer("xcp", cp, vars.ep)
-        else:
-            current_bcs_n = self.tests["bcs"](vars.t)["current n"]
-            current_bcs_p = self.tests["bcs"](vars.t)["current p"]
-            jn = self.tests["sources"](vars.t)["current n"]
-            jp = self.tests["sources"](vars.t)["current p"]
-        dendt = self.electrolyte.current_conservation(
-            "xcn", cn, vars.en, jn, current_bcs_n
-        )
-        depdt = self.electrolyte.current_conservation(
-            "xcp", cp, vars.ep, jp, current_bcs_p
-        )
-
-        return np.concatenate([dendt, depdt])
diff --git a/pybamm/models/reaction_diffusion.py b/pybamm/models/reaction_diffusion.py
index 940b310b..e465032d 100644
--- a/pybamm/models/reaction_diffusion.py
+++ b/pybamm/models/reaction_diffusion.py
@@ -9,45 +9,33 @@ import pybamm
 class ReactionDiffusionModel(pybamm.BaseModel):
     """Reaction-diffusion model.
 
-    Parameters
+    Attributes
     ----------
-    tests : dict, optional
-        A dictionary for testing the convergence of the numerical solution:
-            * {} (default): We are not running in test mode, use built-ins.
-            * {'inits': dict of initial conditions,
-               'bcs': dict of boundary conditions,
-               'sources': dict of source terms
-               }: To be used for testing convergence to an exact solution.
+
+    rhs: dict
+        A dictionary that maps expressions (variables) to expressions that represent
+        the rhs
+    initial_conditions: dict
+        A dictionary that maps expressions (variables) to expressions that represent
+        the initial conditions
+    boundary_conditions: dict
+        A dictionary that maps expressions (variables) to expressions that represent
+        the boundary conditions
+    variables: dict
+        A dictionary that maps strings to expressions that represent
+        the useful variables
+
     """
 
-    def __init__(self, tests={}):
-        super().__init__(tests)
-        self.name = "Reaction Diffusion"
-
-    @property
-    def submodels(self):
-        pdes = {
-            "c": pybamm.electrolyte.StefanMaxwellDiffusion(
-                self.param.electrolyte, self.operators.x, self.mesh.x, self.tests
-            )
-        }
-
-        reactions = {
-            "neg": pybamm.interface.HomogeneousReaction(
-                self.param.neg_reactions, self.mesh.xn
-            ),
-            "pos": pybamm.interface.HomogeneousReaction(
-                self.param.pos_reactions, self.mesh.xp
-            ),
-        }
-
-        # {
-        #     "neg": pybamm.interface.ButlerVolmer(
-        #         self.param.neg_reactions, self.mesh.xn
-        #     ),
-        #     "pos": pybamm.interface.ButlerVolmer(
-        #         self.param.pos_reactions, self.mesh.xp
-        #     ),
-        # }
-
-        return {"pdes": pdes, "reactions": reactions}
+    def __init__(self):
+        super().__init__()
+        # Load reaction flux from submodels
+        G = pybamm.interface.homogeneous_reaction()
+        # Load diffusion model from submodels
+        diffusion_model = pybamm.electrolyte.StefanMaxwellDiffusion(G)
+
+        # Create own model from diffusion model
+        self.update(diffusion_model)
+
+        # Overwrite default solver for faster solution
+        self.default_solver = pybamm.ScipySolver(method="BDF")
diff --git a/pybamm/models/submodels/electrolyte.py b/pybamm/models/submodels/electrolyte.py
index 39609690..b9595bc9 100644
--- a/pybamm/models/submodels/electrolyte.py
+++ b/pybamm/models/submodels/electrolyte.py
@@ -10,6 +10,12 @@ class StefanMaxwellDiffusion(pybamm.BaseModel):
     """A class that generates the expression tree for Stefan-Maxwell Diffusion in the
     electrolyte.
 
+    Parameters
+    ----------
+    G : :class:`pybamm.Symbol`
+        An expression tree that represents the concentration flux at the
+        electrode-electrolyte interface
+
     Attributes
     ----------
 
diff --git a/pybamm/models/submodels/interface.py b/pybamm/models/submodels/interface.py
index 0579ded8..bb709551 100644
--- a/pybamm/models/submodels/interface.py
+++ b/pybamm/models/submodels/interface.py
@@ -7,6 +7,10 @@ import pybamm
 
 
 def homogeneous_reaction():
+    """
+    Homogeneous reaction at the electrode-electrolyte interface
+    """
+
     current_neg = (
         pybamm.Scalar(1, domain=["negative electrode"]) / pybamm.standard_parameters.ln
     )


Test patch:
diff --git a/tests/test_models/standard_model_tests.py b/tests/test_models/standard_model_tests.py
index 190ff701..37b985ca 100644
--- a/tests/test_models/standard_model_tests.py
+++ b/tests/test_models/standard_model_tests.py
@@ -1,40 +1,43 @@
 #
 # Standard basic tests for any model
 #
-import pybamm
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
 
 
 class StandardModelTest(object):
     def __init__(self, model):
         self.model = model
-        self.param = None
-        self.mesh = None
-
-    def test_processing_parameters(self, param_str="LCO"):
-        if param_str == "LCO":
-            self.param = pybamm.ParameterValues(
-                "input/parameters/lithium-ion/parameters/LCO.csv"
-            )
+        # Set defaults
+        self.param = model.default_parameter_values
+        self.disc = model.default_discretisation
+        self.solver = model.default_solver
 
+    def test_processing_parameters(self, param=None):
+        # Overwrite parameters if given
+        if param is not None:
+            self.param = param
         self.param.process_model(self.model)
         # Model should still be well-posed after processing
         self.model.check_well_posedness()
 
-    def test_processing_disc(self, disc_str="Finite Volume"):
-        if disc_str == "Finite Volume":
-            self.mesh = pybamm.FiniteVolumeMacroMesh(self.param, 2)
-            disc = pybamm.FiniteVolumeDiscretisation(self.mesh)
-        disc.process_model(self.model)
+    def test_processing_disc(self, disc=None):
+        # Overwrite discretisation if given
+        if disc is not None:
+            self.disc = disc
+        self.disc.process_model(self.model)
         # Model should still be well-posed after processing
         self.model.check_well_posedness()
 
-    def test_solving(self):
-        solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
-        t_eval = self.mesh["time"]
-        solver.solve(self.model, t_eval)
+    def test_solving(self, solver=None):
+        # Overwrite solver if given
+        if solver is not None:
+            self.solver = solver
+        t_eval = self.disc.mesh["time"]
+        self.solver.solve(self.model, t_eval)
 
-    def test_all(self):
+    def test_all(self, param=None, disc=None, solver=None):
         self.model.check_well_posedness()
-        self.test_processing_parameters()
-        self.test_processing_disc()
-        self.test_solving()
+        self.test_processing_parameters(param)
+        self.test_processing_disc(disc)
+        self.test_solving(solver)
diff --git a/tests/test_models/test_base_model.py b/tests/test_models/test_base_model.py
index 0071f336..72faae8a 100644
--- a/tests/test_models/test_base_model.py
+++ b/tests/test_models/test_base_model.py
@@ -1,6 +1,8 @@
 #
 # Tests for the base model class
 #
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
 import pybamm
 
 import unittest
@@ -85,6 +87,67 @@ class TestBaseModel(unittest.TestCase):
         self.assertEqual(model[key], rhs[key])
         self.assertEqual(model[key], model.rhs[key])
 
+    def test_update(self):
+        # model
+        model = pybamm.BaseModel()
+        c = pybamm.Variable("c", domain=["whole cell"])
+        rhs = {c: 5 * pybamm.div(pybamm.grad(c)) - 1}
+        initial_conditions = {c: 1}
+        boundary_conditions = {c: {"left": 0, "right": 0}}
+        variables = {"c": c}
+        model.rhs = rhs
+        model.initial_conditions = initial_conditions
+        model.boundary_conditions = boundary_conditions
+        model.variables = variables
+
+        # update with submodel
+        submodel = pybamm.BaseModel()
+        d = pybamm.Variable("d", domain=["whole cell"])
+        submodel.rhs = {
+            d: 5 * pybamm.div(pybamm.grad(c)) + pybamm.div(pybamm.grad(d)) - 1
+        }
+        submodel.initial_conditions = {d: 3}
+        submodel.boundary_conditions = {d: {"left": 4, "right": 7}}
+        submodel.variables = {"d": d}
+        model.update(submodel)
+
+        # check
+        self.assertEqual(model.rhs[d], submodel.rhs[d])
+        self.assertEqual(model.initial_conditions[d], submodel.initial_conditions[d])
+        self.assertEqual(model.boundary_conditions[d], submodel.boundary_conditions[d])
+        self.assertEqual(model.variables["d"], submodel.variables["d"])
+        self.assertEqual(model.rhs[c], rhs[c])
+        self.assertEqual(model.initial_conditions[c], initial_conditions[c])
+        self.assertEqual(model.boundary_conditions[c], boundary_conditions[c])
+        self.assertEqual(model.variables["c"], variables["c"])
+
+        # update with conflicting submodel
+        submodel2 = pybamm.BaseModel()
+        submodel2.rhs = {d: pybamm.div(pybamm.grad(d)) - 1}
+        with self.assertRaises(AssertionError) as error:
+            model.update(submodel2)
+        self.assertIsInstance(error.exception.args[0], pybamm.ModelError)
+
+        # update with multiple submodels
+        submodel1 = submodel  # copy submodel from previous test
+        submodel2 = pybamm.BaseModel()
+        e = pybamm.Variable("e", domain=["whole cell"])
+        submodel2.rhs = {
+            e: 5 * pybamm.div(pybamm.grad(d)) + pybamm.div(pybamm.grad(e)) - 1
+        }
+        submodel2.initial_conditions = {e: 3}
+        submodel2.boundary_conditions = {e: {"left": 4, "right": 7}}
+
+        model = pybamm.BaseModel()
+        model.update(submodel1, submodel2)
+
+        self.assertEqual(model.rhs[d], submodel1.rhs[d])
+        self.assertEqual(model.initial_conditions[d], submodel1.initial_conditions[d])
+        self.assertEqual(model.boundary_conditions[d], submodel1.boundary_conditions[d])
+        self.assertEqual(model.rhs[e], submodel2.rhs[e])
+        self.assertEqual(model.initial_conditions[e], submodel2.initial_conditions[e])
+        self.assertEqual(model.boundary_conditions[e], submodel2.boundary_conditions[e])
+
     def test_check_well_posedness(self):
         # Well-posed model - Dirichlet
         model = pybamm.BaseModel()
diff --git a/tests/test_models/test_electrolyte.py b/tests/test_models/test_electrolyte.py
index 87e2fbfa..02a50604 100644
--- a/tests/test_models/test_electrolyte.py
+++ b/tests/test_models/test_electrolyte.py
@@ -1,5 +1,5 @@
 #
-# Tests for the Binary Operator classes
+# Tests for the electrolyte submodels
 #
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
diff --git a/tests/test_models/test_reaction_diffusion.py b/tests/test_models/test_reaction_diffusion.py
new file mode 100644
index 00000000..3bee36a7
--- /dev/null
+++ b/tests/test_models/test_reaction_diffusion.py
@@ -0,0 +1,26 @@
+#
+# Tests for the Reaction diffusion model
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+import pybamm
+import tests
+
+import unittest
+
+
+class TestReactionDiffusionModel(unittest.TestCase):
+    def test_basic_processing(self):
+        model = pybamm.ReactionDiffusionModel()
+
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()