Repository: tinosulzer/PyBaMM

Issue:
Implement Single Particle Model (SPM)
**Summary**
Implement the SPM. We require: 
-  a `StandardParticle` class which implements spherical diffusion on a spherically symmetric particle taking a boundary condition `G` for the flux on the particle surface. 
- an `SPM` class which calls `homogeneous_reactions` to create `G` and then calls `StandardParticle` twice to create a negative electrode particle and a positive electrode particle. 
- The SPM class should also contain an expression for the voltage in terms of the other fundamental variables (which are just c_n and c_p in this case). This will be contained within the `variables` dictionary of the model


Discussion:

Patch:
diff --git a/input/parameters/lithium-ion/parameters/LCO.csv b/input/parameters/lithium-ion/parameters/LCO.csv
index efcf8c3c..141cb2ed 100644
--- a/input/parameters/lithium-ion/parameters/LCO.csv
+++ b/input/parameters/lithium-ion/parameters/LCO.csv
@@ -7,25 +7,26 @@ F,96487,C.mol-1,Faraday’s constant,-,
 T_ref,298.15,K,Reference Temperature,25C,
 ,,,,,
 # Geometric,,,,,
-Ln,1.00E-04,m,Half-width of negative electrode,Scott Moura,
-Ls,2.50E-05,m,Width of separator,Scott Moura,
-Lp,1.00E-04,m,Half-width of positive electrode,Scott Moura,
-H,1.37E-01,m,Electrode height,KOKAM SLPB78205130H,Not needed for 1D
-W,2.07E-01,m,Electrode depth,KOKAM SLPB78205130H,Not needed for 1D
-L_ccn,2.50E-05,m,Negative current collector thickness,Scott Moura,
-L_ccp,2.50E-05,m,Positive current collector thickness,Scott Moura,
+Ln,0.0001,m,Half-width of negative electrode,Scott Moura,
+Ls,2.5E-05,m,Width of separator,Scott Moura,
+Lp,0.0001,m,Half-width of positive electrode,Scott Moura,
+H,0.137,m,Electrode height,KOKAM SLPB78205130H,Not needed for 1D
+W,0.207,m,Electrode depth,KOKAM SLPB78205130H,Not needed for 1D
+L_ccn,2.5E-05,m,Negative current collector thickness,Scott Moura,
+L_ccp,2.5E-05,m,Positive current collector thickness,Scott Moura,
 ,,,,,
 # Electrical,,,,,
 voltage_low_cut,3.105,V,Lower voltage cut-off,-,
 voltage_high_cut,4.7,V,Upper voltage cut-off,-,
 I_typ,24,A/m^2,Typical current density,-,
+Phi_typ,1,V,Typical voltage drop across cell,-,
 ,,,,,
 # Microstructure,,,,,
 epsilon_n,0.3,-,Electrolyte volume fraction in negative electrode region,Scott Moura,
 epsilon_s,1,-,Electrolyte volume fraction in separator region,Scott Moura,
 epsilon_p,0.3,,Electrolyte volume fraction in positive electrode region,Scott Moura,
-R_n,1.00E-05,m,Negative electrode particle radius,Scott Moura,
-R_p,1.00E-05,m,Positive electrode particle radius,Scott Moura,
+R_n,1E-05,m,Negative electrode particle radius,Scott Moura,
+R_p,1E-05,m,Positive electrode particle radius,Scott Moura,
 a_n,90000,m-1,Negative electrode surface area density, Scott Moura,
 a_p,90000,m-1,Positive electrode surface area density, Scott Moura,
 b,1.5,-,Bruggeman coefficient,Scott Moura,
@@ -40,15 +41,15 @@ sigma_n,100,S.m-1,Conductivity of graphite,Scott Moura,
 sigma_p,10,S.m-1,Conductivity of Lithium Cobalt Oxide,Scott Moura,
 cn_max,24983.2619938437,mol.m-3,Maximum lithium concentration in negative electrode ,Scott Moura,
 cp_max,51217.9257309275,mol.m-3,Maximum lithium concentration in positive electrode ,Scott Moura,
-D_n,3.90E-14,m2/s,Negative electrode diffusion coefficient,Scott Moura,
-D_p,1.00E-13,m2/s,Positive electrode diffusion coefficient,Scott Moura,
+Dn_typ,3.9E-14,m2/s,Negative electrode diffusion coefficient,Scott Moura,
+Dp_typ,1E-13,m2/s,Positive electrode diffusion coefficient,Scott Moura,
 ,,,,,
 # Interfacial reactions,,,,,
