Repository: pybamm-team/PyBaMM
Pull request #542 (2019-07-19 14:33:58)

Thesis discharge results
# Description

Scripts for generating results
Hopefully these can be adapted to help with #349 

## Type of change

- [x] New feature (non-breaking change which adds functionality)

# Key checklist:

- [x] No style issues: `$ flake8`
- [x] All tests pass: `$ python run-tests.py --unit`
- [x] The documentation builds: `$ cd docs` and then `$ make clean; make html`

You can run all three at once, using `$ python run-tests.py --quick`.

## Further checks: 

- [x] Code is commented, particularly in hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works
- [x] Any dependent changes have been merged and published in downstream modules



Patch:
diff --git a/docs/index.rst b/docs/index.rst
index bc67cfd11..ca6bcbfd9 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -32,6 +32,7 @@ Contents
     source/spatial_methods/index
     source/solvers/index
     source/processed_variable
+    source/util
 
 Examples
 ========
diff --git a/docs/source/util.rst b/docs/source/util.rst
new file mode 100644
index 000000000..db8ab248c
--- /dev/null
+++ b/docs/source/util.rst
@@ -0,0 +1,13 @@
+Utility functions
+=================
+
+.. autofunction:: pybamm.get_infinite_nested_dict
+
+.. autofunction:: pybamm.load_function
+
+.. autofunction:: pybamm.rmse
+
+.. autofunction:: pybamm.root_dir
+
+.. autoclass:: pybamm.Timer
+  :members:
diff --git a/examples/notebooks/models/lead-acid.ipynb b/examples/notebooks/models/lead-acid.ipynb
index 438c4e2e7..20455e1fe 100644
--- a/examples/notebooks/models/lead-acid.ipynb
+++ b/examples/notebooks/models/lead-acid.ipynb
@@ -321,7 +321,7 @@
     "    t, y = solutions[model].t, solutions[model].y\n",
     "    time = pybamm.ProcessedVariable(model.variables[\"Time [h]\"], t, y)(t)\n",
     "    voltage = pybamm.ProcessedVariable(model.variables[\"Terminal voltage [V]\"], t, y)\n",
-    "    plt.plot(time, voltage(t) * 6, lw=2, label=model.name)\n",
+    "    plt.plot(time, voltage(t), lw=2, label=model.name)\n",
     "plt.xlabel(\"Time [h]\", fontsize=15)\n",
     "plt.ylabel(\"Terminal voltage [V]\", fontsize=15)\n",
     "plt.legend(fontsize=15)\n",
diff --git a/examples/scripts/compare_lead_acid.py b/examples/scripts/compare_lead_acid.py
index bbc2be40d..973c8e579 100644
--- a/examples/scripts/compare_lead_acid.py
+++ b/examples/scripts/compare_lead_acid.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pybamm
 
-pybamm.set_logging_level("INFO")
+pybamm.set_logging_level("DEBUG")
 
 # load models
 models = [
@@ -19,11 +19,14 @@ geometry = models[-1].default_geometry
 param = models[0].default_parameter_values
 param.update(
     {
-        "Typical current [A]": 20,
+        # "Bruggeman coefficient": 0.001,
+        # "Current function": pybamm.GetConstantCurrent(current=0),
+        "Typical current [A]": -20,
         "Initial State of Charge": 1,
         "Typical electrolyte concentration [mol.m-3]": 5600,
         "Negative electrode reference exchange-current density [A.m-2]": 0.08,
         "Positive electrode reference exchange-current density [A.m-2]": 0.006,
+        "Positive electrode reference exchange-current density (oxygen) [A.m-2]": 1e-22,
     }
 )
 for model in models:
@@ -42,7 +45,7 @@ for model in models:
 
 # solve model
 solutions = [None] * len(models)
-t_eval = np.linspace(0, 1, 100)
+t_eval = np.linspace(0, 2, 100)
 for i, model in enumerate(models):
     solution = model.default_solver.solve(model, t_eval)
     solutions[i] = solution
@@ -53,10 +56,10 @@ output_variables = [
         "Average negative electrode interfacial current density [A.m-2]",
         "Average positive electrode interfacial current density [A.m-2]",
     ],
-    "Average negative electrode surface potential difference [V]",
-    "Average positive electrode surface potential difference [V]",
-    "Electrolyte concentration",
-    "Electrolyte flux",
+    "Average negative electrode reaction overpotential [V]",
+    "Average positive electrode reaction overpotential [V]",
+    "State of Charge",
+    "Oxygen concentration [mol.m-3]",
     "Terminal voltage [V]",
 ]
 plot = pybamm.QuickPlot(models, mesh, solutions, output_variables)
diff --git a/input/parameters/lead-acid/default.csv b/input/parameters/lead-acid/default.csv
index b414add87..5f6d2951c 100644
--- a/input/parameters/lead-acid/default.csv
+++ b/input/parameters/lead-acid/default.csv
@@ -17,8 +17,8 @@ Positive tab centre z-coordinate [m],0.114,Tab at top,
 Cell capacity [A.h],17,Manufacturer,
 Number of electrodes connected in parallel to make a cell,8,Manufacturer,
 Number of cells connected in series to make a battery,6,Manufacturer,
-Lower voltage cut-off [V],1.75,,10.5V across 6-cell battery
-Upper voltage cut-off [V],2.42,,14.5V across 6-cell battery
+Lower voltage cut-off [V],1.73,,(just under) 10.5V across 6-cell battery
+Upper voltage cut-off [V],2.44,,(just over) 14.5V across 6-cell battery
 Initial State of Charge,1,-,
 ,,,
 # Electrolyte properties,,,
diff --git a/input/parameters/lithium-ion/mcmb2528_lif6-in-ecdmc_lico2_parameters_Dualfoil.csv b/input/parameters/lithium-ion/mcmb2528_lif6-in-ecdmc_lico2_parameters_Dualfoil.csv
index 22ad7fb6c..7994323d1 100644
--- a/input/parameters/lithium-ion/mcmb2528_lif6-in-ecdmc_lico2_parameters_Dualfoil.csv
+++ b/input/parameters/lithium-ion/mcmb2528_lif6-in-ecdmc_lico2_parameters_Dualfoil.csv
@@ -23,6 +23,7 @@ Positive tab centre z-coordinate [m],0.137,,Need to actual value for KOKAM cell
 ,,,
 # Electrical,,,
 Number of electrodes connected in parallel to make a cell,1,-,
+Number of cells connected in series to make a battery,1,-,
 Lower voltage cut-off [V],3.105,-,
 Upper voltage cut-off [V],4.7,-,
 ,,,
@@ -66,11 +67,11 @@ Charge transfer coefficient,0.5,Scott Moura FastDFN,
 Double-layer capacity [F.m-2],0.2,,
 ,,,
 # Density,,,
-Negative current collector density [kg.m-3],8954,, 
-Negative electrode density [kg.m-3],1657,, 
-Separator density [kg.m-3],397,, 
-Positive electrode density [kg.m-3],3262,, 
-Positive current collector density [kg.m-3],2707,, 
+Negative current collector density [kg.m-3],8954,,
+Negative electrode density [kg.m-3],1657,,
+Separator density [kg.m-3],397,,
+Positive electrode density [kg.m-3],3262,,
+Positive current collector density [kg.m-3],2707,,
 ,,,
 # Specific heat capacity,,,
 Negative current collector specific heat capacity [J.kg-1.K-1],385,,
@@ -80,15 +81,15 @@ Negative electrode specific heat capacity [J.kg-1.K-1],700,,
 Positive current collector specific heat capacity [J.kg-1.K-1],897,,
 ,,,
 # Thermal conductivity,,,
-Negative current collector thermal conductivity [W.m-1.K-1],401,, 
-Negative electrode thermal conductivity [W.m-1.K-1],1.7,, 
-Separator thermal conductivity [W.m-1.K-1],0.16,, 
-Positive electrode thermal conductivity [W.m-1.K-1],2.1,, 
-Positive current collector thermal conductivity [W.m-1.K-1],237,, 
+Negative current collector thermal conductivity [W.m-1.K-1],401,,
+Negative electrode thermal conductivity [W.m-1.K-1],1.7,,
+Separator thermal conductivity [W.m-1.K-1],0.16,,
+Positive electrode thermal conductivity [W.m-1.K-1],2.1,,
+Positive current collector thermal conductivity [W.m-1.K-1],237,,
 ,,,
 # Thermal parameters,,,
 Heat transfer coefficient [W.m-2.K-1],10,,
-Typical temperature variation [K],2.4,, 
+Typical temperature variation [K],2.4,,
 Lumped effective thermal density [J.K-1.m-3],1811600,,
 Effective thermal conductivity [W.m-1.K-1],59.3964,,
 ,,,
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index a51de97bb..c4178f7d5 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -41,9 +41,8 @@ ABSOLUTE_PATH = os.path.join(os.path.split(script_path)[0], "..")
 #
 # Utility classes and methods
 #
-from .util import root_dir
 from .util import Timer
-from .util import load_function
+from .util import root_dir, load_function, rmse, get_infinite_nested_dict
 from .logger import logger, set_logging_level
 from .settings import settings
 
@@ -235,7 +234,7 @@ from .parameters.standard_current_functions.get_current_data import GetCurrentDa
 # other
 #
 from .processed_variable import post_process_variables, ProcessedVariable
-from .quick_plot import QuickPlot
+from .quick_plot import QuickPlot, ax_min, ax_max
 
 #
 # Remove any imported modules, so we don't expose them as part of pybamm
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index f90098f29..90b8396fd 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -97,6 +97,11 @@ class Discretisation(object):
         else:
             # create a blank model so that original model is unchanged
             model_disc = pybamm.BaseModel()
