Repository: tinosulzer/PyBaMM

Issue:
Interfacial reactions
Sub-issue from #120 for the Butler-Volmer interfacial reaction functions

Discussion:

Patch:
diff --git a/docs/source/expression_tree/concatenations.rst b/docs/source/expression_tree/concatenations.rst
index 1665c443..affb5ead 100644
--- a/docs/source/expression_tree/concatenations.rst
+++ b/docs/source/expression_tree/concatenations.rst
@@ -10,4 +10,4 @@ Concatenations
 .. autoclass:: pybamm.DomainConcatenation
   :members:
 
-.
+.. autofunction:: pybamm.piecewise_constant
diff --git a/input/parameters/lead-acid/electrolyte_diffusivity.py b/input/parameters/lead-acid/electrolyte_diffusivity.py
new file mode 100644
index 00000000..d7d9f3c5
--- /dev/null
+++ b/input/parameters/lead-acid/electrolyte_diffusivity.py
@@ -0,0 +1,193 @@
+#
+# Electrolyte diffusivity
+#
+
+
+def electrolyte_diffusivity(c):
+    """
+    Dimensional effective Fickian diffusivity in the electrolyte [m2.s-1].
+
+    Parameters
+    ----------
+    c
+        Electrolyte concentration [mol.m-3]
+    """
+    return (1.75 + 260e-6 * c) * 1e-9
+
+
+#
+# def kappa_hat(param, c):
+#     """Dimensional effective conductivity in the electrolyte [S.m-1]"""
+#     return c * np.exp(6.23 - 1.34e-4 * c - 1.61e-8 * c ** 2) * 1e-4
+#
+#
+# def kappa_eff(param, c, eps):
+#     """Dimensionless molar conductivity in the electrolyte"""
+#     kappa_scale = (
+#         param._raw["F"] ** 2
+#         * param._raw["cmax"]
+#         * D_hat(param._raw["cmax"])
+#         / (param._raw["R"] * param._raw["T_ref"])
+#     )
+#     return kappa_hat(param, c * param._raw["cmax"]) / kappa_scale * (eps ** 1.5)
+#
+#
+# def chi_hat(c):
+#     """Dimensional Darken thermodynamic factor in the electrolyte [-]"""
+#     return 0.49 + 4.1e-4 * c
+#
+#
+# def chi(param, c):
+#     """Dimensionless Darken thermodynamic factor in the electrolyte"""
+#     # Excluded volume fraction
+#     alpha = (2 * param._raw["Vw"] - param._raw["Ve"]) * param.scales["conc"]
+#     chi_scale = 1 / (2 * (1 - param._raw["tpw"]))
+#     return chi_hat(c * param._raw["cmax"]) / chi_scale / (1 + alpha * c)
+#
+#
+# # def curlyK_hat(param, eps):
+# #     """Dimensional permeability [m2]"""
+# #     return eps ** 3 * param._raw["d"] ** 2 / (180 * (1 - eps) ** 2)
+# #
+# #
+# # def curlyK(param, eps):
+# #     """Dimensionless permeability"""
+# #     return curlyK_hat(eps) / param._raw["d"] ** 2
+# #
+#
+#
+# def mu_hat(c):
+#     """Dimensional viscosity of electrolyte [kg.m-1.s-1]"""
+#     return 0.89e-3 + 1.11e-7 * c + 3.29e-11 * c ** 2
+#
+#
+# #
+# # def mu(param, c):
+# #     """Dimensionless viscosity of electrolyte"""
+# #     return mu_hat(c * param._raw["cmax"]) / mu_hat(param._raw["cmax"])
+# #
+# #
+# # def rho_hat(param, c):
+# #     """Dimensional density of electrolyte [kg.m-3]"""
+# #     return (
+# #         Mw
+# #         / param._raw["Vw"]
+# #         * (
+# #             1
+# #             + (
+# #                 param._raw["Me"] * param._raw["Vw"] / param._raw["Mw"]
+# #                 - param._raw["Ve"]
+# #             )
+# #             * c
+# #         )
+# #     )
+# #
+# #
+# # def rho(param, c):
+# #     """Dimensionless density of electrolyte"""
+# #     return rho_hat(c * param._raw["cmax"]) / rho_hat(param._raw["cmax"])
+#
+#
+# def cw_hat(param, c):
+#     """Dimensional solvent concentration [mol.m-3]"""
+#     return (1 - c * param._raw["Ve"]) / param._raw["Vw"]
+#
+#
+# def cw(param, c):
+#     """Dimensionless solvent concentration"""
+#     return cw_hat(param, c * param._raw["cmax"]) / cw_hat(param, param._raw["cmax"])
+#
+#
+# # def dcwdc(param, c):
+# #     """Dimensionless derivative of cw with respect to c"""
+# #     return 0 * c - param._raw["Ve"] / param._raw["Vw"]
+#
+#
+# def m(param, c):
+#     """Dimensional electrolyte molar mass [mol.kg-1]"""
+#     return c * param._raw["Vw"] / ((1 - c * param._raw["Ve"]) * param._raw["Mw"])
+#
+#
+# # def dmdc(param, c):
+# #     """Dimensional derivative of m with respect to c [kg-1]"""
+# #     return param.Vw / ((1 - c * param._raw["Ve"]) ** 2 * param._raw["Mw"])
+#
+#
+# def U_Pb(param, c):
+#     """Dimensionless OCP in the negative electrode"""
+#     m_ = m(param, c * param._raw["cmax"])  # dimensionless
+#     U = (
+#         param._raw["F"]
+#         / (param._raw["R"] * param._raw["T_ref"])
+#         * (
+#             -0.074 * np.log10(m_)
+#             - 0.030 * np.log10(m_) ** 2
+#             - 0.031 * np.log10(m_) ** 3
+#             - 0.012 * np.log10(m_) ** 4
+#         )
+#     )
+#     return U
+#
+#
+# # def U_Pb_hat(param, c):
+# #     """Dimensional OCP in the negative electrode [V]"""
+#     return param._raw["U_Pb_ref"] + param._raw["R"] * param._raw["T_ref"] /param._raw[
+# #         "F"
+# #     ] * U_Pb(c / param._raw["cmax"])
+#
+#
+# # def dUPbdc(param, c):
+# #     """Dimensionless derivative of U_Pb with respect to c"""
+# #     m_ = m(c * param._raw["cmax"])  # dimensionless
+# #     dUdm = (
+# #         param._raw["F"]
+# #         / (param._raw["R"] * param._raw["T_ref"])
+# #         * (
+# #             -0.074 / m_ / np.log(10)
+# #             - 0.030 * 2 * np.log(m_) / (m_ * np.log(10) ** 2)
+# #             - 0.031 * 3 * np.log(m_) ** 2 / m_ / np.log(10) ** 3
+# #             - 0.012 * 4 * np.log(m_) ** 3 / m_ / np.log(10) ** 4
+# #         )
+# #     )
+# #     dmdc = dmdc(c * param._raw["cmax"]) * param._raw["cmax"]  # dimensionless
+# #     return dmdc * dUdm
+#
+#
+# def U_PbO2(param, c):
+#     """Dimensionless OCP in the positive electrode"""
+#     m_ = m(param, c * param._raw["cmax"])
+#     U = (
+#         param._raw["F"]
+#         / (param._raw["R"] * param._raw["T_ref"])
+#         * (
+#             0.074 * np.log10(m_)
+#             + 0.033 * np.log10(m_) ** 2
+#             + 0.043 * np.log10(m_) ** 3
+#             + 0.022 * np.log10(m_) ** 4
+#         )
+#     )
+#     return U
+#
+#
+# # def U_PbO2_hat(param, c):
+# #     """Dimensional OCP in the positive electrode [V]"""
+# #     return param._raw["U_PbO2_ref"] + param._raw["R"] * param._raw[
+# #         "T_ref"
+# #     ] / param._raw["F"] * U_PbO2(c / param._raw["cmax"])
+# #
+# #
+# # def dUPbO2dc(param, c):
+# #     """Dimensionless derivative of U_PbO2 with respect to c"""
+# #     m_ = m(c * param._raw["cmax"])  # dimensionless
+# #     dUdm = (
+# #         param._raw["F"]
+# #         / (param._raw["R"] * param._raw["T_ref"])
+# #         * (
+# #             0.074 / m_ / np.log(10)
+# #             + 0.033 * 2 * np.log(m_) / (m_ * np.log(10) ** 2)
+# #             + 0.043 * 3 * np.log(m_) ** 2 / m_ / np.log(10) ** 3
+# #             + 0.022 * 4 * np.log(m_) ** 3 / m_ / np.log(10) ** 4
+# #         )
+# #     )
+# #     dmdc = dmdc(c * param._raw["cmax"]) * param._raw["cmax"]  # dimensionless
+# #     return dmdc * dUdm
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 91bff3ab..7d9c9a77 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -63,6 +63,7 @@ from .expression_tree.concatenations import (
     Concatenation,
     NumpyModelConcatenation,
     DomainConcatenation,
+    piecewise_constant,
 )
 from .expression_tree.array import Array
 from .expression_tree.matrix import Matrix
