Repository: pybamm-team/PyBaMM (python)

Issue:
Implement DFN as ODEs
**Summary**
Implement the ODE version of the DFN using the capacitance trick.

 **Motivation**
This may offer a faster solution method than solving the DAEs. We can compare the two formulations and determine which is best


Pull request description:
Issue 247 dfn as odes
# Description

Implement ODE formulation of the DFN

Fixes #247 

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/master/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [x] New feature (non-breaking change which adds functionality)
- [ ] Optimization (back-end change that speeds up the code)
- [ ] Bug fix (non-breaking change which fixes an issue)


# Key checklist:

- [ ] No style issues: `$ flake8`
- [ ] All tests pass: `$ python run-tests.py --unit`
- [ ] 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:

- [ ] Code is commented, particularly in hard-to-understand areas
- [ ] Tests added that prove fix is effective or that feature works


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6f2e22249..29a5f3a55 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Features
 
 -   Changed rootfinding algorithm to CasADi, scipy.optimize.root still accessible as an option ([#844](https://github.com/pybamm-team/PyBaMM/pull/844))
+-   Added capacitance effects to lithium-ion models ([#842](https://github.com/pybamm-team/PyBaMM/pull/842))
 -   Added NCA parameter set ([#824](https://github.com/pybamm-team/PyBaMM/pull/824))
 -   Added functionality to `Solution` that automatically gets `t_eval` from the data when simulating drive cycles and performs checks to ensure the output has the required resolution to accurately capture the input current ([#819](https://github.com/pybamm-team/PyBaMM/pull/819))
 -   Added options to export a solution to matlab or csv ([#811](https://github.com/pybamm-team/PyBaMM/pull/811))
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 035e7cb0d..efda022d1 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -224,10 +224,10 @@ class BaseBatteryModel(pybamm.BaseModel):
 
         # Options that are incompatible with models
         if isinstance(self, pybamm.lithium_ion.BaseModel):
-            if options["surface form"] is not False:
-                raise pybamm.OptionError(
-                    "surface form not implemented for lithium-ion models"
-                )
+            # if options["surface form"] is not False:
+            #     raise pybamm.OptionError(
+            #         "surface form not implemented for lithium-ion models"
+            #     )
             if options["convection"] is True:
                 raise pybamm.OptionError(
                     "convection not implemented for lithium-ion models"
diff --git a/pybamm/models/full_battery_models/lithium_ion/dfn.py b/pybamm/models/full_battery_models/lithium_ion/dfn.py
index 0c4bc7c60..582629c92 100644
--- a/pybamm/models/full_battery_models/lithium_ion/dfn.py
+++ b/pybamm/models/full_battery_models/lithium_ion/dfn.py
@@ -84,24 +84,40 @@ class DFN(BaseModel):
 
     def set_solid_submodel(self):
 
-        self.submodels["negative electrode"] = pybamm.electrode.ohm.Full(
-            self.param, "Negative", self.reactions
-        )
-        self.submodels["positive electrode"] = pybamm.electrode.ohm.Full(
-            self.param, "Positive", self.reactions
-        )
+        if self.options["surface form"] is False:
+            submod_n = pybamm.electrode.ohm.Full(self.param, "Negative", self.reactions)
+            submod_p = pybamm.electrode.ohm.Full(self.param, "Positive", self.reactions)
+        else:
+            submod_n = pybamm.electrode.ohm.SurfaceForm(self.param, "Negative")
+            submod_p = pybamm.electrode.ohm.SurfaceForm(self.param, "Positive")
+
+        self.submodels["negative electrode"] = submod_n
+        self.submodels["positive electrode"] = submod_p
 
     def set_electrolyte_submodel(self):
 
         electrolyte = pybamm.electrolyte.stefan_maxwell
+        surf_form = electrolyte.conductivity.surface_potential_form
 
-        self.submodels["electrolyte conductivity"] = electrolyte.conductivity.Full(
-            self.param, self.reactions
-        )
         self.submodels["electrolyte diffusion"] = electrolyte.diffusion.Full(
             self.param, self.reactions
         )
 
+        if self.options["surface form"] is False:
+            self.submodels["electrolyte conductivity"] = electrolyte.conductivity.Full(
+                self.param, self.reactions
+            )
+        elif self.options["surface form"] == "differential":
+            for domain in ["Negative", "Separator", "Positive"]:
+                self.submodels[
+                    domain.lower() + " electrolyte conductivity"
+                ] = surf_form.FullDifferential(self.param, domain, self.reactions)
+        elif self.options["surface form"] == "algebraic":
+            for domain in ["Negative", "Separator", "Positive"]:
+                self.submodels[
+                    domain.lower() + " electrolyte conductivity"
+                ] = surf_form.FullAlgebraic(self.param, domain, self.reactions)
+
     @property
     def default_geometry(self):
         dimensionality = self.options["dimensionality"]
diff --git a/pybamm/models/full_battery_models/lithium_ion/spm.py b/pybamm/models/full_battery_models/lithium_ion/spm.py
index 97e55a278..023b64bc7 100644
--- a/pybamm/models/full_battery_models/lithium_ion/spm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/spm.py
@@ -32,6 +32,7 @@ class SPM(BaseModel):
     def __init__(self, options=None, name="Single Particle Model", build=True):
         super().__init__(options, name)
 
+        self.set_reactions()
         self.set_external_circuit_submodel()
         self.set_porosity_submodel()
         self.set_tortuosity_submodels()
@@ -57,12 +58,21 @@ class SPM(BaseModel):
 
     def set_interfacial_submodel(self):
 
-        self.submodels[
-            "negative interface"
-        ] = pybamm.interface.lithium_ion.InverseButlerVolmer(self.param, "Negative")
-        self.submodels[
-            "positive interface"
-        ] = pybamm.interface.lithium_ion.InverseButlerVolmer(self.param, "Positive")
+        if self.options["surface form"] is False:
+            self.submodels[
+                "negative interface"
+            ] = pybamm.interface.lithium_ion.InverseButlerVolmer(self.param, "Negative")
+            self.submodels[
+                "positive interface"
+            ] = pybamm.interface.lithium_ion.InverseButlerVolmer(self.param, "Positive")
+        else:
+            self.submodels[
+                "negative interface"
+            ] = pybamm.interface.lithium_ion.ButlerVolmer(self.param, "Negative")
+
+            self.submodels[
+                "positive interface"
+            ] = pybamm.interface.lithium_ion.ButlerVolmer(self.param, "Positive")
 
     def set_particle_submodel(self):
 
@@ -96,10 +106,26 @@ class SPM(BaseModel):
     def set_electrolyte_submodel(self):
 
         electrolyte = pybamm.electrolyte.stefan_maxwell
+        surf_form = electrolyte.conductivity.surface_potential_form
 
-        self.submodels[
-            "electrolyte conductivity"
-        ] = electrolyte.conductivity.LeadingOrder(self.param)
+        if self.options["surface form"] is False:
+            self.submodels[
+                "leading-order electrolyte conductivity"
+            ] = electrolyte.conductivity.LeadingOrder(self.param)
+
+        elif self.options["surface form"] == "differential":
+            for domain in ["Negative", "Separator", "Positive"]:
+                self.submodels[
+                    "leading-order " + domain.lower() + " electrolyte conductivity"
+                ] = surf_form.LeadingOrderDifferential(
+                    self.param, domain, self.reactions
+                )
+
+        elif self.options["surface form"] == "algebraic":
+            for domain in ["Negative", "Separator", "Positive"]:
+                self.submodels[
+                    "leading-order " + domain.lower() + " electrolyte conductivity"
+                ] = surf_form.LeadingOrderAlgebraic(self.param, domain, self.reactions)
         self.submodels[
             "electrolyte diffusion"
         ] = electrolyte.diffusion.ConstantConcentration(self.param)
diff --git a/pybamm/models/submodels/interface/lithium_ion.py b/pybamm/models/submodels/interface/lithium_ion.py
index e995a9eb7..f4cc6b9d8 100644
--- a/pybamm/models/submodels/interface/lithium_ion.py
+++ b/pybamm/models/submodels/interface/lithium_ion.py
@@ -3,6 +3,7 @@
 #
 from .base_interface import BaseInterface
 from . import inverse_kinetics, kinetics
+import pybamm
 
 
 class BaseInterfaceLithiumIon(BaseInterface):
@@ -43,6 +44,16 @@ class BaseInterfaceLithiumIon(BaseInterface):
         c_e = variables[self.domain + " electrolyte concentration"]
         T = variables[self.domain + " electrode temperature"]
 
+        # If variable was broadcast, take only the orphan
+        if (
+            isinstance(c_s_surf, pybamm.Broadcast)
+            and isinstance(c_e, pybamm.Broadcast)
+            and isinstance(T, pybamm.Broadcast)
+        ):
+            c_s_surf = c_s_surf.orphans[0]
+            c_e = c_e.orphans[0]
+            T = T.orphans[0]
+
         if self.domain == "Negative":
             prefactor = self.param.m_n(T) / self.param.C_r_n
 
@@ -75,6 +86,11 @@ class BaseInterfaceLithiumIon(BaseInterface):
         c_s_surf = variables[self.domain + " particle surface concentration"]
         T = variables[self.domain + " electrode temperature"]
 
+        # If variable was broadcast, take only the orphan
+        if isinstance(c_s_surf, pybamm.Broadcast) and isinstance(T, pybamm.Broadcast):
+            c_s_surf = c_s_surf.orphans[0]
+            T = T.orphans[0]
+
         if self.domain == "Negative":
             ocp = self.param.U_n(c_s_surf, T)
             dUdT = self.param.dUdT_n(c_s_surf)
diff --git a/pybamm/parameters/print_parameters.py b/pybamm/parameters/print_parameters.py
index 381c3f369..a78b9f31f 100644
--- a/pybamm/parameters/print_parameters.py
+++ b/pybamm/parameters/print_parameters.py
@@ -67,7 +67,10 @@ def print_parameters(parameters, parameter_values, output_file=None):
                 proc_symbol = parameter_values.process_symbol(symbol)
                 if not (
                     callable(proc_symbol)
-                    or isinstance(proc_symbol, pybamm.Concatenation)
+                    or any(
+                        isinstance(x, (pybamm.Concatenation, pybamm.Broadcast))
+                        for x in proc_symbol.pre_order()
+                    )
                 ):
                     evaluated_parameters[name].append(proc_symbol.evaluate(t=0))
 
diff --git a/pybamm/parameters/standard_parameters_lithium_ion.py b/pybamm/parameters/standard_parameters_lithium_ion.py
index bdf1457ed..5818fbe45 100644
--- a/pybamm/parameters/standard_parameters_lithium_ion.py
+++ b/pybamm/parameters/standard_parameters_lithium_ion.py
@@ -92,7 +92,12 @@ b_s_p = pybamm.geometric_parameters.b_s_p
 # Electrochemical reactions
 ne_n = pybamm.Parameter("Negative electrode electrons in reaction")
 ne_p = pybamm.Parameter("Positive electrode electrons in reaction")
-C_dl_dimensional = pybamm.Parameter("Double-layer capacity [F.m-2]")
+C_dl_n_dimensional = pybamm.Parameter(
+    "Negative electrode double-layer capacity [F.m-2]"
+)
+C_dl_p_dimensional = pybamm.Parameter(
+    "Positive electrode double-layer capacity [F.m-2]"
+)
 
 
 # Initial conditions
@@ -283,10 +288,15 @@ centre_y_tab_p = pybamm.geometric_parameters.centre_y_tab_p
 centre_z_tab_p = pybamm.geometric_parameters.centre_z_tab_p
 
 # Microscale geometry
-var = pybamm.standard_spatial_vars
-epsilon_n = pybamm.FunctionParameter("Negative electrode porosity", var.x_n)
-epsilon_s = pybamm.FunctionParameter("Separator porosity", var.x_s)
-epsilon_p = pybamm.FunctionParameter("Positive electrode porosity", var.x_p)
+epsilon_n = pybamm.FunctionParameter(
+    "Negative electrode porosity", pybamm.standard_spatial_vars.x_n
+)
+epsilon_s = pybamm.FunctionParameter(
+    "Separator porosity", pybamm.standard_spatial_vars.x_s
+)
+epsilon_p = pybamm.FunctionParameter(
+    "Positive electrode porosity", pybamm.standard_spatial_vars.x_p
+)
 epsilon = pybamm.Concatenation(epsilon_n, epsilon_s, epsilon_p)
 epsilon_s_n = pybamm.Parameter("Negative electrode active material volume fraction")
 epsilon_s_p = pybamm.Parameter("Positive electrode active material volume fraction")
@@ -313,7 +323,7 @@ alpha_prime = alpha / delta
 
 # Electrolyte Properties
 t_plus = pybamm.Parameter("Cation transference number")
-beta_surf = 0
+beta_surf = pybamm.Scalar(0)
 s = 1 - t_plus
 
 
@@ -327,10 +337,10 @@ def chi(c_e):
 
 # Electrochemical Reactions
 C_dl_n = (
-    C_dl_dimensional * potential_scale / interfacial_current_scale_n / tau_discharge
+    C_dl_n_dimensional * potential_scale / interfacial_current_scale_n / tau_discharge
 )
 C_dl_p = (
-    C_dl_dimensional * potential_scale / interfacial_current_scale_p / tau_discharge
+    C_dl_p_dimensional * potential_scale / interfacial_current_scale_p / tau_discharge
 )
 
 # Electrical
diff --git a/pybamm/processed_variable.py b/pybamm/processed_variable.py
index 96f97f18e..23d74a461 100644
--- a/pybamm/processed_variable.py
+++ b/pybamm/processed_variable.py
@@ -149,22 +149,24 @@ class ProcessedVariable(object):
         self.entries = entries
         self.dimensions = 2
         if self.domain[0] in ["negative particle", "positive particle"]:
-            self.spatial_var_name = "r"
+            self.first_dimension = "r"
             self.r_sol = space
         elif self.domain[0] in [
             "negative electrode",
             "separator",
             "positive electrode",
         ]:
-            self.spatial_var_name = "x"
+            self.first_dimension = "x"
             self.x_sol = space
         elif self.domain == ["current collector"]:
-            self.spatial_var_name = "z"
+            self.first_dimension = "z"
             self.z_sol = space
         else:
-            self.spatial_var_name = "x"
+            self.first_dimension = "x"
             self.x_sol = space
 
+        self.first_dim_pts = space
+
         # set up interpolation
         # note that the order of 't' and 'space' is the reverse of what you'd expect
 
@@ -242,6 +244,8 @@ class ProcessedVariable(object):
         # assign attributes for reference
         self.entries = entries
         self.dimensions = 3
+        self.first_dim_pts = first_dim_pts
+        self.second_dim_pts = second_dim_pts
 
         # set up interpolation
         self._interpolation_function = interp.RegularGridInterpolator(
@@ -335,7 +339,7 @@ class ProcessedVariable(object):
 
     def call_2D(self, t, x, r, z):
         "Evaluate a 2D variable"
-        spatial_var = eval_dimension_name(self.spatial_var_name, x, r, None, z)
+        spatial_var = eval_dimension_name(self.first_dimension, x, r, None, z)
         return self._interpolation_function(t, spatial_var)
 
     def call_3D(self, t, x, r, y, z):
diff --git a/pybamm/quick_plot.py b/pybamm/quick_plot.py
index 4d44fb8d7..9c49794e2 100644
--- a/pybamm/quick_plot.py
+++ b/pybamm/quick_plot.py
@@ -195,8 +195,8 @@ class QuickPlot(object):
 
             # Set the x variable for any two-dimensional variables
             if first_variable.dimensions == 2:
-                spatial_variable_key = first_variable.spatial_var_name
-                spatial_variable_value = first_variable.mesh[0].edges
+                spatial_variable_key = first_variable.first_dimension
+                spatial_variable_value = first_variable.first_dim_pts
                 self.spatial_variable[key] = (
                     spatial_variable_key,
                     spatial_variable_value,


Test patch:
diff --git a/tests/integration/test_models/standard_output_comparison.py b/tests/integration/test_models/standard_output_comparison.py
index e4aa69e0f..b36ac1a7b 100644
--- a/tests/integration/test_models/standard_output_comparison.py
+++ b/tests/integration/test_models/standard_output_comparison.py
@@ -51,9 +51,9 @@ class StandardOutputComparison(object):
         self.run_test_class(VariablesComparison, skip_first_timestep)
 
         if self.chemistry == "Lithium-ion":
-            self.run_test_class(ParticleConcentrationComparison)
+            self.run_test_class(ParticleConcentrationComparison, skip_first_timestep)
         elif self.chemistry == "Lead-acid":
-            self.run_test_class(PorosityComparison)
+            self.run_test_class(PorosityComparison, skip_first_timestep)
 
 
 class BaseOutputComparison(object):
@@ -67,13 +67,14 @@ class BaseOutputComparison(object):
         model_variables = [solution[var] for solution in self.solutions]
         var0 = model_variables[0]
 
-        if var0.mesh is None:
-            x = None
-        else:
-            x = var0.mesh[0].nodes
+        spatial_pts = {}
+        if var0.dimensions >= 2:
+            spatial_pts[var0.first_dimension] = var0.first_dim_pts
+        if var0.dimensions >= 3:
+            spatial_pts[var0.second_dimension] = var0.second_dim_pts
 
         # Calculate tolerance based on the value of var0
-        maxvar0 = np.max(abs(var0(self.t, x)))
+        maxvar0 = np.max(abs(var0(self.t, **spatial_pts)))
         if maxvar0 < 1e-14:
             decimal = -int(np.log10(tol))
         else:
@@ -82,7 +83,7 @@ class BaseOutputComparison(object):
         for model_var in model_variables[1:]:
             np.testing.assert_equal(var0.dimensions, model_var.dimensions)
             np.testing.assert_array_almost_equal(
-                model_var(self.t, x), var0(self.t, x), decimal
+                model_var(self.t, **spatial_pts), var0(self.t, **spatial_pts), decimal
             )
 
 
@@ -123,7 +124,6 @@ class VariablesComparison(BaseOutputComparison):
         self.compare("X-averaged negative electrode open circuit potential")
         self.compare("X-averaged positive electrode open circuit potential")
         self.compare("Terminal voltage")
-        self.compare("X-averaged electrolyte overpotential")
         self.compare("X-averaged solid phase ohmic losses")
         self.compare("Negative electrode reaction overpotential")
         self.compare("Positive electrode reaction overpotential")
diff --git a/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_compare_outputs.py b/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_compare_outputs.py
index a3700d62d..30b2e5980 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_compare_outputs.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lead_acid/test_compare_outputs.py
@@ -48,10 +48,9 @@ class TestCompareOutputs(unittest.TestCase):
         comparison = StandardOutputComparison(solutions)
         comparison.test_averages()
 
-    def test_compare_outputs_capacitance(self):
+    def test_compare_outputs_surface_form(self):
         """
-        Check that the leading-order model solution converges linearly in C_e to the
-        full model solution
+        Check that the models agree with the different surface forms
         """
         # load models
         options = [
@@ -59,6 +58,7 @@ class TestCompareOutputs(unittest.TestCase):
         ]
         model_combos = [
             ([pybamm.lead_acid.LOQS(opt) for opt in options]),
+            ([pybamm.lead_acid.Composite(opt) for opt in options]),
             ([pybamm.lead_acid.Full(opt) for opt in options]),
         ]
 
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs.py
new file mode 100644
index 000000000..973579d62
--- /dev/null
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs.py
@@ -0,0 +1,61 @@
+#
+# Tests for the surface formulation
+#
+import pybamm
+import numpy as np
+import unittest
+from tests import StandardOutputComparison
+
+
+class TestCompareOutputs(unittest.TestCase):
+    def test_compare_outputs_surface_form(self):
+        # load models
+        options = [
+            {"surface form": cap} for cap in [False, "differential", "algebraic"]
+        ]
+        model_combos = [
+            ([pybamm.lithium_ion.SPM(opt) for opt in options]),
+            ([pybamm.lithium_ion.SPMe(opt) for opt in options]),
+            ([pybamm.lithium_ion.DFN(opt) for opt in options]),
+        ]
+
+        for models in model_combos:
+            # load parameter values (same for all models)
+            param = models[0].default_parameter_values
+            param.update({"Current function [A]": 1})
+            for model in models:
+                param.process_model(model)
+
+            # set mesh
+            var = pybamm.standard_spatial_vars
+            var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 5, var.r_p: 5}
+
+            # discretise models
+            discs = {}
+            for model in models:
+                geometry = model.default_geometry
+                param.process_geometry(geometry)
+                mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
+                disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+                disc.process_model(model)
+                discs[model] = disc
+
+            # solve model
+            solutions = []
+            t_eval = np.linspace(0, 3600, 100)
+            for model in models:
+                solution = pybamm.CasadiSolver().solve(model, t_eval)
+                solutions.append(solution)
+
+            # compare outputs
+            comparison = StandardOutputComparison(solutions)
+            comparison.test_all(skip_first_timestep=True)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index 8082aae7d..be6eb491f 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -95,6 +95,18 @@ class TestDFN(unittest.TestCase):
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.DFN(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.DFN(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
     def test_particle_distribution_in_x(self):
         model = pybamm.lithium_ion.DFN()
         param = model.default_parameter_values
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
index cbed675bf..9a933dade 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
@@ -103,6 +103,18 @@ class TestSPM(unittest.TestCase):
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.SPM(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.SPM(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
index bf4866984..0709b9668 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
@@ -95,6 +95,18 @@ class TestSPMe(unittest.TestCase):
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.SPMe(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.SPMe(options)
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py
index 901f43f1c..511b4c354 100644
--- a/tests/unit/test_experiments/test_simulation_with_experiment.py
+++ b/tests/unit/test_experiments/test_simulation_with_experiment.py
@@ -13,7 +13,7 @@ class TestSimulationExperiment(unittest.TestCase):
                 "Charge at 1 A until 4.1 V",
                 "Hold at 4.1 V until 50 mA",
                 "Discharge at 2 W for 1 hour",
-            ],
+            ]
         )
         model = pybamm.lithium_ion.DFN()
         sim = pybamm.Simulation(model, experiment=experiment)
@@ -48,7 +48,7 @@ class TestSimulationExperiment(unittest.TestCase):
         self.assertEqual(sim._experiment_inputs[3]["Voltage cut-off [V]"], -1e10)
 
         self.assertEqual(
-            sim._experiment_times, [3600, 7 * 24 * 3600, 7 * 24 * 3600, 3600],
+            sim._experiment_times, [3600, 7 * 24 * 3600, 7 * 24 * 3600, 3600]
         )
 
         self.assertIn(
@@ -75,11 +75,11 @@ class TestSimulationExperiment(unittest.TestCase):
                 "Charge at 1 A until 4.1 V",
                 "Hold at 4.1 V until C/2",
                 "Discharge at 2 W for 1 hour",
-            ],
+            ]
         )
         model = pybamm.lithium_ion.SPM()
         sim = pybamm.Simulation(model, experiment=experiment)
-        sim.solve()
+        sim.solve(solver=pybamm.CasadiSolver())
         self.assertEqual(sim._solution.termination, "final time")
 
     def test_run_experiment_breaks_early(self):
@@ -89,7 +89,7 @@ class TestSimulationExperiment(unittest.TestCase):
         pybamm.set_logging_level("ERROR")
         # giving the time, should get ignored
         t_eval = [0, 1]
-        sim.solve(t_eval)
+        sim.solve(t_eval, solver=pybamm.CasadiSolver())
         pybamm.set_logging_level("WARNING")
         self.assertIn("event", sim._solution.termination)
 
diff --git a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_composite.py b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_composite.py
index 183fe3cc4..500a00d37 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_composite.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_composite.py
@@ -54,11 +54,6 @@ class TestLeadAcidCompositeMultiDimensional(unittest.TestCase):
 
 
 class TestLeadAcidCompositeWithSideReactions(unittest.TestCase):
-    def test_well_posed_differential(self):
-        options = {"surface form": "differential", "side reactions": ["oxygen"]}
-        model = pybamm.lead_acid.Composite(options)
-        model.check_well_posedness()
-
     def test_well_posed_algebraic(self):
         options = {"surface form": "algebraic", "side reactions": ["oxygen"]}
         model = pybamm.lead_acid.Composite(options)
@@ -69,10 +64,6 @@ class TestLeadAcidCompositeWithSideReactions(unittest.TestCase):
 
 
 class TestLeadAcidCompositeExtended(unittest.TestCase):
-    def test_well_posed(self):
-        model = pybamm.lead_acid.CompositeExtended()
-        model.check_well_posedness()
-
     def test_well_posed_differential_side_reactions(self):
         options = {"surface form": "differential", "side reactions": ["oxygen"]}
         model = pybamm.lead_acid.CompositeExtended(options)
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
index 3ac309316..9e5bda3d5 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
@@ -7,8 +7,6 @@ import unittest
 
 class TestBaseLithiumIonModel(unittest.TestCase):
     def test_incompatible_options(self):
-        with self.assertRaisesRegex(pybamm.OptionError, "surface form not implemented"):
-            pybamm.lithium_ion.BaseModel({"surface form": "differential"})
         with self.assertRaisesRegex(pybamm.OptionError, "convection not implemented"):
             pybamm.lithium_ion.BaseModel({"convection": True})
 
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index 4fcabbccb..f6517319d 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -164,6 +164,16 @@ class TestDFN(unittest.TestCase):
         model = pybamm.lithium_ion.DFN(options)
         model.check_well_posedness()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.DFN(options)
+        model.check_well_posedness()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.DFN(options)
+        model.check_well_posedness()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
index 87da7b98b..9d75316b7 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
@@ -177,6 +177,16 @@ class TestSPM(unittest.TestCase):
         model = pybamm.lithium_ion.SPM(options)
         model.check_well_posedness()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.SPM(options)
+        model.check_well_posedness()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.SPM(options)
+        model.check_well_posedness()
+
 
 class TestSPMExternalCircuits(unittest.TestCase):
     def test_well_posed_voltage(self):
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
index fb1e48dac..9020d70b5 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
@@ -164,6 +164,16 @@ class TestSPMe(unittest.TestCase):
         model = pybamm.lithium_ion.SPMe(options)
         model.check_well_posedness()
 
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        model = pybamm.lithium_ion.SPMe(options)
+        model.check_well_posedness()
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        model = pybamm.lithium_ion.SPMe(options)
+        model.check_well_posedness()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")