From f4470d758e410cbc1f8be2d60cfe0372011aaa44 Mon Sep 17 00:00:00 2001
From: tom <t.g.tranter@gmail.com>
Date: Wed, 18 Sep 2019 11:30:12 +0100
Subject: [PATCH 01/12] add r_average

---
 pybamm/expression_tree/unary_operators.py | 35 +++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index fd3be4c755..910dd8c3d6 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -984,3 +984,38 @@ def boundary_value(symbol, side):
     # Otherwise, calculate boundary value
     else:
         return BoundaryValue(symbol, side)
+
+
+def r_average(symbol):
+    """convenience function for creating an average in the r-direction
+
+    Parameters
+    ----------
+    symbol : :class:`pybamm.Symbol`
+        The function to be averaged
+
+    Returns
+    -------
+    :class:`Symbol`
+        the new averaged symbol
+    """
+    # Symbol must have domain [] or ["current collector"]
+    if "particle" not in symbol.domain and symbol.domain is not []:
+        raise pybamm.DomainError(
+            """r-average only implemented in the 'particle' domain,
+            but symbol has domains {}""".format(
+                symbol.domain
+            )
+        )
+    # If symbol doesn't have a domain, its average value is itself
+    elif symbol.domain is []:
+        new_symbol = symbol.new_copy()
+        new_symbol.parent = None
+        return new_symbol
+    else:
+        r = pybamm.SpatialVariable("r", [symbol.domain])
+        if symbol.domain == ["negative particle"]:
+            l_r = pybamm.geometric_parameters.R_n
+        elif symbol.domain == ["positive particle"]:
+            l_r = pybamm.geometric_parameters.R_p
+        return Integral(symbol, r) / l_r

From b2bba6794d64c75964642d48d51df81caf2af7d9 Mon Sep 17 00:00:00 2001
From: tom <t.g.tranter@gmail.com>
Date: Thu, 19 Sep 2019 11:51:38 +0100
Subject: [PATCH 02/12] update r_average and add test

---
 pybamm/__init__.py                            |  1 +
 pybamm/expression_tree/unary_operators.py     | 17 ++++++------
 .../submodels/particle/base_particle.py       |  8 +++++-
 .../test_unary_operators.py                   | 27 ++++++++++++++++++-
 4 files changed, 43 insertions(+), 10 deletions(-)

diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index b72c8cebc8..e468950c8f 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -108,6 +108,7 @@ def version(formatted=False):
     z_average,
     yz_average,
     boundary_value,
+    r_average
 )
 from .expression_tree.functions import *
 from .expression_tree.parameter import Parameter, FunctionParameter
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 910dd8c3d6..2f75aa3428 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -999,8 +999,9 @@ def r_average(symbol):
     :class:`Symbol`
         the new averaged symbol
     """
+    ok_domains = [['positive particle'], ['negative particle'], []]
     # Symbol must have domain [] or ["current collector"]
-    if "particle" not in symbol.domain and symbol.domain is not []:
+    if symbol.domain not in ok_domains:
         raise pybamm.DomainError(
             """r-average only implemented in the 'particle' domain,
             but symbol has domains {}""".format(
@@ -1008,14 +1009,14 @@ def r_average(symbol):
             )
         )
     # If symbol doesn't have a domain, its average value is itself
-    elif symbol.domain is []:
+    elif symbol.domain == []:
         new_symbol = symbol.new_copy()
         new_symbol.parent = None
         return new_symbol
+    # If symbol is a Broadcast, its average value is its child
+    elif isinstance(symbol, pybamm.Broadcast):
+        return symbol.orphans[0]
     else:
-        r = pybamm.SpatialVariable("r", [symbol.domain])
-        if symbol.domain == ["negative particle"]:
-            l_r = pybamm.geometric_parameters.R_n
-        elif symbol.domain == ["positive particle"]:
-            l_r = pybamm.geometric_parameters.R_p
-        return Integral(symbol, r) / l_r
+        r = pybamm.SpatialVariable("r", symbol.domain)
+        v = pybamm.Broadcast(pybamm.Scalar(1), symbol.domain)
+        return Integral(symbol, r) / Integral(v, r)
diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py
index e454405b55..05b2bf8209 100644
--- a/pybamm/models/submodels/particle/base_particle.py
+++ b/pybamm/models/submodels/particle/base_particle.py
@@ -31,7 +31,9 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             c_scale = self.param.c_n_max
         elif self.domain == "Positive":
             c_scale = self.param.c_p_max
-
+        rad = pybamm.SpatialVariable("r", [self.domain.lower() + " particle"])
+#        c_s_r_av = pybamm.r_average(pybamm.x_average(c_s))
+        c_s_r_av = pybamm.r_average(c_s_xav)
         variables = {
             self.domain + " particle concentration": c_s,
             self.domain + " particle concentration [mol.m-3]": c_s * c_scale,
@@ -39,6 +41,10 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             "X-averaged "
             + self.domain.lower()
             + " particle concentration [mol.m-3]": c_s_xav * c_scale,
+            "RX-averaged " + self.domain.lower() + " particle concentration": c_s_r_av,
+            "RX-averaged "
+            + self.domain.lower()
+            + " particle concentration [mol.m-3]": c_s_r_av * c_scale,
             self.domain + " particle surface concentration": c_s_surf,
             self.domain
             + " particle surface concentration [mol.m-3]": c_scale * c_s_surf,
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 5308eebe3f..9542257de9 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -227,7 +227,7 @@ def test_boundary_value(self):
         self.assertEqual(boundary_a_tert.domain, ["current collector"])
         self.assertEqual(boundary_a_tert.auxiliary_domains, {"secondary": ["bla"]})
 
-    def test_average(self):
+    def test_x_average(self):
         a = pybamm.Scalar(1)
         average_a = pybamm.x_average(a)
         self.assertEqual(average_a.id, a.id)
@@ -262,6 +262,31 @@ def test_average(self):
         with self.assertRaises(pybamm.DomainError):
             pybamm.x_average(a)
 
+    def test_r_average(self):
+        a = pybamm.Scalar(1)
+        average_a = pybamm.r_average(a)
+        self.assertEqual(average_a.id, a.id)
+
+        average_broad_a = pybamm.r_average(pybamm.Broadcast(a, ["negative particle"]))
+        self.assertEqual(average_broad_a.evaluate(), np.array([1]))
+
+        for domain in [
+            ["negative particle"],
+            ["positive particle"]
+        ]:
+            a = pybamm.Symbol("a", domain=domain)
+            r = pybamm.SpatialVariable("r", domain)
+            av_a = pybamm.r_average(a)
+            self.assertIsInstance(av_a, pybamm.Division)
+            self.assertIsInstance(av_a.children[0], pybamm.Integral)
+            self.assertEqual(av_a.children[0].integration_variable[0].domain, r.domain)
+            # electrode domains go to current collector when averaged
+            self.assertEqual(av_a.domain, [])
+
+        a = pybamm.Symbol("a", domain="bad domain")
+        with self.assertRaises(pybamm.DomainError):
+            pybamm.x_average(a)
+
     def test_yz_average(self):
         a = pybamm.Scalar(1)
         z_average_a = pybamm.z_average(a)

From 3ac749437e305311a3aadd5a068ec2fb2bff00b8 Mon Sep 17 00:00:00 2001
From: tom <t.g.tranter@gmail.com>
Date: Fri, 20 Sep 2019 14:34:50 +0100
Subject: [PATCH 03/12] add example SOC script

---
 examples/scripts/SPMe_SOC.py | 78 ++++++++++++++++++++++++++++++++++++
 1 file changed, 78 insertions(+)
 create mode 100644 examples/scripts/SPMe_SOC.py

diff --git a/examples/scripts/SPMe_SOC.py b/examples/scripts/SPMe_SOC.py
new file mode 100644
index 0000000000..24ed89ed13
--- /dev/null
+++ b/examples/scripts/SPMe_SOC.py
@@ -0,0 +1,78 @@
+import pybamm
+import numpy as np
+import matplotlib.pyplot as plt
+plt.close('all')
+pybamm.set_logging_level("INFO")
+
+# load model
+model = pybamm.lithium_ion.SPMe()
+
+# create geometry
+geometry = model.default_geometry
+
+# load parameter values and process model and geometry
+param = model.default_parameter_values
+
+param.update({'Negative electrode surface area density [m-1]': 180000.0,
+              'Positive electrode surface area density [m-1]': 150000.0})
+
+
+max_neg = param['Maximum concentration in negative electrode [mol.m-3]']
+max_pos = param['Maximum concentration in positive electrode [mol.m-3]']
+a_neg = param['Negative electrode surface area density [m-1]']
+a_pos = param['Positive electrode surface area density [m-1]']
+l_neg = param['Negative electrode thickness [m]']
+l_pos = param['Positive electrode thickness [m]']
+por_neg = param['Negative electrode porosity']
+por_pos = param['Positive electrode porosity']
+
+param.process_model(model)
+param.process_geometry(geometry)
+
+s_var = pybamm.standard_spatial_vars
+var_pts = {s_var.x_n: 5, s_var.x_s: 5,
+           s_var.x_p: 5, s_var.r_n: 5,
+           s_var.r_p: 10}
+# set mesh
+mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
+
+# discretise model
+disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+disc.process_model(model)
+
+# solve model
+t_eval = np.linspace(0, 0.2, 100)
+sol = model.default_solver.solve(model, t_eval)
+
+# plot
+out_vars = ['RX-averaged positive particle concentration [mol.m-3]',
+            'RX-averaged negative particle concentration [mol.m-3]',
+            'X-averaged electrolyte concentration [mol.m-3]',
+            'X-averaged positive electrode porosity',
+            'X-averaged negative electrode porosity',
+            'X-averaged separator porosity']
+plot = pybamm.QuickPlot(model, mesh, sol, output_variables=out_vars)
+plot.dynamic_plot()
+keys = list(model.variables.keys())
+keys.sort()
+
+xppc = pybamm.ProcessedVariable(model.variables[out_vars[0]],
+                                sol.t, sol.y, mesh=mesh)
+xnpc = pybamm.ProcessedVariable(model.variables[out_vars[1]],
+                                sol.t, sol.y, mesh=mesh)
+xec = pybamm.ProcessedVariable(model.variables[out_vars[2]],
+                               sol.t, sol.y, mesh=mesh)
+rp = np.linspace(0, 1.0, 11)
+
+plt.figure()
+rplt = 0.0
+plt.plot(np.ones(len(sol.t))*max_neg*por_neg, 'r--', label='Max Neg')
+plt.plot(xnpc(sol.t, r=rplt)*por_neg, 'r', label='Neg Li')
+plt.plot(np.ones(len(sol.t))*max_pos*por_pos, 'b--', label='Max Pos')
+plt.plot(xppc(sol.t, r=rplt)*por_pos, 'b', label='Pos Li')
+plt.plot(xec(sol.t, r=rplt), 'g', label='Elec Li')
+tot = (xnpc(sol.t, r=rplt) * por_neg +
+       xppc(sol.t, r=rplt) * por_pos +
+       xec(sol.t, r=rplt))
+plt.plot(tot, 'k-', label='Total Li')
+plt.legend()

From 4817075d4cd3abb29f3420408d114f795759a4b3 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Sat, 21 Sep 2019 11:18:57 +0100
Subject: [PATCH 04/12] #623 edit r_average to accept electrode domains
 (returns unchanged copy)

---
 pybamm/expression_tree/unary_operators.py      | 18 ++++--------------
 .../particle/fickian/fickian_many_particles.py |  2 ++
 2 files changed, 6 insertions(+), 14 deletions(-)

diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 8c502c3780..5739069831 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -595,9 +595,7 @@ def __init__(self, child, region="entire"):
     def set_id(self):
         """ See :meth:`pybamm.Symbol.set_id()` """
         self._id = hash(
-            (self.__class__, self.name)
-            + (self.children[0].id,)
-            + tuple(self.domain)
+            (self.__class__, self.name) + (self.children[0].id,) + tuple(self.domain)
         )
 
     def _unary_simplify(self, simplified_child):
@@ -798,6 +796,7 @@ def grad_squared(expression):
 
     return Gradient_Squared(expression)
 
+
 #
 # Method to call SurfaceValue
 #
@@ -1009,17 +1008,8 @@ def r_average(symbol):
     :class:`Symbol`
         the new averaged symbol
     """