-m_n,2.00E-05,Am-2(m3/mol)1.5,Reaction rate in negative electrode,Scott Moura,Be careful how we implement BV
-m_p,6.00E-07,Am-2(m3/mol)1.5,Reaction rate in positive electrode,Scott Moura,Be careful how we implement BV
+m_n_dim,2E-05,Am-2(m3/mol)1.5,Reaction rate in negative electrode,Scott Moura,Be careful how we implement BV
+m_p_dim,6E-07,Am-2(m3/mol)1.5,Reaction rate in positive electrode,Scott Moura,Be careful how we implement BV
 alpha,0.5,,Charge transfer coefficient,Scott Moura,
 ,,,,,
 # Initial Conditions,,,,,
-ce0, 1000,mol.m-3, Initial lithium ion concentration in electrolyte,,
-cn0, 1000,mol.m-3, Initial lithium concentration in positive electrode,,
-cp0, 1000,mol.m-3, Initial lithium concentration in negative electrode,,
\ No newline at end of file
+ce0,1000,mol.m-3, Initial lithium ion concentration in electrolyte,,
+cn0,1000,mol.m-3, Initial lithium concentration in positive electrode,,
+cp0,1000,mol.m-3, Initial lithium concentration in negative electrode,,
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 98c5b202..81627d99 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -99,6 +99,7 @@ from .models.base_model import BaseModel
 from .models.reaction_diffusion import ReactionDiffusionModel
 from .models.simple_ode_model import SimpleODEModel
 from .models import lead_acid
+from .models import li_ion
 
 #
 # Submodel classes
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 18cff23e..8e1d0694 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -215,7 +215,7 @@ class Discretisation(object):
         elif isinstance(symbol, pybamm.SurfaceValue):
             child = symbol.children[0]
             discretised_child = self.process_symbol(child)
-            return self._spatial_methods[symbol.domain[0]].surface_value(
+            return self._spatial_methods[child.domain[0]].surface_value(
                 discretised_child
             )
 
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index 8a4940ed..db5acced 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -67,8 +67,8 @@ class NumpyModelConcatenation(pybamm.Symbol):
     def __init__(self, *children):
         children = list(children)
         for i, child in enumerate(children):
-            if isinstance(child, pybamm.Scalar):
-                children[i] = pybamm.Vector(np.array([child.value]))
+            if child.evaluates_to_number():
+                children[i] = pybamm.Vector(np.array([child.evaluate()]))
         super().__init__("model concatenation", children, domain=[])
 
     def evaluate(self, t=None, y=None):
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 78d4ef01..5b850c5e 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -339,6 +339,10 @@ class Symbol(anytree.NodeMixin):
             # return false if NotImplementedError is raised
             # (there is a e.g. Parameter, Variable, ... in the tree)
             return False
+        except TypeError:
+            # return false if TypeError is raised
+            # (there is a e.g. StateVector in the tree)
+            return False
 
     def has_spatial_derivatives(self):
         """Returns True if equation has spatial derivatives (grad or div)."""
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index abe2babb..2a4d87b8 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -219,6 +219,10 @@ class SurfaceValue(SpatialOperator):
     def __init__(self, child):
         super().__init__("surf", child)
 
+        # Domain of SurfaceValue must be ([]) so that expressions can be formed
+        # of surface values of variables in different domains
+        self.domain = []
+
 
 #
 # Methods to call Gradient and Divergence
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index b81d42ac..b09a5337 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -188,26 +188,17 @@ class BaseModel(object):
             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")
 
-            # check and then update dicts
-            self.check_and_combine_dict(self._rhs, submodel.rhs)
-            self.check_and_combine_dict(self._algebraic, submodel.algebraic)
-            self.check_and_combine_dict(
-                self._initial_conditions, submodel.initial_conditions
-            )
-            self.check_and_combine_dict(
-                self._boundary_conditions, submodel.boundary_conditions
-            )
-            self._variables.update(submodel.variables)  # keys are strings so no check
-
-    def check_and_combine_dict(self, dict1, dict2):
-        # check that the key ids are distinct
-        ids1 = set(x.id for x in dict1.keys())
-        ids2 = set(x.id for x in dict2.keys())
-        assert len(ids1.intersection(ids2)) == 0, pybamm.ModelError(
-            "Submodel incompatible: duplicate variables"
-        )
-        dict1.update(dict2)
+            # 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):
         """
diff --git a/pybamm/models/li_ion/__init__.py b/pybamm/models/li_ion/__init__.py
new file mode 100644
index 00000000..0bff98b8
--- /dev/null
+++ b/pybamm/models/li_ion/__init__.py
@@ -0,0 +1,4 @@
+#
+# Root of the li-ion models module.
+#
+from .single_particle_model import SPM
diff --git a/pybamm/models/li_ion/single_particle_model.py b/pybamm/models/li_ion/single_particle_model.py
new file mode 100644
index 00000000..445b2cc9
--- /dev/null
+++ b/pybamm/models/li_ion/single_particle_model.py
@@ -0,0 +1,116 @@
+#
+# Li-ion single particle model
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+import pybamm
+
+
+class SPM(pybamm.BaseModel):
+    """Single Particle Model for li-ion.
+
+    Attributes
+    ----------
+
+    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):
+        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")
+
+        # Parameters
+        ln = pybamm.standard_parameters.ln
+        lp = pybamm.standard_parameters.lp
+        gamma_n = pybamm.standard_parameters.gamma_n
+        gamma_p = pybamm.standard_parameters.gamma_p
+        beta_n = pybamm.standard_parameters.beta_n
+        beta_p = pybamm.standard_parameters.beta_p
+        C_hat_p = pybamm.standard_parameters.C_hat_p
+        m_n = pybamm.standard_parameters.m_n
+        m_p = pybamm.standard_parameters.m_p
+        D_n = pybamm.standard_parameters.D_n
+        D_p = pybamm.standard_parameters.D_p
+        U_n = pybamm.standard_parameters.U_n
+        U_p = pybamm.standard_parameters.U_p
+        Lambda = pybamm.standard_parameters.Lambda
+
+        # Initial conditions
+        cn_init = pybamm.standard_parameters.cn0
+        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)
+        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},
+        }
+
+        # Initial conditions
+        self.initial_conditions = {
+            cn: cn_init,
+            cp: cp_init,
+        }
+
+        # Variables
+        cn_surf = pybamm.surf(cn)
+        cp_surf = pybamm.surf(cp)
+        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)))
+
+        self.variables = {
+            "cn": cn,
+            "cp": cp,
+            "cn_surf": cn_surf,
+            "cp_surf": cp_surf,
+            "V": V,
+        }
diff --git a/pybamm/parameters/standard_parameters.py b/pybamm/parameters/standard_parameters.py
index d7117f46..8ef94c66 100644
--- a/pybamm/parameters/standard_parameters.py
+++ b/pybamm/parameters/standard_parameters.py
@@ -44,6 +44,8 @@ sigma_n, sigma_p
     Electrical conductivities of the negative and positive electrode
 cn_max, cp_max
     Maximum lithium concentration in the negative and positive electrode
