Repository: tinosulzer/PyBaMM

Issue:
Refactor StefanMaxwellDiffusion class
This issue must be completed prior to refactoring the reaction diffusion model #17

StefanMaxwellDiffusion should inherit from BaseModel so that it is a complete model in its own right. It will contain a dictionary with the set of equations in the electrolyte implemented in symbolic form employing the expression tree. 

Will looks something like:

```
class StefanMaxwellDiffusion(BaseModel): 
        def __init__(self): 
            self._rhs = {c: "epression"}
            self._initial_conditions = .... etc 

           self._variable = {"c": c, ...} 
```

EDIT. 
The keys should be of class Symbol expect for the _variable dictionary.


Discussion:

Patch:
diff --git a/input/parameters/lithium-ion/parameters/LCO.csv b/input/parameters/lithium-ion/parameters/LCO.csv
index 63964db9..efcf8c3c 100644
--- a/input/parameters/lithium-ion/parameters/LCO.csv
+++ b/input/parameters/lithium-ion/parameters/LCO.csv
@@ -16,8 +16,9 @@ L_ccn,2.50E-05,m,Negative current collector thickness,Scott Moura,
 L_ccp,2.50E-05,m,Positive current collector thickness,Scott Moura,
 ,,,,,
 # Electrical,,,,,
-voltage_low_cut,3.105,V,Voltage cutoff point for a 6-cell battery,-,
-voltage_high_cut,4.7,V,Number of electrodes connected in parallel to make a cell,-,
+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,-,
 ,,,,,
 # Microstructure,,,,,
 epsilon_n,0.3,-,Electrolyte volume fraction in negative electrode region,Scott Moura,
@@ -25,14 +26,14 @@ 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,
-a_n,90000,m-1,Negative electrode capacity,Scott Moura,
-a_p,90000,m-1,Positive electrode capacity,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,
 ,,,,,
 # Electrolyte physical properties,,,,,
-c_typ,1000,mol.m-3,Typical lithium ion concentration in electrolyte,,
-tpw,0.4,-,Cation transference number,,
-D_e,5.34E-10,m2/s,Typical lithium ion diffusivity,Scott Moura,
+ce_typ,1000,mol.m-3,Typical lithium ion concentration in electrolyte,,
+t_plus,0.4,-,Cation transference number,,
+De_typ,5.34E-10,m2/s,Typical lithium ion diffusivity,Scott Moura,
 ,,,,,
 # Electrode physical properties,,,,,
 sigma_n,100,S.m-1,Conductivity of graphite,Scott Moura,
@@ -46,3 +47,8 @@ D_p,1.00E-13,m2/s,Positive electrode diffusion coefficient,Scott Moura,
 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
 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
diff --git a/pybamm/discretisations/base_discretisation.py b/pybamm/discretisations/base_discretisation.py
index 04def0f9..f262cd54 100644
--- a/pybamm/discretisations/base_discretisation.py
+++ b/pybamm/discretisations/base_discretisation.py
@@ -6,6 +6,7 @@ from __future__ import print_function, unicode_literals
 import pybamm
 
 import numpy as np
+import numbers
 import copy
 
 
@@ -89,10 +90,15 @@ class BaseDiscretisation(object):
 
         """
         for variable, equation in initial_conditions.items():
-            discretised_ic = self.process_symbol(equation, variable.domain)
-            # Turn any scalars into vectors
-            if isinstance(discretised_ic, pybamm.Scalar):
-                discretised_ic = self.scalar_to_vector(discretised_ic, variable.domain)
+            discretised_ic = self.process_symbol(equation, variable.domain).evaluate()
+
+            if isinstance(discretised_ic, numbers.Number):
+                discretised_ic = discretised_ic * self.vector_of_ones(variable.domain)
+            else:
+                raise NotImplementedError(
+                    "Currently only accepts scalar initial conditions"
+                )
+
             initial_conditions[variable] = discretised_ic
 
         # Concatenate and evaluate initial conditions
@@ -221,15 +227,16 @@ class BaseDiscretisation(object):
         """
         raise NotImplementedError
 
-    def scalar_to_vector(self, scalar, domain):
+    def vector_of_ones(self, domain):
         """
-        Convert a Scalar to a uniform Vector of size given by mesh,
-        with same value as Scalar.
+        Returns a Vector of ones of the size given by the mesh.
         """
