Repository: tinosulzer/PyBaMM

Issue:
Allow custom current input into models
## Summary

Most battery models are "current-controlled", i.e. you feed in a current (function of time) and read the resulting voltage.

We should have functionality to allow the user to input their own current function into the model. 

Requires #68 

## Suggested change

1. Pass `current_scale` and `current_function` as parameters to the model (`current_scale` is needed as it is used in the non-dimensionalisation)
2. In `standard_parameters.py`, define `current` as a function that takes in `current_scale` and `current_function` and returns a `pybamm.Function` object:
```python3
def current(current_scale, current_function):
    t = pybamm.t
    return pybamm.Function(t, current_scale*current_function)
```
3. In the model class, call `current` as:
```python3
def __init__(self, current_scale, current_function):
    ...
    current = pybamm.standard_parameters(current_scale, current_function)
    ...
```
We could also have a `parameters.standard_current_functions.py` so that the user can choose from some pre-defined current functions, including interpolation of time series (e.g. from experimental data)

Discussion:

Patch:
diff --git a/docs/source/expression_tree/concatenations.rst b/docs/source/expression_tree/concatenations.rst
index affb5ead..b0bbd364 100644
--- a/docs/source/expression_tree/concatenations.rst
+++ b/docs/source/expression_tree/concatenations.rst
@@ -10,4 +10,5 @@ Concatenations
 .. autoclass:: pybamm.DomainConcatenation
   :members:
 
-.. autofunction:: pybamm.piecewise_constant
+.. autoclass:: pybamm.PiecewiseConstant
+  :members:
diff --git a/input/parameters/lead-acid/default.csv b/input/parameters/lead-acid/default.csv
index d0a2eb5c..bcb23d07 100644
--- a/input/parameters/lead-acid/default.csv
+++ b/input/parameters/lead-acid/default.csv
@@ -16,8 +16,8 @@ W,0.065,m,Electrode depth,Manufacturer,
 # Electrical,,,,,
 Q,17,A.h,Cell capacity,Manufacturer,
 voltage_cutoff_circuit,10.5,V,Voltage cutoff point for a 6-cell battery,-,
-n_electrodes_parallel,8,-,Number of electrodes connected in parallel to make a cell,Manufacturer,
-n_cells_series,6,-,Number of cells connected in series to make a battery,Manufacturer,
+number of electrodes in parallel,8,-,Number of electrodes connected in parallel to make a cell,Manufacturer,
+number of cells in series,6,-,Number of cells connected in series to make a battery,Manufacturer,
 q_init,1,-,Initial State of Charge,-,
 ,,,,,
 # Reactions,,,,,
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 81627d99..50eae1d8 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -63,7 +63,7 @@ from .expression_tree.concatenations import (
     Concatenation,
     NumpyModelConcatenation,
     DomainConcatenation,
-    piecewise_constant,
+    PiecewiseConstant,
 )
 from .expression_tree.array import Array
 from .expression_tree.matrix import Matrix
@@ -111,6 +111,7 @@ from .models.submodels import electrolyte, interface
 #
 from .meshes.meshes import KNOWN_DOMAINS  # need this for importing standard parameters
 from .parameters.parameter_values import ParameterValues
+from .parameters import standard_current_functions
 from .parameters import standard_parameters
 from .parameters import standard_parameters_lead_acid  # calls standard_parameters
 
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 8e1d0694..2e2f14ed 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -222,6 +222,10 @@ class Discretisation(object):
         elif isinstance(symbol, pybamm.BinaryOperator):
             return self.process_binary_operators(symbol)
 
+        elif isinstance(symbol, pybamm.Function):
+            new_child = self.process_symbol(symbol.children[0])
+            return pybamm.Function(symbol.func, new_child)
+
         elif isinstance(symbol, pybamm.UnaryOperator):
             new_child = self.process_symbol(symbol.children[0])
             return symbol.__class__(new_child)
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index db5acced..c17c7c4d 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -6,7 +6,7 @@ from __future__ import print_function, unicode_literals
 import pybamm
 
 import numpy as np