+Dn_typ, Dp_typ
+    Typical diffusivitites in the solid electrode material
 
 Electrochemical Reactions
 --------------------------
@@ -56,6 +58,8 @@ voltage_low_cut, voltage_high_cut
     Low and high voltage cut-offs
 I_typ
     Typical current density
+Phi_typ
+    Typical voltage drop across the cell
 
 Initial Conditions
 -------------------
@@ -101,15 +105,18 @@ sigma_n = Parameter("sigma_n")  # Conductivity in negative electrode
 sigma_p = Parameter("sigma_p")  # Conductivity in positive electrode
 cn_max = Parameter("cn_max")  # Max concentration in negative electrode
 cp_max = Parameter("cp_max")  # Max concentration in positive electrode
+Dn_typ = Parameter("Dn_typ")  # Typical negative particle diffusivity
+Dp_typ = Parameter("Dp_typ")  # Typical positive particle diffusivity
 
 # Electrochemical Reactions
-m_n = Parameter("m_n")  # Reaction rate in negative electrode
-m_p = Parameter("m_p")  # Reaction rate in positive electrode
+m_n_dim = Parameter("m_n_dim")  # Reaction rate in negative electrode
+m_p_dim = Parameter("m_p_dim")  # Reaction rate in positive electrode
 
 # Electrical
 voltage_low_cut = Parameter("voltage_low_cut")  # Lower voltage cut-off
 voltage_high_cut = Parameter("voltage_high_cut")  # Upper voltage cut-off
 I_typ = Parameter("I_typ")  # Typical current density
+Phi_typ = Parameter("Phi_typ")  # Typical voltage drop
 
 # Initial Conditions
 ce0_dimensional = Parameter("ce0")  # Initial li ion concentration in electrolyte
@@ -134,10 +141,92 @@ b = Parameter("b")  # Bruggeman coefficient
 beta_n = a_n * R_n
 beta_p = a_p * R_p
 