+
         mesh_points = np.array([])
+
         for dom in domain:
             mesh_points = np.concatenate([mesh_points, self.mesh[dom].nodes])
-        return pybamm.Vector(scalar.value * np.ones_like(mesh_points))
+        return pybamm.Vector(np.ones_like(mesh_points))
 
     def concatenate(self, *symbols):
         return pybamm.NumpyConcatenation(*symbols)
diff --git a/pybamm/models/core.py b/pybamm/models/core.py
index b023c0ac..084f3200 100644
--- a/pybamm/models/core.py
+++ b/pybamm/models/core.py
@@ -11,13 +11,13 @@ class BaseModel(object):
     Attributes
     ----------
 
-    _rhs: dict
+    rhs: dict
         A dictionary that maps expressions (variables) to expressions that represent
         the rhs
-    _initial_conditions: dict
+    initial_conditions: dict
         A dictionary that maps expressions (variables) to expressions that represent
         the initial conditions
-    _boundary_conditions: dict
+    boundary_conditions: dict
         A dictionary that maps expressions (variables) to expressions that represent
         the boundary conditions
     """
diff --git a/pybamm/models/submodels/electrolyte.py b/pybamm/models/submodels/electrolyte.py
index 9eb26b65..b41b5e50 100644
--- a/pybamm/models/submodels/electrolyte.py
+++ b/pybamm/models/submodels/electrolyte.py
@@ -3,180 +3,54 @@
 #
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
-
-import numpy as np
-
-
-class StefanMaxwellDiffusion(object):
-    """Equations for the electrolyte."""
-
-    def __init__(self, subparam, suboperators, submesh, tests):
-        """
-        Assign simulation-specific objects as attributes.
-
-        Parameters
-        ----------
-        subparam : :class:`pybamm.ParameterValues.Parameters` subclass instance
-            The parameters of the simulation
-        suboperators : :class:`pybamm.operators.Operators` subclass instance
-            The spatial operators.
-        submesh : :class:`pybamm.mesh.Mesh` subclass instance
-            The spatial and temporal discretisation.
-        tests : dict
-            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.
-        """
-        # Apply attributes
-        self.subparam = subparam
-        self.suboperators = suboperators
-        self.submesh = submesh
-        self.tests = tests
-
-    def initial_conditions(self):
-        """Calculates initial conditions for variables in the electrolyte.
-
-        Returns
-        -------
-        array_like
-            The initial conditions
-        """
-        if not self.tests:
-            return self.subparam.c0 * np.ones_like(self.submesh.nodes)
-        else:
-            return self.tests["inits"]["concentration"]
-
-    def pdes_rhs(self, vars):
-        """Conservation of cations.
-
-        Parameters
-        ----------
-        vars : array_like, shape (n,)
-            The variables of the simulation
-        flux_bcs : 2-tuple of array_like, shape (1,)
-            Flux at the boundaries (Neumann BC).
-
-        Returns
-        -------
-        dcdt : array_like, shape (n,)
-            The time derivative of c.
-
-        """
-        if not self.tests:
-            flux_bcs = self.bcs()
-            j = vars.j
-        else:
-            flux_bcs = self.tests["bcs"](vars.t)["concentration"]
-            j = self.tests["sources"](vars.t)["concentration"]
-
-        # Calculate internal flux
-        N_internal = -self.suboperators.grad(vars.c)
-
-        # Add boundary conditions (Neumann)
-        flux_bc_left, flux_bc_right = flux_bcs
-        N = np.concatenate([flux_bc_left, N_internal, flux_bc_right])
-
-        # Calculate time derivative
-        dcdt = -self.suboperators.div(N) + self.subparam.s * j
-
-        return dcdt
-
-    def bcs(self):
-        """Flux boundary conditions for the cation conservation equation.
-
-        Returns
-        -------
-        2-tuple of array_like, shape(1,)
-            The boundary conditions.
-
-        """
-        flux_bc_left = np.array([0])
-        flux_bc_right = np.array([0])
-        return (flux_bc_left, flux_bc_right)
-
-    # def current_conservation(self, domain, c, e, j, current_bcs):
-    #     """The 1D diffusion equation.
-    #
-    #     Parameters
-    #     ----------
-    #     subparam : pybamm.ParameterValues.Parameter() instance
-    #         The parameters of the simulation
-    #     variables : 2-tuple (c, e) of array_like, shape (n,)
-    #         The concentration, and potential difference.
-    #     operators : pybamm.operators.Operators() instance
-    #         The spatial operators.
-    #     j : array_like, shape (n,)
-    #         Interfacial current density.
-    #     current_bcs : 2-tuple of array_like, shape (1,)
-    #         Flux at the boundaries (Neumann BC).
-    #
-    #     Returns
-    #     -------
-    #     dedt : array_like, shape (n,)
-    #         The time derivative of the potential.
-    #
-    #     """
-    #     # Calculate current density
-    #     i = self.macinnes(domain, c, e, current_bcs)
-    #
-    #     # Calculate time derivative
-    #     if domain == "xcn":
-    #         gamma_dl = self.subparam.gamma_dl_n
-    #     elif domain == "xcp":
-    #         gamma_dl = self.subparam.gamma_dl_p
-    #
-    #     dedt = 1 / gamma_dl * (self.operators[domain].div(i) - j)
-    #
-    #     return dedt
-    #
-    # def macinnes(self, domain, c, e, current_bcs):
-    #     """MacInnes equation for the electrolyte current density.
-    #
-    #     Parameters
-    #     ----------
-    #     domain : string
-    #         The domain in which to calculate the electrolyte current density.
-    #     c : array_like, shape (n,)
-    #         The electrolyte concentration.
-    #     e : array_like, shape (n,)
-    #         The potential difference.
-    #     current_bcs : 2-tuple of array_like, shape (1,)
-    #         Flux at the boundaries (Neumann BC).
-    #
-    #     Returns
-    #     -------
-    #     i : array_like, shape (n+1,)
-    #         The current density.
-    #     """
-    #     operators = self.operators[domain]
-    #     kappa_over_c = 1
-    #     kappa = 1
-    #
-    #     # Calculate inner current
-    #     i_inner = kappa_over_c * operators.grad(c) + kappa * operators.grad(e)
-    #
-    #     # Add boundary conditions
-    #     lbc, rbc = current_bcs
-    #     i = np.concatenate([lbc, i_inner, rbc])
-    #
-    #     return i
-    #
-    # def bcs_current(self, domain, t):
-    #     """Boundary conditions for the current conservation equation.
-    #
-    #     Returns
-    #     -------
-    #     2-tuple of array_like, shape(1,)
-    #         The boundary conditions.
-    #
-    #     """
-    #     if domain == "xcn":
-    #         current_bc_left = np.array([0])
-    #         current_bc_right = np.array([self.subparam.icell(t)])
-    #     elif domain == "xcp":
-    #         current_bc_left = np.array([self.subparam.icell(t)])
-    #         current_bc_right = np.array([0])
-    #     return (current_bc_left, current_bc_right)
+import pybamm
+
+
+class StefanMaxwellDiffusion(pybamm.BaseModel):
+    """A class that generates the expression tree for Stefan-Maxwell Diffusion in the
+    electrolyte.
+
+    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
+
+    *Extends:* :class:`BaseModel`
+    """
+
+    def __init__(self, G):
+        super().__init__()
+
+        epsilon = pybamm.standard_parameters.epsilon_s  # make issue for spatially
+        # dependent parameters
+        b = pybamm.standard_parameters.b
+        delta = pybamm.standard_parameters.delta
+        nu = pybamm.standard_parameters.nu
+        t_plus = pybamm.standard_parameters.t_plus
+        ce0 = pybamm.standard_parameters.ce0
+
+        electrolyte_domain = ["whole cell"]
+
+        c_e = pybamm.Variable("c_e", domain=electrolyte_domain)
+
+        N_e = -(epsilon ** b) * pybamm.Gradient(c_e)
+
+        self.rhs = {
+            c_e: -pybamm.Divergence(N_e) / delta / epsilon + nu * (1 - t_plus) * G
+        }
+        self.initial_conditions = {c_e: ce0}
+        self.boundary_conditions = {
+            N_e: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(0)}
+        }
+        self.variables = {"c_e": c_e, "N_e": N_e}
diff --git a/pybamm/parameters/standard_parameters.py b/pybamm/parameters/standard_parameters.py
index 82f7912e..efa65d97 100644
--- a/pybamm/parameters/standard_parameters.py
+++ b/pybamm/parameters/standard_parameters.py
@@ -4,29 +4,134 @@
 """
 Standard Parameters for battery models
 
-Geometric
----------
+Physical Constants
+------------------
+R
+    Ideal gas constant
+F
+    Faraday's constant
+T
+    Reference temperature
+
+Macroscale Geometry
+-------------------
 
 Ln, Ls, Lp
-    The widths of the negative electrode, separator and positive electrode respectively
-L
+    The widths of the negative electrode, separator and positive electrode, respectively
+Lx
     The width of a single cell
 ln, ls, lp
     The dimesionless widths of the negative electrode, separator and positive
     electrode respectively
 
+Microscale Geometry
+-------------------
+R_n, R_p
+    Negative and positive particle radii
+a_n, a_p
+    Negative and positive electrode surface area densities
+
+Electrolyte Properties
+----------------------
+ce_typ
+    Typical lithium ion concentration in electrolyte
+De_typ
+    Typical lithium ion diffusivity in the electrolyte
+
+Electrode Properties
+--------------------
+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
+
+Electrochemical Reactions
+--------------------------
+m_n, m_p
+    Reaction rates in negative and positive electrode regions
+
+Electrical
+----------
+voltage_low_cut, voltage_high_cut
+    Low and high voltage cut-offs
+I_typ
+    Typical current density
+
+Initial Conditions
+-------------------
+ce0_dimensional
+    Initial lithium ion concentration in the electrolyte
+cn0_dimensional, cp0_dimensional
+    Initial lithium concentration in the negative and positive electrodes
 """
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 from pybamm import Parameter
 