-    ok_domains = [['positive particle'], ['negative particle'], []]
-    # Symbol must have domain [] or ["current collector"]
-    if symbol.domain not in ok_domains:
-        raise pybamm.DomainError(
-            """r-average only implemented in the 'particle' domain,
-            but symbol has domains {}""".format(
-                symbol.domain
-            )
-        )
-    # If symbol doesn't have a domain, its average value is itself
-    elif symbol.domain == []:
+    # If symbol doesn't have a particle domain, its r-averaged value is itself
+    if symbol.domain not in [["positive particle"], ["negative particle"]]:
         new_symbol = symbol.new_copy()
         new_symbol.parent = None
         return new_symbol
diff --git a/pybamm/models/submodels/particle/fickian/fickian_many_particles.py b/pybamm/models/submodels/particle/fickian/fickian_many_particles.py
index 91ce9a325a..159a0a710f 100644
--- a/pybamm/models/submodels/particle/fickian/fickian_many_particles.py
+++ b/pybamm/models/submodels/particle/fickian/fickian_many_particles.py
@@ -30,6 +30,8 @@ def get_fundamental_variables(self):
         elif self.domain == "Positive":
             c_s = pybamm.standard_variables.c_s_p
 
+        # TODO: implement c_s_xav for Fickian many particles (tricky because this
+        # requires averaging a secondary domain)
         variables = self._get_standard_concentration_variables(c_s, c_s)
 
         return variables

From 803b419142ddb3e086d2e9ee87ee09e147263f0a Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Sat, 21 Sep 2019 16:42:22 +0100
Subject: [PATCH 05/12] #623 flake8

---
 examples/scripts/SPMe_SOC.py                  | 70 +++++++++----------
 .../submodels/particle/base_particle.py       |  3 +-
 2 files changed, 36 insertions(+), 37 deletions(-)

diff --git a/examples/scripts/SPMe_SOC.py b/examples/scripts/SPMe_SOC.py
index 24ed89ed13..2a70b33ecc 100644
--- a/examples/scripts/SPMe_SOC.py
+++ b/examples/scripts/SPMe_SOC.py
@@ -1,7 +1,8 @@
 import pybamm
 import numpy as np
 import matplotlib.pyplot as plt
-plt.close('all')
+
+plt.close("all")
 pybamm.set_logging_level("INFO")
 
 # load model
@@ -13,26 +14,28 @@
 # load parameter values and process model and geometry
 param = model.default_parameter_values
 
-param.update({'Negative electrode surface area density [m-1]': 180000.0,
-              'Positive electrode surface area density [m-1]': 150000.0})
+param.update(
+    {
+        "Negative electrode surface area density [m-1]": 180000.0,
+        "Positive electrode surface area density [m-1]": 150000.0,
+    }
+)
 
 
-max_neg = param['Maximum concentration in negative electrode [mol.m-3]']
-max_pos = param['Maximum concentration in positive electrode [mol.m-3]']
-a_neg = param['Negative electrode surface area density [m-1]']
-a_pos = param['Positive electrode surface area density [m-1]']
-l_neg = param['Negative electrode thickness [m]']
-l_pos = param['Positive electrode thickness [m]']
-por_neg = param['Negative electrode porosity']
-por_pos = param['Positive electrode porosity']
+max_neg = param["Maximum concentration in negative electrode [mol.m-3]"]
+max_pos = param["Maximum concentration in positive electrode [mol.m-3]"]
+a_neg = param["Negative electrode surface area density [m-1]"]
+a_pos = param["Positive electrode surface area density [m-1]"]
+l_neg = param["Negative electrode thickness [m]"]
+l_pos = param["Positive electrode thickness [m]"]
+por_neg = param["Negative electrode porosity"]
+por_pos = param["Positive electrode porosity"]
 
 param.process_model(model)
 param.process_geometry(geometry)
 
 s_var = pybamm.standard_spatial_vars
-var_pts = {s_var.x_n: 5, s_var.x_s: 5,
-           s_var.x_p: 5, s_var.r_n: 5,
-           s_var.r_p: 10}
+var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5, s_var.r_n: 5, s_var.r_p: 10}
 # set mesh
 mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
 
@@ -45,34 +48,31 @@
 sol = model.default_solver.solve(model, t_eval)
 
 # plot
-out_vars = ['RX-averaged positive particle concentration [mol.m-3]',
-            'RX-averaged negative particle concentration [mol.m-3]',
-            'X-averaged electrolyte concentration [mol.m-3]',
-            'X-averaged positive electrode porosity',
-            'X-averaged negative electrode porosity',
-            'X-averaged separator porosity']
+out_vars = [
+    "RX-averaged positive particle concentration [mol.m-3]",
+    "RX-averaged negative particle concentration [mol.m-3]",
+    "X-averaged electrolyte concentration [mol.m-3]",
+    "X-averaged positive electrode porosity",
+    "X-averaged negative electrode porosity",
+    "X-averaged separator porosity",
+]
 plot = pybamm.QuickPlot(model, mesh, sol, output_variables=out_vars)
 plot.dynamic_plot()
 keys = list(model.variables.keys())
 keys.sort()
 
