diff --git a/CHANGELOG.md b/CHANGELOG.md
index baaa8cfd13..f9c492050e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added parameter set for an A123 LFP cell ([#1209](https://github.com/pybamm-team/PyBaMM/pull/1209))
 -   Added variables related to equivalent circuit models ([#1204](https://github.com/pybamm-team/PyBaMM/pull/1204))
 -   Added an example script to check conservation of lithium ([#1186](https://github.com/pybamm-team/PyBaMM/pull/1186))
 -   Added `erf` and `erfc` functions ([#1184](https://github.com/pybamm-team/PyBaMM/pull/1184))
diff --git a/examples/scripts/compare_spectral_volume.py b/examples/scripts/compare_spectral_volume.py
index 5e5e984903..537289938a 100644
--- a/examples/scripts/compare_spectral_volume.py
+++ b/examples/scripts/compare_spectral_volume.py
@@ -8,8 +8,10 @@
 
 # load model
 # don't use new_copy
-models = [pybamm.lithium_ion.DFN(name="Finite Volume"),
-          pybamm.lithium_ion.DFN(name="Spectral Volume")]
+models = [
+    pybamm.lithium_ion.DFN(name="Finite Volume"),
+    pybamm.lithium_ion.DFN(name="Spectral Volume"),
+]
 
 # create geometry
 geometries = [m.default_geometry for m in models]
@@ -24,32 +26,31 @@
 var = pybamm.standard_spatial_vars
 var_pts = {var.x_n: 1, var.x_s: 1, var.x_p: 1, var.r_n: 1, var.r_p: 1}
 # the Finite Volume method also works on spectral meshes
-meshes = [pybamm.Mesh(
-    geometry,
-    {
-        "negative particle": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
-        ),
-        "positive particle": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
-        ),
-        "negative electrode": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
-        ),
-        "separator": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
-        ),
-        "positive electrode": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
-        ),
-        "current collector": pybamm.SubMesh0D,
-    },
-    var_pts) for geometry in geometries]
+meshes = [
+    pybamm.Mesh(
+        geometry,
+        {
+            "negative particle": pybamm.MeshGenerator(
+                pybamm.SpectralVolume1DSubMesh, {"order": order}
+            ),
+            "positive particle": pybamm.MeshGenerator(
+                pybamm.SpectralVolume1DSubMesh, {"order": order}
+            ),
+            "negative electrode": pybamm.MeshGenerator(
+                pybamm.SpectralVolume1DSubMesh, {"order": order}
+            ),
+            "separator": pybamm.MeshGenerator(
+                pybamm.SpectralVolume1DSubMesh, {"order": order}
+            ),
+            "positive electrode": pybamm.MeshGenerator(
+                pybamm.SpectralVolume1DSubMesh, {"order": order}
+            ),
+            "current collector": pybamm.SubMesh0D,
+        },
+        var_pts,
+    )
+    for geometry in geometries
+]
 
 # discretise model
 disc_fv = pybamm.Discretisation(meshes[0], models[0].default_spatial_methods)
@@ -61,8 +62,8 @@
         "negative electrode": pybamm.SpectralVolume(order=order),
         "separator": pybamm.SpectralVolume(order=order),
         "positive electrode": pybamm.SpectralVolume(order=order),
-        "current collector": pybamm.ZeroDimensionalSpatialMethod()
-    }
+        "current collector": pybamm.ZeroDimensionalSpatialMethod(),
+    },
 )
 
 disc_fv.process_model(models[0])
diff --git a/pybamm/CITATIONS.txt b/pybamm/CITATIONS.txt
index 206f9cfc94..afff80c0a3 100644
--- a/pybamm/CITATIONS.txt
+++ b/pybamm/CITATIONS.txt
@@ -78,12 +78,15 @@ year={2019},
 publisher={The Electrochemical Society}
 }
 
