Issue 247 dfn as odes

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))

pybamm/models/full_battery_models/base_battery_model.py

@@ -224,10 +224,10 @@ def options(self, extra_options):
 
         # 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"

pybamm/models/full_battery_models/lithium_ion/dfn.py

@@ -84,24 +84,40 @@ def set_particle_submodel(self):
 
     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"]

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 @@ def set_convection_submodel(self):
 
     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 @@ def set_positive_electrode_submodel(self):
     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)

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 @@ def _get_exchange_current_density(self, variables):
         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 @@ def _get_open_circuit_potential(self, variables):
         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)

pybamm/parameters/standard_parameters_lithium_ion.py

@@ -92,7 +92,12 @@
 # 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
@@ -327,10 +332,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

pybamm/processed_variable.py

@@ -149,22 +149,24 @@ def initialise_2D(self):
         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 @@ def initialise_3D(self):
         # 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 @@ def __call__(self, t=None, x=None, r=None, y=None, z=None, warn=True):
 
     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):

pybamm/quick_plot.py

@@ -195,8 +195,8 @@ def set_output_variables(self, output_variables, solutions):
 
             # 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,

tests/integration/test_models/standard_output_comparison.py

@@ -51,9 +51,9 @@ def test_all(self, skip_first_timestep=False):
         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 @@ def compare(self, var, tol=1e-2):
         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 @@ def compare(self, var, tol=1e-2):
         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 @@ def test_all(self):
         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")

tests/integration/test_models/test_full_battery_models/test_lead_acid/test_compare_outputs.py

@@ -48,17 +48,17 @@ def test_compare_averages_asymptotics(self):
         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 = [
             {"surface form": cap} for cap in [False, "differential", "algebraic"]
         ]
         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]),
         ]
 

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()