@@ -105,8 +106,8 @@ from .models.submodels import electrolyte, interface
 #
 # Parameters class and methods
 #
+from .meshes.meshes import KNOWN_DOMAINS  # need this for importing standard parameters
 from .parameters.parameter_values import ParameterValues
-from .parameters import functions_lead_acid
 from .parameters import standard_parameters
 from .parameters import standard_parameters_lead_acid  # calls standard_parameters
 
@@ -128,7 +129,6 @@ from .discretisations.finite_volume_discretisations import (
     FiniteVolumeDiscretisation,
     NodeToEdge,
 )
-from .meshes.meshes import KNOWN_DOMAINS
 from .meshes.meshes import Mesh
 from .meshes.submeshes import SubMesh1D, Uniform1DSubMesh
 
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index 6f99e939..3cd0b3e6 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -177,3 +177,25 @@ class DomainConcatenation(Concatenation):
                 vector[self._slices[dom]] = child_vector[slices[dom]]
 
         return vector
+
+
+def piecewise_constant(neg_value, sep_value, pos_value):
+    """Concatenate three values that don't have a domain into a piecewise constant
+    concatenation. This is useful when we don't want to assign a domain to the inputs
+
+    Parameters
+    ----------
+    neg_value, sep_value, pos_value : :class:`numbers.Number` or :class:`pybamm.Symbol`
+        The constant values to be concatenated
+
+    Returns
+    -------
+    :class:`pybamm.Concantenation`
+        The piecewise constant concatenation
+    """
+    neg_value_with_domain = neg_value * pybamm.Scalar(1, domain=["negative electrode"])
+    sep_value_with_domain = sep_value * pybamm.Scalar(1, domain=["separator"])
+    pos_value_with_domain = pos_value * pybamm.Scalar(1, domain=["positive electrode"])
+    return pybamm.Concatenation(
+        neg_value_with_domain, sep_value_with_domain, pos_value_with_domain
+    )
diff --git a/pybamm/models/submodels/interface.py b/pybamm/models/submodels/interface.py
index bb709551..5a32ac42 100644
--- a/pybamm/models/submodels/interface.py
+++ b/pybamm/models/submodels/interface.py
@@ -19,3 +19,121 @@ def homogeneous_reaction():
         -pybamm.Scalar(1, domain=["positive electrode"]) / pybamm.standard_parameters.lp
     )
     return pybamm.Concatenation(current_neg, current_sep, current_pos)