-import numbers
+import warnings
 
 
 class Concatenation(pybamm.Symbol):
@@ -106,7 +106,7 @@ class DomainConcatenation(Concatenation):
         children = list(children)
 
         for i, child in enumerate(children):
-            if child.is_constant():
+            if child.is_constant() or child.evaluates_to_number():
                 children[i] = self.process_node_for_concatenate(child, mesh)
 
         # Allow the base class to sort the domains into the correct order
@@ -144,25 +144,22 @@ class DomainConcatenation(Concatenation):
             the sub-expression to process (node.is_constant() is true)
 
         """
-
-        # node must be constant
-        value = node.evaluate()
-
-        # correct size of vector should be number of points in the domains
-        subvector_size = sum([mesh[dom].npts for dom in node.domain])
-
-        # check if its a scalar, if so convert to vector
-        if isinstance(value, numbers.Number):
-            value = np.full(subvector_size, value)
-
-        # check it is the right size
-        if value.size != subvector_size:
-            raise ValueError(
-                "Error: expression evaluated to a vector of incorrect length"
-            )
-
-        # convert to a Vector node
-        return pybamm.Vector(value, domain=node.domain)
+        try:
+            node_size = node.size
+        except AttributeError:
+            node_size = 0
+
+        if node_size > 1:
+            if node_size == sum([mesh[dom].npts for dom in node.domain]):
+                return node
+            else:
+                raise ValueError(
+                    "Error: expression evaluated to a vector of incorrect length"
+                )
+        else:
+            warnings.warn("Implementation not robust")
+            npts = {dom: submesh.npts for dom, submesh in mesh.items()}
+            return pybamm.NumpyBroadcast(node, node.domain, npts)
 
     def evaluate(self, t=None, y=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
@@ -179,23 +176,29 @@ class DomainConcatenation(Concatenation):
         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
+class PiecewiseConstant(Concatenation):
+    """Piecewise constant concatenation of three symbols.
+    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
+    neg_value: :class:`numbers.Number` or :class:`pybamm.Symbol`
+        The value in the negative electrode
+    sep_value: :class:`numbers.Number` or :class:`pybamm.Symbol`
+        The value in the separator
+    pos_value: :class:`numbers.Number` or :class:`pybamm.Symbol`
+        The value in the positive electrode
 
-    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
-    )
+
+    def __init__(self, neg_value, sep_value, pos_value):
+        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"]
+        )
+        super().__init__(
+            neg_value_with_domain, sep_value_with_domain, pos_value_with_domain
+        )
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 5b850c5e..8a7eeff8 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -334,7 +334,7 @@ class Symbol(anytree.NodeMixin):
         """
         try:
             # return true if node evaluates to a number
-            return isinstance(self.evaluate(), numbers.Number)
+            return isinstance(self.evaluate(t=0), numbers.Number)
         except NotImplementedError:
             # return false if NotImplementedError is raised
             # (there is a e.g. Parameter, Variable, ... in the tree)
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 2a4d87b8..cf42c215 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -147,9 +147,12 @@ class Broadcast(SpatialOperator):
     """
 
     def __init__(self, child, domain, name=None):