-xppc = pybamm.ProcessedVariable(model.variables[out_vars[0]],
-                                sol.t, sol.y, mesh=mesh)
-xnpc = pybamm.ProcessedVariable(model.variables[out_vars[1]],
-                                sol.t, sol.y, mesh=mesh)
-xec = pybamm.ProcessedVariable(model.variables[out_vars[2]],
-                               sol.t, sol.y, mesh=mesh)
+xppc = pybamm.ProcessedVariable(model.variables[out_vars[0]], sol.t, sol.y, mesh=mesh)
+xnpc = pybamm.ProcessedVariable(model.variables[out_vars[1]], sol.t, sol.y, mesh=mesh)
+xec = pybamm.ProcessedVariable(model.variables[out_vars[2]], sol.t, sol.y, mesh=mesh)
 rp = np.linspace(0, 1.0, 11)
 
 plt.figure()
 rplt = 0.0
-plt.plot(np.ones(len(sol.t))*max_neg*por_neg, 'r--', label='Max Neg')
-plt.plot(xnpc(sol.t, r=rplt)*por_neg, 'r', label='Neg Li')
-plt.plot(np.ones(len(sol.t))*max_pos*por_pos, 'b--', label='Max Pos')
-plt.plot(xppc(sol.t, r=rplt)*por_pos, 'b', label='Pos Li')
-plt.plot(xec(sol.t, r=rplt), 'g', label='Elec Li')
-tot = (xnpc(sol.t, r=rplt) * por_neg +
-       xppc(sol.t, r=rplt) * por_pos +
-       xec(sol.t, r=rplt))
-plt.plot(tot, 'k-', label='Total Li')
+plt.plot(np.ones(len(sol.t)) * max_neg * por_neg, "r--", label="Max Neg")
+plt.plot(xnpc(sol.t, r=rplt) * por_neg, "r", label="Neg Li")
+plt.plot(np.ones(len(sol.t)) * max_pos * por_pos, "b--", label="Max Pos")
+plt.plot(xppc(sol.t, r=rplt) * por_pos, "b", label="Pos Li")
+plt.plot(xec(sol.t, r=rplt), "g", label="Elec Li")
+tot = xnpc(sol.t, r=rplt) * por_neg + xppc(sol.t, r=rplt) * por_pos + xec(sol.t, r=rplt)
+plt.plot(tot, "k-", label="Total Li")
 plt.legend()
diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py
index 05b2bf8209..043accf70c 100644
--- a/pybamm/models/submodels/particle/base_particle.py
+++ b/pybamm/models/submodels/particle/base_particle.py
@@ -31,8 +31,7 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             c_scale = self.param.c_n_max
         elif self.domain == "Positive":
             c_scale = self.param.c_p_max