-# Lengths
+# --------------------------------------------------------------------------------------
+"""Dimensional Parameters"""
+
+# Physical Constants
+R = Parameter("R")
+F = Parameter("F")
+T = Parameter("T_ref")
+
+# Macroscale Geometry
 Ln = Parameter("Ln")
 Ls = Parameter("Ls")
 Lp = Parameter("Lp")
-L = Ln + Ls + Lp
+Lx = Ln + Ls + Lp
+
+# Microscale Geometry
+R_n = Parameter("R_n")
+R_p = Parameter("R_p")
+a_n = Parameter("a_n")
+a_p = Parameter("a_p")
+
+# Electrolyte Properties
+ce_typ = Parameter("ce_typ")  # Typical lithium ion concentration in electrolyte
+De_typ = Parameter("De_typ")  # Typical electrolyte diffusivity
+
+# Electrode Properties
+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
+
+# Electrochemical Reactions
+m_n = Parameter("m_n")  # Reaction rate in negative electrode
+m_p = Parameter("m_p")  # 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
+
+# Initial Conditions
+ce0_dimensional = Parameter("ce0")  # Initial li ion concentration in electrolyte
+cn0_dimensional = Parameter("cn0")  # Initial li concentration in neg electrode
+cp0_dimensional = Parameter("cp0")  # Initial li concentration in pos electrode
+
+# --------------------------------------------------------------------------------------
+"""Dimensionless Parameters"""
+
+# Macroscale Geometry
+ln = Ln / Lx
+ls = Ls / Lx
+lp = Lp / Lx
+
+# Microscale Geometry
+epsilon_n = Parameter("epsilon_n")  # Electrolyte volume fraction in neg electrode
+epsilon_s = Parameter("epsilon_s")  # Electrolyte volume fraction in separator
+epsilon_p = Parameter("epsilon_p")  # Electrolyte volume fraction in pos electrode
+b = Parameter("b")  # Bruggeman coefficient
+beta_n = a_n * R_n
+beta_p = a_p * R_p
+
+# 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
 