-        if child.domain != []:
+        if child.domain not in [[], domain]:
             raise pybamm.DomainError(
-                """Domain of a broadcasted child must be [] but is '{}'""".format(
+                """
+                Domain of a broadcasted child must be []
+                or same as 'domain' but is '{}'
+                """.format(
                     child.domain
                 )
             )
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index b81d42ac..d34e6a0d 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -6,6 +6,7 @@ from __future__ import print_function, unicode_literals
 import pybamm
 
 import numbers
+import os
 
 
 class BaseModel(object):
@@ -50,7 +51,17 @@ class BaseModel(object):
 
         # Default parameter values, discretisation and solver
         self.default_parameter_values = pybamm.ParameterValues(
-            "input/parameters/lithium-ion/parameters/LCO.csv"
+            "input/parameters/lithium-ion/parameters/LCO.csv",
+            {
+                "I_typ": 1,
+                "current function": os.path.join(
+                    os.getcwd(),
+                    "pybamm",
+                    "parameters",
+                    "standard_current_functions",
+                    "constant_current.py",
+                ),
+            },
         )
 
         self.default_geometry = pybamm.Geometry1DMacro()
diff --git a/pybamm/models/lead_acid/loqs.py b/pybamm/models/lead_acid/loqs.py
index aabb3bba..74047e0f 100644
--- a/pybamm/models/lead_acid/loqs.py
+++ b/pybamm/models/lead_acid/loqs.py
@@ -5,6 +5,8 @@ from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 import pybamm
 
+import os
+
 
 class LOQS(pybamm.BaseModel):
     """Leading-Order Quasi-Static model for lead-acid.
@@ -37,7 +39,9 @@ class LOQS(pybamm.BaseModel):
         epss = pybamm.Variable("epss", domain=[])
         epsp = pybamm.Variable("epsp", domain=[])
 
-        # Parameters
+        # Current function
+        icell = pybamm.standard_parameters_lead_acid.current_with_time
+        # Parameters and functions
         ln = pybamm.standard_parameters.ln
         ls = pybamm.standard_parameters.ls
         lp = pybamm.standard_parameters.lp
@@ -56,14 +60,14 @@ class LOQS(pybamm.BaseModel):
         epsp_init = pybamm.standard_parameters_lead_acid.epsp_init
 
         # ODEs
-        jn = pybamm.standard_parameters_lead_acid.icell / ln
-        jp = -pybamm.standard_parameters_lead_acid.icell / lp
+        jn = icell / ln
+        jp = -icell / lp
         depsndt = -beta_surf_n * jn
         depspdt = -beta_surf_p * jp
         dcdt = (
             1
             / (ln * epsn + ls * epss + lp * epsp)
-            * ((sn - sp) - c * (ln * depsndt + lp * depspdt))
+            * ((sn - sp) * icell - c * (ln * depsndt + lp * depspdt))
         )
         self.rhs = {c: dcdt, epsn: depsndt, epss: pybamm.Scalar(0), epsp: depspdt}
         # Initial conditions
@@ -94,5 +98,15 @@ class LOQS(pybamm.BaseModel):
 
         # Overwrite default parameter values
         self.default_parameter_values = pybamm.ParameterValues(
-            "input/parameters/lead-acid/default.csv", {"current scale": 1}
+            "input/parameters/lead-acid/default.csv",
+            {
+                "current scale": 1,
+                "current function": os.path.join(
+                    os.getcwd(),
+                    "pybamm",
+                    "parameters",
+                    "standard_current_functions",
+                    "constant_current.py",
+                ),
+            },
         )
diff --git a/pybamm/models/reaction_diffusion.py b/pybamm/models/reaction_diffusion.py
index e465032d..a0784cbb 100644
--- a/pybamm/models/reaction_diffusion.py
+++ b/pybamm/models/reaction_diffusion.py
@@ -29,8 +29,10 @@ class ReactionDiffusionModel(pybamm.BaseModel):
 
     def __init__(self):
         super().__init__()
+        # Define current
+        current = pybamm.standard_parameters.current_with_time
         # Load reaction flux from submodels
-        G = pybamm.interface.homogeneous_reaction()
+        G = pybamm.interface.homogeneous_reaction(current)
         # Load diffusion model from submodels
         diffusion_model = pybamm.electrolyte.StefanMaxwellDiffusion(G)
 
diff --git a/pybamm/models/submodels/interface.py b/pybamm/models/submodels/interface.py
index 5a32ac42..77dd3226 100644
--- a/pybamm/models/submodels/interface.py
+++ b/pybamm/models/submodels/interface.py
@@ -6,19 +6,14 @@ from __future__ import print_function, unicode_literals
 import pybamm
 
 
-def homogeneous_reaction():
+def homogeneous_reaction(current):
     """
     Homogeneous reaction at the electrode-electrolyte interface
     """
-
-    current_neg = (
-        pybamm.Scalar(1, domain=["negative electrode"]) / pybamm.standard_parameters.ln
-    )
-    current_sep = pybamm.Scalar(0, domain=["separator"])
-    current_pos = (
-        -pybamm.Scalar(1, domain=["positive electrode"]) / pybamm.standard_parameters.lp
-    )
-    return pybamm.Concatenation(current_neg, current_sep, current_pos)
+    # hack to make the concatenation work. Concatenation needs some work
+    current_neg = current / pybamm.standard_parameters.ln
+    current_pos = -current / pybamm.standard_parameters.lp
+    return pybamm.PiecewiseConstant(current_neg, 0, current_pos)
 
 
 def butler_volmer_lead_acid(c, phi, domain=None):
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index a32f1482..568f1b60 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -161,6 +161,10 @@ class ParameterValues(dict):
             new_child = self.process_symbol(symbol.children[0])
             return pybamm.Broadcast(new_child, symbol.domain)
 
+        elif isinstance(symbol, pybamm.Function):
+            new_child = self.process_symbol(symbol.children[0])
+            return pybamm.Function(symbol.func, new_child)
+
         elif isinstance(symbol, pybamm.UnaryOperator):
             new_child = self.process_symbol(symbol.children[0])
             return symbol.__class__(new_child)
diff --git a/pybamm/parameters/standard_current_functions/__init__.py b/pybamm/parameters/standard_current_functions/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/pybamm/parameters/standard_current_functions/constant_current.py b/pybamm/parameters/standard_current_functions/constant_current.py
new file mode 100644
index 00000000..4320d0a0
--- /dev/null
+++ b/pybamm/parameters/standard_current_functions/constant_current.py
@@ -0,0 +1,4 @@
+def constant_current(t):
+    """Constant current function"""
+    # output has to have same shape and type as t
+    return 0 * t + 1
diff --git a/pybamm/parameters/standard_parameters.py b/pybamm/parameters/standard_parameters.py
index 8ef94c66..b335cb90 100644
--- a/pybamm/parameters/standard_parameters.py
+++ b/pybamm/parameters/standard_parameters.py
@@ -2,7 +2,7 @@
 # Standard parameters for battery models
 #
 """
-Standard Parameters for battery models
+Standard pybamm.Parameters for battery models
 
 Physical Constants
 ------------------
@@ -70,58 +70,60 @@ cn0_dimensional, cp0_dimensional
 """
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
-from pybamm import Parameter
+import pybamm
 
 # --------------------------------------------------------------------------------------
 """Dimensional Parameters"""
 
 # Physical Constants
-R = Parameter("R")
-F = Parameter("F")
-T = Parameter("T_ref")
+R = pybamm.Parameter("R")
+F = pybamm.Parameter("F")
+T = pybamm.Parameter("T_ref")
 
 # Macroscale Geometry
-Ln = Parameter("Ln")
-Ls = Parameter("Ls")
-Lp = Parameter("Lp")
+Ln = pybamm.Parameter("Ln")
+Ls = pybamm.Parameter("Ls")
+Lp = pybamm.Parameter("Lp")
 Lx = Ln + Ls + Lp
 
 # 3D Geometry
-Ly = Parameter("Ly")
-Lz = Parameter("Lz")
+Ly = pybamm.Parameter("Ly")
+Lz = pybamm.Parameter("Lz")
 
 # Microscale Geometry
-R_n = Parameter("R_n")
-R_p = Parameter("R_p")
-a_n = Parameter("a_n")
-a_p = Parameter("a_p")
+R_n = pybamm.Parameter("R_n")
+R_p = pybamm.Parameter("R_p")
+a_n = pybamm.Parameter("a_n")
+a_p = pybamm.Parameter("a_p")
 
 # Electrolyte Properties
-ce_typ = Parameter("ce_typ")  # Typical lithium ion concentration in electrolyte
-De_typ = Parameter("De_typ")  # Typical electrolyte diffusivity
+ce_typ = pybamm.Parameter("ce_typ")  # Typical lithium ion concentration in electrolyte
+De_typ = pybamm.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
-Dn_typ = Parameter("Dn_typ")  # Typical negative particle diffusivity
-Dp_typ = Parameter("Dp_typ")  # Typical positive particle diffusivity
+sigma_n = pybamm.Parameter("sigma_n")  # Conductivity in negative electrode
+sigma_p = pybamm.Parameter("sigma_p")  # Conductivity in positive electrode
+cn_max = pybamm.Parameter("cn_max")  # Max concentration in negative electrode
+cp_max = pybamm.Parameter("cp_max")  # Max concentration in positive electrode
+Dn_typ = pybamm.Parameter("Dn_typ")  # Typical negative particle diffusivity
+Dp_typ = pybamm.Parameter("Dp_typ")  # Typical positive particle diffusivity
 
 # Electrochemical Reactions
-m_n_dim = Parameter("m_n_dim")  # Reaction rate in negative electrode
-m_p_dim = Parameter("m_p_dim")  # Reaction rate in positive electrode
+m_n_dim = pybamm.Parameter("m_n_dim")  # Reaction rate in negative electrode
+m_p_dim = pybamm.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
+voltage_low_cut = pybamm.Parameter("voltage_low_cut")  # Lower voltage cut-off
+voltage_high_cut = pybamm.Parameter("voltage_high_cut")  # Upper voltage cut-off
+I_typ = pybamm.Parameter("I_typ")  # Typical current density
+Phi_typ = pybamm.Parameter("Phi_typ")  # Typical voltage drop
+current_with_time = pybamm.FunctionParameter("current function", pybamm.t)
+dimensional_current_with_time = I_typ * current_with_time
 
 # 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
+ce0_dimensional = pybamm.Parameter("ce0")  # Initial li ion concentration in electrolyte
+cn0_dimensional = pybamm.Parameter("cn0")  # Initial li concentration in neg electrode
+cp0_dimensional = pybamm.Parameter("cp0")  # Initial li concentration in pos electrode
 
 # --------------------------------------------------------------------------------------
 """Dimensionless Parameters"""
@@ -134,10 +136,14 @@ ly = Ly / Lz
 lz = Lz / Lz
 
 # 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
+epsilon_n = pybamm.Parameter(
+    "epsilon_n"
+)  # Electrolyte volume fraction in neg electrode
+epsilon_s = pybamm.Parameter("epsilon_s")  # Electrolyte volume fraction in separator
+epsilon_p = pybamm.Parameter(
+    "epsilon_p"
+)  # Electrolyte volume fraction in pos electrode
+b = pybamm.Parameter("b")  # Bruggeman coefficient
 beta_n = a_n * R_n
 beta_p = a_p * R_p
 
@@ -165,7 +171,7 @@ m_n = tau_d / tau_r_n
 m_p = tau_d / tau_r_n
 
 # Electrolyte Properties
-t_plus = Parameter("t_plus")  # cation transference number
+t_plus = pybamm.Parameter("t_plus")  # cation transference number
 delta = (Lx ** 2 / De_typ) * (I_typ / (F * cn_max * Lx))
 nu = cn_max / ce_typ
 
diff --git a/pybamm/parameters/standard_parameters_lead_acid.py b/pybamm/parameters/standard_parameters_lead_acid.py
index 455ba324..c80a0d3d 100644
--- a/pybamm/parameters/standard_parameters_lead_acid.py
+++ b/pybamm/parameters/standard_parameters_lead_acid.py
@@ -30,10 +30,13 @@ H = pybamm.Parameter("H")
 W = pybamm.Parameter("W")
 A_cs = H * W  # Area of the current collectors [m2]
 Ibar = pybamm.Parameter("current scale")
-ibar = Ibar / (8 * A_cs)  # Specified scale for the current [A.m-2]
+n_electrodes_parallel = pybamm.Parameter("number of electrodes in parallel")
+ibar = Ibar / (n_electrodes_parallel * A_cs)  # Specified scale for the current [A.m-2]
 Q = 17  # Capacity [Ah]
 Crate = Ibar / Q  # C-rate [-]
-icell = Ibar
+current_with_time = pybamm.FunctionParameter("current function", pybamm.t)
+dimensional_current_with_time = ibar * current_with_time
+
 
 # Microstructure
 An = pybamm.Parameter(
@@ -91,6 +94,22 @@ 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]
 
+# --------------------------------------------------------------------------------------
+"""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
+
 # --------------------------------------------------------------------------------------
 """Functions"""
 
@@ -105,7 +124,7 @@ def D_dim(c):
 
 def D(c):
     "Dimensionless electrolyte diffusivity"
-    return D_dim(c * cmax) / D_dim(cmax)
+    return D_dim(c * concentration_scale) / D_dim(concentration_scale)
 
 
 def U_Pb(c):
@@ -118,22 +137,6 @@ def U_PbO2(c):
     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"""
 
@@ -207,6 +210,6 @@ 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)
+s = pybamm.PiecewiseConstant(sn, 0, sp)
+beta_surf = pybamm.PiecewiseConstant(beta_surf_n, 0, beta_surf_p)
+eps_init = pybamm.PiecewiseConstant(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 21cca07d..3f929e58 100644
--- a/tests/test_expression_tree/test_concatenations.py
+++ b/tests/test_expression_tree/test_concatenations.py
@@ -144,15 +144,21 @@ class TestConcatenations(unittest.TestCase):
         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)
+        a = pybamm.Scalar(1)
+        b = pybamm.Scalar(2)
+        c = pybamm.Scalar(3)
+        conc = pybamm.PiecewiseConstant(a, b, c)
 
         self.assertIsInstance(conc, pybamm.Concatenation)
         self.assertEqual(
             conc.domain, ["negative electrode", "separator", "positive electrode"]
         )
+        self.assertEqual(conc.children[0].domain, ["negative electrode"])
+        self.assertEqual(conc.children[1].domain, ["separator"])
+        self.assertEqual(conc.children[2].domain, ["positive electrode"])
+        self.assertEqual(conc.children[0].evaluate(), 1)
+        self.assertEqual(conc.children[1].evaluate(), 2)
+        self.assertEqual(conc.children[2].evaluate(), 3)
 
 
 if __name__ == "__main__":
diff --git a/tests/test_expression_tree/test_symbol.py b/tests/test_expression_tree/test_symbol.py
index 85b145ed..e35c0749 100644
--- a/tests/test_expression_tree/test_symbol.py
+++ b/tests/test_expression_tree/test_symbol.py
@@ -158,6 +158,10 @@ class TestSymbol(unittest.TestCase):
         a = pybamm.Matrix(np.ones((4, 6)))
         self.assertFalse(a.evaluates_to_number())
 
+        # Time variable returns true
+        a = 3 * pybamm.t + 2
+        self.assertTrue(a.evaluates_to_number())
+
     def test_symbol_repr(self):
         """
         test that __repr___ returns the string
diff --git a/tests/test_models/test_interface.py b/tests/test_models/test_interface.py
index 2b36824d..73c3862c 100644
--- a/tests/test_models/test_interface.py
+++ b/tests/test_models/test_interface.py
@@ -15,7 +15,7 @@ class TestHomogeneousReaction(unittest.TestCase):
         param = pybamm.ParameterValues(
             "input/parameters/lithium-ion/parameters/LCO.csv"
         )
-        rxn = pybamm.interface.homogeneous_reaction()
+        rxn = pybamm.interface.homogeneous_reaction(1)
 
         processed_rxn = param.process_symbol(rxn)
 
@@ -38,7 +38,6 @@ class TestHomogeneousReaction(unittest.TestCase):
                 ]
             )
         )