-@article{Mohtat2020,
-author = {Mohtat, Peyman and Siegel, Jason B and Stefanopoulou, Anna G},
-title = {{High C-rate Differential Expansion and Voltage Model for Li-ion Batteries}},
-journal = {Submitted for publication},
-year = {2019},
-publisher={ECSarXiv}
+@article{mohtat2020differential,
+  title={Differential Expansion and Voltage Model for Li-ion Batteries at Practical Charging Rates},
+  author={Mohtat, Peyman and Lee, Suhak and Sulzer, Valentin and Siegel, Jason B and Stefanopoulou, Anna G},
+  journal={Journal of The Electrochemical Society},
+  volume={167},
+  number={11},
+  pages={110561},
+  year={2020},
+  publisher={IOP Publishing}
 }
 
 @article{scikits-odes,
@@ -204,3 +207,25 @@ primaryClass={physics.app-ph},
   year={2005},
   publisher={IOP Publishing}
 }
+
+@article{lain2019design,
+  title={Design Strategies for High Power vs. High Energy Lithium Ion Cells},
+  author={Lain, Michael J and Brandon, James and Kendrick, Emma},
+  journal={Batteries},
+  volume={5},
+  number={4},
+  pages={64},
+  year={2019},
+  publisher={Multidisciplinary Digital Publishing Institute}
+}
+
+@article{prada2013simplified,
+  title={A simplified electrochemical and thermal aging model of LiFePO4-graphite Li-ion batteries: power and capacity fade simulations},
+  author={Prada, Eric and Di Domenico, D and Creff, Y and Bernard, J and Sauvant-Moynot, Val{\'e}rie and Huet, Fran{\c{c}}ois},
+  journal={Journal of The Electrochemical Society},
+  volume={160},
+  number={4},
+  pages={A616},
+  year={2013},
+  publisher={IOP Publishing}
+}
\ No newline at end of file
diff --git a/pybamm/expression_tree/operations/evaluate.py b/pybamm/expression_tree/operations/evaluate.py
index 7aa14e60ac..cab74588cc 100644
--- a/pybamm/expression_tree/operations/evaluate.py
+++ b/pybamm/expression_tree/operations/evaluate.py
@@ -9,10 +9,12 @@
 
 import numbers
 from platform import system
+
 if system() != "Windows":
     import jax
 
     from jax.config import config
+
     config.update("jax_enable_x64", True)
 
 
@@ -95,18 +97,21 @@ def find_symbols(symbol, constant_symbols, variable_symbols, to_dense=False):
             dummy_eval_left = symbol.children[0].evaluate_for_shape()
             dummy_eval_right = symbol.children[1].evaluate_for_shape()
             if not to_dense and scipy.sparse.issparse(dummy_eval_left):
-                symbol_str = "{0}.multiply({1})"\
-                    .format(children_vars[0], children_vars[1])
+                symbol_str = "{0}.multiply({1})".format(
+                    children_vars[0], children_vars[1]
+                )
             elif not to_dense and scipy.sparse.issparse(dummy_eval_right):
-                symbol_str = "{1}.multiply({0})"\
-                    .format(children_vars[0], children_vars[1])
+                symbol_str = "{1}.multiply({0})".format(
+                    children_vars[0], children_vars[1]
+                )
             else:
                 symbol_str = "{0} * {1}".format(children_vars[0], children_vars[1])
         elif isinstance(symbol, pybamm.Division):
             dummy_eval_left = symbol.children[0].evaluate_for_shape()
             if not to_dense and scipy.sparse.issparse(dummy_eval_left):
-                symbol_str = "{0}.multiply(1/{1})"\
-                    .format(children_vars[0], children_vars[1])
+                symbol_str = "{0}.multiply(1/{1})".format(
+                    children_vars[0], children_vars[1]
+                )
             else:
                 symbol_str = "{0} / {1}".format(children_vars[0], children_vars[1])
 
@@ -114,11 +119,13 @@ def find_symbols(symbol, constant_symbols, variable_symbols, to_dense=False):
             dummy_eval_left = symbol.children[0].evaluate_for_shape()
             dummy_eval_right = symbol.children[1].evaluate_for_shape()
             if not to_dense and scipy.sparse.issparse(dummy_eval_left):
-                symbol_str = "{0}.multiply({1})"\
-                    .format(children_vars[0], children_vars[1])
+                symbol_str = "{0}.multiply({1})".format(
+                    children_vars[0], children_vars[1]
+                )
             elif not to_dense and scipy.sparse.issparse(dummy_eval_right):
-                symbol_str = "{1}.multiply({0})"\
-                    .format(children_vars[0], children_vars[1])
+                symbol_str = "{1}.multiply({0})".format(
+                    children_vars[0], children_vars[1]
+                )
             else:
                 symbol_str = "{0} * {1}".format(children_vars[0], children_vars[1])
 
@@ -294,18 +301,20 @@ def __init__(self, symbol):
         # extract constants in generated function
         for i, symbol_id in enumerate(constants.keys()):
             const_name = id_to_python_variable(symbol_id, True)
-            python_str = '{} = constants[{}]\n'.format(const_name, i) + python_str
+            python_str = "{} = constants[{}]\n".format(const_name, i) + python_str
 
         # constants passed in as an ordered dict, convert to list
         self._constants = list(constants.values())
 
         # indent code
-        python_str = '   ' + python_str
-        python_str = python_str.replace('\n', '\n   ')
+        python_str = "   " + python_str
+        python_str = python_str.replace("\n", "\n   ")
 
         # add function def to first line
-        python_str = 'def evaluate(constants, t=None, y=None, '\
-            'y_dot=None, inputs=None, known_evals=None):\n' + python_str
+        python_str = (
+            "def evaluate(constants, t=None, y=None, "
+            "y_dot=None, inputs=None, known_evals=None):\n" + python_str
+        )
 
         # calculate the final variable that will output the result of calling `evaluate`
         # on `symbol`
@@ -315,21 +324,18 @@ def __init__(self, symbol):
 
         # add return line
         if symbol.is_constant() and isinstance(result_value, numbers.Number):
-            python_str = python_str + '\n   return ' + str(result_value)
+            python_str = python_str + "\n   return " + str(result_value)
         else:
-            python_str = python_str + '\n   return ' + result_var
+            python_str = python_str + "\n   return " + result_var
 
         # store a copy of examine_jaxpr
-        python_str = python_str + \
-            '\nself._evaluate = evaluate'
+        python_str = python_str + "\nself._evaluate = evaluate"
 
         self._python_str = python_str
         self._symbol = symbol
 
         # compile and run the generated python code,
-        compiled_function = compile(
-            python_str, result_var, "exec"
-        )
+        compiled_function = compile(python_str, result_var, "exec")
         exec(compiled_function)
 
     def evaluate(self, t=None, y=None, y_dot=None, inputs=None, known_evals=None):
@@ -377,7 +383,7 @@ def __init__(self, symbol):
         constants, python_str = pybamm.to_python(symbol, debug=False, to_dense=True)
 
         # replace numpy function calls to jax numpy calls
-        python_str = python_str.replace('np.', 'jax.numpy.')
+        python_str = python_str.replace("np.", "jax.numpy.")
 
         # convert all numpy constants to device vectors
         for symbol_id in constants:
@@ -387,18 +393,20 @@ def __init__(self, symbol):
         # extract constants in generated function
         for i, symbol_id in enumerate(constants.keys()):
             const_name = id_to_python_variable(symbol_id, True)
-            python_str = '{} = constants[{}]\n'.format(const_name, i) + python_str
+            python_str = "{} = constants[{}]\n".format(const_name, i) + python_str
 
         # constants passed in as an ordered dict, convert to list
         self._constants = list(constants.values())
 
         # indent code
-        python_str = '   ' + python_str
-        python_str = python_str.replace('\n', '\n   ')
+        python_str = "   " + python_str
+        python_str = python_str.replace("\n", "\n   ")
 
         # add function def to first line
-        python_str = 'def evaluate_jax(constants, t=None, y=None, '\
-            'y_dot=None, inputs=None, known_evals=None):\n' + python_str
+        python_str = (
+            "def evaluate_jax(constants, t=None, y=None, "
+            "y_dot=None, inputs=None, known_evals=None):\n" + python_str
+        )
 
         # calculate the final variable that will output the result of calling `evaluate`
         # on `symbol`
@@ -408,18 +416,15 @@ def __init__(self, symbol):
 
         # add return line
         if symbol.is_constant() and isinstance(result_value, numbers.Number):
-            python_str = python_str + '\n   return ' + str(result_value)
+            python_str = python_str + "\n   return " + str(result_value)
         else:
-            python_str = python_str + '\n   return ' + result_var
+            python_str = python_str + "\n   return " + result_var
 
         # store a copy of examine_jaxpr
-        python_str = python_str + \
-            '\nself._evaluate_jax = evaluate_jax'
+        python_str = python_str + "\nself._evaluate_jax = evaluate_jax"
 
         # compile and run the generated python code,
-        compiled_function = compile(
-            python_str, result_var, "exec"
-        )
+        compiled_function = compile(python_str, result_var, "exec")
         exec(compiled_function)
 
         self._jit_evaluate = jax.jit(self._evaluate_jax, static_argnums=(0, 4, 5))
diff --git a/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/LFP_electrolyte_exchange_current_density_kashkooli2017.py b/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/LFP_electrolyte_exchange_current_density_kashkooli2017.py
new file mode 100644
index 0000000000..083c352806
--- /dev/null
+++ b/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/LFP_electrolyte_exchange_current_density_kashkooli2017.py
@@ -0,0 +1,37 @@
+from pybamm import exp, constants, Parameter
+
+
+def LFP_electrolyte_exchange_current_density_kashkooli2017(c_e, c_s_surf, T):  # , 1
+    """
+    Exchange-current density for Butler-Volmer reactions between LFP and electrolyte
+
+    References
+    ----------
+    .. [1] Kashkooli, A. G., Amirfazli, A., Farhad, S., Lee, D. U., Felicelli, S., Park,
+    H. W., ... & Chen, Z. (2017). Representative volume element model of lithium-ion
+    battery electrodes based on X-ray nano-tomography. Journal of Applied
+    Electrochemistry, 47(3), 281-293.
+
+    Parameters
+    ----------
+    c_e : :class:`pybamm.Symbol`
+        Electrolyte concentration [mol.m-3]
+    c_s_surf : :class:`pybamm.Symbol`
+        Particle concentration [mol.m-3]
+    T : :class:`pybamm.Symbol`
+        Temperature [K]
+
+    Returns
+    -------
+    :class:`pybamm.Symbol`
+        Exchange-current density [A.m-2]
+    """
+
+    m_ref = 6 * 10 ** (-7)  # (A/m2)(mol/m3)**1.5 - includes ref concentrations
+    E_r = 39570
+    arrhenius = exp(E_r / constants.R * (1 / 298.15 - 1 / T))
+    c_p_max = Parameter("Maximum concentration in positive electrode [mol.m-3]")
+
+    return (
+        m_ref * arrhenius * c_e ** 0.5 * c_s_surf ** 0.5 * (c_p_max - c_s_surf) ** 0.5
+    )
diff --git a/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/LFP_ocp_ashfar2017.py b/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/LFP_ocp_ashfar2017.py
new file mode 100644
index 0000000000..9695acfe4b
--- /dev/null
+++ b/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/LFP_ocp_ashfar2017.py
@@ -0,0 +1,24 @@
+from pybamm import exp
+
+
+def LFP_ocp_ashfar2017(sto):
+    """
+    Open-circuit potential for LFP
+
+    References
+    ----------
+    .. [1] Afshar, S., Morris, K., & Khajepour, A. (2017). Efficient electrochemical
+    model for lithium-ion cells. arXiv preprint arXiv:1709.03970.
+
+    Parameters
+    ----------
+    sto : :class:`pybamm.Symbol`
+       Stochiometry of material (li-fraction)
+
+    """
+
+    c1 = -150 * sto
+    c2 = -30 * (1 - sto)
+    k = 3.4077 - 0.020269 * sto + 0.5 * exp(c1) - 0.9 * exp(c2)
+
+    return k
diff --git a/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/README.md b/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/README.md
new file mode 100644
index 0000000000..681185b06c
--- /dev/null
+++ b/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/README.md
@@ -0,0 +1,7 @@
+# Lithium Cobalt Oxide cathode parameters
+
+Parameters for an LFP cathode, from the paper
+
+> Prada, E., Di Domenico, D., Creff, Y., Bernard, J., Sauvant-Moynot, V., & Huet, F. (2013). A simplified electrochemical and thermal aging model of LiFePO4-graphite Li-ion batteries: power and capacity fade simulations. [Journal of The Electrochemical Society](https://doi.org/10.1149/2.053304jes), 160(4), A616.
+
+and references therein. The functions used for OCP and exchange-current density are from separate references (documented within the functions), to provide better fit to data
diff --git a/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/parameters.csv b/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/parameters.csv
new file mode 100644
index 0000000000..c500cc53e8
--- /dev/null
+++ b/pybamm/input/parameters/lithium-ion/cathodes/LFP_Prada2013/parameters.csv
@@ -0,0 +1,35 @@
+Name [units],Value,Reference,Notes
+# Empty rows and rows starting with # will be ignored,,,
+,,,
+# Electrode properties,,,
+Positive electrode conductivity [S.m-1],0.33795074,"Calculation, consistant with Hosseinzadeh 2018",
+Maximum concentration in positive electrode [mol.m-3],22806,Prada Thermal Aging,
+Positive electrode diffusivity [m2.s-1],5.90E-18,Prada Fast Charging,
+Positive electrode OCP [V],[function]LFP_ocp_ashfar2017,,
+,,,
+# Microstructure,,,
+Positive electrode porosity,0.12728395,Calculation minimized to Severson,
+Positive electrode active material volume fraction,0.28485556,Calculation minimized to Severson,
+Positive particle radius [m],1.00E-08,Calculation minimized to Severson,
+Positive particle distribution in x,1,,
+Positive electrode Bruggeman coefficient (electrode),1.5,,
+Positive electrode Bruggeman coefficient (electrolyte),1.5,,
+,,,
+# Interfacial reactions,,,
+Positive electrode cation signed stoichiometry,-1,,
+Positive electrode electrons in reaction,1,,
+Reference OCP vs SHE in the positive electrode [V],,adjust,
+Positive electrode charge transfer coefficient,0.5,adjust,
+Positive electrode double-layer capacity [F.m-2],0.2,adjust,
+,,,
+# Density,,,
+Positive electrode density [kg.m-3],2341.17,default,
+,,,
+,,,
+# Thermal parameters,,,
+Positive electrode specific heat capacity [J.kg-1.K-1],1100,,
+Positive electrode thermal conductivity [W.m-1.K-1],2.1,,
+Positive electrode OCP entropic change [V.K-1],0,,
+,,,
+# exchange-current density,,,
+Positive electrode exchange-current density [A.m-2],[function]LFP_electrolyte_exchange_current_density_kashkooli2017,,
diff --git a/pybamm/input/parameters/lithium-ion/cathodes/lico2_Ramadass2004/lico2_ocp_Ramadass2004.py b/pybamm/input/parameters/lithium-ion/cathodes/lico2_Ramadass2004/lico2_ocp_Ramadass2004.py
index 7a5e34a622..8541121a84 100644
--- a/pybamm/input/parameters/lithium-ion/cathodes/lico2_Ramadass2004/lico2_ocp_Ramadass2004.py
+++ b/pybamm/input/parameters/lithium-ion/cathodes/lico2_Ramadass2004/lico2_ocp_Ramadass2004.py
@@ -1,4 +1,3 @@
-
 def lico2_ocp_Ramadass2004(sto):
     """
     Lithium Cobalt Oxide (LiCO2) Open Circuit Potential (OCP) as a a function of the
@@ -21,12 +20,20 @@ def lico2_ocp_Ramadass2004(sto):
     stretch = 1.13
     sto = stretch * sto
 
-    u_eq = ((- 4.656 + 88.669 * (sto ** 2)
-             - 401.119 * (sto ** 4) + 342.909 * (sto ** 6)
-             - 462.471 * (sto ** 8) + 433.434 * (sto ** 10)) / (
-                 - 1 + 18.933 * (sto ** 2) - 79.532 * (sto ** 4)
-                 + 37.311 * (sto ** 6) - 73.083 * (sto ** 8)
-                 + 95.96 * (sto ** 10))
-            )
+    u_eq = (
+        -4.656
+        + 88.669 * (sto ** 2)
+        - 401.119 * (sto ** 4)
+        + 342.909 * (sto ** 6)
+        - 462.471 * (sto ** 8)
+        + 433.434 * (sto ** 10)
+    ) / (
+        -1
+        + 18.933 * (sto ** 2)
+        - 79.532 * (sto ** 4)
+        + 37.311 * (sto ** 6)
+        - 73.083 * (sto ** 8)
+        + 95.96 * (sto ** 10)
+    )
 
     return u_eq
diff --git a/pybamm/input/parameters/lithium-ion/cells/A123_Lain2019/README.md b/pybamm/input/parameters/lithium-ion/cells/A123_Lain2019/README.md
new file mode 100644
index 0000000000..8c1bc70405
--- /dev/null
+++ b/pybamm/input/parameters/lithium-ion/cells/A123_Lain2019/README.md
@@ -0,0 +1,5 @@
+# Pouch cell parameters
+
+Parameters for an A123 LFP cell, from the paper
+
+> Lain, M. J., Brandon, J., & Kendrick, E. (2019). Design Strategies for High Power vs. High Energy Lithium Ion Cells. [Batteries](https://doi.org/10.3390/batteries5040064), 5(4), 64.
diff --git a/pybamm/input/parameters/lithium-ion/cells/A123_Lain2019/parameters.csv b/pybamm/input/parameters/lithium-ion/cells/A123_Lain2019/parameters.csv
new file mode 100644
index 0000000000..3056894a56
--- /dev/null
+++ b/pybamm/input/parameters/lithium-ion/cells/A123_Lain2019/parameters.csv
@@ -0,0 +1,37 @@
+Name [units],Value,Reference,Notes
+# Empty rows and rows starting with ‘#’ will be ignored,,,
+,,,
+# Macroscale geometry,,,
+Negative current collector thickness [m],1.00E-05,Lain 2019,
+Negative electrode thickness [m],3.60E-05,Lain 2019,
+Separator thickness [m],1.80E-05,Lain 2019,
+Positive electrode thickness [m],8.10E-05,Lain 2019,
+Positive current collector thickness [m],1.90E-05,Lain 2019,
+Electrode height [m],6.49E-02,A123 18650 datasheet,Not needed for 1D
+Electrode width [m],1.78,fit to Severson 2019,Not needed for 1D
+Negative tab width [m],0.04,default,Not needed for 1D
+Negative tab centre y-coordinate [m],0.06,default,Not needed for 1D
+Negative tab centre z-coordinate [m],0.137,default,Not needed for 1D
+Positive tab width [m],0.04,default,Not needed for 1D
+Positive tab centre y-coordinate [m],0.147,default,Not needed for 1D
+Positive tab centre z-coordinate [m],0.137,default,Not needed for 1D
+,,,
+# Current collector properties ,,,
+Negative current collector conductivity [S.m-1],58411000,CRC Handbook,copper
+Positive current collector conductivity [S.m-1],36914000,CRC Handbook,aluminium
+,,,
+# Density,,,
+Negative current collector density [kg.m-3],8960,CRC Handbook,copper
+Positive current collector density [kg.m-3],2700,CRC Handbook,aluminium
+,,,
+# Specific heat capacity,,,
+Negative current collector specific heat capacity [J.kg-1.K-1],385,CRC Handbook,copper
+Positive current collector specific heat capacity [J.kg-1.K-1],897,CRC Handbook,aluminium
+,,,
+# Thermal conductivity,,,
+Negative current collector thermal conductivity [W.m-1.K-1],401,CRC Handbook,copper
+Positive current collector thermal conductivity [W.m-1.K-1],237,CRC Handbook,aluminium
+,,,
+# Electrical,,,
+Cell capacity [A.h],1.1,Lain 2019,
+Typical current [A],30,Lain 2019,
diff --git a/pybamm/input/parameters/lithium-ion/cells/UMBL_Mohtat2020/parameters.csv b/pybamm/input/parameters/lithium-ion/cells/UMBL_Mohtat2020/parameters.csv
index bfe0cb39ba..7f2492965d 100644
--- a/pybamm/input/parameters/lithium-ion/cells/UMBL_Mohtat2020/parameters.csv
+++ b/pybamm/input/parameters/lithium-ion/cells/UMBL_Mohtat2020/parameters.csv
@@ -10,7 +10,7 @@ Positive current collector thickness [m],2.5E-05,Scott Moura FastDFN,no info fro
 Electrode height [m],1,KOKAM SLPB78205130H,Not needed for 1D
 Electrode width [m],0.2050,KOKAM SLPB78205130H,Not needed for 1D
 Cell cooling surface area [m2],0.41,,pouch
-Cell volume [m3]3.92E-5,,pouch
+Cell volume [m3],3.92E-5,,pouch
 ,,,
 # Current collector properties ,,,
 Negative current collector conductivity [S.m-1],59600000,LIONSIMBA,carbon
diff --git a/pybamm/input/parameters/lithium-ion/electrolytes/lipf6_Nyman2008/parameters.csv b/pybamm/input/parameters/lithium-ion/electrolytes/lipf6_Nyman2008/parameters.csv
index c6e77e42a6..238c6fb4c8 100644
--- a/pybamm/input/parameters/lithium-ion/electrolytes/lipf6_Nyman2008/parameters.csv
+++ b/pybamm/input/parameters/lithium-ion/electrolytes/lipf6_Nyman2008/parameters.csv
@@ -1,9 +1,9 @@
-Name [units],Value,Reference,Notes
-# Empty rows and rows starting with ‘#’ will be ignored,,,
-,,,
-# Electrolyte properties,,,
-Typical electrolyte concentration [mol.m-3],1000,Chen 2020,
-Cation transference number,0.2594,Chen 2020,
-1 + dlnf/dlnc,1,,
-Electrolyte diffusivity [m2.s-1],[function]electrolyte_diffusivity_Nyman2008,Nyman 2008," "
-Electrolyte conductivity [S.m-1],[function]electrolyte_conductivity_Nyman2008,Nyman 2008," "
\ No newline at end of file
+Name [units],Value,Reference,Notes
+# Empty rows and rows starting with ‘#’ will be ignored,,,
+,,,
+# Electrolyte properties,,,
+Typical electrolyte concentration [mol.m-3],1000,Chen 2020,
+Cation transference number,0.2594,Chen 2020,
+1 + dlnf/dlnc,1,,
+Electrolyte diffusivity [m2.s-1],[function]electrolyte_diffusivity_Nyman2008,Nyman 2008, 
+Electrolyte conductivity [S.m-1],[function]electrolyte_conductivity_Nyman2008,Nyman 2008, 
diff --git a/pybamm/input/parameters/lithium-ion/experiments/1C_discharge_from_full_Chen2020/parameters.csv b/pybamm/input/parameters/lithium-ion/experiments/1C_discharge_from_full_Chen2020/parameters.csv
index 685cdbf661..e9dd8657c3 100644
--- a/pybamm/input/parameters/lithium-ion/experiments/1C_discharge_from_full_Chen2020/parameters.csv
+++ b/pybamm/input/parameters/lithium-ion/experiments/1C_discharge_from_full_Chen2020/parameters.csv
@@ -16,4 +16,4 @@ Upper voltage cut-off [V],4.4,,
 Initial concentration in negative electrode [mol.m-3],29866,Chen 2020,
 Initial concentration in positive electrode [mol.m-3],17038,Chen 2020,
 Initial concentration in electrolyte [mol.m-3],1000,Chen 2020,
-Initial temperature [K],298.15,,
\ No newline at end of file
+Initial temperature [K],298.15,,
diff --git a/pybamm/input/parameters/lithium-ion/experiments/4C_discharge_from_full_Prada2013/README.md b/pybamm/input/parameters/lithium-ion/experiments/4C_discharge_from_full_Prada2013/README.md
new file mode 100644
index 0000000000..997e2e9233
--- /dev/null
+++ b/pybamm/input/parameters/lithium-ion/experiments/4C_discharge_from_full_Prada2013/README.md
@@ -0,0 +1,7 @@
+# 1C discharge from full
+
+Discharge lithium-ion battery from full charge at 4C, using the initial conditions from the paper
+
+> Prada, E., Di Domenico, D., Creff, Y., Bernard, J., Sauvant-Moynot, V., & Huet, F. (2013). A simplified electrochemical and thermal aging model of LiFePO4-graphite Li-ion batteries: power and capacity fade simulations. [Journal of The Electrochemical Society](https://doi.org/10.1149/2.053304jes), 160(4), A616.
+
+and references therein.
diff --git a/pybamm/input/parameters/lithium-ion/experiments/4C_discharge_from_full_Prada2013/parameters.csv b/pybamm/input/parameters/lithium-ion/experiments/4C_discharge_from_full_Prada2013/parameters.csv
new file mode 100644
index 0000000000..b6abadce83
--- /dev/null
+++ b/pybamm/input/parameters/lithium-ion/experiments/4C_discharge_from_full_Prada2013/parameters.csv
@@ -0,0 +1,21 @@
+Name [units],Value,Reference,Notes
+# Empty rows and rows starting with ‘#’ will be ignored,,,
+,,,
+# Temperature,,,
+Reference temperature [K],298.15,25C,
+Heat transfer coefficient [W.m-2.K-1],10,,
+Ambient temperature [K],298.15,,
+,,,
+,,,
+# 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],2,,
+Upper voltage cut-off [V],4.4,,
+Current function [A],4.4,,
+,,,
+# Initial conditions,,,
+Initial concentration in negative electrode [mol.m-3],28831.45783,Minimized to Severson Data,
+Initial concentration in positive electrode [mol.m-3],35.3766672,Minimized to Severson Data,
+Initial concentration in electrolyte [mol.m-3],1200,Lain,
+Initial temperature [K],298.15,,
diff --git a/pybamm/meshes/one_dimensional_submeshes.py b/pybamm/meshes/one_dimensional_submeshes.py
index ae12a3f2f0..f4f0bdd66f 100644
--- a/pybamm/meshes/one_dimensional_submeshes.py
+++ b/pybamm/meshes/one_dimensional_submeshes.py
@@ -340,8 +340,7 @@ def __init__(self, lims, npts, edges=None, order=2):
 
         # default: Spectral Volumes of equal size
         if edges is None:
-            edges = np.linspace(spatial_lims["min"], spatial_lims["max"],
-                                npts + 1)
+            edges = np.linspace(spatial_lims["min"], spatial_lims["max"], npts + 1)
         # check that npts + 1 equals number of user-supplied edges
         elif (npts + 1) != len(edges):
             raise pybamm.GeometryError(
@@ -369,16 +368,30 @@ def __init__(self, lims, npts, edges=None, order=2):
 
         coord_sys = spatial_var.coord_sys
 
-        cv_edges = np.array([edges[0]] + [
-            x
-            for (a, b) in zip(edges[:-1], edges[1:])
-            for x in np.flip(
-                a + 0.5 * (b - a) * (1 + np.sin(np.pi * np.array(
-                    [((order + 1) - 1 - 2 * i) / (2 * (order + 1) - 2)
-                     for i in range(order + 1)]
-                )))
-            )[1:]
-        ])
+        cv_edges = np.array(
+            [edges[0]]
+            + [
+                x
+                for (a, b) in zip(edges[:-1], edges[1:])
+                for x in np.flip(
+                    a
+                    + 0.5
+                    * (b - a)
+                    * (
+                        1
+                        + np.sin(
+                            np.pi
+                            * np.array(
+                                [
+                                    ((order + 1) - 1 - 2 * i) / (2 * (order + 1) - 2)
+                                    for i in range(order + 1)
+                                ]
+                            )
+                        )
+                    )
+                )[1:]
+            ]
+        )
 
         self.sv_edges = edges
         self.sv_nodes = (edges[:-1] + edges[1:]) / 2
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py
index cb98a2dec9..050f472c5c 100644
--- a/pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py
@@ -33,7 +33,9 @@ class BasicDFNHalfCell(BaseModel):
     """
 
     def __init__(
-        self, name="Doyle-Fuller-Newman half cell model", options=None,
+        self,
+        name="Doyle-Fuller-Newman half cell model",
+        options=None,
     ):
         super().__init__({}, name)
         pybamm.citations.register("marquis2019asymptotic")
@@ -438,8 +440,8 @@ def __init__(
             "Negative particle concentration": c_s_n,
             "Negative particle surface concentration [mol.m-3]": param.c_n_max
             * c_s_surf_n,
-            "X-averaged negative particle surface concentration [mol.m-3]":
-            param.c_n_max * c_s_surf_n_av,
+            "X-averaged negative particle surface concentration "
+            "[mol.m-3]": param.c_n_max * c_s_surf_n_av,
             "Negative particle concentration [mol.m-3]": param.c_n_max * c_s_n,
             "Electrolyte concentration": c_e,
             "Electrolyte concentration [mol.m-3]": param.c_e_typ * c_e,
@@ -448,8 +450,8 @@ def __init__(
             "Positive particle concentration": c_s_p,
             "Positive particle surface concentration [mol.m-3]": param.c_p_max
             * c_s_surf_p,
-            "X-averaged positive particle surface concentration [mol.m-3]":
-            param.c_p_max * c_s_surf_p_av,
+            "X-averaged positive particle surface concentration "
+            "[mol.m-3]": param.c_p_max * c_s_surf_p_av,
             "Positive particle concentration [mol.m-3]": param.c_p_max * c_s_p,
             "Current [A]": I,
             "Negative electrode potential": phi_s_n,
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
index 8e027d0e76..ecaa68837c 100644
--- a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
@@ -170,4 +170,3 @@ def __init__(self, name="Single Particle Model"):
 
     def new_copy(self, build=False):
         return pybamm.BaseModel.new_copy(self)
-
diff --git a/pybamm/parameters/parameter_sets.py b/pybamm/parameters/parameter_sets.py
index d14bb97ee2..7eb90bbd3c 100644
--- a/pybamm/parameters/parameter_sets.py
+++ b/pybamm/parameters/parameter_sets.py
@@ -10,42 +10,56 @@
 --------------------------
     * Chen2020 :
         C.-H. Chen, F. Brosa Planella, K. O’Regan, D. Gastol, W. D. Widanage, and E.
-        Kendrick. “Development of Experimental Techniques for Parameterization of
-        Multi-scale Lithium-ion Battery Models.” Journal of the Electrochemical Society,
+        Kendrick. "Development of Experimental Techniques for Parameterization of
+        Multi-scale Lithium-ion Battery Models." Journal of the Electrochemical Society,
         167(8), 080534 (2020).
     * Ecker2015 :
         M. Ecker, T. K. D. Tran, P. Dechent, S. Käbitz, A. Warnecke, and D. U. Sauer.
-        “Parameterization of a Physico-Chemical Model of a Lithium-Ion Battery. I.
-        Determination of Parameters.” Journal of the Electrochemical Society, 162(9),
+        "Parameterization of a Physico-Chemical Model of a Lithium-Ion Battery. I.
+        Determination of Parameters." Journal of the Electrochemical Society, 162(9),
         A1836-A1848 (2015).
     * Marquis2019 :
-        S. G. Marquis, V. Sulzer, R. Timms, C. P. Please and S. J. Chapman. “An
-        asymptotic derivation of a single particle model with electrolyte.” Journal of
+        S. G. Marquis, V. Sulzer, R. Timms, C. P. Please and S. J. Chapman. "An
+        asymptotic derivation of a single particle model with electrolyte." Journal of
         the Electrochemical Society, 166(15), A3693–A3706 (2019).
     * Mohtat2020 :
-        Submitted for publication.
+        P. Mohtat, S. Lee, V. Sulzer, J. B. Siegel & A. G. Stefanopoulou. "Differential
+        Expansion and Voltage Model for Li-ion Batteries at Practical Charging Rates".
+        Journal of The Electrochemical Society, 167(11), 110561 (2020).
     * NCA_Kim2011 :
         G. H. Kim, K. Smith, K. J. Lee, S. Santhanagopalan, and A. Pesaran.
-        “Multi-domain modeling of lithium-ion batteries encompassing multi-physics in
-        varied length scales.” Journal of The Electrochemical Society, 158(8), A955-A969
+        "Multi-domain modeling of lithium-ion batteries encompassing multi-physics in
+        varied length scales." Journal of The Electrochemical Society, 158(8), A955-A969
         (2011).
-    * Ramadass 2004 :
-        P. Ramadass, B. Haran, P. M. Gomadam, R. White, and B. N. Popov. “Development
-        of First Principles Capacity Fade Model for Li-Ion Cells.” Journal of the
+    * Ramadass2004 :
+        P. Ramadass, B. Haran, P. M. Gomadam, R. White, and B. N. Popov. "Development
+        of First Principles Capacity Fade Model for Li-Ion Cells." Journal of the
         Electrochemical Society, 151(2), A196-A203 (2004).
+    * Prada2013
+        C.-H. Chen, F. Brosa Planella, K. O’Regan, D. Gastol, W. D. Widanage, and E.
+        Kendrick. "Development of Experimental Techniques for Parameterization of
+        Multi-scale Lithium-ion Battery Models." Journal of the Electrochemical Society,
+        167(8), 080534 (2020).
+        E. Prada, D., Di Domenico, Y., Creff, J., Bernard, V., Sauvant-Moynot, and F.
+        Huet. "A simplified electrochemical and thermal aging model of LiFePO4-graphite
+        Li-ion batteries: power and capacity fade simulations". Journal of The
+        Electrochemical Society, 160(4), A616 (2013).
+        M. J. Lain, J., Brandon, and E. Kendrick. "Design Strategies for High Power vs.
+        High Energy Lithium Ion Cells". Batteries, 5(4), 64 (2019).
 
 Lead-acid parameter sets
 --------------------------
     * Sulzer2019 :
-        V. Sulzer, S. J. Chapman, C. P. Please, D. A. Howey, and C. W. Monroe, “Faster
+        V. Sulzer, S. J. Chapman, C. P. Please, D. A. Howey, and C. W. Monroe, "Faster
         lead-acid battery simulations from porous-electrode theory: Part I. Physical
-        model.” Journal of the Electrochemical Society, 166(12), 2363 (2019).
+        model." Journal of the Electrochemical Society, 166(12), 2363 (2019).
 
 """
 
 #
 # Lithium-ion
 #
+
 NCA_Kim2011 = {
     "chemistry": "lithium-ion",
     "cell": "Kim2011",
@@ -103,8 +117,9 @@
     "electrolyte": "LiPF6_Mohtat2020",
     "experiment": "1C_charge_from_empty_Mohtat2020",
     "sei": "example",
-    "citation": "Mohtat2020",
+    "citation": "mohtat2020differential",
 }
+
 Ramadass2004 = {
     "chemistry": "lithium-ion",
     "cell": "sony_Ramadass2004",
@@ -117,6 +132,17 @@
     "citation": "marquis2019asymptotic",
 }
 
+Prada2013 = {
+    "chemistry": "lithium-ion",
+    "cell": "A123_Lain2019",
+    "anode": "graphite_Chen2020",
+    "separator": "separator_Chen2020",
+    "cathode": "LFP_Prada2013",
+    "electrolyte": "lipf6_Nyman2008",
+    "experiment": "4C_discharge_from_full_Prada2013",
+    "citation": ["Chen2020", "lain2019design", "prada2013simplified"],
+}
+
 #
 # Lead-acid
 #
diff --git a/pybamm/solvers/jax_bdf_solver.py b/pybamm/solvers/jax_bdf_solver.py
index d9024117b9..2600e71f24 100644
--- a/pybamm/solvers/jax_bdf_solver.py
+++ b/pybamm/solvers/jax_bdf_solver.py
@@ -82,9 +82,7 @@ def fun_bind_inputs(y, t):
     t0 = t_eval[0]
     h0 = t_eval[1] - t0
 
-    stepper = _bdf_init(
-        fun_bind_inputs, jac_bind_inputs, mass, t0, y0, h0, rtol, atol
-    )
+    stepper = _bdf_init(fun_bind_inputs, jac_bind_inputs, mass, t0, y0, h0, rtol, atol)
     i = 0
     y_out = jnp.empty((len(t_eval), len(y0)), dtype=y0.dtype)
 
@@ -101,30 +99,49 @@ def body_fun(state):
 
         def for_body(j, y_out):
             t = t_eval[j]
-            y_out = jax.ops.index_update(y_out, jax.ops.index[j, :],
-                                         _bdf_interpolate(stepper, t))
+            y_out = jax.ops.index_update(
+                y_out, jax.ops.index[j, :], _bdf_interpolate(stepper, t)
+            )
             return y_out
 
         y_out = jax.lax.fori_loop(i, index, for_body, y_out)
         return [stepper, t_eval, index, y_out]
 
-    stepper, t_eval, i, y_out = jax.lax.while_loop(cond_fun, body_fun,
-                                                   init_state)
+    stepper, t_eval, i, y_out = jax.lax.while_loop(cond_fun, body_fun, init_state)
     return y_out
 
 
 BDFInternalStates = [
-    't', 'atol', 'rtol', 'M', 'newton_tol', 'order', 'h', 'n_equal_steps', 'D',
-    'y0', 'scale_y0', 'kappa', 'gamma', 'alpha', 'c', 'error_const', 'J', 'LU', 'U',
-    'psi', 'n_function_evals', 'n_jacobian_evals', 'n_lu_decompositions', 'n_steps',
-    'consistent_y0_failed'
+    "t",
+    "atol",
+    "rtol",
+    "M",
+    "newton_tol",
+    "order",
+    "h",
+    "n_equal_steps",
+    "D",
+    "y0",
+    "scale_y0",
+    "kappa",
+    "gamma",
+    "alpha",
+    "c",
+    "error_const",
+    "J",
+    "LU",
+    "U",
+    "psi",
+    "n_function_evals",
+    "n_jacobian_evals",
+    "n_lu_decompositions",
+    "n_steps",
+    "consistent_y0_failed",
 ]
-BDFState = collections.namedtuple('BDFState', BDFInternalStates)
+BDFState = collections.namedtuple("BDFState", BDFInternalStates)
 
 jax.tree_util.register_pytree_node(
-    BDFState,
-    lambda xs: (tuple(xs), None),
-    lambda _, xs: BDFState(*xs)
+    BDFState, lambda xs: (tuple(xs), None), lambda _, xs: BDFState(*xs)
 )
 
 
@@ -161,56 +178,56 @@ def _bdf_init(fun, jac, mass, t0, y0, h0, rtol, atol):
     """
 
     state = {}
-    state['t'] = t0
-    state['atol'] = atol
-    state['rtol'] = rtol
-    state['M'] = mass
+    state["t"] = t0
+    state["atol"] = atol
+    state["rtol"] = rtol
+    state["M"] = mass
     EPS = jnp.finfo(y0.dtype).eps
-    state['newton_tol'] = jnp.maximum(10 * EPS / rtol, jnp.minimum(0.03, rtol ** 0.5))
+    state["newton_tol"] = jnp.maximum(10 * EPS / rtol, jnp.minimum(0.03, rtol ** 0.5))
 
     scale_y0 = atol + rtol * jnp.abs(y0)
     y0, not_converged = _select_initial_conditions(
-        fun, mass, t0, y0, state['newton_tol'], scale_y0
+        fun, mass, t0, y0, state["newton_tol"], scale_y0
     )
-    state['consistent_y0_failed'] = not_converged
+    state["consistent_y0_failed"] = not_converged
 
     f0 = fun(y0, t0)
     order = 1
-    state['order'] = order
-    state['h'] = _select_initial_step(atol, rtol, fun, t0, y0, f0, h0)
-    state['n_equal_steps'] = 0
+    state["order"] = order
+    state["h"] = _select_initial_step(atol, rtol, fun, t0, y0, f0, h0)
+    state["n_equal_steps"] = 0
     D = jnp.empty((MAX_ORDER + 1, len(y0)), dtype=y0.dtype)
     D = jax.ops.index_update(D, jax.ops.index[0, :], y0)
-    D = jax.ops.index_update(D, jax.ops.index[1, :], f0 * state['h'])
-    state['D'] = D
-    state['y0'] = y0
-    state['scale_y0'] = scale_y0
+    D = jax.ops.index_update(D, jax.ops.index[1, :], f0 * state["h"])
+    state["D"] = D
+    state["y0"] = y0
+    state["scale_y0"] = scale_y0
 
     # kappa values for difference orders, taken from Table 1 of [1]
     kappa = jnp.array([0, -0.1850, -1 / 9, -0.0823, -0.0415, 0])
     gamma = jnp.hstack((0, jnp.cumsum(1 / jnp.arange(1, MAX_ORDER + 1))))
     alpha = 1.0 / ((1 - kappa) * gamma)
-    c = state['h'] * alpha[order]
+    c = state["h"] * alpha[order]
     error_const = kappa * gamma + 1 / jnp.arange(1, MAX_ORDER + 2)
 
-    state['kappa'] = kappa
-    state['gamma'] = gamma
-    state['alpha'] = alpha
-    state['c'] = c
-    state['error_const'] = error_const
+    state["kappa"] = kappa
+    state["gamma"] = gamma
+    state["alpha"] = alpha
+    state["c"] = c
+    state["error_const"] = error_const
 
     J = jac(y0, t0)
-    state['J'] = J
+    state["J"] = J
 
-    state['LU'] = jax.scipy.linalg.lu_factor(state['M'] - c * J)
+    state["LU"] = jax.scipy.linalg.lu_factor(state["M"] - c * J)
 
-    state['U'] = _compute_R(order, 1)
-    state['psi'] = None
+    state["U"] = _compute_R(order, 1)
+    state["psi"] = None
 
-    state['n_function_evals'] = 2
-    state['n_jacobian_evals'] = 1
-    state['n_lu_decompositions'] = 1
-    state['n_steps'] = 0
+    state["n_function_evals"] = 2
+    state["n_jacobian_evals"] = 1
+    state["n_lu_decompositions"] = 1
+    state["n_steps"] = 0
 
     tuple_state = BDFState(*[state[k] for k in BDFInternalStates])
     y0, scale_y0 = _predict(tuple_state, D)
@@ -232,8 +249,7 @@ def _compute_R(order, factor):
     I = jnp.arange(1, MAX_ORDER + 1).reshape(-1, 1)
     J = jnp.arange(1, MAX_ORDER + 1)
     M = jnp.empty((MAX_ORDER + 1, MAX_ORDER + 1))
-    M = jax.ops.index_update(M, jax.ops.index[1:, 1:],
-                             (I - 1 - factor * J) / I)
+    M = jax.ops.index_update(M, jax.ops.index[1:, 1:], (I - 1 - factor * J) / I)
     M = jax.ops.index_update(M, jax.ops.index[0], 1)
     R = jnp.cumprod(M, axis=0)
 
@@ -242,7 +258,7 @@ def _compute_R(order, factor):
 
 def _select_initial_conditions(fun, M, t0, y0, tol, scale_y0):
     # identify algebraic variables as zeros on diagonal
-    algebraic_variables = jnp.diag(M == 0.)
+    algebraic_variables = jnp.diag(M == 0.0)
 
     # if all differentiable variables then return y0 (can use normal python if since M
     # is static)
@@ -283,24 +299,24 @@ def while_body(while_state):
         k, converged, dy_norm_old, d, y_a = while_state
         f_eval = fun_a(y_a)
         dy = jax.scipy.linalg.lu_solve(LU, f_eval)
-        dy_norm = jnp.sqrt(jnp.mean((dy / scale_y0_a)**2))
+        dy_norm = jnp.sqrt(jnp.mean((dy / scale_y0_a) ** 2))
         rate = dy_norm / dy_norm_old
 
         d += dy
         y_a = y0_a + d
 
         # if converged then break out of iteration early
-        pred = dy_norm_old >= 0.
+        pred = dy_norm_old >= 0.0
         pred *= rate / (1 - rate) * dy_norm < tol
-        converged = (dy_norm == 0.) + pred
+        converged = (dy_norm == 0.0) + pred
 
         dy_norm_old = dy_norm
 
         return [k + 1, converged, dy_norm_old, d, y_a]
 
-    k, converged, dy_norm_old, d, y_a = jax.lax.while_loop(while_cond,
-                                                           while_body,
-                                                           while_state)
+    k, converged, dy_norm_old, d, y_a = jax.lax.while_loop(
+        while_cond, while_body, while_state
+    )
     y_tilde = jax.ops.index_update(y0, algebraic_variables, y_a)
 
     return y_tilde, converged
@@ -322,7 +338,7 @@ def _select_initial_step(atol, rtol, fun, t0, y0, f0, h0):
     scale = atol + jnp.abs(y0) * rtol
     y1 = y0 + h0 * f0
     f1 = fun(y1, t0 + h0)
-    d2 = jnp.sqrt(jnp.mean(((f1 - f0) / scale)**2))
+    d2 = jnp.sqrt(jnp.mean(((f1 - f0) / scale) ** 2))
     order = 1
     h1 = h0 * d2 ** (-1 / (order + 1))
     return jnp.minimum(100 * h0, h1)
@@ -351,10 +367,7 @@ def _update_psi(state, D):
     subGamma = jnp.where(orders > 0, jnp.where(orders <= order, state.gamma, 0), 0)
     orders = jnp.repeat(orders.reshape(-1, 1), n, axis=1)
     subD = jnp.where(orders > 0, jnp.where(orders <= order, D, 0), 0)
-    psi = jnp.dot(
-        subD.T,
-        subGamma
-    ) * state.alpha[order]
+    psi = jnp.dot(subD.T, subGamma) * state.alpha[order]
     return psi
 
 
@@ -373,10 +386,8 @@ def _update_difference_for_next_step(state, d):
     """
     order = state.order
     D = state.D
-    D = jax.ops.index_update(D, jax.ops.index[order + 2],
-                             d - D[order + 1])
-    D = jax.ops.index_update(D, jax.ops.index[order + 1],
-                             d)
+    D = jax.ops.index_update(D, jax.ops.index[order + 2], d - D[order + 1])
+    D = jax.ops.index_update(D, jax.ops.index[order + 1], d)
     i = order
     while_state = [i, D]
 
@@ -386,8 +397,7 @@ def while_cond(while_state):
 
     def while_body(while_state):
         i, D = while_state
-        D = jax.ops.index_add(D, jax.ops.index[i],
-                              D[i + 1])
+        D = jax.ops.index_add(D, jax.ops.index[i], D[i + 1])
         i -= 1
         return [i, D]
 
@@ -426,8 +436,9 @@ def _update_step_size(state, factor):
     J = jnp.arange(0, MAX_ORDER + 1)
 
     # only update order+1, order+1 entries of D
-    RU = jnp.where(jnp.logical_and(I <= order, J <= order),
-                   RU, jnp.identity(MAX_ORDER + 1))
+    RU = jnp.where(
+        jnp.logical_and(I <= order, J <= order), RU, jnp.identity(MAX_ORDER + 1)
+    )
     D = state.D
     D = jnp.dot(RU.T, D)
     # D = jax.ops.index_update(D, jax.ops.index[:order + 1],
@@ -439,9 +450,9 @@ def _update_step_size(state, factor):
     # update y0 (D has changed)
     y0, scale_y0 = _predict(state, D)
 
-    return state._replace(n_equal_steps=n_equal_steps,
-                          h=h, c=c,
-                          D=D, psi=psi, y0=y0, scale_y0=scale_y0)
+    return state._replace(
+        n_equal_steps=n_equal_steps, h=h, c=c, D=D, psi=psi, y0=y0, scale_y0=scale_y0
+    )
 
 
 def _update_jacobian(state, jac):
@@ -453,8 +464,12 @@ def _update_jacobian(state, jac):
     n_jacobian_evals = state.n_jacobian_evals + 1
     LU = jax.scipy.linalg.lu_factor(state.M - state.c * J)
     n_lu_decompositions = state.n_lu_decompositions + 1
-    return state._replace(J=J, n_jacobian_evals=n_jacobian_evals, LU=LU,
-                          n_lu_decompositions=n_lu_decompositions)
+    return state._replace(
+        J=J,
+        n_jacobian_evals=n_jacobian_evals,
+        LU=LU,
+        n_lu_decompositions=n_lu_decompositions,
+    )
 
 
 def _newton_iteration(state, fun):
@@ -485,7 +500,7 @@ def while_body(while_state):
         n_function_evals += 1
         b = c * f_eval - M @ (psi + d)
         dy = jax.scipy.linalg.lu_solve(LU, b)
-        dy_norm = jnp.sqrt(jnp.mean((dy / scale_y0)**2))
+        dy_norm = jnp.sqrt(jnp.mean((dy / scale_y0) ** 2))
         rate = dy_norm / dy_norm_old
 
         # if iteration is not going to converge in NEWTON_MAXITER
@@ -499,23 +514,22 @@ def while_body(while_state):
         y = y0 + d
 
         # if converged then break out of iteration early
-        pred = dy_norm_old >= 0.
+        pred = dy_norm_old >= 0.0
         pred *= rate / (1 - rate) * dy_norm < tol
-        converged = (dy_norm == 0.) + pred
+        converged = (dy_norm == 0.0) + pred
 
         dy_norm_old = dy_norm
 
         return [k + 1, converged, dy_norm_old, d, y, n_function_evals]
 
-    k, converged, dy_norm_old, d, y, n_function_evals = \
-        jax.lax.while_loop(while_cond,
-                           while_body,
-                           while_state)
+    k, converged, dy_norm_old, d, y, n_function_evals = jax.lax.while_loop(
+        while_cond, while_body, while_state
+    )
     return converged, k, y, d, state._replace(n_function_evals=n_function_evals)
 
 
 def rms_norm(arg):
-    return jnp.sqrt(jnp.mean(arg**2))
+    return jnp.sqrt(jnp.mean(arg ** 2))
 
 
 def _prepare_next_step(state, d):
@@ -534,21 +548,18 @@ def _prepare_next_step_order_change(state, d, y, n_iter):
     scale_y = state.atol + state.rtol * jnp.abs(y)
     error = state.error_const[order] * d
     error_norm = rms_norm(error / scale_y)
-    safety = 0.9 * (2 * NEWTON_MAXITER + 1) / (2 * NEWTON_MAXITER
-                                               + n_iter)
+    safety = 0.9 * (2 * NEWTON_MAXITER + 1) / (2 * NEWTON_MAXITER + n_iter)
 
     # similar to the optimal step size factor we calculated above for the current
     # order k, we need to calculate the optimal step size factors for orders
     # k-1 and k+1. To do this, we note that the error = C_k * D^{k+1} y_n
     error_m_norm = jnp.where(
-        order > 1,
-        rms_norm(state.error_const[order - 1] * D[order] / scale_y),
-        jnp.inf
+        order > 1, rms_norm(state.error_const[order - 1] * D[order] / scale_y), jnp.inf
     )
     error_p_norm = jnp.where(
         order < MAX_ORDER,
         rms_norm(state.error_const[order + 1] * D[order + 2] / scale_y),
-        jnp.inf
+        jnp.inf,
     )
 
     error_norms = jnp.array([error_m_norm, error_norm, error_p_norm])
@@ -566,7 +577,7 @@ def _prepare_next_step_order_change(state, d, y, n_iter):
 
 
 def _bdf_step(state, fun, jac):
-    #print('bdf_step', state.t, state.h)
+    # print('bdf_step', state.t, state.h)
     # we will try and use the old jacobian unless convergence of newton iteration
     # fails
     updated_jacobian = False
@@ -596,7 +607,7 @@ def while_body(while_state):
         state = tree_multimap(
             partial(jnp.where, not_converged * updated_jacobian),
             _update_step_size_and_lu(state, 0.3),
-            state
+            state,
         )
 
         # if not_converged * updated_jacobian:
@@ -608,11 +619,10 @@ def while_body(while_state):
         # again
         (state, updated_jacobian) = tree_multimap(
             partial(
-                jnp.where,
-                not_converged * (updated_jacobian == False)  # noqa: E712
+                jnp.where, not_converged * (updated_jacobian == False)  # noqa: E712
             ),
             (_update_jacobian(state, jac), True),
-            (state, False + updated_jacobian)
+            (state, False + updated_jacobian),
         )
 
         safety = 0.9 * (2 * NEWTON_MAXITER + 1) / (2 * NEWTON_MAXITER + n_iter)
@@ -626,26 +636,24 @@ def while_body(while_state):
         error_norm = rms_norm(error / scale_y)
 
         # calculate optimal step size factor as per eq 2.46 of [2]
-        factor = jnp.maximum(MIN_FACTOR,
-                             safety *
-                             error_norm ** (-1 / (state.order + 1)))
+        factor = jnp.maximum(
+            MIN_FACTOR, safety * error_norm ** (-1 / (state.order + 1))
+        )
 
         # if converged * (error_norm > 1):
         #    print('converged, but error is too large',error_norm, factor, d, scale_y)
 
         (state, step_accepted) = tree_multimap(
-            partial(
-                jnp.where,
-                converged * (error_norm > 1)  # noqa: E712
-            ),
+            partial(jnp.where, converged * (error_norm > 1)),  # noqa: E712
             (_update_step_size_and_lu(state, factor), False),
-            (state, converged)
+            (state, converged),
         )
 
         return [state, step_accepted, updated_jacobian, y, d, n_iter]
 
-    state, step_accepted, updated_jacobian, y, d, n_iter = \
-        jax.lax.while_loop(while_cond, while_body, while_state)
+    state, step_accepted, updated_jacobian, y, d, n_iter = jax.lax.while_loop(
+        while_cond, while_body, while_state
+    )
 
     # take the accepted step
     n_steps = state.n_steps + 1
@@ -660,7 +668,7 @@ def while_body(while_state):
     state = tree_multimap(
         partial(jnp.where, n_equal_steps < state.order + 1),
         _prepare_next_step(state, d),
-        _prepare_next_step_order_change(state, d, y, n_iter)
+        _prepare_next_step_order_change(state, d, y, n_iter),
     )
 
     return state
@@ -692,9 +700,9 @@ def while_body(while_state):
         j += 1
         return [j, time_factor, order_summation]
 
-    j, time_factor, order_summation = jax.lax.while_loop(while_cond,
-                                                         while_body,
-                                                         while_state)
+    j, time_factor, order_summation = jax.lax.while_loop(
+        while_cond, while_body, while_state
+    )
     return order_summation
 
 
@@ -708,7 +716,7 @@ def block_fun(i, j, Ai, Aj):
                     Ai.shape[0] if Ai.ndim > 1 else 1,
                     Aj.shape[1] if Aj.ndim > 1 else 1,
                 ),
-                dtype=Ai.dtype
+                dtype=Ai.dtype,
             )
 
     blocks = [
@@ -718,6 +726,7 @@ def block_fun(i, j, Ai, Aj):
 
     return jnp.block(blocks)
 
+
 # NOTE: the code below (except the docstring on jax_bdf_integrate and other minor
 # edits), has been modified from the JAX library at https://github.com/google/jax.
 # The main difference is the addition of support for semi-explicit dae index 1 problems
@@ -785,10 +794,13 @@ def jax_bdf_integrate(func, y0, t_eval, *args, rtol=1e-6, atol=1e-6, mass=None):
            fundamental algorithms for scientific computing in Python.
            Nature methods, 17(3), 261-272.
     """
+
     def _check_arg(arg):
         if not isinstance(arg, core.Tracer) and not core.valid_jaxtype(arg):
-            msg = ("The contents of odeint *args must be arrays or scalars, but got "
-                   "\n{}.")
+            msg = (
+                "The contents of odeint *args must be arrays or scalars, but got "
+                "\n{}."
+            )
         raise TypeError(msg.format(arg))
 
     flat_args, in_tree = tree_flatten((y0, t_eval[0], *args))
@@ -871,7 +883,7 @@ def aug_dynamics(augmented_state, t, *args):
 
         return (-y_dot, y_bar_dot, *rest)
 
-    algebraic_variables = jnp.diag(mass) == 0.
+    algebraic_variables = jnp.diag(mass) == 0.0
     differentiable_variables = algebraic_variables == False  # noqa: E712
     mass_is_I = (mass == jnp.eye(mass.shape[0])).all()
     is_dae = jnp.any(algebraic_variables)
@@ -894,9 +906,9 @@ def initialise(g0, y0, t0):
             g0_a = g0[algebraic_variables]
             invJ_aa = jax.scipy.linalg.lu_solve(LU, g0_a)
             y_bar = jax.ops.index_update(
-                g0, differentiable_variables,
-                jax.scipy.linalg.lu_solve(LU_invM_dd,
-                                          g0_a - J_ad @ invJ_aa)
+                g0,
+                differentiable_variables,
+                jax.scipy.linalg.lu_solve(LU_invM_dd, g0_a - J_ad @ invJ_aa),
             )
         else:
             y_bar = jax.scipy.linalg.lu_solve(LU_invM_dd, g0)
@@ -904,12 +916,12 @@ def initialise(g0, y0, t0):
 
     y_bar = initialise(g[-1], ys[-1], ts[-1])
     ts_bar = []
-    t0_bar = 0.
+    t0_bar = 0.0
 
     def arg_to_identity(arg):
         return onp.identity(arg.shape[0] if arg.ndim > 0 else 1, dtype=arg.dtype)
 
-    aug_mass = (mass, mass, jnp.array(1.), tree_map(arg_to_identity, args))
+    aug_mass = (mass, mass, jnp.array(1.0), tree_map(arg_to_identity, args))
 
     def scan_fun(carry, i):
         y_bar, t0_bar, args_bar = carry
@@ -918,10 +930,14 @@ def scan_fun(carry, i):
         t0_bar = t0_bar - t_bar
         # Run augmented system backwards to previous observation
         _, y_bar, t0_bar, args_bar = jax_bdf_integrate(
-            aug_dynamics, (ys[i], y_bar, t0_bar, args_bar),
+            aug_dynamics,
+            (ys[i], y_bar, t0_bar, args_bar),
             jnp.array([-ts[i], -ts[i - 1]]),
-            *args, mass=aug_mass,
-            rtol=rtol, atol=atol)
+            *args,
+            mass=aug_mass,
+            rtol=rtol,
+            atol=atol
+        )
         y_bar, t0_bar, args_bar = tree_map(op.itemgetter(1), (y_bar, t0_bar, args_bar))
         # Add gradient from current output
         y_bar = y_bar + initialise(g[i - 1], ys[i - 1], ts[i - 1])
@@ -929,7 +945,8 @@ def scan_fun(carry, i):
 
     init_carry = (y_bar, t0_bar, tree_map(jnp.zeros_like, args))
     (y_bar, t0_bar, args_bar), rev_ts_bar = jax.lax.scan(
-        scan_fun, init_carry, jnp.arange(len(ts) - 1, 0, -1))
+        scan_fun, init_carry, jnp.arange(len(ts) - 1, 0, -1)
+    )
     ts_bar = jnp.concatenate([jnp.array([t0_bar]), rev_ts_bar[::-1]])
     return (y_bar, ts_bar, *args_bar)
 
@@ -951,8 +968,7 @@ def closure_convert(fun, in_tree, in_avals):
     # differentiating. As a proxy for that, we hoist consts with float dtype.
     # TODO(mattjj): revise this approach
     (closure_consts, hoisted_consts), merge = partition_list(
-        lambda c: dtypes.issubdtype(dtypes.dtype(c), jnp.inexact),
-        consts
+        lambda c: dtypes.issubdtype(dtypes.dtype(c), jnp.inexact), consts
     )
     num_consts = len(hoisted_consts)
 
@@ -973,6 +989,7 @@ def partition_list(choice, lst):
     def merge(l1, l2):
         i1, i2 = iter(l1), iter(l2)
         return [next(i2 if snd else i1) for snd in which]
+
     return out, merge
 
 
diff --git a/pybamm/solvers/jax_solver.py b/pybamm/solvers/jax_solver.py
index db912d0da1..77595b0a7a 100644
--- a/pybamm/solvers/jax_solver.py
+++ b/pybamm/solvers/jax_solver.py
@@ -45,16 +45,17 @@ class JaxSolver(pybamm.BaseSolver):
         for details.
     """
 
-    def __init__(self, method='RK45', root_method=None,
-                 rtol=1e-6, atol=1e-6, extra_options=None):
+    def __init__(
+        self, method="RK45", root_method=None, rtol=1e-6, atol=1e-6, extra_options=None
+    ):
         # note: bdf solver itself calculates consistent initial conditions so can set
         # root_method to none, allow user to override this behavior
         super().__init__(method, rtol, atol, root_method=root_method)
-        method_options = ['RK45', 'BDF']
+        method_options = ["RK45", "BDF"]
         if method not in method_options:
-            raise ValueError('method must be one of {}'.format(method_options))
+            raise ValueError("method must be one of {}".format(method_options))
         self.ode_solver = False
-        if method == 'RK45':
+        if method == "RK45":
             self.ode_solver = True
         self.extra_options = extra_options or {}
         self.name = "JAX solver ({})".format(method)
@@ -80,8 +81,9 @@ def get_solve(self, model, t_eval):
         """
         if model not in self._cached_solves:
             if model not in self.models_set_up:
-                raise RuntimeError("Model is not set up for solving, run"
-                                   "`solver.solve(model)` first")
+                raise RuntimeError(
+                    "Model is not set up for solving, run" "`solver.solve(model)` first"
+                )
 
             self._cached_solves[model] = self.create_solve(model, t_eval)
 
@@ -106,32 +108,35 @@ def create_solve(self, model, t_eval):
 
         """
         if model.convert_to_format != "jax":
-            raise RuntimeError("Model must be converted to JAX to use this solver"
-                               " (i.e. `model.convert_to_format = 'jax')")
+            raise RuntimeError(
+                "Model must be converted to JAX to use this solver"
+                " (i.e. `model.convert_to_format = 'jax')"
+            )
 
         if model.terminate_events_eval:
-            raise RuntimeError("Terminate events not supported for this solver."
-                               " Model has the following events:"
-                               " {}.\nYou can remove events using `model.events = []`."
-                               " It might be useful to first solve the model using a"
-                               " different solver to obtain the time of the event, then"
-                               " re-solve using no events and a fixed"
-                               " end-time".format(model.events))
+            raise RuntimeError(
+                "Terminate events not supported for this solver."
+                " Model has the following events:"
+                " {}.\nYou can remove events using `model.events = []`."
+                " It might be useful to first solve the model using a"
+                " different solver to obtain the time of the event, then"
+                " re-solve using no events and a fixed"
+                " end-time".format(model.events)
+            )
 
         # Initial conditions, make sure they are an 0D array
         y0 = jnp.array(model.y0).reshape(-1)
         mass = None
-        if self.method == 'BDF':
+        if self.method == "BDF":
             mass = model.mass_matrix.entries.toarray()
 
         def rhs_ode(y, t, inputs):
-            return model.rhs_eval(t, y, inputs),
+            return (model.rhs_eval(t, y, inputs),)
 
         def rhs_dae(y, t, inputs):
-            return jnp.concatenate([
-                model.rhs_eval(t, y, inputs),
-                model.algebraic_eval(t, y, inputs),
-            ])
+            return jnp.concatenate(
+                [model.rhs_eval(t, y, inputs), model.algebraic_eval(t, y, inputs)]
+            )
 
         def solve_model_rk45(inputs):
             y = odeint(
@@ -158,7 +163,7 @@ def solve_model_bdf(inputs):
             )
             return jnp.transpose(y)
 
-        if self.method == 'RK45':
+        if self.method == "RK45":
             return jax.jit(solve_model_rk45)
         else:
             return jax.jit(solve_model_bdf)
@@ -194,5 +199,4 @@ def _integrate(self, model, t_eval, inputs=None):
         termination = "final time"
         t_event = None
         y_event = onp.array(None)
-        return pybamm.Solution(t_eval, y,
-                               t_event, y_event, termination)
+        return pybamm.Solution(t_eval, y, t_event, y_event, termination)
diff --git a/scripts/install_KLU_Sundials.py b/scripts/install_KLU_Sundials.py
index 539e41a268..92f3f27e5e 100755
--- a/scripts/install_KLU_Sundials.py
+++ b/scripts/install_KLU_Sundials.py
@@ -121,5 +121,3 @@ def download_extract_library(url, download_dir):
 print("-" * 10, "Building the sundials", "-" * 40)
 make_cmd = ["make", "install"]
 subprocess.run(make_cmd, cwd=build_dir)
-
-
diff --git a/tests/integration/test_models/standard_model_tests.py b/tests/integration/test_models/standard_model_tests.py
index f1176493a5..48267a8bbf 100644
--- a/tests/integration/test_models/standard_model_tests.py
+++ b/tests/integration/test_models/standard_model_tests.py
@@ -128,8 +128,9 @@ def __init__(self, model, parameter_values=None, disc=None):
 
         self.model = model
 
-    def evaluate_model(self, simplify=False, use_known_evals=False,
-                       to_python=False, to_jax=False):
+    def evaluate_model(
+        self, simplify=False, use_known_evals=False, to_python=False, to_jax=False
+    ):
         result = np.empty((0, 1))
         for eqn in [self.model.concatenated_rhs, self.model.concatenated_algebraic]:
             if simplify:
diff --git a/tests/integration/test_spatial_methods/test_spectral_volume.py b/tests/integration/test_spatial_methods/test_spectral_volume.py
index 699e5f3afe..09b10d6426 100644
--- a/tests/integration/test_spatial_methods/test_spectral_volume.py
+++ b/tests/integration/test_spatial_methods/test_spectral_volume.py
@@ -8,8 +8,7 @@
 
 
 def get_mesh_for_testing(
-    xpts=None, rpts=10, ypts=15, zpts=15, geometry=None, cc_submesh=None,
-    order=2
+    xpts=None, rpts=10, ypts=15, zpts=15, geometry=None, cc_submesh=None, order=2
 ):
     param = pybamm.ParameterValues(
         values={
@@ -33,22 +32,19 @@ def get_mesh_for_testing(
 
     submesh_types = {
         "negative electrode": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
+            pybamm.SpectralVolume1DSubMesh, {"order": order}
+        ),
+        "separator": pybamm.MeshGenerator(
+            pybamm.SpectralVolume1DSubMesh, {"order": order}
         ),
-        "separator": pybamm.MeshGenerator(pybamm.SpectralVolume1DSubMesh,
-                                          {"order": order}),
         "positive electrode": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
+            pybamm.SpectralVolume1DSubMesh, {"order": order}
         ),
         "negative particle": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
+            pybamm.SpectralVolume1DSubMesh, {"order": order}
         ),
         "positive particle": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
+            pybamm.SpectralVolume1DSubMesh, {"order": order}
         ),
         "current collector": pybamm.MeshGenerator(pybamm.SubMesh0D),
     }
diff --git a/tests/unit/test_meshes/test_one_dimensional_submesh.py b/tests/unit/test_meshes/test_one_dimensional_submesh.py
index 8cb084e786..59d8013419 100644
--- a/tests/unit/test_meshes/test_one_dimensional_submesh.py
+++ b/tests/unit/test_meshes/test_one_dimensional_submesh.py
@@ -277,8 +277,7 @@ class TestSpectralVolume1DSubMesh(unittest.TestCase):
     def test_exceptions(self):
         edges = np.array([0, 0.3, 1])
         submesh_params = {"edges": edges}
-        mesh = pybamm.MeshGenerator(pybamm.SpectralVolume1DSubMesh,
-                                    submesh_params)
+        mesh = pybamm.MeshGenerator(pybamm.SpectralVolume1DSubMesh, submesh_params)
 
         x_n = pybamm.standard_spatial_vars.x_n
 
@@ -306,8 +305,7 @@ def test_mesh_creation_no_parameters(self):
         )
 
         geometry = {
-            "negative particle": {r: {"min": pybamm.Scalar(0),
-                                      "max": pybamm.Scalar(1)}}
+            "negative particle": {r: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}
         }
 
         edges = np.array([0, 0.3, 1])
@@ -329,16 +327,17 @@ def test_mesh_creation_no_parameters(self):
 
         # check number of edges and nodes
         self.assertEqual(len(mesh["negative particle"].sv_nodes), var_pts[r])
-        self.assertEqual(len(mesh["negative particle"].nodes),
-                         order * var_pts[r])
+        self.assertEqual(len(mesh["negative particle"].nodes), order * var_pts[r])
         self.assertEqual(
             len(mesh["negative particle"].edges),
             len(mesh["negative particle"].nodes) + 1,
         )
 
         # check Chebyshev subdivision locations
-        for (a, b) in zip(mesh["negative particle"].edges.tolist(),
-                          [0, 0.075, 0.225, 0.3, 0.475, 0.825, 1]):
+        for (a, b) in zip(
+            mesh["negative particle"].edges.tolist(),
+            [0, 0.075, 0.225, 0.3, 0.475, 0.825, 1],
+        ):
             self.assertAlmostEqual(a, b)
 
         # test uniform submesh creation
@@ -352,8 +351,10 @@ def test_mesh_creation_no_parameters(self):
 
         # create mesh
         mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
-        for (a, b) in zip(mesh["negative particle"].edges.tolist(),
-                          [0.0, 0.125, 0.375, 0.5, 0.625, 0.875, 1.0]):
+        for (a, b) in zip(
+            mesh["negative particle"].edges.tolist(),
+            [0.0, 0.125, 0.375, 0.5, 0.625, 0.875, 1.0],
+        ):
             self.assertAlmostEqual(a, b)
 
 
diff --git a/tests/unit/test_parameters/test_parameter_sets/test_LFP_Prada2013.py b/tests/unit/test_parameters/test_parameter_sets/test_LFP_Prada2013.py
new file mode 100644
index 0000000000..e8af896605
--- /dev/null
+++ b/tests/unit/test_parameters/test_parameter_sets/test_LFP_Prada2013.py
@@ -0,0 +1,42 @@
+#
+# Tests for LFP Prada 2013 parameter set
+#
+import pybamm
+import unittest
+
+
+class TestLFPPrada2013(unittest.TestCase):
+    def test_load_params(self):
+        cathode = pybamm.ParameterValues({}).read_parameters_csv(
+            pybamm.get_parameters_filepath(
+                "input/parameters/lithium-ion/cathodes/LFP_Prada2013/parameters.csv"
+            )
+        )
+        self.assertEqual(cathode["Positive electrode porosity"], "0.12728395")
+
+        cell = pybamm.ParameterValues({}).read_parameters_csv(
+            pybamm.get_parameters_filepath(
+                "input/parameters/lithium-ion/cells/A123_Lain2019/parameters.csv"
+            )
+        )
+        self.assertAlmostEqual(cell["Negative current collector thickness [m]"], 1e-5)
+
+    def test_standard_lithium_parameters(self):
+
+        chemistry = pybamm.parameter_sets.Prada2013
+        parameter_values = pybamm.ParameterValues(chemistry=chemistry)
+
+        model = pybamm.lithium_ion.DFN()
+        sim = pybamm.Simulation(model, parameter_values=parameter_values)
+        sim.set_parameters()
+        sim.build()
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()
diff --git a/tests/unit/test_solvers/test_jax_bdf_solver.py b/tests/unit/test_solvers/test_jax_bdf_solver.py
index 1ec810c09d..d336956c2f 100644
--- a/tests/unit/test_solvers/test_jax_bdf_solver.py
+++ b/tests/unit/test_solvers/test_jax_bdf_solver.py
@@ -5,6 +5,7 @@
 import time
 import numpy as np
 from platform import system
+
 if system() != "Windows":
     import jax
 
@@ -40,8 +41,7 @@ def fun(y, t):
         t1 = time.perf_counter() - t0
 
         # test accuracy
-        np.testing.assert_allclose(y[:, 0], np.exp(0.1 * t_eval),
-                                   rtol=1e-6, atol=1e-6)
+        np.testing.assert_allclose(y[:, 0], np.exp(0.1 * t_eval), rtol=1e-6, atol=1e-6)
 
         t0 = time.perf_counter()
         y = pybamm.jax_bdf_integrate(fun, y0, t_eval, rtol=1e-8, atol=1e-8)
@@ -51,23 +51,21 @@ def fun(y, t):
         self.assertLess(t2, t1)
 
         # test second run is accurate
-        np.testing.assert_allclose(y[:, 0], np.exp(0.1 * t_eval),
-                                   rtol=1e-6, atol=1e-6)
+        np.testing.assert_allclose(y[:, 0], np.exp(0.1 * t_eval), rtol=1e-6, atol=1e-6)
 
     def test_mass_matrix(self):
         # Solve
         t_eval = np.linspace(0.0, 1.0, 80)
 
         def fun(y, t):
-            return jax.numpy.stack([
-                0.1 * y[0],
-                y[1] - 2.0 * y[0],
-            ])
+            return jax.numpy.stack([0.1 * y[0], y[1] - 2.0 * y[0]])
 
-        mass = jax.numpy.array([
-            [2.0, 0.0],
-            [0.0, 0.0],
-        ])
+        mass = jax.numpy.array(
+            [
+                [2.0, 0.0],
+                [0.0, 0.0],
+            ]
+        )
 
         # give some bad initial conditions, solver should calculate correct ones using
         # this as a guess
@@ -79,10 +77,8 @@ def fun(y, t):
 
         # test accuracy
         soln = np.exp(0.05 * t_eval)
-        np.testing.assert_allclose(y[:, 0], soln,
-                                   rtol=1e-7, atol=1e-7)
-        np.testing.assert_allclose(y[:, 1], 2.0 * soln,
-                                   rtol=1e-7, atol=1e-7)
+        np.testing.assert_allclose(y[:, 0], soln, rtol=1e-7, atol=1e-7)
+        np.testing.assert_allclose(y[:, 1], 2.0 * soln, rtol=1e-7, atol=1e-7)
 
         t0 = time.perf_counter()
         y = pybamm.jax_bdf_integrate(fun, y0, t_eval, mass=mass, rtol=1e-8, atol=1e-8)
@@ -92,8 +88,7 @@ def fun(y, t):
         self.assertLess(t2, t1)
 
         # test second run is accurate
-        np.testing.assert_allclose(y[:, 0], np.exp(0.05 * t_eval),
-                                   rtol=1e-7, atol=1e-7)
+        np.testing.assert_allclose(y[:, 0], np.exp(0.05 * t_eval), rtol=1e-7, atol=1e-7)
 
     def test_solver_sensitivities(self):
         # Create model
@@ -125,9 +120,9 @@ def fun(y, t, inputs):
         @jax.jit
         def solve_bdf(rate):
             return jax.numpy.sum(
-                pybamm.jax_bdf_integrate(fun, y0, t_eval,
-                                         {'rate': rate},
-                                         rtol=1e-9, atol=1e-9)
+                pybamm.jax_bdf_integrate(
+                    fun, y0, t_eval, {"rate": rate}, rtol=1e-9, atol=1e-9
+                )
             )
 
         # check answers with finite difference
@@ -145,15 +140,19 @@ def test_mass_matrix_with_sensitivities(self):
         t_eval = np.linspace(0.0, 1.0, 80)
 
         def fun(y, t, inputs):
-            return jax.numpy.stack([
-                inputs['rate'] * y[0],
-                y[1] - 2.0 * y[0],
-            ])
+            return jax.numpy.stack(
+                [
+                    inputs["rate"] * y[0],
+                    y[1] - 2.0 * y[0],
+                ]
+            )
 
-        mass = jax.numpy.array([
-            [2.0, 0.0],
-            [0.0, 0.0],
-        ])
+        mass = jax.numpy.array(
+            [
+                [2.0, 0.0],
+                [0.0, 0.0],
+            ]
+        )
 
         y0 = jax.numpy.array([1.0, 2.0])
 
@@ -164,10 +163,9 @@ def fun(y, t, inputs):
         @jax.jit
         def solve_bdf(rate):
             return jax.numpy.sum(
-                pybamm.jax_bdf_integrate(fun, y0, t_eval,
-                                         {'rate': rate},
-                                         mass=mass,
-                                         rtol=1e-9, atol=1e-9)
+                pybamm.jax_bdf_integrate(
+                    fun, y0, t_eval, {"rate": rate}, mass=mass, rtol=1e-9, atol=1e-9
+                )
             )
 
         # check answers with finite difference
@@ -203,8 +201,9 @@ def test_solver_with_inputs(self):
         def fun(y, t, inputs):
             return rhs.evaluate(t=t, y=y, inputs=inputs).reshape(-1)
 
-        y = pybamm.jax_bdf_integrate(fun, y0, t_eval, {
-            "rate": 0.1}, rtol=1e-9, atol=1e-9)
+        y = pybamm.jax_bdf_integrate(
+            fun, y0, t_eval, {"rate": 0.1}, rtol=1e-9, atol=1e-9
+        )
 
         np.testing.assert_allclose(y[:, 0].reshape(-1), np.exp(-0.1 * t_eval))
 
diff --git a/tests/unit/test_solvers/test_jax_solver.py b/tests/unit/test_solvers/test_jax_solver.py
index d30f3f76b6..3c6f727583 100644
--- a/tests/unit/test_solvers/test_jax_solver.py
+++ b/tests/unit/test_solvers/test_jax_solver.py
@@ -5,6 +5,7 @@
 import time
 import numpy as np
 from platform import system
+
 if system() != "Windows":
     import jax
 
@@ -27,18 +28,17 @@ def test_model_solver(self):
         disc = pybamm.Discretisation(mesh, spatial_methods)
         disc.process_model(model)
 
-        for method in ['RK45', 'BDF']:
+        for method in ["RK45", "BDF"]:
             # Solve
-            solver = pybamm.JaxSolver(
-                method=method, rtol=1e-8, atol=1e-8
-            )
+            solver = pybamm.JaxSolver(method=method, rtol=1e-8, atol=1e-8)
             t_eval = np.linspace(0, 1, 80)
             t0 = time.perf_counter()
             solution = solver.solve(model, t_eval)
             t_first_solve = time.perf_counter() - t0
             np.testing.assert_array_equal(solution.t, t_eval)
-            np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t),
-                                       rtol=1e-6, atol=1e-6)
+            np.testing.assert_allclose(
+                solution.y[0], np.exp(0.1 * solution.t), rtol=1e-6, atol=1e-6
+            )
 
             # Test time
             self.assertEqual(
@@ -73,19 +73,15 @@ def test_semi_explicit_model(self):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.JaxSolver(
-            method='BDF', rtol=1e-8, atol=1e-8
-        )
+        solver = pybamm.JaxSolver(method="BDF", rtol=1e-8, atol=1e-8)
         t_eval = np.linspace(0, 1, 80)
         t0 = time.perf_counter()
         solution = solver.solve(model, t_eval)
         t_first_solve = time.perf_counter() - t0
         np.testing.assert_array_equal(solution.t, t_eval)
         soln = np.exp(0.1 * solution.t)
-        np.testing.assert_allclose(solution.y[0], soln,
-                                   rtol=1e-7, atol=1e-7)
-        np.testing.assert_allclose(solution.y[-1], 2 * soln,
-                                   rtol=1e-7, atol=1e-7)
+        np.testing.assert_allclose(solution.y[0], soln, rtol=1e-7, atol=1e-7)
+        np.testing.assert_allclose(solution.y[-1], 2 * soln, rtol=1e-7, atol=1e-7)
 
         # Test time
         self.assertEqual(
@@ -116,23 +112,21 @@ def test_solver_sensitivities(self):
         disc = pybamm.Discretisation(mesh, spatial_methods)
         disc.process_model(model)
 
-        for method in ['RK45', 'BDF']:
+        for method in ["RK45", "BDF"]:
             # Solve
-            solver = pybamm.JaxSolver(
-                method=method, rtol=1e-8, atol=1e-8
-            )
+            solver = pybamm.JaxSolver(method=method, rtol=1e-8, atol=1e-8)
             t_eval = np.linspace(0, 1, 80)
 
             h = 0.0001
             rate = 0.1
 
             # need to solve the model once to get it set up by the base solver
-            solver.solve(model, t_eval, inputs={'rate': rate})
+            solver.solve(model, t_eval, inputs={"rate": rate})
             solve = solver.get_solve(model, t_eval)
 
             # create a dummy "model" where we calculate the sum of the time series
             def solve_model(rate):
-                return jax.numpy.sum(solve({'rate': rate}))
+                return jax.numpy.sum(solve({"rate": rate}))
 
             # check answers with finite difference
             eval_plus = solve_model(rate + h)
@@ -216,15 +210,17 @@ def test_model_solver_with_inputs(self):
         solution = solver.solve(model, t_eval, inputs={"rate": 0.1})
         t_first_solve = time.perf_counter() - t0
 
-        np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t),
-                                   rtol=1e-6, atol=1e-6)
+        np.testing.assert_allclose(
+            solution.y[0], np.exp(-0.1 * solution.t), rtol=1e-6, atol=1e-6
+        )
 
         t0 = time.perf_counter()
         solution = solver.solve(model, t_eval, inputs={"rate": 0.2})
         t_second_solve = time.perf_counter() - t0
 
-        np.testing.assert_allclose(solution.y[0], np.exp(-0.2 * solution.t),
-                                   rtol=1e-6, atol=1e-6)
+        np.testing.assert_allclose(
+            solution.y[0], np.exp(-0.2 * solution.t), rtol=1e-6, atol=1e-6
+        )
 
         self.assertLess(t_second_solve, t_first_solve)
 
@@ -247,7 +243,7 @@ def test_get_solve(self):
 
         # test that another method string gives error
         with self.assertRaises(ValueError):
-            solver = pybamm.JaxSolver(method='not_real')
+            solver = pybamm.JaxSolver(method="not_real")
 
         # Solve
         solver = pybamm.JaxSolver(rtol=1e-8, atol=1e-8)
@@ -260,13 +256,11 @@ def test_get_solve(self):
         solver = solver.get_solve(model, t_eval)
         y = solver({"rate": 0.1})
 
-        np.testing.assert_allclose(y[0], np.exp(-0.1 * t_eval),
-                                   rtol=1e-6, atol=1e-6)
+        np.testing.assert_allclose(y[0], np.exp(-0.1 * t_eval), rtol=1e-6, atol=1e-6)
 
         y = solver({"rate": 0.2})
 
-        np.testing.assert_allclose(y[0], np.exp(-0.2 * t_eval),
-                                   rtol=1e-6, atol=1e-6)
+        np.testing.assert_allclose(y[0], np.exp(-0.2 * t_eval), rtol=1e-6, atol=1e-6)
 
 
 if __name__ == "__main__":
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index bc8d5a22c1..58eb86fe73 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -353,10 +353,7 @@ def test_model_solver_dae_multiple_nonsmooth_python(self):
         ]
         for discontinuity in discontinuities:
             model.events.append(
-                pybamm.Event(
-                    "nonsmooth rate",
-                    pybamm.Scalar(discontinuity),
-                )
+                pybamm.Event("nonsmooth rate", pybamm.Scalar(discontinuity),)
             )
         disc = get_discretisation_for_testing()
         disc.process_model(model)
@@ -760,10 +757,7 @@ def test_model_step_nonsmooth_events(self):
         ]
         for discontinuity in discontinuities:
             model.events.append(
-                pybamm.Event(
-                    "nonsmooth rate",
-                    pybamm.Scalar(discontinuity),
-                )
+                pybamm.Event("nonsmooth rate", pybamm.Scalar(discontinuity),)
             )
         disc = get_discretisation_for_testing()
         disc.process_model(model)
@@ -781,12 +775,8 @@ def test_model_step_nonsmooth_events(self):
         np.testing.assert_array_less(step_solution.y[-1], 1.2)
         var1_soln = (step_solution.t % a) ** 2 / 2 + a ** 2 / 2 * (step_solution.t // a)
         var2_soln = 2 * var1_soln
-        np.testing.assert_array_almost_equal(
-            step_solution.y[0], var1_soln, decimal=5
-        )
-        np.testing.assert_array_almost_equal(
-            step_solution.y[-1], var2_soln, decimal=5
-        )
+        np.testing.assert_array_almost_equal(step_solution.y[0], var1_soln, decimal=5)
+        np.testing.assert_array_almost_equal(step_solution.y[-1], var2_soln, decimal=5)
 
     def test_model_solver_dae_nonsmooth(self):
         whole_cell = ["negative electrode", "separator", "positive electrode"]
diff --git a/tests/unit/test_spatial_methods/test_spectral_volume.py b/tests/unit/test_spatial_methods/test_spectral_volume.py
index 425bc0fabc..0e9d89920c 100644
--- a/tests/unit/test_spatial_methods/test_spectral_volume.py
+++ b/tests/unit/test_spatial_methods/test_spectral_volume.py
@@ -8,8 +8,7 @@
 
 
 def get_mesh_for_testing(
-    xpts=None, rpts=10, ypts=15, zpts=15, geometry=None, cc_submesh=None,
-    order=2
+    xpts=None, rpts=10, ypts=15, zpts=15, geometry=None, cc_submesh=None, order=2
 ):
     param = pybamm.ParameterValues(
         values={
@@ -33,22 +32,19 @@ def get_mesh_for_testing(
 
     submesh_types = {
         "negative electrode": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
+            pybamm.SpectralVolume1DSubMesh, {"order": order}
+        ),
+        "separator": pybamm.MeshGenerator(
+            pybamm.SpectralVolume1DSubMesh, {"order": order}
         ),
-        "separator": pybamm.MeshGenerator(pybamm.SpectralVolume1DSubMesh,
-                                          {"order": order}),
         "positive electrode": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
+            pybamm.SpectralVolume1DSubMesh, {"order": order}
         ),
         "negative particle": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
+            pybamm.SpectralVolume1DSubMesh, {"order": order}
         ),
         "positive particle": pybamm.MeshGenerator(
-            pybamm.SpectralVolume1DSubMesh,
-            {"order": order}
+            pybamm.SpectralVolume1DSubMesh, {"order": order}
         ),
         "current collector": pybamm.MeshGenerator(pybamm.SubMesh0D),
     }