+# Discharge timescale
+tau_d = (F * cn_max * Lx) / I_typ
+
+# Particle diffusion timescales
+tau_n = R_n ** 2 / Dn_typ
+tau_p = R_n ** 2 / Dp_typ
+
+# Reaction timescales
+tau_r_n = F / (m_n_dim * a_n * ce_typ ** 0.5)
+tau_r_p = F / (m_p_dim * a_p * ce_typ ** 0.5)
+
+# Scaled maximum concentration in positive particle
+# Note: C_hat_n == 1, so not needed
+C_hat_p = cp_max / cn_max
+
+# Ratio of discharge and solid diffusion timescales
+gamma_n = tau_d / tau_n
+gamma_p = tau_d / tau_p
+
+# Reaction properties
+m_n = tau_d / tau_r_n
+m_p = tau_d / tau_r_n
+
 # Electrolyte Properties
 t_plus = Parameter("t_plus")  # cation transference number
 delta = (Lx ** 2 / De_typ) * (I_typ / (F * cn_max * Lx))
 nu = cn_max / ce_typ
 
+# Ratio of typical to thermal voltage
+Lambda = Phi_typ / (R * T / F)
+
 # Initial conditions
 ce0 = ce0_dimensional / ce_typ
+cn0 = cn0_dimensional / cn_max
+cp0 = cp0_dimensional / cp_max
+
+# -----------------------------------------------------------------------------
+"""Functions"""
+
+
+def D_n_dim(c):
+    "Dimensional diffusivity in negative particle"
+    return Dn_typ
+
+
+def D_n(c):
+    "Dimensionless negative particle diffusivity"
+    return D_n_dim(c) / Dn_typ
+
+
+def D_p_dim(c):
+    "Dimensional diffusivity in positive particle"
+    return Dp_typ
+
+
+def D_p(c):
+    "Dimensionless positive particle diffusivity"
+    return D_p_dim(c) / Dp_typ
+
+
+def U_n(c):
+    "Dimensionless open circuit potential in the negative electrode"
+    #  out = (0.194 + 1.5 * np.exp(-120.0 * c)
+    #       + 0.0351 * np.tanh((c - 0.286) / 0.083)
+    #       - 0.0045 * np.tanh((c - 0.849) / 0.119)
+    #       - 0.035 * np.tanh((c - 0.9233) / 0.05)
+    #       - 0.0147 * np.tanh((c - 0.5) / 0.034)
+    #       - 0.102 * np.tanh((c - 0.194) / 0.142)
+    #       - 0.022 * np.tanh((c - 0.9) / 0.0164)
+    #       - 0.011 * np.tanh((c - 0.124) / 0.0226)
+    #       + 0.0155 * np.tanh((c - 0.105) / 0.029))
+    # Set constant until functions implemented correctly
+    out = 0.2230
+    return out / Phi_typ
+
+
+def U_p(c):
+    "Dimensionless open circuit potential in the positive electrode"
+    # stretch = 1.062
+    # sto = stretch * c
+    # out = (2.16216 + 0.07645 * np.tanh(30.834 - 54.4806 * sto)
+    #       + 2.1581 * np.tanh(52.294 - 50.294 * sto)
+    #       - 0.14169 * np.tanh(11.0923 - 19.8543 * sto)
+    #       + 0.2051 * np.tanh(1.4684 - 5.4888 * sto)
+    #       + 0.2531 * np.tanh((-sto + 0.56478) / 0.1316)
+    #       - 0.02167 * np.tanh((sto - 0.525) / 0.006))
+    # Set constant until functions implemented correctly
+    out = 4.1212
+    return out / Phi_typ


Test patch:
diff --git a/tests/test_models/test_li_ion.py b/tests/test_models/test_li_ion.py
new file mode 100644
index 00000000..d14c23c8
--- /dev/null
+++ b/tests/test_models/test_li_ion.py
@@ -0,0 +1,39 @@
+#
+# Tests for the li-ion models
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+import pybamm
+import tests
+
+import unittest
+import numpy as np
+
+
+class TestLiIonSPM(unittest.TestCase):
+    def test_basic_processing(self):
+        model = pybamm.li_ion.SPM()
+        modeltest = tests.StandardModelTest(model)
+
+        modeltest.test_all()
+
+
+def test_surface_concentrartion(self):
+    model = pybamm.li_ion.SPM()
+    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
+
+    # check surface concentration decreases in negative particle and
+    # increases in positive particle for discharge
+    np.testing.assert_array_less(
+        model.variables["cn_surf"].evaluate(T, Y)[:, 1:],
+        model.variables["cn_surf"].evaluate(T, Y)[:, :-1],
+    )
+    np.testing.assert_array_less(
+        model.variables["cp_surf"].evaluate(T, Y)[:, :-1],
+        model.variables["cp_surf"].evaluate(T, Y)[:, 1:],
+    )