-        self.assertIsInstance(processed_rxn.children[1], pybamm.Scalar)
         self.assertEqual(processed_rxn.children[0].domain, ["negative electrode"])
         self.assertEqual(processed_rxn.children[1].domain, ["separator"])
         self.assertEqual(processed_rxn.children[2].domain, ["positive electrode"])
@@ -47,6 +46,7 @@ class TestHomogeneousReaction(unittest.TestCase):
         ln = param.process_symbol(pybamm.standard_parameters.ln)
         lp = param.process_symbol(pybamm.standard_parameters.lp)
         self.assertEqual(processed_rxn.children[0].evaluate() * ln.evaluate(), 1)
+        self.assertEqual(processed_rxn.children[1].evaluate(), 0)
         self.assertEqual(processed_rxn.children[2].evaluate() * lp.evaluate(), -1)
 
     def test_discretisation(self):
@@ -56,7 +56,7 @@ class TestHomogeneousReaction(unittest.TestCase):
         defaults = shared.TestDefaults1DMacro()
         disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
 
-        rxn = pybamm.interface.homogeneous_reaction()
+        rxn = pybamm.interface.homogeneous_reaction(1)
 
         param_rxn = param.process_symbol(rxn)
         processed_rxn = disc.process_symbol(param_rxn)