-        rad = pybamm.SpatialVariable("r", [self.domain.lower() + " particle"])
-#        c_s_r_av = pybamm.r_average(pybamm.x_average(c_s))
+        # c_s_r_av = pybamm.r_average(pybamm.x_average(c_s))
         c_s_r_av = pybamm.r_average(c_s_xav)
         variables = {
             self.domain + " particle concentration": c_s,

From e15f9aa2e42fcb8ec7ae28a0e857190fdae60c20 Mon Sep 17 00:00:00 2001
From: tom <t.g.tranter@gmail.com>
Date: Mon, 23 Sep 2019 14:37:16 +0100
Subject: [PATCH 06/12] add new variables for SOC calc

---
 examples/scripts/SPMe_SOC.py                      | 15 ++++++++-------
 pybamm/models/submodels/particle/base_particle.py |  6 ++++++
 2 files changed, 14 insertions(+), 7 deletions(-)

diff --git a/examples/scripts/SPMe_SOC.py b/examples/scripts/SPMe_SOC.py
index 2a70b33ecc..644ba66101 100644
--- a/examples/scripts/SPMe_SOC.py
+++ b/examples/scripts/SPMe_SOC.py
@@ -16,7 +16,7 @@
 
 param.update(
     {
-        "Negative electrode surface area density [m-1]": 180000.0,
+        "Negative electrode surface area density [m-1]": 150000.0,
         "Positive electrode surface area density [m-1]": 150000.0,
     }
 )
@@ -49,12 +49,13 @@
 
 # plot
 out_vars = [
+    "Negative active lithium",
+    "Positive active lithium",
+    "X-averaged electrolyte concentration [mol.m-3]",
+    "Negative active volume density",
+    "Positive active volume density",
     "RX-averaged positive particle concentration [mol.m-3]",
     "RX-averaged negative particle concentration [mol.m-3]",
-    "X-averaged electrolyte concentration [mol.m-3]",
-    "X-averaged positive electrode porosity",
-    "X-averaged negative electrode porosity",
-    "X-averaged separator porosity",
 ]
 plot = pybamm.QuickPlot(model, mesh, sol, output_variables=out_vars)
 plot.dynamic_plot()
@@ -68,9 +69,9 @@
 
 plt.figure()
 rplt = 0.0
-plt.plot(np.ones(len(sol.t)) * max_neg * por_neg, "r--", label="Max Neg")
+#plt.plot(np.ones(len(sol.t)) * max_neg * por_neg, "r--", label="Max Neg")
 plt.plot(xnpc(sol.t, r=rplt) * por_neg, "r", label="Neg Li")
-plt.plot(np.ones(len(sol.t)) * max_pos * por_pos, "b--", label="Max Pos")
+#plt.plot(np.ones(len(sol.t)) * max_pos * por_pos, "b--", label="Max Pos")
 plt.plot(xppc(sol.t, r=rplt) * por_pos, "b", label="Pos Li")
 plt.plot(xec(sol.t, r=rplt), "g", label="Elec Li")
 tot = xnpc(sol.t, r=rplt) * por_neg + xppc(sol.t, r=rplt) * por_pos + xec(sol.t, r=rplt)
diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py
index 043accf70c..0ea45e3998 100644
--- a/pybamm/models/submodels/particle/base_particle.py
+++ b/pybamm/models/submodels/particle/base_particle.py
@@ -26,13 +26,17 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
         c_s_surf = pybamm.surf(c_s, set_domain=True)
 
         c_s_surf_av = pybamm.x_average(c_s_surf)
+        geo_param = pybamm.geometric_parameters
 
         if self.domain == "Negative":
             c_scale = self.param.c_n_max
+            active_volume = geo_param.a_n_dim * geo_param.R_n / 3
         elif self.domain == "Positive":
             c_scale = self.param.c_p_max
+            active_volume = geo_param.a_p_dim * geo_param.R_p / 3
         # c_s_r_av = pybamm.r_average(pybamm.x_average(c_s))
         c_s_r_av = pybamm.r_average(c_s_xav)
+
         variables = {
             self.domain + " particle concentration": c_s,
             self.domain + " particle concentration [mol.m-3]": c_s * c_scale,
@@ -53,6 +57,8 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             "X-averaged "
             + self.domain.lower()
             + " particle surface concentration [mol.m-3]": c_scale * c_s_surf_av,
+            self.domain + " active volume density": active_volume,
+            self.domain + " active lithium": active_volume * c_s_r_av * c_scale
         }
 
         return variables

From e935dd5334fa49013ee3d1cbcac3207aa9f928fb Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Fri, 27 Sep 2019 15:26:31 +0100
Subject: [PATCH 07/12] output more variables and investigate SOC

---
 examples/scripts/SPMe_SOC.py                  | 196 ++++++++++++------
 .../submodels/particle/base_particle.py       |  14 +-
 2 files changed, 137 insertions(+), 73 deletions(-)

diff --git a/examples/scripts/SPMe_SOC.py b/examples/scripts/SPMe_SOC.py
index 644ba66101..b8ba5668bd 100644
--- a/examples/scripts/SPMe_SOC.py
+++ b/examples/scripts/SPMe_SOC.py
@@ -5,75 +5,139 @@
 plt.close("all")
 pybamm.set_logging_level("INFO")
 
-# load model
-model = pybamm.lithium_ion.SPMe()
+# Dimensions
+h = 0.137
+w = 0.207
+A = h * w
+l_n = 1e-3
+l_p = 1e-3
+l_s = 2.5e-5
+l1d = (l_n + l_p + l_s)
+vol = h * w * l1d # m^3
+vol_cm3 = vol * 1e6
 
-# create geometry
-geometry = model.default_geometry
+tot_cap = 0.0
+tot_time = 0.0
 
-# load parameter values and process model and geometry
-param = model.default_parameter_values
+I_mag = 1.0
+for I_app in [-1.0, 1.0]:
+    I_app*=I_mag
+    # load model
+    model = pybamm.lithium_ion.SPMe()
+    # create geometry
+    geometry = model.default_geometry
+    # load parameter values and process model and geometry
+    param = model.default_parameter_values
+    #I_app = 1.0 # [A]
+    param.update(
+        {
+#            "Electrode height [m]": h,
+#            "Electrode width [m]": w,
+            "Negative electrode thickness [m]": l_n,
+            "Positive electrode thickness [m]": l_p,
+            "Separator thickness [m]": l_s,
+            "Lower voltage cut-off [V]": 3.105,
+            "Upper voltage cut-off [V]": 4.7,
+            "Maximum concentration in negative electrode [mol.m-3]": 24983.261993843702,
+            "Maximum concentration in positive electrode [mol.m-3]": 51217.9257309275,
+            "Negative electrode surface area density [m-1]": 180000.0,
+            "Positive electrode surface area density [m-1]": 150000.0,
+            "Typical current [A]": I_app,
+        }
+    )
 
-param.update(
-    {
-        "Negative electrode surface area density [m-1]": 150000.0,
-        "Positive electrode surface area density [m-1]": 150000.0,
-    }
-)
+#    max_neg = param["Maximum concentration in negative electrode [mol.m-3]"]
+#    max_pos = param["Maximum concentration in positive electrode [mol.m-3]"]
+#    init_neg = param["Initial concentration in negative electrode [mol.m-3]"]
+#    init_pos = param["Initial concentration in positive electrode [mol.m-3]"]
+#    a_neg = param["Negative electrode surface area density [m-1]"]
+#    a_pos = param["Positive electrode surface area density [m-1]"]
+#    l_neg = param["Negative electrode thickness [m]"]
+#    l_pos = param["Positive electrode thickness [m]"]
+#    l_sep = param["Separator thickness [m]"]
+#    por_neg = param["Negative electrode porosity"]
+#    por_pos = param["Positive electrode porosity"]
+#    por_sep = param["Separator porosity"]
+    
+    param.process_model(model)
+    param.process_geometry(geometry)
+    s_var = pybamm.standard_spatial_vars
+    var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5, s_var.r_n: 5, s_var.r_p: 10}
+    # set mesh
+    mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
+    # discretise model
+    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+    disc.process_model(model)
+    # solve model
+    t_eval = np.linspace(0, 0.2, 100)
+    sol = model.default_solver.solve(model, t_eval)
+    
+    # plot
+    vars = [
+        "Electrolyte concentration",
+        "Positive electrode volume-averaged concentration [mol.m-3]",
+        "Negative electrode volume-averaged concentration [mol.m-3]",
+        "Positive electrode average extent of lithiation",
+        "Negative electrode average extent of lithiation",
+        "X-averaged positive electrolyte concentration [mol.m-3]",
+        "X-averaged negative electrolyte concentration [mol.m-3]",
+        "X-averaged separator electrolyte concentration [mol.m-3]"
+    ]
+    plot = pybamm.QuickPlot(model, mesh, sol, output_variables=vars)
+    plot.dynamic_plot()
+#    keys = list(model.variables.keys())
+#    keys.sort()
+    
+#    avol_pos = pybamm.ProcessedVariable(model.variables[vars[0]], sol.t, sol.y, mesh=mesh)
+#    avol_neg = pybamm.ProcessedVariable(model.variables[vars[1]], sol.t, sol.y, mesh=mesh)
+#    xppc = pybamm.ProcessedVariable(model.variables[vars[2]], sol.t, sol.y, mesh=mesh)
+#    xnpc = pybamm.ProcessedVariable(model.variables[vars[3]], sol.t, sol.y, mesh=mesh)
+    xpext = pybamm.ProcessedVariable(model.variables[vars[3]], sol.t, sol.y, mesh=mesh)
+    xnext = pybamm.ProcessedVariable(model.variables[vars[4]], sol.t, sol.y, mesh=mesh)
+#    xppce = pybamm.ProcessedVariable(model.variables[vars[6]], sol.t, sol.y, mesh=mesh)
+#    xnpce = pybamm.ProcessedVariable(model.variables[vars[7]], sol.t, sol.y, mesh=mesh)
+#    xsepe = pybamm.ProcessedVariable(model.variables[vars[8]], sol.t, sol.y, mesh=mesh)
+    time = pybamm.ProcessedVariable(model.variables["Time [h]"], sol.t, sol.y, mesh=mesh)
+#    rp = np.linspace(0, 1.0, 11)
+#        
+#    plt.figure()
+#    rplt = 0.0
+#    plt.plot(avol_neg(sol.t, r=rplt) * max_neg, "r--", label="Max Neg")
+#    plt.plot(avol_neg(sol.t, r=rplt) * init_neg, "r*-", label="Init Neg")
+#    plt.plot(avol_pos(sol.t, r=rplt) * max_pos, "b--", label="Max Pos")
+#    plt.plot(avol_pos(sol.t, r=rplt) * init_pos, "b*-", label="Init Pos")
+#    plt.plot(xnpc(sol.t, r=rplt), "r", label="Neg Li")
+#    plt.plot(xppc(sol.t, r=rplt), "b", label="Pos Li")
+#    
+#    pos_electrolyte_li = xppce(sol.t)*l_pos*por_pos
+#    neg_electrolyte_li = xnpce(sol.t)*l_neg*por_neg
+#    sep_electrolyte_li = xsepe(sol.t)*l_sep*por_sep
+#    tot_electrolyte_li = pos_electrolyte_li+neg_electrolyte_li+sep_electrolyte_li
+#    plt.plot(tot_electrolyte_li, "g", label="Elec Li")
+#    tot_li = xnpc(sol.t, r=rplt) + xppc(sol.t, r=rplt) + tot_electrolyte_li
+#    plt.plot(tot_li, "k-", label="Total Li")
+#    plt.legend()
+    
+    # Coulomb counting
+    time_hours = time(sol.t)
+    dc_time = np.around(time_hours[-1], 3)
+    cap = np.absolute(I_app * 1000 * dc_time) # mAh
+    plt.figure()
+    plt.plot(np.absolute(I_app)*1000*time_hours, xnext(sol.t), 'r-', label='Negative')
+    plt.plot(np.absolute(I_app)*1000*time_hours, xpext(sol.t), 'b-', label='Positive')
+    plt.xlabel('Capacity [mAh]')
+    plt.ylabel('Extent of Lithiation')
+    plt.legend()
+    if I_app < 0.0:
+        plt.title('Charge')
+    else:
+        plt.title('Discharge')
+    print('Applied Current', I_app, 'A', 'Time', dc_time, 'hrs', 'Capacity', cap, 'mAh')
+    tot_cap += cap
+    tot_time += dc_time
 
+print('Total Charge/Discharge Time', tot_time, 'hrs')
+print('Total Capacity', np.around(tot_cap, 3), 'mAh')
+print('Total Capacity', np.around(tot_cap, 3)/vol_cm3, 'mAh.cm-3')
 
-max_neg = param["Maximum concentration in negative electrode [mol.m-3]"]
-max_pos = param["Maximum concentration in positive electrode [mol.m-3]"]
-a_neg = param["Negative electrode surface area density [m-1]"]
-a_pos = param["Positive electrode surface area density [m-1]"]
-l_neg = param["Negative electrode thickness [m]"]
-l_pos = param["Positive electrode thickness [m]"]
-por_neg = param["Negative electrode porosity"]
-por_pos = param["Positive electrode porosity"]
 
-param.process_model(model)
-param.process_geometry(geometry)
-
-s_var = pybamm.standard_spatial_vars
-var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5, s_var.r_n: 5, s_var.r_p: 10}
-# set mesh
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-# solve model
-t_eval = np.linspace(0, 0.2, 100)
-sol = model.default_solver.solve(model, t_eval)
-
-# plot
-out_vars = [
-    "Negative active lithium",
-    "Positive active lithium",
-    "X-averaged electrolyte concentration [mol.m-3]",
-    "Negative active volume density",
-    "Positive active volume density",
-    "RX-averaged positive particle concentration [mol.m-3]",
-    "RX-averaged negative particle concentration [mol.m-3]",
-]
-plot = pybamm.QuickPlot(model, mesh, sol, output_variables=out_vars)
-plot.dynamic_plot()
-keys = list(model.variables.keys())
-keys.sort()
-
-xppc = pybamm.ProcessedVariable(model.variables[out_vars[0]], sol.t, sol.y, mesh=mesh)
-xnpc = pybamm.ProcessedVariable(model.variables[out_vars[1]], sol.t, sol.y, mesh=mesh)
-xec = pybamm.ProcessedVariable(model.variables[out_vars[2]], sol.t, sol.y, mesh=mesh)
-rp = np.linspace(0, 1.0, 11)
-
-plt.figure()
-rplt = 0.0
-#plt.plot(np.ones(len(sol.t)) * max_neg * por_neg, "r--", label="Max Neg")
-plt.plot(xnpc(sol.t, r=rplt) * por_neg, "r", label="Neg Li")
-#plt.plot(np.ones(len(sol.t)) * max_pos * por_pos, "b--", label="Max Pos")
-plt.plot(xppc(sol.t, r=rplt) * por_pos, "b", label="Pos Li")
-plt.plot(xec(sol.t, r=rplt), "g", label="Elec Li")
-tot = xnpc(sol.t, r=rplt) * por_neg + xppc(sol.t, r=rplt) * por_pos + xec(sol.t, r=rplt)
-plt.plot(tot, "k-", label="Total Li")
-plt.legend()
diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py
index 0ea45e3998..e0f1b055e1 100644
--- a/pybamm/models/submodels/particle/base_particle.py
+++ b/pybamm/models/submodels/particle/base_particle.py
@@ -34,7 +34,6 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
         elif self.domain == "Positive":
             c_scale = self.param.c_p_max
             active_volume = geo_param.a_p_dim * geo_param.R_p / 3