-# Dimensionless lengths
-ln = Ln / L
-ls = Ls / L
-lp = Lp / L
+# Initial conditions
+ce0 = ce0_dimensional / ce_typ


Test patch:
diff --git a/tests/test_discretisations/test_base_discretisations.py b/tests/test_discretisations/test_base_discretisations.py
index e2f873d5..e4378d8e 100644
--- a/tests/test_discretisations/test_base_discretisations.py
+++ b/tests/test_discretisations/test_base_discretisations.py
@@ -270,13 +270,12 @@ class TestDiscretise(unittest.TestCase):
         )
         np.testing.assert_array_equal(y0, model.rhs.evaluate(None, y0))
 
-    def test_scalar_to_vector(self):
-        a = pybamm.Scalar(5)
+    def test_vector_of_ones(self):
         mesh = MeshForTesting()
         disc = pybamm.BaseDiscretisation(mesh)
-        a_vec = disc.scalar_to_vector(a, ["whole cell"])
-        self.assertEqual(a_vec.evaluate(None)[0], a.value)
-        self.assertEqual(a_vec.shape, mesh["whole cell"].nodes.shape)
+        vec = disc.vector_of_ones(["whole cell"])
+        self.assertEqual(vec.evaluate(None)[0], 1)
+        self.assertEqual(vec.shape, mesh["whole cell"].nodes.shape)
 
     def test_concatenation(self):
         a = pybamm.Symbol("a")