+            model_disc.name = model.name
+            model_disc.options = model.options
+            model_disc.use_jacobian = model.use_jacobian
+            model_disc.use_simplify = model.use_simplify
+            model_disc.use_to_python = model.use_to_python
 
         model_disc.bcs = self.bcs
 
@@ -123,7 +128,7 @@ class Discretisation(object):
         model_disc.events = processed_events
 
         # Create mass matrix
-        model_disc.mass_matrix = self.create_mass_matrix(model)
+        model_disc.mass_matrix = self.create_mass_matrix(model_disc)
 
         # Check that resulting model makes sense
         self.check_model(model_disc)
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 86b9f3611..ec9e1d7af 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -40,6 +40,7 @@ class BaseModel(object):
 
     def __init__(self, name="Unnamed model"):
         self.name = name
+        self.options = {}
 
         # Initialise empty model
         self._rhs = {}
@@ -192,6 +193,14 @@ class BaseModel(object):
     def set_of_parameters(self):
         return self._set_of_parameters
 
+    @property
+    def options(self):
+        return self._options
+
+    @options.setter
+    def options(self, options):
+        self._options = options
+
     def __getitem__(self, key):
         return self.rhs[key]
 
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 73a03f9be..00f287cab 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -355,6 +355,9 @@ class BaseBatteryModel(pybamm.BaseModel):
         pybamm.logger.debug("Setting voltage variables")
         self.set_voltage_variables()
 
+        pybamm.logger.debug("Setting SoC variables")
+        self.set_soc_variables()
+
         self._built = True
 
     def set_thermal_submodel(self):
@@ -469,9 +472,39 @@ class BaseBatteryModel(pybamm.BaseModel):
             }
         )
 
+        # Battery-wide variables
+        eta_e_av_dim = self.variables.get("Average electrolyte ohmic losses [V]", 0)
+        eta_c_av_dim = self.variables.get("Average concentration overpotential [V]", 0)
+        num_cells = pybamm.Parameter(
+            "Number of cells connected in series to make a battery"
+        )
+
+        self.variables.update(
+            {
+                "Average battery open circuit voltage [V]": ocv_av_dim * num_cells,
+                "Measured battery open circuit voltage [V]": ocv_dim * num_cells,
+                "Average battery reaction overpotential [V]": eta_r_av_dim * num_cells,
+                "Average battery solid phase ohmic losses [V]": delta_phi_s_av_dim
+                * num_cells,
+                "Average battery electrolyte ohmic losses [V]": eta_e_av_dim
+                * num_cells,
+                "Average battery concentration overpotential [V]": eta_c_av_dim
+                * num_cells,
+                "Battery voltage [V]": V_dim * num_cells,
+            }
+        )
+
         # Cut-off voltage
         voltage = self.variables["Terminal voltage"]
         self.events["Minimum voltage"] = voltage - self.param.voltage_low_cut