+
+
+def butler_volmer_lead_acid(c, phi, domain=None):
+    """
+    Butler-Volmer reactions for lead-acid chemistry
+
+    .. math::
+        j = j_0(c) * \\sinh(\\phi - U(c)),
+
+        \\text{where} \\phi = \\Phi_\\text{s} - \\Phi
+
+    Parameters
+    ----------
+    c : :class:`pybamm.Symbol`
+        The electrolyte concentration
+    phi : :class:`pybamm.Symbol`
+        The difference betweent the solid potential and electrolyte potential
+    domain : iterable of strings
+        The domain in which to calculate the interfacial current density. Default is
+        None, in which case the domain is calculated based on c and phi or defaults to
+        the domain spanning the whole cell
+
+    Returns
+    -------
+    :class:`pybamm.Symbol`
+        The interfacial current density in the appropriate domain
+    """
+    if domain is None:
+        # raise error if no domain can be found
+        if c.domain == [] and phi.domain == []:
+            raise ValueError(
+                "domain cannot be None if c.domain and phi.domain are empty"
+            )
+        # otherwise read domain from c and phi, making sure they are consistent with
+        # each other
+        else:
+            if phi.domain == c.domain or c.domain == [] or phi.domain == []:
+                domain = c.domain
+            else:
+                raise pybamm.DomainError(
+                    "c and phi must have the same (or empty) domain"
+                )
+
+    # Get the current densities based on domain
+    if domain == ["negative electrode"]:
+        j0n = exchange_current_density(c, ["negative electrode"])
+        etan = phi - pybamm.standard_parameters_lead_acid.U_Pb(c)
+        return j0n * etan  # Function(etan, np.sinh)
+    elif domain == ["positive electrode"]:
+        j0p = exchange_current_density(c, ["positive electrode"])
+        etap = phi - pybamm.standard_parameters_lead_acid.U_PbO2(c)
+        return j0p * etap  # Function(etap, np.sinh)
+    # To get current density across the whole domain, unpack and call this function
+    # again in the subdomains, then concatenate
+    elif domain == ["negative electrode", "separator", "positive electrode"]:
+        # Unpack c
+        if all([isinstance(var, pybamm.Concatenation) for var in [c, phi]]):
+            cn, cs, cp = c.orphans
+            phin, phis, phip = phi.orphans
+        else:
+            raise ValueError(
+                "c and phi must both be Concatenations, not '{}' and '{}'".format(
+                    type(c), type(phi)
+                )
+            )
+        # Negative electrode
+        current_neg = butler_volmer_lead_acid(cn, phin, domain=["negative electrode"])
+        # Separator
+        current_sep = pybamm.Scalar(0, domain=["separator"])
+        # Positive electrode
+        current_pos = butler_volmer_lead_acid(cp, phip, domain=["positive electrode"])
+        # Concatenate
+        return pybamm.Concatenation(current_neg, current_sep, current_pos)
+    else:
+        raise pybamm.DomainError("domain '{}' not recognised".format(domain))
+
+
+def exchange_current_density(c, domain=None):
+    """The exchange current-density as a function of concentration
+
+    Parameters
+    ----------
+    c : :class:`pybamm.Variable`
+        A Variable representing the concentration
+    domain : string
+        Which domain to calculate the exchange current density in ("negative electrode"
+        or "positive electrode"). Default is None, in which case the domain is c.domain
+
+    Returns
+    -------
+    :class:`pybamm.Symbol`
+        The exchange-current density
+    """
+    if domain is None:
+        # read domain from c if it exists
+        if c.domain != []:
+            domain = c.domain
+        # otherwise raise error
+        else:
+            raise ValueError("domain cannot be None if c.domain is empty")
+
+    if domain == ["negative electrode"]:
+        # concentration domain should be empty or the same as domain
+        if c.domain not in [["negative electrode"], []]:
+            raise pybamm.DomainError("""concentration and domain do not match""")
+        iota_ref_n = pybamm.standard_parameters_lead_acid.iota_ref_n
+        return iota_ref_n * c
+    elif domain == ["positive electrode"]:
+        # concentration domain should be empty or the same as domain
+        if c.domain not in [["positive electrode"], []]:
+            raise pybamm.DomainError("""concentration and domain do not match""")
+        iota_ref_p = pybamm.standard_parameters_lead_acid.iota_ref_p
+        Ve = pybamm.standard_parameters_lead_acid.Ve
+        Vw = pybamm.standard_parameters_lead_acid.Vw
+        cw = (1 - c * Ve) / Vw
+        return iota_ref_p * c ** 2 * cw
+    else:
+        raise pybamm.DomainError("domain '{}' not recognised".format(domain))
diff --git a/pybamm/parameters/functions_lead_acid.py b/pybamm/parameters/functions_lead_acid.py
deleted file mode 100644
index 7a5c0638..00000000
--- a/pybamm/parameters/functions_lead_acid.py
+++ /dev/null
@@ -1,196 +0,0 @@
-#
-# Parameter functions for the lead-acid models
-#
-from __future__ import absolute_import, division
-from __future__ import print_function, unicode_literals
-
-import numpy as np
-
-
-def D_hat(c):
-    """
-    Dimensional effective Fickian diffusivity in the electrolyte [m2.s-1].
-    """
-    return (1.75 + 260e-6 * c) * 1e-9
-
-
-def D_eff(param, c, eps):
-    """Dimensionless effective Fickian diffusivity in the electrolyte."""
-    return D_hat(c * param._raw["cmax"]) / D_hat(param._raw["cmax"]) * (eps ** 1.5)
-
-
-def kappa_hat(param, c):
-    """Dimensional effective conductivity in the electrolyte [S.m-1]"""
-    return c * np.exp(6.23 - 1.34e-4 * c - 1.61e-8 * c ** 2) * 1e-4
-
-
-def kappa_eff(param, c, eps):
-    """Dimensionless molar conductivity in the electrolyte"""
-    kappa_scale = (
-        param._raw["F"] ** 2
-        * param._raw["cmax"]
-        * D_hat(param._raw["cmax"])
-        / (param._raw["R"] * param._raw["T_ref"])
-    )
-    return kappa_hat(param, c * param._raw["cmax"]) / kappa_scale * (eps ** 1.5)
-
-
-def chi_hat(c):
-    """Dimensional Darken thermodynamic factor in the electrolyte [-]"""
-    return 0.49 + 4.1e-4 * c
-
-
-def chi(param, c):
-    """Dimensionless Darken thermodynamic factor in the electrolyte"""
-    # Excluded volume fraction
-    alpha = (2 * param._raw["Vw"] - param._raw["Ve"]) * param.scales["conc"]
-    chi_scale = 1 / (2 * (1 - param._raw["tpw"]))
-    return chi_hat(c * param._raw["cmax"]) / chi_scale / (1 + alpha * c)
-
-
-# def curlyK_hat(param, eps):
-#     """Dimensional permeability [m2]"""
-#     return eps ** 3 * param._raw["d"] ** 2 / (180 * (1 - eps) ** 2)
-#
-#
-# def curlyK(param, eps):
-#     """Dimensionless permeability"""
-#     return curlyK_hat(eps) / param._raw["d"] ** 2
-#
-
-
-def mu_hat(c):
-    """Dimensional viscosity of electrolyte [kg.m-1.s-1]"""
-    return 0.89e-3 + 1.11e-7 * c + 3.29e-11 * c ** 2
-
-
-#
-# def mu(param, c):
-#     """Dimensionless viscosity of electrolyte"""
-#     return mu_hat(c * param._raw["cmax"]) / mu_hat(param._raw["cmax"])
-#
-#
-# def rho_hat(param, c):
-#     """Dimensional density of electrolyte [kg.m-3]"""
-#     return (
-#         Mw
-#         / param._raw["Vw"]
-#         * (
-#             1
-#             + (
-#                 param._raw["Me"] * param._raw["Vw"] / param._raw["Mw"]
-#                 - param._raw["Ve"]
-#             )
-#             * c
-#         )
-#     )
-#
-#
-# def rho(param, c):
-#     """Dimensionless density of electrolyte"""
-#     return rho_hat(c * param._raw["cmax"]) / rho_hat(param._raw["cmax"])
-
-
-def cw_hat(param, c):
-    """Dimensional solvent concentration [mol.m-3]"""
-    return (1 - c * param._raw["Ve"]) / param._raw["Vw"]
-
-
-def cw(param, c):
-    """Dimensionless solvent concentration"""
-    return cw_hat(param, c * param._raw["cmax"]) / cw_hat(param, param._raw["cmax"])
-
-
-# def dcwdc(param, c):
-#     """Dimensionless derivative of cw with respect to c"""
-#     return 0 * c - param._raw["Ve"] / param._raw["Vw"]
-
-
-def m(param, c):
-    """Dimensional electrolyte molar mass [mol.kg-1]"""
-    return c * param._raw["Vw"] / ((1 - c * param._raw["Ve"]) * param._raw["Mw"])
-
-
-# def dmdc(param, c):
-#     """Dimensional derivative of m with respect to c [kg-1]"""
-#     return param.Vw / ((1 - c * param._raw["Ve"]) ** 2 * param._raw["Mw"])
-
-
-def U_Pb(param, c):
-    """Dimensionless OCP in the negative electrode"""
-    m_ = m(param, c * param._raw["cmax"])  # dimensionless
-    U = (
-        param._raw["F"]
-        / (param._raw["R"] * param._raw["T_ref"])
-        * (
-            -0.074 * np.log10(m_)
-            - 0.030 * np.log10(m_) ** 2
-            - 0.031 * np.log10(m_) ** 3
-            - 0.012 * np.log10(m_) ** 4
-        )
-    )
-    return U
-
-
-# def U_Pb_hat(param, c):
-#     """Dimensional OCP in the negative electrode [V]"""
-#     return param._raw["U_Pb_ref"] + param._raw["R"] * param._raw["T_ref"] /param._raw[
-#         "F"
-#     ] * U_Pb(c / param._raw["cmax"])
-
-
-# def dUPbdc(param, c):
-#     """Dimensionless derivative of U_Pb with respect to c"""
-#     m_ = m(c * param._raw["cmax"])  # dimensionless
-#     dUdm = (
-#         param._raw["F"]
-#         / (param._raw["R"] * param._raw["T_ref"])
-#         * (
-#             -0.074 / m_ / np.log(10)
-#             - 0.030 * 2 * np.log(m_) / (m_ * np.log(10) ** 2)
-#             - 0.031 * 3 * np.log(m_) ** 2 / m_ / np.log(10) ** 3
-#             - 0.012 * 4 * np.log(m_) ** 3 / m_ / np.log(10) ** 4
-#         )
-#     )
-#     dmdc = dmdc(c * param._raw["cmax"]) * param._raw["cmax"]  # dimensionless
-#     return dmdc * dUdm
-
-
-def U_PbO2(param, c):
-    """Dimensionless OCP in the positive electrode"""
-    m_ = m(param, c * param._raw["cmax"])
-    U = (
-        param._raw["F"]
-        / (param._raw["R"] * param._raw["T_ref"])
-        * (
-            0.074 * np.log10(m_)
-            + 0.033 * np.log10(m_) ** 2
-            + 0.043 * np.log10(m_) ** 3
-            + 0.022 * np.log10(m_) ** 4
-        )
-    )
-    return U
-
-
-# def U_PbO2_hat(param, c):
-#     """Dimensional OCP in the positive electrode [V]"""
-#     return param._raw["U_PbO2_ref"] + param._raw["R"] * param._raw[
-#         "T_ref"
-#     ] / param._raw["F"] * U_PbO2(c / param._raw["cmax"])
-#
-#
-# def dUPbO2dc(param, c):
-#     """Dimensionless derivative of U_PbO2 with respect to c"""
-#     m_ = m(c * param._raw["cmax"])  # dimensionless
-#     dUdm = (
-#         param._raw["F"]
-#         / (param._raw["R"] * param._raw["T_ref"])
-#         * (
-#             0.074 / m_ / np.log(10)
-#             + 0.033 * 2 * np.log(m_) / (m_ * np.log(10) ** 2)
-#             + 0.043 * 3 * np.log(m_) ** 2 / m_ / np.log(10) ** 3
-#             + 0.022 * 4 * np.log(m_) ** 3 / m_ / np.log(10) ** 4
-#         )
-#     )
-#     dmdc = dmdc(c * param._raw["cmax"]) * param._raw["cmax"]  # dimensionless
-#     return dmdc * dUdm
diff --git a/pybamm/parameters/standard_parameters_lead_acid.py b/pybamm/parameters/standard_parameters_lead_acid.py
index 83bd5b5a..455ba324 100644
--- a/pybamm/parameters/standard_parameters_lead_acid.py
+++ b/pybamm/parameters/standard_parameters_lead_acid.py
@@ -91,18 +91,54 @@ Cdl = pybamm.Parameter("Cdl")  # Double-layer capacity [F.m-2]
 U_Pb_ref = pybamm.Parameter("U_Pb_ref")  # Reference OCP in the lead [V]
 U_PbO2_ref = pybamm.Parameter("U_PbO2_ref")  # Reference OCP in the lead dioxide [V]
 