-        # c_s_r_av = pybamm.r_average(pybamm.x_average(c_s))
         c_s_r_av = pybamm.r_average(c_s_xav)
 
         variables = {
@@ -44,10 +43,10 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             "X-averaged "
             + self.domain.lower()
             + " particle concentration [mol.m-3]": c_s_xav * c_scale,
-            "RX-averaged " + self.domain.lower() + " particle concentration": c_s_r_av,
-            "RX-averaged "
-            + self.domain.lower()
-            + " particle concentration [mol.m-3]": c_s_r_av * c_scale,
+#            "RX-averaged " + self.domain.lower() + " particle concentration": c_s_r_av,
+#            "RX-averaged "
+#            + self.domain.lower()
+#            + " particle concentration [mol.m-3]": c_s_r_av * c_scale,
             self.domain + " particle surface concentration": c_s_surf,
             self.domain
             + " particle surface concentration [mol.m-3]": c_scale * c_s_surf,
@@ -57,8 +56,9 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             "X-averaged "
             + self.domain.lower()
             + " particle surface concentration [mol.m-3]": c_scale * c_s_surf_av,
-            self.domain + " active volume density": active_volume,
-            self.domain + " active lithium": active_volume * c_s_r_av * c_scale
+            self.domain + " electrode active volume fraction": active_volume,
+            self.domain + " electrode volume-averaged concentration [mol.m-3]": active_volume * c_s_r_av * c_scale,
+            self.domain + " electrode average extent of lithiation": c_s_r_av
         }
 
         return variables

From ddef3c3c2ccd0518f0ea84983f8aa6e1e136b503 Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Fri, 27 Sep 2019 15:39:37 +0100
Subject: [PATCH 08/12] correct the dimensions

---
 examples/scripts/SPMe_SOC.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/examples/scripts/SPMe_SOC.py b/examples/scripts/SPMe_SOC.py
index b8ba5668bd..9738686e64 100644
--- a/examples/scripts/SPMe_SOC.py
+++ b/examples/scripts/SPMe_SOC.py
@@ -5,12 +5,14 @@
 plt.close("all")
 pybamm.set_logging_level("INFO")
 
+factor = 6.38
+
 # Dimensions
 h = 0.137
-w = 0.207
+w = 0.207/factor
 A = h * w
-l_n = 1e-3
-l_p = 1e-3
+l_n = 1e-4
+l_p = 1e-4
 l_s = 2.5e-5
 l1d = (l_n + l_p + l_s)
 vol = h * w * l1d # m^3
@@ -19,7 +21,7 @@
 tot_cap = 0.0
 tot_time = 0.0
 
-I_mag = 1.0
+I_mag = 862e-3
 for I_app in [-1.0, 1.0]:
     I_app*=I_mag
     # load model
@@ -83,8 +85,8 @@
         "X-averaged negative electrolyte concentration [mol.m-3]",
         "X-averaged separator electrolyte concentration [mol.m-3]"
     ]
-    plot = pybamm.QuickPlot(model, mesh, sol, output_variables=vars)
-    plot.dynamic_plot()
+#    plot = pybamm.QuickPlot(model, mesh, sol, output_variables=vars)
+#    plot.dynamic_plot()
 #    keys = list(model.variables.keys())
 #    keys.sort()
     

From c5b5d6f0ac74dc3475b75eacf823189114fd1fbd Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Fri, 27 Sep 2019 16:26:42 +0100
Subject: [PATCH 09/12] tidy up example

---
 examples/scripts/SPMe_SOC.py | 119 +++++++++++------------------------
 1 file changed, 37 insertions(+), 82 deletions(-)

diff --git a/examples/scripts/SPMe_SOC.py b/examples/scripts/SPMe_SOC.py
index 9738686e64..f06fe89ac4 100644
--- a/examples/scripts/SPMe_SOC.py
+++ b/examples/scripts/SPMe_SOC.py
@@ -15,56 +15,46 @@
 l_p = 1e-4
 l_s = 2.5e-5
 l1d = (l_n + l_p + l_s)
-vol = h * w * l1d # m^3
+vol = h * w * l1d
 vol_cm3 = vol * 1e6
 
 tot_cap = 0.0
 tot_time = 0.0
-
-I_mag = 862e-3
-for I_app in [-1.0, 1.0]:
-    I_app*=I_mag
+fig, axes = plt.subplots(1, 2, sharey=True)
+I_mag = 1.01/factor
+for enum, I_app in enumerate([-1.0, 1.0]):
+    I_app *= I_mag
     # load model
     model = pybamm.lithium_ion.SPMe()
     # create geometry
     geometry = model.default_geometry
     # load parameter values and process model and geometry
     param = model.default_parameter_values
-    #I_app = 1.0 # [A]
+
     param.update(
         {
-#            "Electrode height [m]": h,
-#            "Electrode width [m]": w,
+            "Electrode height [m]": h,
+            "Electrode width [m]": w,
             "Negative electrode thickness [m]": l_n,
             "Positive electrode thickness [m]": l_p,
             "Separator thickness [m]": l_s,
             "Lower voltage cut-off [V]": 3.105,
             "Upper voltage cut-off [V]": 4.7,
-            "Maximum concentration in negative electrode [mol.m-3]": 24983.261993843702,
-            "Maximum concentration in positive electrode [mol.m-3]": 51217.9257309275,
+            "Maximum concentration in negative electrode [mol.m-3]": 25000,
+            "Maximum concentration in positive electrode [mol.m-3]": 50000,
+            "Initial concentration in negative electrode [mol.m-3]": 12500,
+            "Initial concentration in positive electrode [mol.m-3]": 25000,
             "Negative electrode surface area density [m-1]": 180000.0,
             "Positive electrode surface area density [m-1]": 150000.0,
             "Typical current [A]": I_app,
         }
     )
 
-#    max_neg = param["Maximum concentration in negative electrode [mol.m-3]"]
-#    max_pos = param["Maximum concentration in positive electrode [mol.m-3]"]
-#    init_neg = param["Initial concentration in negative electrode [mol.m-3]"]
-#    init_pos = param["Initial concentration in positive electrode [mol.m-3]"]
-#    a_neg = param["Negative electrode surface area density [m-1]"]
-#    a_pos = param["Positive electrode surface area density [m-1]"]
-#    l_neg = param["Negative electrode thickness [m]"]
-#    l_pos = param["Positive electrode thickness [m]"]
-#    l_sep = param["Separator thickness [m]"]
-#    por_neg = param["Negative electrode porosity"]
-#    por_pos = param["Positive electrode porosity"]
-#    por_sep = param["Separator porosity"]
-    
     param.process_model(model)
     param.process_geometry(geometry)
     s_var = pybamm.standard_spatial_vars
-    var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5, s_var.r_n: 5, s_var.r_p: 10}
+    var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5,
+               s_var.r_n: 5, s_var.r_p: 10}
     # set mesh
     mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
     # discretise model
@@ -73,73 +63,38 @@
     # solve model
     t_eval = np.linspace(0, 0.2, 100)
     sol = model.default_solver.solve(model, t_eval)