diff --git a/tests/test_parameters/test_standard_parameters.py b/tests/test_parameters/test_standard_parameters.py
index d45ae623..d264a7f0 100644
--- a/tests/test_parameters/test_standard_parameters.py
+++ b/tests/test_parameters/test_standard_parameters.py
@@ -3,6 +3,7 @@
 #
 import pybamm
 
+import os
 import unittest
 
 
@@ -30,6 +31,31 @@ class TestStandardParameters(unittest.TestCase):
         lp_eval = parameter_values.process_symbol(lp)
         self.assertAlmostEqual((ln_eval + ls_eval + lp_eval).evaluate(), 1)
 
+    def test_functions(self):
+        # create current functions
+        current = pybamm.standard_parameters.dimensional_current_with_time
+        dimensionless_current = pybamm.standard_parameters.current_with_time
+
+        # process
+        parameter_values = pybamm.ParameterValues(
+            {
+                "I_typ": 2,
+                "current function": os.path.join(
+                    os.getcwd(),
+                    "pybamm",
+                    "parameters",
+                    "standard_current_functions",
+                    "constant_current.py",
+                ),
+            }
+        )
+        current_eval = parameter_values.process_symbol(current)
+        dimensionless_current_eval = parameter_values.process_symbol(
+            dimensionless_current
+        )
+        self.assertEqual(current_eval.evaluate(t=3), 2)
+        self.assertEqual(dimensionless_current_eval.evaluate(t=3), 1)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/test_parameters/test_standard_parameters_lead_acid.py b/tests/test_parameters/test_standard_parameters_lead_acid.py
index 39363eb1..6ac13b08 100644
--- a/tests/test_parameters/test_standard_parameters_lead_acid.py
+++ b/tests/test_parameters/test_standard_parameters_lead_acid.py
@@ -6,6 +6,7 @@ from __future__ import print_function, unicode_literals
 import pybamm
 from tests import shared
 