-# Functions
-D_dim = pybamm.Parameter("epsn_max")
+# --------------------------------------------------------------------------------------
+"""Functions"""
+
 rho_dim = pybamm.Parameter("epsn_max")
 mu_dim = pybamm.Parameter("epsn_max")
 
+
+def D_dim(c):
+    "Dimensional electrolyte diffusivity"
+    return 1  # pybamm.FunctionParameter("electrolyte_diffusivity", c)
+
+
+def D(c):
+    "Dimensionless electrolyte diffusivity"
+    return D_dim(c * cmax) / D_dim(cmax)
+
+
+def U_Pb(c):
+    "Dimensional open-circuit potential for lead electrode"
+    return U_Pb_ref
+
+
+def U_PbO2(c):
+    "Dimensionless open-circuit potential for lead-dioxide electrode"
+    return U_PbO2_ref
+
+
+# --------------------------------------------------------------------------------------
+"""Scales"""
+Lx = pybamm.standard_parameters.Lx
+F = pybamm.standard_parameters.F
+R = pybamm.standard_parameters.R
+T = pybamm.standard_parameters.T
+
+concentration_scale = cmax
+length_scale = Lx
+time_scale = cmax * F * Lx / ibar  # discharge time scale
+potential_scale = R * T / F
+current_scale = ibar
+interfacial_current_scale_neg = ibar / (An * Lx)
+interfacial_current_scale_pos = ibar / (Ap * Lx)
+velocity_scale = ibar / (cmax * F)  # Reaction velocity scale
+
 # --------------------------------------------------------------------------------------
 """Dimensionless Parameters"""
 
 Cd = (
-    (pybamm.standard_parameters.Lx ** 2)
-    / D_dim
-    / (cmax * pybamm.standard_parameters.F * pybamm.standard_parameters.Lx / ibar)
+    (Lx ** 2) / D_dim(cmax) / time_scale
 )  # Diffusional C-rate: diffusion timescale/discharge timescale
 alpha = (nu * Vw - Ve) * cmax  # Excluded volume fraction
 sn = (
@@ -112,24 +148,16 @@ sp = (
     -(s_plus_p + ne_p * pybamm.standard_parameters.t_plus) / ne_p
 )  # Dimensionless rection rate (pos)
 iota_s_n = (
-    sigma_eff_n
-    * pybamm.standard_parameters.R
-    * pybamm.standard_parameters.T
-    / (pybamm.standard_parameters.F * pybamm.standard_parameters.Lx)
-    / ibar
+    sigma_eff_n * potential_scale / current_scale
 )  # Dimensionless lead conductivity
 iota_s_p = (
-    sigma_eff_p
-    * pybamm.standard_parameters.R
-    * pybamm.standard_parameters.T
-    / (pybamm.standard_parameters.F * pybamm.standard_parameters.Lx)
-    / ibar
+    sigma_eff_p * potential_scale / current_scale
 )  # Dimensionless lead dioxide conductivity
-iota_ref_n = jref_n / (
-    ibar / (An * pybamm.standard_parameters.Lx)
+iota_ref_n = (
+    jref_n / interfacial_current_scale_neg
 )  # Dimensionless exchange-current density (neg)
-iota_ref_p = jref_p / (
-    ibar / (Ap * pybamm.standard_parameters.Lx)
+iota_ref_p = (
+    jref_p / interfacial_current_scale_pos
 )  # Dimensionless exchange-current density (pos)
 beta_surf_n = -cmax * DeltaVsurf_n / ne_n  # Molar volume change (lead)
 beta_surf_p = -cmax * DeltaVsurf_p / ne_p  # Molar volume change (lead dioxide)
@@ -144,36 +172,21 @@ omega_c = (
     cmax * Me / rho_dim * (pybamm.standard_parameters.t_plus + Mn / Me)
 )  # Migrative kinematic relationship coefficient
 gamma_dl_n = (
-    Cdl
-    * pybamm.standard_parameters.R
-    * pybamm.standard_parameters.T
-    * An
-    * pybamm.standard_parameters.Lx
-    / (pybamm.standard_parameters.F * ibar)
-    / (cmax * pybamm.standard_parameters.F * pybamm.standard_parameters.Lx / ibar)
+    Cdl * potential_scale / interfacial_current_scale_neg / time_scale
 )  # Dimensionless double-layer capacity (neg)
 gamma_dl_p = (
-    Cdl
-    * pybamm.standard_parameters.R
-    * pybamm.standard_parameters.T
-    * Ap
-    * pybamm.standard_parameters.Lx
-    / (pybamm.standard_parameters.F * ibar)
-    / (cmax * pybamm.standard_parameters.F * pybamm.standard_parameters.Lx / ibar)
+    Cdl * potential_scale / interfacial_current_scale_pos / time_scale
 )  # Dimensionless double-layer capacity (pos)
 voltage_cutoff = (
-    pybamm.standard_parameters.F
-    / (pybamm.standard_parameters.R * pybamm.standard_parameters.T)
+    1
+    / potential_scale
     * (pybamm.standard_parameters.voltage_low_cut - (U_PbO2_ref - U_Pb_ref))
 )  # Dimensionless voltage cut-off
-U_rxn = ibar / (cmax * pybamm.standard_parameters.F)  # Reaction velocity scale
 pi_os = (
-    mu_dim
-    * U_rxn
-    * pybamm.standard_parameters.Lx
-    / (d ** 2 * pybamm.standard_parameters.R * pybamm.standard_parameters.T * cmax)
+    mu_dim * velocity_scale * Lx / (d ** 2 * R * T * cmax)
 )  # Ratio of viscous pressure scale to osmotic pressure scale
 
+
 # Initial conditions
 q_init = pybamm.Parameter("q_init")  # Initial SOC [-]
 qmax = (
@@ -182,7 +195,7 @@ qmax = (
         + pybamm.standard_parameters.Ls * epss_max
         + pybamm.standard_parameters.Lp * epsp_max
     )
-    / pybamm.standard_parameters.Lx
+    / Lx
     / (sp - sn)
 )  # Dimensionless max capacity
 epsDelta_n = beta_surf_n / pybamm.standard_parameters.Ln * qmax
@@ -191,3 +204,9 @@ c_init = q_init
 epsn_init = epsn_max - epsDelta_n * (1 - q_init)  # Initial pororsity (neg) [-]
 epss_init = epss_max  # Initial pororsity (sep) [-]
 epsp_init = epsp_max - epsDelta_p * (1 - q_init)  # Initial pororsity (pos) [-]
+
+
+# Concatenated symbols
+s = pybamm.piecewise_constant(sn, 0, sp)
+beta_surf = pybamm.piecewise_constant(beta_surf_n, 0, beta_surf_p)
+eps_init = pybamm.piecewise_constant(epsn_init, epss_init, epsp_init)


Test patch:
diff --git a/tests/test_expression_tree/test_concatenations.py b/tests/test_expression_tree/test_concatenations.py
index 6043a316..1512db99 100644
--- a/tests/test_expression_tree/test_concatenations.py
+++ b/tests/test_expression_tree/test_concatenations.py
@@ -143,6 +143,17 @@ class TestConcatenations(unittest.TestCase):
         self.assertEqual(c.id, c_new.id)
         self.assertEqual(conc.id, pybamm.Concatenation(a_new, b_new, c_new).id)
 
+    def test_piecewise_constant(self):
+        a = pybamm.Scalar("a")
+        b = pybamm.Scalar("b")
+        c = pybamm.Scalar("c")
+        conc = pybamm.piecewise_constant(a, b, c)
+
+        self.assertIsInstance(conc, pybamm.Concatenation)
+        self.assertEqual(
+            conc.domain, ["negative electrode", "separator", "positive electrode"]
+        )
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/test_models/test_interface.py b/tests/test_models/test_interface.py
index 5c3a1db2..b83e2574 100644
--- a/tests/test_models/test_interface.py
+++ b/tests/test_models/test_interface.py
@@ -4,9 +4,10 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 import pybamm
-import tests.shared as shared
+from tests import shared
 
 import unittest
+import numpy as np
 
 
 class TestHomogeneousReaction(unittest.TestCase):
@@ -44,18 +45,20 @@ class TestHomogeneousReaction(unittest.TestCase):
 
         # test values
         ln = param.process_symbol(pybamm.standard_parameters.ln)
-        ln * 1 / ln
         lp = param.process_symbol(pybamm.standard_parameters.lp)
         self.assertEqual(processed_rxn.children[0].evaluate() * ln.evaluate(), 1)
         self.assertEqual(processed_rxn.children[2].evaluate() * lp.evaluate(), -1)
 
     def test_discretisation(self):
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        param = pybamm.ParameterValues(
+            "input/parameters/lithium-ion/parameters/LCO.csv"
+        )
+        mesh = shared.TestDefaults1DMacro().mesh
+        disc = pybamm.BaseDiscretisation(mesh)
 
         rxn = pybamm.interface.homogeneous_reaction()
 
-        param_rxn = defaults.param.process_symbol(rxn)
+        param_rxn = param.process_symbol(rxn)
         processed_rxn = disc.process_symbol(param_rxn)
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
@@ -66,6 +69,229 @@ class TestHomogeneousReaction(unittest.TestCase):
         self.assertEqual(processed_rxn.shape, combined_submeshes.nodes.shape)
 
 
+class TestButlerVolmerLeadAcid(unittest.TestCase):
+    def setUp(self):
+        self.cn = pybamm.Variable("concentration", domain=["negative electrode"])
+        self.cs = pybamm.Variable("concentration", domain=["separator"])
+        self.cp = pybamm.Variable("concentration", domain=["positive electrode"])
+        self.c = pybamm.Concatenation(self.cn, self.cs, self.cp)
+        self.phin = pybamm.Variable(
+            "potential difference", domain=["negative electrode"]
+        )
+        self.phis = pybamm.Variable("potential difference", domain=["separator"])
+        self.phip = pybamm.Variable(
+            "potential difference", domain=["positive electrode"]
+        )
+        self.phi = pybamm.Concatenation(self.phin, self.phis, self.phip)
+
+    def tearDown(self):
+        del self.cn
+        del self.cs
+        del self.cp
+        del self.c
+        del self.phin
+        del self.phis
+        del self.phip
+        del self.phi
+
+    def test_creation(self):
+        # negative electrode passes, returns Multiplication
+        bv = pybamm.interface.butler_volmer_lead_acid(self.cn, self.phin)
+        self.assertIsInstance(bv, pybamm.Multiplication)
+        self.assertEqual(bv.domain, ["negative electrode"])
+
+        # positive electrode passes, returns Multiplication
+        bv = pybamm.interface.butler_volmer_lead_acid(self.cp, self.phip)
+        self.assertIsInstance(bv, pybamm.Multiplication)
+        self.assertEqual(bv.domain, ["positive electrode"])
+
+        # whole cell domain passes, returns concatenation
+        bv = pybamm.interface.butler_volmer_lead_acid(self.c, self.phi)
+        self.assertIsInstance(bv, pybamm.Concatenation)
+        self.assertEqual(
+            bv.domain, ["negative electrode", "separator", "positive electrode"]
+        )
+
+        # c and phi without domain, domain gets input
+        c = pybamm.Variable("concentration", domain=[])
+        phi = pybamm.Variable("potential", domain=[])
+        bv = pybamm.interface.butler_volmer_lead_acid(
+            c, phi, domain=["negative electrode"]
+        )
+        self.assertIsInstance(bv, pybamm.Multiplication)
+        self.assertEqual(bv.domain, [])
+
+    def test_failure(self):
+        # domain that doesn't exist
+        with self.assertRaises(pybamm.DomainError):
+            pybamm.interface.butler_volmer_lead_acid(
+                None, None, domain=["not a domain"]
+            )
+        # no domain
+        c = pybamm.Variable("concentration", domain=[])
+        phi = pybamm.Variable("potential", domain=[])
+        with self.assertRaises(ValueError):
+            pybamm.interface.butler_volmer_lead_acid(c, phi)
+        # mismatched domains
+        with self.assertRaises(pybamm.DomainError):
+            pybamm.interface.butler_volmer_lead_acid(self.cn, self.phip)
+        # c, phi not concatenations with domain as whole domain
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        c = pybamm.Variable("concentration", domain=whole_cell)
+        phi = pybamm.Variable("potential", domain=whole_cell)
+        with self.assertRaises(ValueError):
+            pybamm.interface.butler_volmer_lead_acid(c, phi, whole_cell)
+
+    def test_set_parameters(self):
+        bv = pybamm.interface.butler_volmer_lead_acid(self.c, self.phi)
+        param = pybamm.ParameterValues(
+            "input/parameters/lead-acid/default.csv", {"current scale": 1}
+        )
+        proc_bv = param.process_symbol(bv)
+        [self.assertNotIsInstance(x, pybamm.Parameter) for x in proc_bv.pre_order()]
+
+    def test_discretisation(self):
+        bv_n = pybamm.interface.butler_volmer_lead_acid(self.cn, self.phin)
+        bv_p = pybamm.interface.butler_volmer_lead_acid(self.cp, self.phip)
+
+        # process parameters
+        param = pybamm.ParameterValues(
+            "input/parameters/lead-acid/default.csv", {"current scale": 1}
+        )
+        param_bv_n = param.process_symbol(bv_n)
+        param_bv_p = param.process_symbol(bv_p)
+
+        # discretise
+        mesh = shared.TestDefaults1DMacro().mesh
+        disc = pybamm.BaseDiscretisation(mesh)
+        y_slices = disc.get_variable_slices([self.cn, self.cp, self.phin, self.phip])
+        processed_bv_n = disc.process_symbol(param_bv_n, y_slices)
+        processed_bv_p = disc.process_symbol(param_bv_p, y_slices)
+
+        submesh = np.concatenate(
+            [mesh["negative electrode"].nodes, mesh["positive electrode"].nodes]
+        )
+        y = np.concatenate([submesh ** 2, submesh ** 3, -submesh])
+
+        # should evaluate to vectors with the right shape
+        self.assertEqual(
+            processed_bv_n.evaluate(None, y).shape,
+            mesh["negative electrode"].nodes.shape,
+        )
+        self.assertEqual(
+            processed_bv_p.evaluate(None, y).shape,
+            mesh["positive electrode"].nodes.shape,
+        )
+
+    @unittest.skip("needs smarter domain concatenation")
+    def test_discretisation_whole(self):
+        bv_whole = pybamm.interface.butler_volmer_lead_acid(self.c, self.phi)
+
+        # process parameters
+        param = pybamm.ParameterValues(
+            "input/parameters/lead-acid/default.csv", {"current scale": 1}
+        )
+        param_bv_whole = param.process_symbol(bv_whole)
+
+        # discretise
+        mesh = shared.TestDefaults1DMacro().mesh
+        disc = pybamm.BaseDiscretisation(mesh)
+        y_slices = disc.get_variable_slices([self.cn, self.cp, self.phin, self.phip])
+        processed_bv_whole = disc.process_symbol(param_bv_whole, y_slices)
+
+        # test
+        submesh = np.concatenate(
+            [mesh["negative electrode"].nodes, mesh["positive electrode"].nodes]
+        )
+        y = np.concatenate([submesh ** 2, submesh ** 3, -submesh])
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        combined_submeshes = disc.mesh.combine_submeshes(*whole_cell)
+        self.assertEqual(
+            processed_bv_whole.evaluate(None, y).shape, combined_submeshes.nodes.shape
+        )
+
+
+class TestExchangeCurrentDensity(unittest.TestCase):
+    def setUp(self):
+        self.cn = pybamm.Variable("concentration", domain=["negative electrode"])
+        self.cs = pybamm.Variable("concentration", domain=["separator"])
+        self.cp = pybamm.Variable("concentration", domain=["positive electrode"])
+        self.c = pybamm.Concatenation(self.cn, self.cs, self.cp)
+
+    def tearDown(self):
+        del self.cn
+        del self.cs
+        del self.cp
+        del self.c
+
+    def test_creation(self):
+        # Concentration without domain passes
+        c = pybamm.Variable("c")
+        pybamm.interface.exchange_current_density(c, ["negative electrode"])
+        pybamm.interface.exchange_current_density(c, ["positive electrode"])
+
+        # Concentration with correct domain passes
+        j0n = pybamm.interface.exchange_current_density(self.cn)
+        j0p = pybamm.interface.exchange_current_density(self.cp)
+        self.assertEqual(j0n.domain, ["negative electrode"])
+        self.assertEqual(j0p.domain, ["positive electrode"])
+
+        # Concentration with wrong domain fails
+        with self.assertRaises(pybamm.DomainError):
+            pybamm.interface.exchange_current_density(self.cp, ["negative electrode"])
+        with self.assertRaises(pybamm.DomainError):
+            pybamm.interface.exchange_current_density(self.cn, ["positive electrode"])
+        c = pybamm.Variable("concentration", domain=[])
+        with self.assertRaises(ValueError):
+            pybamm.interface.exchange_current_density(c)
+
+    def test_failure(self):
+        with self.assertRaises(pybamm.DomainError):
+            pybamm.interface.exchange_current_density(None, ["not a domain"])
+
+    def test_set_parameters(self):
+        param = pybamm.ParameterValues(
+            "input/parameters/lead-acid/default.csv", {"current scale": 1}
+        )
+        j0n = pybamm.interface.exchange_current_density(self.cn)
+        j0p = pybamm.interface.exchange_current_density(self.cp)
+        proc_j0n = param.process_symbol(j0n)
+        proc_j0p = param.process_symbol(j0p)
+        [self.assertNotIsInstance(x, pybamm.Parameter) for x in proc_j0n.pre_order()]
+        [self.assertNotIsInstance(x, pybamm.Parameter) for x in proc_j0p.pre_order()]
+
+    def test_discretisation(self):
+        # create exchange-current densities
+        j0n = pybamm.interface.exchange_current_density(self.cn)
+        j0p = pybamm.interface.exchange_current_density(self.cp)
+
+        # process parameters
+        param = pybamm.ParameterValues(
+            "input/parameters/lead-acid/default.csv", {"current scale": 1}
+        )
+        param_j0n = param.process_symbol(j0n)
+        param_j0p = param.process_symbol(j0p)
+
+        # discretise
+        mesh = shared.TestDefaults1DMacro().mesh
+        disc = pybamm.BaseDiscretisation(mesh)
+        y_slices = disc.get_variable_slices([self.cn, self.cp])
+        processed_j0n = disc.process_symbol(param_j0n, y_slices)
+        processed_j0p = disc.process_symbol(param_j0p, y_slices)
+
+        submesh = np.concatenate(
+            [mesh["negative electrode"].nodes, mesh["positive electrode"].nodes]
+        )
+        y = submesh ** 2
+        # should evaluate to vectors with the right shape
+        self.assertEqual(
+            processed_j0n.evaluate(y=y).shape, mesh["negative electrode"].nodes.shape
+        )
+        self.assertEqual(
+            processed_j0p.evaluate(y=y).shape, mesh["positive electrode"].nodes.shape
+        )
+
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys
diff --git a/tests/test_parameters/test_standard_parameters_lead_acid.py b/tests/test_parameters/test_standard_parameters_lead_acid.py
index 6b513d59..c680b557 100644
--- a/tests/test_parameters/test_standard_parameters_lead_acid.py
+++ b/tests/test_parameters/test_standard_parameters_lead_acid.py
@@ -4,6 +4,7 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 import pybamm
+from tests import shared
 
 import unittest
 
@@ -47,6 +48,30 @@ class TestStandardParametersLeadAcid(unittest.TestCase):
         # Excluded volume fraction should be less than 0.1
         self.assertLess(abs(param_eval["alpha"]), 1e-1)
 
+    def test_concatenated_parameters(self):
+        # create
+        s = pybamm.standard_parameters_lead_acid.s
+        self.assertIsInstance(s, pybamm.Concatenation)
+        self.assertEqual(
+            s.domain, ["negative electrode", "separator", "positive electrode"]
+        )
+
+        # process parameters and discretise
+        parameter_values = pybamm.ParameterValues(
+            "input/parameters/lead-acid/default.csv", {"current scale": 1}
+        )
+        mesh = shared.TestDefaults1DMacro().mesh
+        disc = pybamm.BaseDiscretisation(mesh)
+        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"
+        )
+        self.assertEqual(processed_s.shape, combined_submeshes.nodes.shape)
+
     @unittest.skip("lead acid functions not yet implemented")
     def test_functions_lead_acid(self):
         # Tests on how the parameters interact