-    
-    # plot
-    vars = [
-        "Electrolyte concentration",
-        "Positive electrode volume-averaged concentration [mol.m-3]",
-        "Negative electrode volume-averaged concentration [mol.m-3]",
-        "Positive electrode average extent of lithiation",
-        "Negative electrode average extent of lithiation",
-        "X-averaged positive electrolyte concentration [mol.m-3]",
-        "X-averaged negative electrolyte concentration [mol.m-3]",
-        "X-averaged separator electrolyte concentration [mol.m-3]"
-    ]
-#    plot = pybamm.QuickPlot(model, mesh, sol, output_variables=vars)
-#    plot.dynamic_plot()
-#    keys = list(model.variables.keys())
-#    keys.sort()
-    
-#    avol_pos = pybamm.ProcessedVariable(model.variables[vars[0]], sol.t, sol.y, mesh=mesh)
-#    avol_neg = pybamm.ProcessedVariable(model.variables[vars[1]], sol.t, sol.y, mesh=mesh)
-#    xppc = pybamm.ProcessedVariable(model.variables[vars[2]], sol.t, sol.y, mesh=mesh)
-#    xnpc = pybamm.ProcessedVariable(model.variables[vars[3]], sol.t, sol.y, mesh=mesh)
-    xpext = pybamm.ProcessedVariable(model.variables[vars[3]], sol.t, sol.y, mesh=mesh)
-    xnext = pybamm.ProcessedVariable(model.variables[vars[4]], sol.t, sol.y, mesh=mesh)
-#    xppce = pybamm.ProcessedVariable(model.variables[vars[6]], sol.t, sol.y, mesh=mesh)
-#    xnpce = pybamm.ProcessedVariable(model.variables[vars[7]], sol.t, sol.y, mesh=mesh)
-#    xsepe = pybamm.ProcessedVariable(model.variables[vars[8]], sol.t, sol.y, mesh=mesh)
-    time = pybamm.ProcessedVariable(model.variables["Time [h]"], sol.t, sol.y, mesh=mesh)
-#    rp = np.linspace(0, 1.0, 11)
-#        
-#    plt.figure()
-#    rplt = 0.0
-#    plt.plot(avol_neg(sol.t, r=rplt) * max_neg, "r--", label="Max Neg")
-#    plt.plot(avol_neg(sol.t, r=rplt) * init_neg, "r*-", label="Init Neg")
-#    plt.plot(avol_pos(sol.t, r=rplt) * max_pos, "b--", label="Max Pos")
-#    plt.plot(avol_pos(sol.t, r=rplt) * init_pos, "b*-", label="Init Pos")
-#    plt.plot(xnpc(sol.t, r=rplt), "r", label="Neg Li")
-#    plt.plot(xppc(sol.t, r=rplt), "b", label="Pos Li")
-#    
-#    pos_electrolyte_li = xppce(sol.t)*l_pos*por_pos
-#    neg_electrolyte_li = xnpce(sol.t)*l_neg*por_neg
-#    sep_electrolyte_li = xsepe(sol.t)*l_sep*por_sep
-#    tot_electrolyte_li = pos_electrolyte_li+neg_electrolyte_li+sep_electrolyte_li
-#    plt.plot(tot_electrolyte_li, "g", label="Elec Li")
-#    tot_li = xnpc(sol.t, r=rplt) + xppc(sol.t, r=rplt) + tot_electrolyte_li
-#    plt.plot(tot_li, "k-", label="Total Li")
-#    plt.legend()
-    
+
+    var = "Positive electrode average extent of lithiation"
+    xpext = pybamm.ProcessedVariable(model.variables[var],
+                                     sol.t, sol.y, mesh=mesh)
+    var = "Negative electrode average extent of lithiation"
+    xnext = pybamm.ProcessedVariable(model.variables[var],
+                                     sol.t, sol.y, mesh=mesh)
+    time = pybamm.ProcessedVariable(model.variables["Time [h]"],
+                                    sol.t, sol.y, mesh=mesh)
+
     # Coulomb counting
     time_hours = time(sol.t)
     dc_time = np.around(time_hours[-1], 3)
-    cap = np.absolute(I_app * 1000 * dc_time) # mAh
-    plt.figure()
-    plt.plot(np.absolute(I_app)*1000*time_hours, xnext(sol.t), 'r-', label='Negative')
-    plt.plot(np.absolute(I_app)*1000*time_hours, xpext(sol.t), 'b-', label='Positive')
-    plt.xlabel('Capacity [mAh]')
-    plt.ylabel('Extent of Lithiation')
+    # Capacity mAh
+    cap = np.absolute(I_app * 1000 * dc_time)
+
+    axes[enum].plot(np.absolute(I_app)*1000*time_hours,
+                    xnext(sol.t), 'r-', label='Negative')
+    axes[enum].plot(np.absolute(I_app)*1000*time_hours,
+                    xpext(sol.t), 'b-', label='Positive')
+    axes[enum].set_xlabel('Capacity [mAh]')
     plt.legend()
     if I_app < 0.0:
-        plt.title('Charge')
+        axes[enum].set_ylabel('Extent of Lithiation')
+        axes[enum].title.set_text('Charge')
     else:
-        plt.title('Discharge')
-    print('Applied Current', I_app, 'A', 'Time', dc_time, 'hrs', 'Capacity', cap, 'mAh')
+        axes[enum].title.set_text('Discharge')
+    print('Applied Current', I_app, 'A', 'Time',
+          dc_time, 'hrs', 'Capacity', cap, 'mAh')
     tot_cap += cap
     tot_time += dc_time
 
 print('Total Charge/Discharge Time', tot_time, 'hrs')
 print('Total Capacity', np.around(tot_cap, 3), 'mAh')
 print('Total Capacity', np.around(tot_cap, 3)/vol_cm3, 'mAh.cm-3')
-
-

From 075ec894a65e69284cd672698ec3a53599b84827 Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Fri, 27 Sep 2019 16:33:43 +0100
Subject: [PATCH 10/12] add surface conc to plot

---
 examples/scripts/SPMe_SOC.py | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/examples/scripts/SPMe_SOC.py b/examples/scripts/SPMe_SOC.py
index f06fe89ac4..8de141ccaf 100644
--- a/examples/scripts/SPMe_SOC.py
+++ b/examples/scripts/SPMe_SOC.py
@@ -70,6 +70,12 @@
     var = "Negative electrode average extent of lithiation"
     xnext = pybamm.ProcessedVariable(model.variables[var],
                                      sol.t, sol.y, mesh=mesh)
+    var = "X-averaged positive particle surface concentration"
+    xpsurf = pybamm.ProcessedVariable(model.variables[var],
+                                      sol.t, sol.y, mesh=mesh)
+    var = "X-averaged negative particle surface concentration"
+    xnsurf = pybamm.ProcessedVariable(model.variables[var],
+                                      sol.t, sol.y, mesh=mesh)
     time = pybamm.ProcessedVariable(model.variables["Time [h]"],
                                     sol.t, sol.y, mesh=mesh)
 
@@ -80,9 +86,13 @@
     cap = np.absolute(I_app * 1000 * dc_time)
 
     axes[enum].plot(np.absolute(I_app)*1000*time_hours,
-                    xnext(sol.t), 'r-', label='Negative')
+                    xnext(sol.t), 'r-', label='Average Neg')
     axes[enum].plot(np.absolute(I_app)*1000*time_hours,
-                    xpext(sol.t), 'b-', label='Positive')
+                    xpext(sol.t), 'b-', label='Average Pos')
+    axes[enum].plot(np.absolute(I_app)*1000*time_hours,
+                    xnsurf(sol.t), 'r--', label='Surface Neg')
+    axes[enum].plot(np.absolute(I_app)*1000*time_hours,
+                    xpsurf(sol.t), 'b--', label='Surface Pos')
     axes[enum].set_xlabel('Capacity [mAh]')
     plt.legend()
     if I_app < 0.0:

From 40524124d042f78769a5e67c32740675c098e18c Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Mon, 30 Sep 2019 11:49:44 +0100
Subject: [PATCH 11/12] flake8

---
 examples/scripts/SPMe_SOC.py                  | 46 +++++++++----------
 .../submodels/particle/base_particle.py       |  9 ++--
 2 files changed, 26 insertions(+), 29 deletions(-)

diff --git a/examples/scripts/SPMe_SOC.py b/examples/scripts/SPMe_SOC.py
index 8de141ccaf..32adb72064 100644
--- a/examples/scripts/SPMe_SOC.py
+++ b/examples/scripts/SPMe_SOC.py
@@ -9,7 +9,7 @@
 
 # Dimensions
 h = 0.137
-w = 0.207/factor
+w = 0.207 / factor
 A = h * w
 l_n = 1e-4
 l_p = 1e-4
@@ -21,7 +21,7 @@
 tot_cap = 0.0
 tot_time = 0.0
 fig, axes = plt.subplots(1, 2, sharey=True)
-I_mag = 1.01/factor
+I_mag = 1.01 / factor
 for enum, I_app in enumerate([-1.0, 1.0]):
     I_app *= I_mag
     # load model