+        self.events["Maximum voltage"] = voltage - self.param.voltage_high_cut
+
+    def set_soc_variables(self):
+        """
+        Set variables relating to the state of charge.
+        This function is overriden by the base battery models
+        """
+        pass
 
     def process_parameters_and_discretise(self, symbol):
         """
diff --git a/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py b/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py
index e6f77dfa7..85485dc17 100644
--- a/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py
+++ b/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py
@@ -109,3 +109,9 @@ class BaseModel(pybamm.BaseBatteryModel):
             }
             self.reactions["main"]["Negative"]["s_ox"] = 0
             self.reactions["main"]["Positive"]["s_ox"] = 0
+
+    def set_soc_variables(self):
+        "Set variables relating to the state of charge."
+        # State of Charge defined as function of dimensionless electrolyte concentration
+        soc = self.variables["Average electrolyte concentration"] * 100
+        self.variables.update({"State of Charge": soc, "Depth of Discharge": 100 - soc})
diff --git a/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py b/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py
index 96db3c8bf..ed47a52a3 100644
--- a/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py
@@ -48,14 +48,19 @@ class BaseElectrolyteDiffusion(pybamm.BaseSubModel):
         variables = {
             "Electrolyte concentration": c_e,
             "Electrolyte concentration [mol.m-3]": c_e_typ * c_e,
+            "Electrolyte concentration [Molar]": c_e_typ * c_e / 1000,
             "Average electrolyte concentration": c_e_av,
             "Average electrolyte concentration [mol.m-3]": c_e_typ * c_e_av,
+            "Average electrolyte concentration [Molar]": c_e_typ * c_e_av / 1000,
             "Negative electrolyte concentration": c_e_n,
             "Negative electrolyte concentration [mol.m-3]": c_e_typ * c_e_n,
+            "Negative electrolyte concentration [Molar]": c_e_typ * c_e_n / 1000,
             "Separator electrolyte concentration": c_e_s,
             "Separator electrolyte concentration [mol.m-3]": c_e_typ * c_e_s,
+            "Separator electrolyte concentration [Molar]": c_e_typ * c_e_s / 1000,
             "Positive electrolyte concentration": c_e_p,
             "Positive electrolyte concentration [mol.m-3]": c_e_typ * c_e_p,
+            "Positive electrolyte concentration [Molar]": c_e_typ * c_e_p / 1000,
             "Average negative electrolyte concentration": c_e_n_av,
             "Average negative electrolyte concentration [mol.m-3]": c_e_typ * c_e_n_av,
             "Average separator electrolyte concentration": c_e_s_av,
diff --git a/pybamm/models/submodels/interface/base_interface.py b/pybamm/models/submodels/interface/base_interface.py
index ec9553aac..6c679320e 100644
--- a/pybamm/models/submodels/interface/base_interface.py
+++ b/pybamm/models/submodels/interface/base_interface.py
@@ -48,10 +48,12 @@ class BaseInterface(pybamm.BaseSubModel):
 
     def _get_standard_interfacial_current_variables(self, j):
 
+        i_typ = self.param.i_typ
+        L_x = self.param.L_x
         if self.domain == "Negative":
-            j_scale = self.param.i_typ / (self.param.a_n_dim * self.param.L_x)
+            j_scale = i_typ / (self.param.a_n_dim * L_x)
         elif self.domain == "Positive":
-            j_scale = self.param.i_typ / (self.param.a_p_dim * self.param.L_x)
+            j_scale = i_typ / (self.param.a_p_dim * L_x)
 
         # Average, and broadcast if necessary
         j_av = pybamm.average(j)
@@ -71,6 +73,7 @@ class BaseInterface(pybamm.BaseSubModel):
             + self.reaction_name
             + " interfacial current density": j_av,
             self.domain
+            + " electrode"
             + self.reaction_name
             + " interfacial current density [A.m-2]": j_scale * j,
             "Average "
@@ -78,16 +81,27 @@ class BaseInterface(pybamm.BaseSubModel):
             + " electrode"
             + self.reaction_name
             + " interfacial current density [A.m-2]": j_scale * j_av,
+            self.domain
+            + " electrode"
+            + self.reaction_name
+            + " interfacial current density per volume [A.m-3]": i_typ / L_x * j,
+            "Average "
+            + self.domain.lower()
+            + " electrode"
+            + self.reaction_name
+            + " interfacial current density per volume [A.m-3]": i_typ / L_x * j_av,
         }
 
         return variables
 
     def _get_standard_total_interfacial_current_variables(self, j_tot_av):
 
+        i_typ = self.param.i_typ
+        L_x = self.param.L_x
         if self.domain == "Negative":
-            j_scale = self.param.i_typ / (self.param.a_n_dim * self.param.L_x)
+            j_scale = i_typ / (self.param.a_n_dim * L_x)
         elif self.domain == "Positive":
-            j_scale = self.param.i_typ / (self.param.a_p_dim * self.param.L_x)
+            j_scale = i_typ / (self.param.a_p_dim * L_x)
 
         variables = {
             "Average "
@@ -97,14 +111,21 @@ class BaseInterface(pybamm.BaseSubModel):
             + self.domain.lower()
             + " electrode total interfacial current density [A.m-2]": j_scale
             * j_tot_av,
+            "Average "
+            + self.domain.lower()
+            + " electrode total interfacial current density per volume [A.m-3]": i_typ
+            / L_x
+            * j_tot_av,
         }
 
         return variables
 
     def _get_standard_whole_cell_interfacial_current_variables(self, variables):
 
-        j_n_scale = self.param.i_typ / (self.param.a_n_dim * self.param.L_x)
-        j_p_scale = self.param.i_typ / (self.param.a_p_dim * self.param.L_x)
+        i_typ = self.param.i_typ
+        L_x = self.param.L_x
+        j_n_scale = i_typ / (self.param.a_n_dim * L_x)
+        j_p_scale = i_typ / (self.param.a_p_dim * L_x)
 
         j_n = variables[
             "Negative electrode" + self.reaction_name + " interfacial current density"
@@ -120,22 +141,27 @@ class BaseInterface(pybamm.BaseSubModel):
             variables = {
                 "Interfacial current density": j,
                 "Interfacial current density [A.m-2]": j_dim,
+                "Interfacial current density per volume [A.m-3]": i_typ / L_x * j,
             }
         else:
             reaction_name = self.reaction_name[1:].capitalize()
             variables = {
                 reaction_name + " interfacial current density": j,
                 reaction_name + " interfacial current density [A.m-2]": j_dim,
+                reaction_name
+                + " interfacial current density per volume [A.m-3]": i_typ / L_x * j,
             }
 
         return variables
 
     def _get_standard_exchange_current_variables(self, j0):
 
+        i_typ = self.param.i_typ
+        L_x = self.param.L_x
         if self.domain == "Negative":
-            j_scale = self.param.i_typ / (self.param.a_n_dim * self.param.L_x)
+            j_scale = i_typ / (self.param.a_n_dim * L_x)
         elif self.domain == "Positive":
-            j_scale = self.param.i_typ / (self.param.a_p_dim * self.param.L_x)
+            j_scale = i_typ / (self.param.a_p_dim * L_x)
 
         # Average, and broadcast if necessary
         j0_av = pybamm.average(j0)
@@ -165,14 +191,25 @@ class BaseInterface(pybamm.BaseSubModel):
             + " electrode"
             + self.reaction_name
             + " exchange current density [A.m-2]": j_scale * j0_av,
+            self.domain
+            + " electrode"
+            + self.reaction_name
+            + " exchange current density per volume [A.m-3]": i_typ / L_x * j0,
+            "Average "
+            + self.domain.lower()
+            + " electrode"
+            + self.reaction_name
+            + " exchange current density per volume [A.m-3]": i_typ / L_x * j0_av,
         }
 
         return variables
 
     def _get_standard_whole_cell_exchange_current_variables(self, variables):
 
-        j_n_scale = self.param.i_typ / (self.param.a_n_dim * self.param.L_x)
-        j_p_scale = self.param.i_typ / (self.param.a_p_dim * self.param.L_x)
+        i_typ = self.param.i_typ
+        L_x = self.param.L_x
+        j_n_scale = i_typ / (self.param.a_n_dim * L_x)
+        j_p_scale = i_typ / (self.param.a_p_dim * L_x)
 
         j0_n = variables[
             "Negative electrode" + self.reaction_name + " exchange current density"
@@ -188,12 +225,15 @@ class BaseInterface(pybamm.BaseSubModel):
             variables = {
                 "Exchange current density": j0,
                 "Exchange current density [A.m-2]": j0_dim,
+                "Exchange current density per volume [A.m-3]": i_typ / L_x * j0,
             }
         else:
             reaction_name = self.reaction_name[1:].capitalize()
             variables = {
                 reaction_name + " exchange current density": j0,
                 reaction_name + " exchange current density [A.m-2]": j0_dim,
+                reaction_name
+                + " exchange current density per volume [A.m-3]": i_typ / L_x * j0,
             }
 
         return variables
diff --git a/pybamm/parameters/standard_parameters_lead_acid.py b/pybamm/parameters/standard_parameters_lead_acid.py
index 31d1b005a..fec141b03 100644
--- a/pybamm/parameters/standard_parameters_lead_acid.py
+++ b/pybamm/parameters/standard_parameters_lead_acid.py
@@ -96,7 +96,6 @@ s_plus_Ox_dim = pybamm.Parameter("Signed stoichiometry of cations (oxygen reacti
 s_w_Ox_dim = pybamm.Parameter("Signed stoichiometry of water (oxygen reaction)")
 s_ox_Ox_dim = pybamm.Parameter("Signed stoichiometry of oxygen (oxygen reaction)")
 ne_Ox = pybamm.Parameter("Electrons in oxygen reaction")
-c_ox_ref = pybamm.Parameter("Reference oxygen molecule concentration [mol.m-3]")
 U_Ox_dim = pybamm.Parameter("Oxygen reference OCP vs SHE [V]")
 # Hydrogen
 j0_n_Hy_ref_dimensional = pybamm.Parameter(
diff --git a/pybamm/util.py b/pybamm/util.py
index da7459e35..094fc8ca7 100644
--- a/pybamm/util.py
+++ b/pybamm/util.py
@@ -5,11 +5,13 @@
 # (see https://github.com/pints-team/pints)
 #
 import importlib
+import numpy as np
 import os
 import sys
 import timeit
 import pathlib
 import pybamm
+from collections import defaultdict
 
 
 def root_dir():
@@ -151,3 +153,33 @@ def load_function(filename):
         )
 
     return getattr(module_object, valid_module)
+
+
+def rmse(x, y):
+    "Calculate the root-mean-square-error between two vectors x and y, ignoring NaNs"
+    # Check lengths
+    if len(x) != len(y):
+        raise ValueError("Vectors must have the same length")
+    return np.sqrt(np.nanmean((x - y) ** 2))
+
+
+def get_infinite_nested_dict():
+    """
+    Return a dictionary that allows infinite nesting without having to define level by
+    level.
+
+    See:
+    https://stackoverflow.com/questions/651794/whats-the-best-way-to-initialize-a-dict-of-dicts-in-python/652226#652226
+
+    Example
+    -------
+    >>> import pybamm
+    >>> d = pybamm.get_infinite_nested_dict()
+    >>> d["a"] = 1
+    >>> d["a"]
+    1
+    >>> d["b"]["c"]["d"] = 2
+    >>> d["b"]["c"] == {"d": 2}
+    True
+    """
+    return defaultdict(get_infinite_nested_dict)
diff --git a/results/2019_08_sulzer_thesis/README.md b/results/2019_08_sulzer_thesis/README.md
index 91c12fb51..46300fae7 100644
--- a/results/2019_08_sulzer_thesis/README.md
+++ b/results/2019_08_sulzer_thesis/README.md
@@ -1,6 +1,7 @@
 # Results: Valentin Sulzer's Thesis
 
-This folder contains the scripts used to generate results for Valentin Sulzers thesis
+This folder contains the scripts used to generate results for Valentin Sulzer's thesis
+The plots were formatted using a formatting file `matplotlibrc` identical to [this one](_matplotlibrc) (but not included in the GitHub repo to avoid clashes with different formatting files).
 
 ## Chapter 2 - Model
 
@@ -11,6 +12,39 @@ This folder contains the scripts used to generate results for Valentin Sulzers t
 
 - [Effect of capacitance](effect_of_capacitance.py): comparison of the one-dimensional model with and without capacitance terms included
     - Voltages
+    - Errors
     - Time taken for full solution
-- [Discharge plots](lead_acid_discharge.py):
-    - Comparison of voltage curves for the full porous electrode model and a hierarchy of simplified models (Leading-Order Quasi-Static and Composite)
+- [Discharge asymptotics results](lead_acid_discharge.py):
+    - Comparison of voltage curves for the full porous electrode model and a hierarchy of simplified models (Leading-Order Quasi-Static, First-Order Quasi-Static and Composite)
+    - Comparison of variable profiles at various times for the full and reduced-order models
+        - Electrolyte concentration
+        - Electrolyte potential
+        - Interfacial current density
+    - Errors compared to full model and time taken to solve each model
+    - Breakdown of voltage into constituent overpotentials
+- Effect of convection (to do):
+    - Voltage at various C-rates with and without convection
+    - Velocity profiles
+    - Increasing the volume changes to see more of an effect
+- [Effect of side reactions](effect_of_side_reactions.py):
+    - Voltage at various C-rates with and without side reactions
+    - Interfacial current densities
+    - Electrolyte concentrations
+- [Charge asymptotics results](lead_acid_charge.py):
+    - Comparison of voltage curves for the full porous electrode model and a hierarchy of simplified models (Leading-Order Quasi-Static, First-Order Quasi-Static and Composite)
+    - Comparison of average interfacial current densities for each of the side reactions for the full porous electrode model and a hierarchy of simplified models (Leading-Order Quasi-Static, First-Order Quasi-Static and Composite)
+    - Comparison of variable profiles at various times for the full and reduced-order models
+        - Electrolyte concentration
+        - Oxygen concentration
+    - Errors compared to full model and time taken to solve each model
+    - Breakdown of voltage into constituent overpotentials
+- [Self-discharge](self_discharge.py):
+    - Self-discharge voltages
+
+## Chapter 4 - Small aspect ratio cells
+
+- 2+1D model
+    - Model and capacitance formulation
+    - Concentrations and potentials as functions of x, y, z
+    - Times taken
+- Further asymptotics
diff --git a/results/2019_08_sulzer_thesis/_matplotlibrc b/results/2019_08_sulzer_thesis/_matplotlibrc
new file mode 100644
index 000000000..7e8738810
--- /dev/null
+++ b/results/2019_08_sulzer_thesis/_matplotlibrc
@@ -0,0 +1,27 @@
+# Set a 2/1 aspect ratio with automatic layout.
+figure.figsize: 6.4, 4.
+
+# Use LaTeX fonts
+font.family: serif
+font.serif: ComputerModern
+text.usetex: true
+
+# Backend options (42: use TrueType fonts)
+pdf.fonttype: 42
+ps.fonttype: 42
+
+# Set large font sizes for paper figures
+axes.labelsize: 11
+axes.titlesize: 11
+xtick.labelsize: 11
+ytick.labelsize: 11
+legend.fontsize: 11
+legend.numpoints: 1
+legend.scatterpoints: 1
+
+# Set lines and ticks according to the `seaborn-paper` style sheet
+grid.linewidth: 0.8
+lines.linewidth: 1.4
+patch.linewidth: 0.24
+lines.markersize: 5.6
+lines.markeredgewidth: 0
diff --git a/results/2019_08_sulzer_thesis/effect_of_capacitance.py b/results/2019_08_sulzer_thesis/effect_of_capacitance.py
index 942912e46..826205675 100644
--- a/results/2019_08_sulzer_thesis/effect_of_capacitance.py
+++ b/results/2019_08_sulzer_thesis/effect_of_capacitance.py
@@ -6,128 +6,153 @@ import matplotlib.pyplot as plt
 import numpy as np
 import pickle
 import pybamm
+import shared_plotting
 from config import OUTPUT_DIR
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
-from shared import model_comparison
+from shared_solutions import model_comparison, convergence_study
 
+save_folder = "results/2019_08_sulzer_thesis/data/capacitance_results/"
 
-def plot_voltages(all_variables, t_eval, Crates):
-    # Only use some Crates
-    all_variables = {k: v for k, v in all_variables.items() if k in Crates}
-    # Plot
-    plt.subplots()
-    n = int(len(all_variables) // np.sqrt(len(all_variables)))
-    m = np.ceil(len(all_variables) / n)
+
+def plot_voltages(all_variables, t_eval):
+    linestyles = ["k-", "b-.", "r--"]
+    _, axes = shared_plotting.plot_voltages(
+        all_variables, t_eval, linestyles=linestyles
+    )
+
+    # Add inset plot
     for k, (Crate, models_variables) in enumerate(all_variables.items()):
-        t_max = max(
-            np.nanmax(var["Time [s]"](t_eval)) for var in models_variables.values()
-        )
-        ax = plt.subplot(n, m, k + 1)
-        plt.axis([0, t_max, 10.5, 13])
-        ax.set_xlabel("Time [s]")
-        if len(Crates) > 1:
-            plt.title("\\textbf{{({})}} {} C".format(chr(97 + k), Crate), y=-0.4)
-
-        # Hide the right and top spines
-        ax.spines["right"].set_visible(False)
-        ax.spines["top"].set_visible(False)
-
-        # Only show ticks on the left and bottom spines
-        ax.yaxis.set_ticks_position("left")
-        ax.xaxis.set_ticks_position("bottom")
-
-        # Add inset plot
-        inset = inset_axes(ax, width="30%", height="30%", loc=3, borderpad=3)
-
-        # Linestyles
-        linestyles = ["k-", "b-.", "r--"]
-        for j, (model_name, variables) in enumerate(models_variables.items()):
-            if k == 0:
-                label = model_name
-            else:
-                label = None
-            if k % m == 0:
-                ax.set_ylabel("Voltage [V]")
-
-            ax.plot(
-                variables["Time [s]"](t_eval),
-                variables["Terminal voltage [V]"](t_eval) * 6,
-                linestyles[j],
-                label=label,
-            )
-            inset.plot(
-                variables["Time [s]"](t_eval[:40]),
-                variables["Terminal voltage [V]"](t_eval[:40]) * 6,
-                linestyles[j],
-            )
-        # plt.legend(loc="upper right")
-    file_name = "capacitance_voltage_comparison.eps".format(Crate)
-    plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000)
+        ax = axes.flat[k]
+        y_min = ax.get_ylim()[0]
+        ax.set_ylim([y_min, 13.6])
+        inset = inset_axes(ax, width="40%", height="40%", loc=1, borderpad=0)
+        for j, variables in enumerate(models_variables.values()):
+            time = variables["Time [s]"](t_eval)
+            capacitance_indices = np.where(time < 50)
+            time = time[capacitance_indices]
+            voltage = variables["Battery voltage [V]"](t_eval)[capacitance_indices]
+            inset.plot(time, voltage, linestyles[j])
+            inset.set_xlabel("Time [s]", fontsize=10)
+            inset.set_xlim([0, 3])
+
+    file_name = "capacitance_voltage_comparison.eps"
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight")
 
 
 def plot_errors(all_variables, t_eval, Crates):
-    def rmse(predictions, targets):
-        return np.sqrt(np.nanmean((predictions - targets) ** 2))
-
+    # Linestyles
+    linestyles = ["k-", "b-.", "r--"]
     # Only use some Crates
     all_variables = {k: v for k, v in all_variables.items() if k in Crates}
     # Plot
-    plt.subplots()
+    fig, ax = plt.subplots(1, 1)
     for k, (Crate, models_variables) in enumerate(all_variables.items()):
-        ax = plt.subplot(1, 1, 1)
-        ax.set_xlabel("Time [s]")
+        ax.set_xlabel("Time [h]")
         ax.set_ylabel("Error [V]")
 
-        # Hide the right and top spines
-        ax.spines["right"].set_visible(False)
-        ax.spines["top"].set_visible(False)
-
-        # Only show ticks on the left and bottom spines
-        ax.yaxis.set_ticks_position("left")
-        ax.xaxis.set_ticks_position("bottom")
-
-        # Linestyles
-        linestyles = ["k-", "b-.", "r--"]
-
         for j, (model, variables) in enumerate(models_variables.items()):
-            options = dict(model[1])
-            if options["surface form"] is False:
+            if model == "direct form":
                 base_model_results = models_variables[model]
                 continue
-            if k == 0:
-                label = model[0]
-            else:
-                label = None
-
             error = np.abs(
-                variables["Terminal voltage [V]"](t_eval)
-                - base_model_results["Terminal voltage [V]"](t_eval)
+                variables["Battery voltage [V]"](t_eval)
+                - base_model_results["Battery voltage [V]"](t_eval)
             )
-            ax.loglog(variables["Time [s]"](t_eval), error, linestyles[j], label=label)
-        # plt.legend(loc="upper right")
+            ax.loglog(variables["Time [h]"](t_eval), error, linestyles[j], label=model)
+        ax.legend(loc="best")
+    fig.tight_layout()
     file_name = "capacitance_errors_voltages.eps".format(Crate)
-    plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000)
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000)
 
 
-if __name__ == "__main__":
-    parser = argparse.ArgumentParser()
-    parser.add_argument("--compute", action="store_true", help="(Re)-compute results.")
-    args = parser.parse_args()
-    t_eval = np.concatenate([np.logspace(-6, -3, 50), np.linspace(0.001, 1, 100)])
-    if args.compute:
-        pybamm.set_logging_level("INFO")
+def discharge_states(compute):
+    savefile = "effect_of_capacitance_data.pickle"
+    if compute:
         models = [
-            pybamm.lead_acid.NewmanTiedemann(),
-            pybamm.lead_acid.NewmanTiedemann({"surface form": "differential"}),
-            pybamm.lead_acid.NewmanTiedemann({"surface form": "algebraic"}),
+            pybamm.lead_acid.NewmanTiedemann(name="direct form"),
+            pybamm.lead_acid.NewmanTiedemann(
+                {"surface form": "differential"},
+                name="capacitance form\n(differential)",
+            ),
+            pybamm.lead_acid.NewmanTiedemann(
+                {"surface form": "algebraic"}, name="capacitance form\n(algebraic)"
+            ),
         ]
-        Crates = [0.1, 0.5, 1, 2]
+        Crates = [0.1, 1, 5]
+        t_eval = np.concatenate(
+            [np.logspace(-6, -3, 50), np.linspace(0.001, 1, 100)[1:]]
+        )
         all_variables, t_eval = model_comparison(models, Crates, t_eval)
-        with open("capacitance_data.pickle", "wb") as f:
+        with open(savefile, "wb") as f:
             data = (all_variables, t_eval)
             pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)
+    else:
+        try:
+            with open(savefile, "rb") as f:
+                (all_variables, t_eval) = pickle.load(f)
+        except FileNotFoundError:
+            raise FileNotFoundError(
+                "Run script with '--compute' first to generate results"
+            )
+    plot_voltages(all_variables, t_eval)
+    plot_errors(all_variables, t_eval, [5])
+
+
+def plot_times(models_times_and_voltages):
+    shared_plotting.plot_times(
+        models_times_and_voltages, Crate=1, linestyles=["k-", "b-.", "r--"]
+    )
+    file_name = "capacitance_solver_times.eps"
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight")
+
+
+def discharge_times_and_errors(compute):
+    savefile = "capacitance_times_and_errors.pickle"
+    if compute:
+        try:
+            with open(savefile, "rb") as f:
+                models_times_and_voltages = pickle.load(f)
+        except FileNotFoundError:
+            models_times_and_voltages = pybamm.get_infinite_nested_dict()
+        models = [
+            pybamm.lead_acid.NewmanTiedemann(name="direct form"),
+            pybamm.lead_acid.NewmanTiedemann(
+                {"surface form": "differential"},
+                name="capacitance form\n(differential)",
+            ),
+            pybamm.lead_acid.NewmanTiedemann(
+                {"surface form": "algebraic"}, name="capacitance form\n(algebraic)"
+            ),
+        ]
+        Crates = [1]
+        all_npts = np.linspace(10, 100, 2)
+        t_eval = np.linspace(0, 1, 100)
+        new_models_times_and_voltages = convergence_study(
+            models, Crates, all_npts, t_eval
+        )
+        models_times_and_voltages.update(new_models_times_and_voltages)
+        with open(savefile, "wb") as f:
+            pickle.dump(models_times_and_voltages, f, pickle.HIGHEST_PROTOCOL)
+    else:
+        try:
+            with open(savefile, "rb") as f:
+                models_times_and_voltages = pickle.load(f)
+        except FileNotFoundError:
+            raise FileNotFoundError(
+                "Run script with '--compute' first to generate results"
+            )
+    # plot_errors(models_times_and_voltages)
+    plot_times(models_times_and_voltages)
 
-    with open("capacitance_data.pickle", "rb") as f:
-        (all_variables, t_eval) = pickle.load(f)
-    plot_voltages(all_variables, t_eval, [1])
-    plot_errors(all_variables, t_eval, [0.1, 0.5, 1, 2])
+
+if __name__ == "__main__":
+    pybamm.set_logging_level("INFO")
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--compute", action="store_true", help="(Re)-compute results.")
+    args = parser.parse_args()
+    # discharge_states(args.compute)
+    discharge_times_and_errors(args.compute)
+    plt.show()
diff --git a/results/2019_08_sulzer_thesis/effect_of_side_reactions.py b/results/2019_08_sulzer_thesis/effect_of_side_reactions.py
new file mode 100644
index 000000000..7b23539e7
--- /dev/null
+++ b/results/2019_08_sulzer_thesis/effect_of_side_reactions.py
@@ -0,0 +1,136 @@
+#
+# Simulations: effect of side reactions for charge of a lead-acid battery
+#
+import argparse
+import matplotlib.pyplot as plt
+import numpy as np
+import pickle
+import pybamm
+import shared_plotting
+from shared_solutions import model_comparison
+
+try:
+    from config import OUTPUT_DIR
+except ImportError:
+    OUTPUT_DIR = None
+
+
+def plot_voltages(all_variables, t_eval):
+    shared_plotting.plot_voltages(all_variables, t_eval)
+    file_name = "side_reactions_voltage_comparison.eps"
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight")
+
+
+def plot_interfacial_currents(models_variables, t_eval):
+    models = list(models_variables.keys())
+    file_name = "side_reactions_interfacial_current_density_comparison.eps"
+    fig, ax = plt.subplots(1, 1)
+    output_vars = [
+        "Average positive electrode interfacial current density",
+        "Average positive electrode oxygen interfacial current density",
+        "Average negative electrode oxygen interfacial current density",
+        "Average negative electrode interfacial current density",
+    ]
+    labels = [
+        "Positive electrode (main)",
+        "Positive electrode (oxygen)",
+        "Negative electrode (oxygen)",
+        "Negative electrode (main)",
+    ]
+    t_max = max(np.nanmax(var["Time [h]"](t_eval)) for var in models_variables.values())
+    ax.set_xlim([0, t_max])
+    ax.set_xlabel("Time [h]")
+    ax.set_ylabel("Interfacial current densities")
+    linestyles = ["--", ":", "-.", "-"]
+    colors = ["k", "r"]
+    plots = {}
+    for j, (model, variables) in enumerate(models_variables.items()):
+        if model in models:
+            for k, var in enumerate(output_vars):
+                plots[(model, k)], = ax.plot(
+                    variables["Time [h]"](t_eval),
+                    variables[var](t_eval),
+                    linestyle=linestyles[k],
+                    color=colors[j],
+                )
+    leg1 = ax.legend(
+        [plots[(model, len(linestyles) - 1)] for model in models],
+        models,
+        loc="center left",
+        bbox_to_anchor=(1, 0.25),
+    )
+    ax.legend(labels, loc="center left", bbox_to_anchor=(1, 0.75))
+    ax.add_artist(leg1)
+    fig.tight_layout()
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000)
+
+
+def plot_variables(all_variables, t_eval):
+    # Set up
+    Crates = [-0.1, -1, -5]
+    times = np.linspace(0, 2, 4)
+    var_file_names = {
+        "Electrolyte concentration [Molar]"
+        + "": "charge_electrolyte_concentration_comparison.eps"
+    }
+    limits_exceptions = {"Electrolyte concentration [Molar]": {"min": 0}}
+    all_variables = {k: v for k, v in all_variables.items() if k in Crates}
+    for var, file_name in var_file_names.items():
+        if var in limits_exceptions:
+            exceptions = limits_exceptions[var]
+        else:
+            exceptions = {}
+        shared_plotting.plot_variable(all_variables, times, var, exceptions)
+        if OUTPUT_DIR is not None:
+            plt.savefig(
+                OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight"
+            )
+
+
+def charge_states(compute):
+    savefile = "effect_of_side_reactions_data.pickle"
+    if compute:
+        models = [
+            pybamm.lead_acid.NewmanTiedemann(
+                {"surface form": "algebraic"}, name="Without oxygen"
+            ),
+            pybamm.lead_acid.NewmanTiedemann(
+                {"surface form": "algebraic", "side reactions": ["oxygen"]},
+                name="With oxygen",
+            ),
+        ]
+        Crates = [-0.1, -1, -5]
+        t_eval = np.linspace(0, 2, 100)
+        extra_parameter_values = {
+            "Positive electrode"
+            + "reference exchange-current density (oxygen) [A.m-2]": 1e-24,
+            "Initial State of Charge": 0.5,
+        }
+        all_variables, t_eval = model_comparison(
+            models, Crates, t_eval, extra_parameter_values=extra_parameter_values
+        )
+        with open(savefile, "wb") as f:
+            data = (all_variables, t_eval)
+            pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)
+    else:
+        try:
+            with open(savefile, "rb") as f:
+                (all_variables, t_eval) = pickle.load(f)
+        except FileNotFoundError:
+            raise FileNotFoundError(
+                "Run script with '--compute' first to generate results"
+            )
+    plot_voltages(all_variables, t_eval)
+    plot_interfacial_currents(all_variables[-1], t_eval)
+    plot_variables(all_variables, t_eval)
+
+
+if __name__ == "__main__":
+    pybamm.set_logging_level("DEBUG")
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--compute", action="store_true", help="(Re)-compute results.")
+    args = parser.parse_args()
+    charge_states(args.compute)
+    plt.show()
diff --git a/results/2019_08_sulzer_thesis/lead_acid_charge.py b/results/2019_08_sulzer_thesis/lead_acid_charge.py
new file mode 100644
index 000000000..b0d2f82e7
--- /dev/null
+++ b/results/2019_08_sulzer_thesis/lead_acid_charge.py
@@ -0,0 +1,146 @@
+#
+# Simulations: charge of a lead-acid battery
+#
+import argparse
+import matplotlib.pyplot as plt
+import numpy as np
+import pickle
+import pybamm
+import shared_plotting
+from shared_solutions import model_comparison
+
+try:
+    from config import OUTPUT_DIR
+except ImportError:
+    OUTPUT_DIR = None
+
+
+def plot_voltages(all_variables, t_eval):
+    shared_plotting.plot_voltages(all_variables, t_eval)
+    file_name = "charge_voltage_comparison.eps"
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight")
+
+
+def plot_interfacial_currents(models_variables, t_eval):
+    models = ["Full", "LOQS"]
+    file_name = "charge_interfacial_current_density_comparison.eps"
+    fig, ax = plt.subplots(1, 1)
+    output_vars = [
+        "Average positive electrode interfacial current density",
+        "Average positive electrode oxygen interfacial current density",
+        "Average negative electrode oxygen interfacial current density",
+        "Average negative electrode interfacial current density",
+    ]
+    labels = [
+        "Positive electrode (main)",
+        "Positive electrode (oxygen)",
+        "Negative electrode (oxygen)",
+        "Negative electrode (main)",
+    ]
+    t_max = max(np.nanmax(var["Time [h]"](t_eval)) for var in models_variables.values())
+    ax.set_xlim([0, t_max])
+    ax.set_xlabel("Time [h]")
+    ax.set_ylabel("Interfacial current densities")
+    linestyles = ["--", ":", "-.", "-"]
+    colors = ["k", "r"]
+    plots = {}
+    for j, (model, variables) in enumerate(models_variables.items()):
+        if model in models:
+            for k, var in enumerate(output_vars):
+                plots[(model, k)], = ax.plot(
+                    variables["Time [h]"](t_eval),
+                    variables[var](t_eval),
+                    linestyle=linestyles[k],
+                    color=colors[j],
+                )
+    leg1 = ax.legend(
+        [plots[("Full", 3)], plots[("LOQS", 3)]],
+        ["Full", "LOQS"],
+        loc="center left",
+        bbox_to_anchor=(1, 0.25),
+    )
+    ax.legend(labels, loc="center left", bbox_to_anchor=(1, 0.75))
+    ax.add_artist(leg1)
+    fig.tight_layout()
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000)
+
+
+def plot_variables(all_variables, t_eval):
+    # Set up
+    Crates = [-0.1, -2, -5]
+    times = np.linspace(0, 2, 4)
+    var_file_names = {
+        "Electrolyte concentration [Molar]"
+        + "": "charge_electrolyte_concentration_comparison.eps",
+        "Oxygen concentration [Molar]": "charge_oxygen_concentration_comparison.eps",
+    }
+    limits_exceptions = {"Electrolyte concentration [Molar]": {"min": 0}}
+    all_variables = {k: v for k, v in all_variables.items() if k in Crates}
+    for var, file_name in var_file_names.items():
+        if var in limits_exceptions:
+            exceptions = limits_exceptions[var]
+        else:
+            exceptions = {}
+        shared_plotting.plot_variable(all_variables, times, var, exceptions)
+        if OUTPUT_DIR is not None:
+            plt.savefig(
+                OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight"
+            )
+
+
+def plot_voltage_breakdown(all_variables, t_eval):
+    Crates = [-0.1, -2, -5]
+    model = "LOQS"
+    shared_plotting.plot_voltage_breakdown(all_variables, t_eval, model, Crates)
+    file_name = "charge_voltage_breakdown.eps"
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight")
+
+
+def charge_states(compute):
+    if compute:
+        models = [
+            pybamm.lead_acid.NewmanTiedemann(
+                {"side reactions": ["oxygen"]}, name="Full"
+            ),
+            pybamm.lead_acid.LOQS(
+                {"surface form": "algebraic", "side reactions": ["oxygen"]}, name="LOQS"
+            ),
+        ]
+        Crates = [-0.1, -0.2, -0.5, -1, -2, -5]
+        t_eval = np.linspace(0, 2, 100)
+        extra_parameter_values = {
+            "Positive electrode"
+            + "reference exchange-current density (oxygen) [A.m-2]": 1e-24,
+            "Initial State of Charge": 0.5,
+        }
+        all_variables, t_eval = model_comparison(
+            models, Crates, t_eval, extra_parameter_values=extra_parameter_values
+        )
+        with open("charge_asymptotics_data.pickle", "wb") as f:
+            data = (all_variables, t_eval)
+            pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)
+    else:
+        try:
+            with open("charge_asymptotics_data.pickle", "rb") as f:
+                (all_variables, t_eval) = pickle.load(f)
+        except FileNotFoundError:
+            raise FileNotFoundError(
+                "Run script with '--compute' first to generate results"
+            )
+    # plot_voltages(all_variables, t_eval)
+    plot_interfacial_currents(all_variables[-1], t_eval)
+    # plot_variables(all_variables, t_eval)
+    # plot_voltage_breakdown(all_variables, t_eval)
+
+
+if __name__ == "__main__":
+    pybamm.set_logging_level("INFO")
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--compute", action="store_true", help="(Re)-compute results.")
+    args = parser.parse_args()
+    charge_states(args.compute)
+    # charge_times_and_errors(args.compute)
+    plt.show()
diff --git a/results/2019_08_sulzer_thesis/lead_acid_discharge.py b/results/2019_08_sulzer_thesis/lead_acid_discharge.py
index 386c2f139..23f931c85 100644
--- a/results/2019_08_sulzer_thesis/lead_acid_discharge.py
+++ b/results/2019_08_sulzer_thesis/lead_acid_discharge.py
@@ -6,68 +6,164 @@ import matplotlib.pyplot as plt
 import numpy as np
 import pickle
 import pybamm
-from config import OUTPUT_DIR
-from shared import model_comparison
+import shared_plotting
+from collections import defaultdict
+from shared_solutions import model_comparison, convergence_study
+
+try:
+    from config import OUTPUT_DIR
+except ImportError:
+    OUTPUT_DIR = None
 
 
 def plot_voltages(all_variables, t_eval):
-    # Plot
-    plt.subplots(figsize=(6, 4))
-    n = int(len(all_variables) // np.sqrt(len(all_variables)))
-    m = np.ceil(len(all_variables) / n)
-    for k, Crate in enumerate(all_variables.keys()):
-        models_variables = all_variables[Crate]
-        t_max = max(
-            np.nanmax(var["Time [h]"](t_eval)) for var in models_variables.values()
-        )
-        ax = plt.subplot(n, m, k + 1)
-        plt.axis([0, t_max, 10.5, 13])
-        plt.xlabel("Time [h]")
-        plt.title("\\textbf{{({})}} {}C".format(chr(97 + k), Crate))
-
-        # Hide the right and top spines
-        ax.spines["right"].set_visible(False)
-        ax.spines["top"].set_visible(False)
-
-        # Only show ticks on the left and bottom spines
-        ax.yaxis.set_ticks_position("left")
-        ax.xaxis.set_ticks_position("bottom")
-        for model, variables in models_variables.items():
-            if k == 0:
-                label = model[0]
-            else:
-                label = None
-            if k % m == 0:
-                plt.ylabel("Voltage [V]")
-
-            plt.plot(
-                variables["Time [h]"](t_eval),
-                variables["Terminal voltage [V]"](t_eval) * 6,
-                label=label,
+    shared_plotting.plot_voltages(all_variables, t_eval)
+    file_name = "discharge_voltage_comparison.eps"
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight")
+
+
+def plot_variables(all_variables, t_eval):
+    # Set up
+    Crates = [0.1, 2, 5]
+    times = np.linspace(0, 0.5, 4)
+    var_file_names = {
+        "Electrolyte concentration [Molar]"
+        + "": "discharge_electrolyte_concentration_comparison.eps",
+        "Electrolyte potential [V]": "discharge_electrolyte_potential_comparison.eps",
+        "Interfacial current density"
+        + "": "discharge_interfacial_current_density_comparison.eps",
+    }
+    limits_exceptions = {"Electrolyte concentration [Molar]": {"min": 0}}
+    all_variables = {k: v for k, v in all_variables.items() if k in Crates}
+    for var, file_name in var_file_names.items():
+        if var in limits_exceptions:
+            exceptions = limits_exceptions[var]
+        else:
+            exceptions = {}
+        shared_plotting.plot_variable(all_variables, times, var, exceptions)
+        if OUTPUT_DIR is not None:
+            plt.savefig(
+                OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight"
             )
-        # plt.legend(loc="upper right")
-    file_name = "discharge_voltage_comparison.eps".format(Crate)
-    plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000)
 
 
-if __name__ == "__main__":
-    parser = argparse.ArgumentParser()
-    parser.add_argument("--compute", action="store_true", help="(Re)-compute results.")
-    args = parser.parse_args()
-    t_eval = np.linspace(0, 1, 100)
-    if args.compute:
-        pybamm.set_logging_level("INFO")
+def plot_voltage_breakdown(all_variables, t_eval):
+    Crates = [0.1, 2, 5]
+    model = "Composite"
+    shared_plotting.plot_voltage_breakdown(all_variables, t_eval, model, Crates)
+    file_name = "discharge_voltage_breakdown.eps"
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight")
+
+
+def discharge_states(compute):
+    savefile = "discharge_asymptotics_data.pickle"
+    if compute:
         models = [
-            pybamm.lead_acid.LOQS(),
-            pybamm.lead_acid.Composite(),
-            pybamm.lead_acid.NewmanTiedemann(),
+            pybamm.lead_acid.NewmanTiedemann(name="Full"),
+            pybamm.lead_acid.LOQS(name="LOQS"),
+            pybamm.lead_acid.FOQS(name="FOQS"),
+            pybamm.lead_acid.Composite(name="Composite"),
         ]
         Crates = [0.1, 0.2, 0.5, 1, 2, 5]
-        all_variables, t_eval = model_comparison(models, Crates, t_eval)
-        with open("discharge_asymptotics_data.pickle", "wb") as f:
+        t_eval = np.linspace(0, 1, 100)
+        extra_parameter_values = {"Bruggeman coefficient": 0.001}
+        all_variables, t_eval = model_comparison(
+            models, Crates, t_eval, extra_parameter_values=extra_parameter_values
+        )
+        with open(savefile, "wb") as f:
             data = (all_variables, t_eval)
             pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)
+    else:
+        try:
+            with open(savefile, "rb") as f:
+                (all_variables, t_eval) = pickle.load(f)
+        except FileNotFoundError:
+            raise FileNotFoundError(
+                "Run script with '--compute' first to generate results"
+            )
+    # plot_voltages(all_variables, t_eval)
+    # plot_variables(all_variables, t_eval)
+    plot_voltage_breakdown(all_variables, t_eval)
+
+
+def plot_errors(models_times_and_voltages):
+    npts = 20
+    linestyles = ["k-", "g--", "r:", "b-."]
+    Crates = defaultdict(list)
+    voltage_errors = defaultdict(list)
+    fig, ax = plt.subplots(1, 1)
+    for i, (model, times_and_voltages) in enumerate(models_times_and_voltages.items()):
+        if model != "Full":
+            for Crate, variables in times_and_voltages[npts].items():
+                Crates[model].append(Crate)
+                full_voltage = models_times_and_voltages["Full"][npts][Crate][
+                    "Battery voltage [V]"
+                ]
+                reduced_voltage = variables["Battery voltage [V]"]
+                voltage_errors[model].append(pybamm.rmse(full_voltage, reduced_voltage))
+            ax.semilogx(
+                Crates[model], voltage_errors[model], linestyles[i], label=model
+            )
+    ax.set_xlabel("C-rate")
+    ax.set_ylabel("RMSE [V]")
+    ax.legend(loc="best")
+    fig.tight_layout()
+    file_name = "discharge_asymptotics_rmse.eps"
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight")
+
+
+def plot_times(models_times_and_voltages):
+    shared_plotting.plot_times(models_times_and_voltages, Crate=1)
+    file_name = "discharge_asymptotics_solver_times.eps"
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight")
+
+
+def discharge_times_and_errors(compute):
+    savefile = "discharge_asymptotics_times_and_errors.pickle"
+    if compute:
+        try:
+            with open(savefile, "rb") as f:
+                models_times_and_voltages = pickle.load(f)
+        except FileNotFoundError:
+            models_times_and_voltages = pybamm.get_infinite_nested_dict()
+        models = [
+            pybamm.lead_acid.NewmanTiedemann(
+                {"surface form": "algebraic"}, name="Full"
+            ),
+            pybamm.lead_acid.LOQS(name="LOQS"),
+            # pybamm.lead_acid.FOQS(name="FOQS"),
+            # pybamm.lead_acid.Composite(name="Composite"),
+        ]
+        Crates = np.linspace(0.01, 5, 2)
+        all_npts = [20]
+        t_eval = np.linspace(0, 1, 100)
+        new_models_times_and_voltages = convergence_study(
+            models, Crates, all_npts, t_eval
+        )
+        models_times_and_voltages.update(new_models_times_and_voltages)
+        with open(savefile, "wb") as f:
+            pickle.dump(models_times_and_voltages, f, pickle.HIGHEST_PROTOCOL)
+    else:
+        try:
+            with open(savefile, "rb") as f:
+                models_times_and_voltages = pickle.load(f)
+        except FileNotFoundError:
+            raise FileNotFoundError(
+                "Run script with '--compute' first to generate results"
+            )
+    plot_errors(models_times_and_voltages)
+    # plot_times(models_times_and_voltages)
 
-    with open("discharge_asymptotics_data.pickle", "rb") as f:
-        (all_variables, t_eval) = pickle.load(f)
-    plot_voltages(all_variables, t_eval)
+
+if __name__ == "__main__":
+    pybamm.set_logging_level("INFO")
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--compute", action="store_true", help="(Re)-compute results.")
+    args = parser.parse_args()
+    discharge_states(args.compute)
+    # discharge_times_and_errors(args.compute)
+    plt.show()
diff --git a/results/2019_08_sulzer_thesis/self_discharge.py b/results/2019_08_sulzer_thesis/self_discharge.py
new file mode 100644
index 000000000..ca9182586
--- /dev/null
+++ b/results/2019_08_sulzer_thesis/self_discharge.py
@@ -0,0 +1,61 @@
+#
+# Simulations: self-discharge
+#
+import argparse
+import matplotlib.pyplot as plt
+import numpy as np
+import pickle
+import pybamm
+import shared_plotting
+from config import OUTPUT_DIR
+from shared_solutions import model_comparison
+
+
+def plot_voltages(all_variables, t_eval):
+    shared_plotting.plot_voltages(all_variables, t_eval)
+    file_name = "sefl_discharge_voltage_comparison.eps"
+    if OUTPUT_DIR is not None:
+        plt.savefig(OUTPUT_DIR + file_name, format="eps", dpi=1000, bbox_inches="tight")
+
+
+def self_discharge_states(compute):
+    save_file = "self_discharge_data.pickle"
+    if compute:
+        models = [
+            pybamm.lead_acid.NewmanTiedemann(name="Full, without oxygen"),
+            pybamm.lead_acid.NewmanTiedemann(
+                {"side reactions": ["oxygen"]}, name="Full, with oxygen"
+            ),
+            pybamm.lead_acid.LOQS(
+                {"surface form": "algebraic", "side reactions": ["oxygen"]},
+                name="LOQS, with oxygen",
+            ),
+        ]
+        extra_parameter_values = {
+            "Current function": pybamm.GetConstantCurrent(current=0)
+        }
+        t_eval = np.linspace(0, 1000, 100)
+        all_variables, t_eval = model_comparison(
+            models, [1], t_eval, extra_parameter_values=extra_parameter_values
+        )
+        with open(save_file, "wb") as f:
+            data = (all_variables, t_eval)
+            pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)
+    else:
+        try:
+            with open(save_file, "rb") as f:
+                (all_variables, t_eval) = pickle.load(f)
+        except FileNotFoundError:
+            raise FileNotFoundError(
+                "Run script with '--compute' first to generate results"
+            )
+    plot_voltages(all_variables, t_eval)
+
+
+if __name__ == "__main__":
+    pybamm.set_logging_level("INFO")
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--compute", action="store_true", help="(Re)-compute results.")
+    args = parser.parse_args()
+    self_discharge_states(args.compute)
+    plt.show()
diff --git a/results/2019_08_sulzer_thesis/shared.py b/results/2019_08_sulzer_thesis/shared.py
deleted file mode 100644
index 470a28212..000000000
--- a/results/2019_08_sulzer_thesis/shared.py
+++ /dev/null
@@ -1,55 +0,0 @@
-#
-# Simulations: discharge of a lead-acid battery
-#
-import pybamm
-
-
-def options_to_tuple(options):
-    bc_options = tuple(options["bc_options"].items())
-    other_options = tuple(
-        {k: v for k, v in options.items() if k != "bc_options"}.items()
-    )
-    return (*bc_options, *other_options)
-
-
-def model_comparison(models, Crates, t_eval):
-    # load parameter values and geometry
-    geometry = models[0].default_geometry
-    param = models[0].default_parameter_values
-
-    # Process parameters (same parameters for all models)
-    for model in models:
-        param.process_model(model)
-    param.process_geometry(geometry)
-
-    # set mesh
-    var = pybamm.standard_spatial_vars
-    var_pts = {var.x_n: 20, var.x_s: 20, var.x_p: 20}
-    mesh = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts)
-
-    # discretise models
-    discs = {}
-    for model in models:
-        disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-        disc.process_model(model)
-        # Store discretisation
-        discs[model] = disc
-
-    # solve model for range of Crates
-    all_variables = {}
-    for Crate in Crates:
-        all_variables[Crate] = {}
-        current = Crate * 17
-        pybamm.logger.info("Setting typical current to {} A".format(current))
-        param.update({"Typical current [A]": current})
-        for model in models:
-            param.update_model(model, discs[model])
-            solver = model.default_solver
-            solution = solver.solve(model, t_eval)
-            all_variables[Crate][
-                (model.name, options_to_tuple(model.options))
-            ] = pybamm.post_process_variables(
-                model.variables, solution.t, solution.y, mesh
-            )
-
-    return all_variables, t_eval
diff --git a/results/2019_08_sulzer_thesis/shared_plotting.py b/results/2019_08_sulzer_thesis/shared_plotting.py
new file mode 100644
index 000000000..a6383bfaf
--- /dev/null
+++ b/results/2019_08_sulzer_thesis/shared_plotting.py
@@ -0,0 +1,229 @@
+#
+# Shared plotting
+#
+import matplotlib.pyplot as plt
+import numpy as np
+import pybamm
+from collections import defaultdict
+
+
+def plot_voltages(all_variables, t_eval, linestyles=None):
+    # Plot
+    linestyles = linestyles or ["k-", "g--", "r:", "b-."]
+    n = int(len(all_variables) // np.sqrt(len(all_variables)))
+    m = int(np.ceil(len(all_variables) / n))
+    fig, axes = plt.subplots(n, m, figsize=(6.4, 4.5))
+    labels = [model for model in [x for x in all_variables.values()][0].keys()]
+    y_min = 0.98 * min(
+        np.nanmin(variables["Battery voltage [V]"](t_eval))
+        for models_variables in all_variables.values()
+        for variables in models_variables.values()
+    )
+    y_max = 1.02 * max(
+        np.nanmax(variables["Battery voltage [V]"](t_eval))
+        for models_variables in all_variables.values()
+        for variables in models_variables.values()
+    )
+    # Strict voltage cut-offs
+    y_min = max(y_min, 10.5)
+    y_max = min(y_max, 14.5)
+    for k, (Crate, models_variables) in enumerate(all_variables.items()):
+        if len(all_variables) == 1:
+            ax = axes
+        else:
+            ax = axes.flat[k]
+        t_max = max(
+            np.nanmax(var["Time [h]"](t_eval)) for var in models_variables.values()
+        )
+        ax.set_xlim([0, t_max])
+        ax.set_ylim([y_min, y_max])
+        ax.set_xlabel("Time [h]")
+        if len(all_variables) > 1:
+            ax.set_title(
+                "\\textbf{{({})}} {}C ($\\mathcal{{C}}_e={}$)".format(
+                    chr(97 + k), abs(Crate), abs(Crate) * 0.6
+                )
+            )
+            # # Hide the right and top spines
+            # ax.spines["right"].set_visible(False)
+            # ax.spines["top"].set_visible(False)
+            #
+            # # Only show ticks on the left and bottom spines
+            # ax.yaxis.set_ticks_position("left")
+            # ax.xaxis.set_ticks_position("bottom")
+        ax.xaxis.set_major_locator(plt.MaxNLocator(3))
+        if k % m == 0:
+            ax.set_ylabel("Voltage [V]")
+        else:
+            ax.set_yticklabels([])
+
+        for j, variables in enumerate(models_variables.values()):
+            ax.plot(
+                variables["Time [h]"](t_eval),
+                variables["Battery voltage [V]"](t_eval),
+                linestyles[j],
+            )
+    if len(all_variables) == 1:
+        leg = ax.legend(labels, loc="best")
+        fig.tight_layout()
+    else:
+        leg = fig.legend(labels, loc="lower center", ncol=len(labels))
+        plt.subplots_adjust(bottom=0.25, right=0.95, hspace=1.1, wspace=0.3)
+    leg.get_frame().set_edgecolor("k")
+    return fig, axes
+
+
+def plot_variable(all_variables, times, variable, limits_exceptions=None):
+    linestyles = ["k-", "g--", "r:", "b-."]
+    n = len(times)
+    m = len(all_variables)
+    Crates = list(all_variables.keys())
+    labels = [model for model in [x for x in all_variables.values()][0].keys()]
+    x = all_variables[Crates[0]][labels[0]]["x"](0, np.linspace(0, 1))[:, 0]
+    x_dim = all_variables[Crates[0]][labels[0]]["x [m]"](0, np.linspace(0, 1))[:, 0]
+
+    fig, axes = plt.subplots(n, m, figsize=(6.4, 6.4))
+
+    # Default limits
+    y_min = pybamm.ax_min(
+        [
+            np.nanmin(variables[variable](times, x))
+            for Crate, models_variables in all_variables.items()
+            for variables in models_variables.values()
+        ]
+    )
+    y_max = pybamm.ax_max(
+        [
+            np.nanmax(variables[variable](times, x))
+            for Crate, models_variables in all_variables.items()
+            for variables in models_variables.values()
+        ]
+    )
+    # Exceptions
+    if "min" in limits_exceptions:
+        y_min = limits_exceptions["min"]
+    if "max" in limits_exceptions:
+        y_max = limits_exceptions["max"]
+
+    # Plot
+    for i, (Crate, models_variables) in enumerate(all_variables.items()):
+        for j, time in enumerate(times):
+            ax = axes[j, i]
+            ax.set_xlim([x_dim[0], x_dim[-1]])
+            ax.set_ylim([y_min, y_max])
+            ax.yaxis.set_major_locator(plt.MaxNLocator(3))
+
+            # Title
+            if j == 0:
+                ax.set_title(
+                    "\\textbf{{({})}} {}C ($\\mathcal{{C}}_e={}$)".format(
+                        chr(97 + i), abs(Crate), abs(Crate) * 0.6
+                    )
+                )
+            # x-axis
+            if j == len(times) - 1:
+                ax.set_xlabel("x [m]")
+            else:
+                ax.set_xticklabels([])
+
+            # y-axis
+            if i == 0:
+                for variables in models_variables.values():
+                    try:
+                        soc = variables["State of Charge"](time)
+                        ax.set_ylabel(
+                            "{}\% SoC".format(int(soc)), rotation=0, labelpad=30
+                        )
+                        ax.yaxis.get_label().set_verticalalignment("center")
+                    except ValueError:
+                        pass
+            else:
+                ax.set_yticklabels([])
+
+            # Plot
+            for j, variables in enumerate(models_variables.values()):
+                ax.plot(x_dim, variables[variable](time, x), linestyles[j])
+    leg = fig.legend(labels, loc="lower center", ncol=len(labels), frameon=True)
+    leg.get_frame().set_edgecolor("k")
+    plt.subplots_adjust(
+        bottom=0.15, top=0.95, left=0.18, right=0.97, hspace=0.08, wspace=0.05
+    )
+
+
+def plot_voltage_breakdown(all_variables, t_eval, model, Crates):
+    n = int(len(Crates) // np.sqrt(len(Crates)))
+    m = int(np.ceil(len(Crates) / n))
+    fig, axes = plt.subplots(n, m, figsize=(6.4, 2.3))
+    labels = ["V", "$V_U$", "$V_k$", "$V_c$", "$V_o$"]
+    overpotentials = [
+        "Average battery reaction overpotential [V]",
+        "Average battery concentration overpotential [V]",
+        "Average battery electrolyte ohmic losses [V]",
+    ]
+    y_min = 0.95 * min(
+        np.nanmin(models_variables[model]["Battery voltage [V]"](t_eval))
+        for models_variables in all_variables.values()
+    )
+    y_max = 1.05 * max(
+        np.nanmax(models_variables[model]["Battery voltage [V]"](t_eval))
+        for models_variables in all_variables.values()
+    )
+    for k, Crate in enumerate(Crates):
+        variables = all_variables[Crate][model]
+        ax = axes.flat[k]
+
+        # Set up
+        t_max = np.nanmax(variables["Time [h]"](t_eval))
+        ax.set_xlim([0, t_max])
+        ax.set_ylim([y_min, y_max])
+        ax.set_xlabel("Time [h]")
+        ax.set_title(
+            "\\textbf{{({})}} {}C ($\\mathcal{{C}}_e={}$)".format(
+                chr(97 + k), abs(Crate), abs(Crate) * 0.6
+            )
+        )
+        ax.xaxis.set_major_locator(plt.MaxNLocator(3))
+        if k % m == 0:
+            ax.set_ylabel("Voltage [V]")
+
+        # Plot
+        # Initialise
+        # for lead-acid we multiply everything by 6 to
+        time = variables["Time [h]"](t_eval)
+        initial_ocv = variables["Average battery open circuit voltage [V]"](0)
+        ocv = variables["Average battery open circuit voltage [V]"](t_eval)
+        ax.fill_between(time, ocv, initial_ocv)
+        top = ocv
+        # Plot
+        for overpotential in overpotentials:
+            bottom = top + variables[overpotential](t_eval)
+            ax.fill_between(time, bottom, top)
+            top = bottom
+        ax.plot(time, variables["Battery voltage [V]"](t_eval), "k--")
+    leg = axes.flat[-1].legend(
+        labels, bbox_to_anchor=(1.05, 0.5), loc="center left", frameon=True
+    )
+    leg.get_frame().set_edgecolor("k")
+    fig.tight_layout()
+
+
+def plot_times(models_times_and_voltages, Crate=1, linestyles=None):
+    linestyles = linestyles or ["k-", "g--", "r:", "b-."]
+    all_npts = defaultdict(list)
+    solver_times = defaultdict(list)
+    fig, ax = plt.subplots(1, 1)
+    for i, (model, times_and_voltages) in enumerate(models_times_and_voltages.items()):
+        for npts in times_and_voltages.keys():
+            try:
+                solver_time = times_and_voltages[npts][Crate][
+                    "solution object"
+                ].solve_time
+            except KeyError:
+                continue
+            all_npts[model].append(npts * 3)
+            solver_times[model].append(solver_time)
+        ax.loglog(all_npts[model], solver_times[model], linestyles[i], label=model)
+    ax.set_xlabel("Number of grid points")
+    ax.set_ylabel("Solver time [s]")
+    ax.legend(loc="best")
+    fig.tight_layout()
diff --git a/results/2019_08_sulzer_thesis/shared_solutions.py b/results/2019_08_sulzer_thesis/shared_solutions.py
new file mode 100644
index 000000000..86fe5911b
--- /dev/null
+++ b/results/2019_08_sulzer_thesis/shared_solutions.py
@@ -0,0 +1,142 @@
+#
+# Simulations
+#
+import pybamm
+
+
+def model_comparison(models, Crates, t_eval, extra_parameter_values=None):
+    " Solve models at a range of Crates "
+    # load parameter values and geometry
+    geometry = models[0].default_geometry
+    extra_parameter_values = extra_parameter_values or {}
+    param = models[0].default_parameter_values
+    param.update(extra_parameter_values)
+
+    # Process parameters (same parameters for all models)
+    for model in models:
+        param.process_model(model)
+    param.process_geometry(geometry)
+
+    # set mesh
+    var = pybamm.standard_spatial_vars
+    var_pts = {var.x_n: 20, var.x_s: 20, var.x_p: 20}
+    mesh = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts)
+
+    # discretise models
+    discs = {}
+    for model in models:
+        disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+        disc.process_model(model)
+        # Store discretisation
+        discs[model] = disc
+
+    # solve model for range of Crates
+    all_variables = {}
+    for Crate in Crates:
+        all_variables[Crate] = {}
+        current = Crate * 17
+        pybamm.logger.info("Setting typical current to {} A".format(current))
+        param.update({"Typical current [A]": current})
+        for model in models:
+            param.update_model(model, discs[model])
+            solution = model.default_solver.solve(model, t_eval)
+            variables = pybamm.post_process_variables(
+                model.variables, solution.t, solution.y, mesh
+            )
+            variables["solution"] = solution
+            all_variables[Crate][model.name] = variables
+
+    return all_variables, t_eval
+
+
+def convergence_study(models, Crates, all_npts, t_eval, extra_parameter_values=None):
+    " Solve models at a range of number of grid points "
+    # load parameter values and geometry
+    geometry = models[0].default_geometry
+    param = models[0].default_parameter_values
+    # Update parameters
+    extra_parameter_values = extra_parameter_values or {}
+    param.update(extra_parameter_values)
+
+    # Process parameters (same parameters for all models)
+    for model in models:
+        param.process_model(model)
+    param.process_geometry(geometry)
+
+    # set mesh
+    var = pybamm.standard_spatial_vars
+
+    # solve model for range of Crates and npts
+    models_times_and_voltages = {model.name: {} for model in models}
+    for npts in all_npts:
+        pybamm.logger.info("Setting number of grid points to {}".format(npts))
+        var_pts = {var.x_n: npts, var.x_s: npts, var.x_p: npts}
+        mesh = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts)
+
+        # discretise models, store discretised model and discretisation
+        models_disc = {}
+        discs = {}
+        for model in models:
+            disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+            models_times_and_voltages[model.name][npts] = {}
+            models_disc[model.name] = disc.process_model(model, inplace=False)
+            discs[model.name] = disc
+
+        # Solve for a range of C-rates
+        for Crate in Crates:
+            current = Crate * 17
+            pybamm.logger.info("Setting typical current to {} A".format(current))
+            param.update({"Typical current [A]": current})
+            for model in models:
+                model_disc = models_disc[model.name]
+                disc = discs[model.name]
+                param.update_model(model_disc, disc)
+                try:
+                    solution = model.default_solver.solve(model_disc, t_eval)
+                except pybamm.SolverError:
+                    continue
+                voltage = pybamm.ProcessedVariable(
+                    model_disc.variables["Battery voltage [V]"], solution.t, solution.y
+                )(t_eval)
+                variables = {
+                    "Battery voltage [V]": voltage,
+                    "solution object": solution,
+                }
+                models_times_and_voltages[model.name][npts][Crate] = variables
+
+    return models_times_and_voltages
+
+
+def simulation(models, t_eval, extra_parameter_values=None, disc_only=False):
+
+    # create geometry
+    geometry = models[-1].default_geometry
+
+    # load parameter values and process models and geometry
+    param = models[0].default_parameter_values
+    extra_parameter_values = extra_parameter_values or {}
+    param.update(extra_parameter_values)
+    for model in models:
+        param.process_model(model)
+    param.process_geometry(geometry)
+
+    # set mesh
+    var = pybamm.standard_spatial_vars
+    var_pts = {var.x_n: 25, var.x_s: 41, var.x_p: 34}
+    mesh = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts)
+
+    # discretise models
+    for model in models:
+        disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+        disc.process_model(model)
+
+    if disc_only:
+        return model, mesh
+
+    # solve model
+    solutions = [None] * len(models)
+    for i, model in enumerate(models):
+        solution = model.default_solver.solve(model, t_eval)
+        solutions[i] = solution
+
+    return models, mesh, solutions


Test patch:
diff --git a/tests/unit/test_util.py b/tests/unit/test_util.py
index 0b4adecbd..465bfcfd2 100644
--- a/tests/unit/test_util.py
+++ b/tests/unit/test_util.py
@@ -1,9 +1,7 @@
 #
-# Tests the Timer class.
-#
-# The code in this file is adapted from Pints
-# (see https://github.com/pints-team/pints)
+# Tests the utility functions.
 #
+import numpy as np
 import os
 import pybamm
 import unittest
@@ -61,6 +59,24 @@ class TestUtil(unittest.TestCase):
         func = pybamm.load_function("process_symbol_test_function.py")
         self.assertEqual(func(3), 369)
 
+    def test_rmse(self):
+        self.assertEqual(pybamm.rmse(np.ones(5), np.zeros(5)), 1)
+        self.assertEqual(pybamm.rmse(2 * np.ones(5), np.zeros(5)), 2)
+        self.assertEqual(pybamm.rmse(2 * np.ones(5), np.ones(5)), 1)
+
+        x = np.array([1, 2, 3, 4, 5])
+        self.assertEqual(pybamm.rmse(x, x), 0)
+
+        with self.assertRaisesRegex(ValueError, "same length"):
+            pybamm.rmse(np.ones(5), np.zeros(3))
+
+    def test_infinite_nested_dict(self):
+        d = pybamm.get_infinite_nested_dict()
+        d[1][2][3] = "x"
+        self.assertEqual(d[1][2][3], "x")
+        d[4][5] = "y"
+        self.assertEqual(d[4][5], "y")
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")