+import os
 import unittest
 
 
@@ -72,6 +73,34 @@ class TestStandardParametersLeadAcid(unittest.TestCase):
         )
         self.assertEqual(processed_s.shape, combined_submeshes.nodes.shape)
 
+    def test_current_functions(self):
+        # create current functions
+        current = pybamm.standard_parameters_lead_acid.dimensional_current_with_time
+        dimensionless_current = pybamm.standard_parameters_lead_acid.current_with_time
+
+        # process
+        parameter_values = pybamm.ParameterValues(
+            {
+                "H": 0.1,
+                "W": 0.1,
+                "number of electrodes in parallel": 8,
+                "current scale": 2,
+                "current function": os.path.join(
+                    os.getcwd(),
+                    "pybamm",
+                    "parameters",
+                    "standard_current_functions",
+                    "constant_current.py",
+                ),
+            }
+        )
+        current_eval = parameter_values.process_symbol(current)
+        dimensionless_current_eval = parameter_values.process_symbol(
+            dimensionless_current
+        )
+        self.assertAlmostEqual(current_eval.evaluate(t=3), 2 / (8 * 0.1 * 0.1))
+        self.assertEqual(dimensionless_current_eval.evaluate(t=3), 1)
+
     @unittest.skip("lead acid functions not yet implemented")
     def test_functions_lead_acid(self):
         # Tests on how the parameters interact