@@ -32,22 +32,21 @@
     param = model.default_parameter_values
 
     param.update(
-        {
-            "Electrode height [m]": h,
-            "Electrode width [m]": w,
-            "Negative electrode thickness [m]": l_n,
-            "Positive electrode thickness [m]": l_p,
-            "Separator thickness [m]": l_s,
-            "Lower voltage cut-off [V]": 3.105,
-            "Upper voltage cut-off [V]": 4.7,
-            "Maximum concentration in negative electrode [mol.m-3]": 25000,
-            "Maximum concentration in positive electrode [mol.m-3]": 50000,
-            "Initial concentration in negative electrode [mol.m-3]": 12500,
-            "Initial concentration in positive electrode [mol.m-3]": 25000,
-            "Negative electrode surface area density [m-1]": 180000.0,
-            "Positive electrode surface area density [m-1]": 150000.0,
-            "Typical current [A]": I_app,
-        }
+        {"Electrode height [m]": h,
+         "Electrode width [m]": w,
+         "Negative electrode thickness [m]": l_n,
+         "Positive electrode thickness [m]": l_p,
+         "Separator thickness [m]": l_s,
+         "Lower voltage cut-off [V]": 3.105,
+         "Upper voltage cut-off [V]": 4.7,
+         "Maximum concentration in negative electrode [mol.m-3]": 25000,
+         "Maximum concentration in positive electrode [mol.m-3]": 50000,
+         "Initial concentration in negative electrode [mol.m-3]": 12500,
+         "Initial concentration in positive electrode [mol.m-3]": 25000,
+         "Negative electrode surface area density [m-1]": 180000.0,
+         "Positive electrode surface area density [m-1]": 150000.0,
+         "Typical current [A]": I_app,
+         }
     )
 
     param.process_model(model)
@@ -84,14 +83,15 @@
     dc_time = np.around(time_hours[-1], 3)
     # Capacity mAh
     cap = np.absolute(I_app * 1000 * dc_time)
+    cap_time = np.absolute(I_app * 1000 * time_hours)
 