diff --git a/tests/test_models/__init__.py b/tests/test_models/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/test_models/test_electrolyte.py b/tests/test_models/test_electrolyte.py
new file mode 100644
index 00000000..8b8a2cd0
--- /dev/null
+++ b/tests/test_models/test_electrolyte.py
@@ -0,0 +1,40 @@
+#
+# Tests for the Binary Operator classes
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+import pybamm
+
+import unittest
+
+
+class TestStefanMaxwellDiffusion(unittest.TestCase):
+    def test_make_tree(self):
+        G = pybamm.Scalar(1)
+        pybamm.electrolyte.StefanMaxwellDiffusion(G)
+
+    def test_processing_parameters(self):
+        G = pybamm.Scalar(1)
+        model = pybamm.electrolyte.StefanMaxwellDiffusion(G)
+
+        param = pybamm.ParameterValues(
+            "input/parameters/lithium-ion/parameters/LCO.csv"
+        )
+
+        param.process_model(model)
+
+    def test_processing_disc(self):
+        G = pybamm.Scalar(1)
+        model = pybamm.electrolyte.StefanMaxwellDiffusion(G)
+
+        param = pybamm.ParameterValues(
+            "input/parameters/lithium-ion/parameters/LCO.csv"
+        )
+
+        param.process_model(model)
+
+        print(model.initial_conditions)
+
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 2)
+        disc = pybamm.FiniteVolumeDiscretisation(mesh)
+        disc.process_model(model)
diff --git a/tests/test_parameters/test_standard_parameters.py b/tests/test_parameters/test_standard_parameters.py
index b5d8cf5c..d45ae623 100644
--- a/tests/test_parameters/test_standard_parameters.py
+++ b/tests/test_parameters/test_standard_parameters.py
@@ -11,7 +11,7 @@ class TestStandardParameters(unittest.TestCase):
         Ln = pybamm.standard_parameters.Ln
         Ls = pybamm.standard_parameters.Ls
         Lp = pybamm.standard_parameters.Lp
-        L = pybamm.standard_parameters.L
+        Lx = pybamm.standard_parameters.Lx
         ln = pybamm.standard_parameters.ln
         ls = pybamm.standard_parameters.ls
         lp = pybamm.standard_parameters.lp
@@ -22,9 +22,9 @@ class TestStandardParameters(unittest.TestCase):
         Ln_eval = parameter_values.process_symbol(Ln)
         Ls_eval = parameter_values.process_symbol(Ls)
         Lp_eval = parameter_values.process_symbol(Lp)
-        L_eval = parameter_values.process_symbol(L)
-        self.assertEqual((Ln_eval + Ls_eval + Lp_eval).evaluate(), L_eval.evaluate())
-        self.assertEqual((Ln_eval + Ls_eval + Lp_eval).id, L_eval.id)
+        Lx_eval = parameter_values.process_symbol(Lx)
+        self.assertEqual((Ln_eval + Ls_eval + Lp_eval).evaluate(), Lx_eval.evaluate())
+        self.assertEqual((Ln_eval + Ls_eval + Lp_eval).id, Lx_eval.id)
         ln_eval = parameter_values.process_symbol(ln)
         ls_eval = parameter_values.process_symbol(ls)
         lp_eval = parameter_values.process_symbol(lp)