-    axes[enum].plot(np.absolute(I_app)*1000*time_hours,
+    axes[enum].plot(cap_time,
                     xnext(sol.t), 'r-', label='Average Neg')
-    axes[enum].plot(np.absolute(I_app)*1000*time_hours,
+    axes[enum].plot(cap_time,
                     xpext(sol.t), 'b-', label='Average Pos')
-    axes[enum].plot(np.absolute(I_app)*1000*time_hours,
+    axes[enum].plot(cap_time,
                     xnsurf(sol.t), 'r--', label='Surface Neg')
-    axes[enum].plot(np.absolute(I_app)*1000*time_hours,
+    axes[enum].plot(cap_time,
                     xpsurf(sol.t), 'b--', label='Surface Pos')
     axes[enum].set_xlabel('Capacity [mAh]')
     plt.legend()
@@ -107,4 +107,4 @@
 
 print('Total Charge/Discharge Time', tot_time, 'hrs')
 print('Total Capacity', np.around(tot_cap, 3), 'mAh')
-print('Total Capacity', np.around(tot_cap, 3)/vol_cm3, 'mAh.cm-3')
+print('Total Capacity', np.around(tot_cap, 3) / vol_cm3, 'mAh.cm-3')
diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py
index e0f1b055e1..491c607f02 100644
--- a/pybamm/models/submodels/particle/base_particle.py
+++ b/pybamm/models/submodels/particle/base_particle.py
@@ -35,7 +35,7 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             c_scale = self.param.c_p_max
             active_volume = geo_param.a_p_dim * geo_param.R_p / 3
         c_s_r_av = pybamm.r_average(c_s_xav)
-
+        c_s_r_av_vol = active_volume * c_s_r_av * c_scale
         variables = {
             self.domain + " particle concentration": c_s,
             self.domain + " particle concentration [mol.m-3]": c_s * c_scale,
@@ -43,10 +43,6 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             "X-averaged "
             + self.domain.lower()
             + " particle concentration [mol.m-3]": c_s_xav * c_scale,
-#            "RX-averaged " + self.domain.lower() + " particle concentration": c_s_r_av,
-#            "RX-averaged "
-#            + self.domain.lower()
-#            + " particle concentration [mol.m-3]": c_s_r_av * c_scale,
             self.domain + " particle surface concentration": c_s_surf,
             self.domain
             + " particle surface concentration [mol.m-3]": c_scale * c_s_surf,
@@ -57,7 +53,8 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             + self.domain.lower()
             + " particle surface concentration [mol.m-3]": c_scale * c_s_surf_av,
             self.domain + " electrode active volume fraction": active_volume,
-            self.domain + " electrode volume-averaged concentration [mol.m-3]": active_volume * c_s_r_av * c_scale,
+            self.domain
+            + " electrode volume-averaged concentration [mol.m-3]": c_s_r_av_vol,
             self.domain + " electrode average extent of lithiation": c_s_r_av
         }
 

From 9a9aace16f5f1de201887890c5dfd2fc379c9029 Mon Sep 17 00:00:00 2001
From: Tom Tranter <t.g.tranter@gmail.com>
Date: Tue, 1 Oct 2019 11:19:53 +0100
Subject: [PATCH 12/12] tidy up SOC and make example more interesting + change
 solve termination to warning

---
 examples/scripts/SPMe_SOC.py                  | 221 ++++++++++--------
 .../submodels/particle/base_particle.py       |   6 +-
 pybamm/solvers/base_solver.py                 |   2 +-
 3 files changed, 126 insertions(+), 103 deletions(-)

diff --git a/examples/scripts/SPMe_SOC.py b/examples/scripts/SPMe_SOC.py
index 32adb72064..411eb50f67 100644
--- a/examples/scripts/SPMe_SOC.py
+++ b/examples/scripts/SPMe_SOC.py
@@ -1,110 +1,131 @@
+#
+#  Example to show the state of charge of a battery using the SPMe model
+#  Initial conditions are specified to start each electrode in 1/2 charged
+#  state. A charge and discharge are performed with current chosen to be
+#  1C rate when electrode dimensions are euqal.
+#  Coulomb counting is performed to calculate the capacity of the
+#  battery within the operating voltage limits and maximum particle concs.
+#  The anode thickenss is varied to highlight the importance of electrode
+#  sizing to enable full lithium utilization
+#
 import pybamm
 import numpy as np
 import matplotlib.pyplot as plt
 
 plt.close("all")
-pybamm.set_logging_level("INFO")
+pybamm.set_logging_level(30)
 
 factor = 6.38
-
-# Dimensions
-h = 0.137
-w = 0.207 / factor
-A = h * w
-l_n = 1e-4
+capacities = []
+specific_capacities = []
 l_p = 1e-4
-l_s = 2.5e-5
-l1d = (l_n + l_p + l_s)
-vol = h * w * l1d
-vol_cm3 = vol * 1e6
-
-tot_cap = 0.0
-tot_time = 0.0
-fig, axes = plt.subplots(1, 2, sharey=True)
-I_mag = 1.01 / factor
-for enum, I_app in enumerate([-1.0, 1.0]):
-    I_app *= I_mag
-    # load model
-    model = pybamm.lithium_ion.SPMe()
-    # create geometry
-    geometry = model.default_geometry
-    # load parameter values and process model and geometry
-    param = model.default_parameter_values
-
-    param.update(
-        {"Electrode height [m]": h,
-         "Electrode width [m]": w,
-         "Negative electrode thickness [m]": l_n,
-         "Positive electrode thickness [m]": l_p,
-         "Separator thickness [m]": l_s,
-         "Lower voltage cut-off [V]": 3.105,
-         "Upper voltage cut-off [V]": 4.7,
-         "Maximum concentration in negative electrode [mol.m-3]": 25000,
-         "Maximum concentration in positive electrode [mol.m-3]": 50000,
-         "Initial concentration in negative electrode [mol.m-3]": 12500,
-         "Initial concentration in positive electrode [mol.m-3]": 25000,
-         "Negative electrode surface area density [m-1]": 180000.0,
-         "Positive electrode surface area density [m-1]": 150000.0,
-         "Typical current [A]": I_app,
-         }
-    )
-
-    param.process_model(model)
-    param.process_geometry(geometry)
-    s_var = pybamm.standard_spatial_vars
-    var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5,
-               s_var.r_n: 5, s_var.r_p: 10}
-    # set mesh
-    mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-    # discretise model
-    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-    disc.process_model(model)
-    # solve model
-    t_eval = np.linspace(0, 0.2, 100)
-    sol = model.default_solver.solve(model, t_eval)
-
-    var = "Positive electrode average extent of lithiation"
-    xpext = pybamm.ProcessedVariable(model.variables[var],
-                                     sol.t, sol.y, mesh=mesh)
-    var = "Negative electrode average extent of lithiation"
-    xnext = pybamm.ProcessedVariable(model.variables[var],
-                                     sol.t, sol.y, mesh=mesh)
-    var = "X-averaged positive particle surface concentration"
-    xpsurf = pybamm.ProcessedVariable(model.variables[var],
-                                      sol.t, sol.y, mesh=mesh)
-    var = "X-averaged negative particle surface concentration"
-    xnsurf = pybamm.ProcessedVariable(model.variables[var],
-                                      sol.t, sol.y, mesh=mesh)
-    time = pybamm.ProcessedVariable(model.variables["Time [h]"],
-                                    sol.t, sol.y, mesh=mesh)
-
-    # Coulomb counting
-    time_hours = time(sol.t)
-    dc_time = np.around(time_hours[-1], 3)
-    # Capacity mAh
-    cap = np.absolute(I_app * 1000 * dc_time)
-    cap_time = np.absolute(I_app * 1000 * time_hours)
+thicknesses = np.linspace(1.0, 2.5, 11) * l_p
+for l_n in thicknesses:
+    e_ratio = np.around(l_n / l_p, 3)
+    # Dimensions
+    h = 0.137
+    w = 0.207 / factor
+    A = h * w
+    l_s = 2.5e-5
+    l1d = (l_n + l_p + l_s)
+    vol = h * w * l1d
+    vol_cm3 = vol * 1e6
+    tot_cap = 0.0
+    tot_time = 0.0
+    fig, axes = plt.subplots(1, 2, sharey=True)
+    I_mag = 1.01 / factor
+    print('*' * 30)
+    for enum, I_app in enumerate([-1.0, 1.0]):
+        I_app *= I_mag
+        # load model
+        model = pybamm.lithium_ion.SPMe()
+        # create geometry
+        geometry = model.default_geometry
+        # load parameter values and process model and geometry
+        param = model.default_parameter_values
+        param.update(
+            {"Electrode height [m]": h,
+             "Electrode width [m]": w,
+             "Negative electrode thickness [m]": l_n,
+             "Positive electrode thickness [m]": l_p,
+             "Separator thickness [m]": l_s,
+             "Lower voltage cut-off [V]": 2.8,
+             "Upper voltage cut-off [V]": 4.7,
+             "Maximum concentration in negative electrode [mol.m-3]": 25000,
+             "Maximum concentration in positive electrode [mol.m-3]": 50000,
+             "Initial concentration in negative electrode [mol.m-3]": 12500,
+             "Initial concentration in positive electrode [mol.m-3]": 25000,
+             "Negative electrode surface area density [m-1]": 180000.0,
+             "Positive electrode surface area density [m-1]": 150000.0,
+             "Typical current [A]": I_app,
+             }
+        )
+        param.process_model(model)
+        param.process_geometry(geometry)
+        s_var = pybamm.standard_spatial_vars
+        var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5,
+                   s_var.r_n: 5, s_var.r_p: 10}
+        # set mesh
+        mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
+        # discretise model
+        disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+        disc.process_model(model)
+        # solve model
+        t_eval = np.linspace(0, 0.2, 100)
+        sol = model.default_solver.solve(model, t_eval)
+        var = "Positive electrode average extent of lithiation"
+        xpext = pybamm.ProcessedVariable(model.variables[var],
+                                         sol.t, sol.y, mesh=mesh)
+        var = "Negative electrode average extent of lithiation"
+        xnext = pybamm.ProcessedVariable(model.variables[var],
+                                         sol.t, sol.y, mesh=mesh)
+        var = "X-averaged positive particle surface concentration"
+        xpsurf = pybamm.ProcessedVariable(model.variables[var],
+                                          sol.t, sol.y, mesh=mesh)
+        var = "X-averaged negative particle surface concentration"
+        xnsurf = pybamm.ProcessedVariable(model.variables[var],
+                                          sol.t, sol.y, mesh=mesh)
+        time = pybamm.ProcessedVariable(model.variables["Time [h]"],
+                                        sol.t, sol.y, mesh=mesh)
+        # Coulomb counting
+        time_hours = time(sol.t)
+        dc_time = np.around(time_hours[-1], 3)
+        # Capacity mAh
+        cap = np.absolute(I_app * 1000 * dc_time)
+        cap_time = np.absolute(I_app * 1000 * time_hours)
+        axes[enum].plot(cap_time,
+                        xnext(sol.t), 'r-', label='Average Neg')
+        axes[enum].plot(cap_time,
+                        xpext(sol.t), 'b-', label='Average Pos')
+        axes[enum].plot(cap_time,
+                        xnsurf(sol.t), 'r--', label='Surface Neg')
+        axes[enum].plot(cap_time,
+                        xpsurf(sol.t), 'b--', label='Surface Pos')
+        axes[enum].set_xlabel('Capacity [mAh]')
+        handles, labels = axes[enum].get_legend_handles_labels()
+        axes[enum].legend(handles, labels)
+        if I_app < 0.0:
+            axes[enum].set_ylabel('Extent of Lithiation, Elecrode Ratio: '
+                                  + str(e_ratio))
+            axes[enum].title.set_text('Charge')
+        else:
+            axes[enum].title.set_text('Discharge')
+        print('Applied Current', I_app, 'A', 'Time',
+              dc_time, 'hrs', 'Capacity', cap, 'mAh')
+        tot_cap += cap
+        tot_time += dc_time
 
-    axes[enum].plot(cap_time,
-                    xnext(sol.t), 'r-', label='Average Neg')
-    axes[enum].plot(cap_time,
-                    xpext(sol.t), 'b-', label='Average Pos')
-    axes[enum].plot(cap_time,
-                    xnsurf(sol.t), 'r--', label='Surface Neg')
-    axes[enum].plot(cap_time,
-                    xpsurf(sol.t), 'b--', label='Surface Pos')
-    axes[enum].set_xlabel('Capacity [mAh]')
-    plt.legend()
-    if I_app < 0.0:
-        axes[enum].set_ylabel('Extent of Lithiation')
-        axes[enum].title.set_text('Charge')
-    else:
-        axes[enum].title.set_text('Discharge')
-    print('Applied Current', I_app, 'A', 'Time',
-          dc_time, 'hrs', 'Capacity', cap, 'mAh')
-    tot_cap += cap
-    tot_time += dc_time
+    print('Anode : Cathode thickness', e_ratio)
+    print('Total Charge/Discharge Time', tot_time, 'hrs')
+    print('Total Capacity', np.around(tot_cap, 3), 'mAh')
+    specific_cap = np.around(tot_cap, 3) / vol_cm3
+    print('Total Capacity', specific_cap, 'mAh.cm-3')
+    capacities.append(tot_cap)
+    specific_capacities.append(specific_cap)
 
-print('Total Charge/Discharge Time', tot_time, 'hrs')
-print('Total Capacity', np.around(tot_cap, 3), 'mAh')
-print('Total Capacity', np.around(tot_cap, 3) / vol_cm3, 'mAh.cm-3')
+fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
+ax1.plot(thicknesses / l_p, capacities)
+ax2.plot(thicknesses / l_p, specific_capacities)
+ax1.set_ylabel('Capacity [mAh]')
+ax2.set_ylabel('Specific Capacity [mAh.cm-3]')
+ax2.set_xlabel('Anode : Cathode thickness')
diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py
index 491c607f02..ab3e5c7f0f 100644
--- a/pybamm/models/submodels/particle/base_particle.py
+++ b/pybamm/models/submodels/particle/base_particle.py
@@ -35,7 +35,7 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             c_scale = self.param.c_p_max
             active_volume = geo_param.a_p_dim * geo_param.R_p / 3
         c_s_r_av = pybamm.r_average(c_s_xav)
-        c_s_r_av_vol = active_volume * c_s_r_av * c_scale
+        c_s_r_av_vol = active_volume * c_s_r_av
         variables = {
             self.domain + " particle concentration": c_s,
             self.domain + " particle concentration [mol.m-3]": c_s * c_scale,
@@ -54,7 +54,9 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             + " particle surface concentration [mol.m-3]": c_scale * c_s_surf_av,
             self.domain + " electrode active volume fraction": active_volume,
             self.domain
-            + " electrode volume-averaged concentration [mol.m-3]": c_s_r_av_vol,
+            + " electrode volume-averaged concentration": c_s_r_av_vol,
+            self.domain + " electrode "
+            + "volume-averaged concentration [mol.m-3]": c_s_r_av_vol * c_scale,
             self.domain + " electrode average extent of lithiation": c_s_r_av
         }
 
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index eb4a93e224..9bd9608044 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -64,7 +64,7 @@ def solve(self, model, t_eval):
         solution.total_time = timer.time() - start_time
         solution.set_up_time = set_up_time
 
-        pybamm.logger.info("Finish solving {} ({})".format(model.name, termination))
+        pybamm.logger.warning("Finish solving {} ({})".format(model.name, termination))
         pybamm.logger.info(
             "Set-up time: {}, Solve time: {}, Total time: {}".format(
                 timer.format(solution.set_up_time),