Add half cell submodels

CHANGELOG.md

@@ -1,3 +1,9 @@
+# [Unreleased](https://github.com/pybamm-team/PyBaMM)
+
+## Features
+
+-   The DFN model can now be used directly (instead of `BasicDFNHalfCell`) to simulate a half-cell ([#1600](https://github.com/pybamm-team/PyBaMM/pull/1600))
+
 # [v21.08](https://github.com/pybamm-team/PyBaMM/tree/v21.08) - 2021-08-26
 
 This release introduces:

docs/source/models/lithium_ion/electrode_soh.rst

@@ -0,0 +1,10 @@
+Electrode SOH models
+====================
+
+.. autoclass:: pybamm.lithium_ion.ElectrodeSOH
+    :members:
+
+.. autoclass:: pybamm.lithium_ion.ElectrodeSOHHalfCell
+    :members:
+
+.. autofunction:: pybamm.lithium_ion.get_initial_stoichiometries

docs/source/models/lithium_ion/index.rst

@@ -10,3 +10,4 @@ Lithium-ion Models
   dfn
   newman_tobias
   yang2017
+  electrode_soh

docs/source/models/submodels/electrode/ohm/index.rst

@@ -8,3 +8,4 @@ Ohmic
   composite_ohm
   full_ohm
   surface_form_ohm
+  li_metal_explicit

docs/source/models/submodels/electrode/ohm/li_metal_explicit.rst

@@ -0,0 +1,5 @@
+Explicit potential drop for lithium metal
+=========================================
+
+.. autoclass:: pybamm.electrode.ohm.LithiumMetalExplicit
+    :members:

docs/tutorials/add-parameter-values.rst

@@ -16,7 +16,7 @@ At the stage of creating a model, we use :class:`pybamm.Parameter` and :class:`p
 We then create a :class:`ParameterValues` class, using a specific set of parameters, to iterate through the model and replace any :class:`pybamm.Parameter` objects with a :class:`pybamm.Scalar` and any :class:`pybamm.FunctionParameter` objects with a :class:`pybamm.Function`.
 
 For an example of how the parameter values work, see the
-`parameter values notebook <https://github.com/pybamm-team/PyBaMM/blob/develop/examples/notebooks/parameter-values.ipynb>`_.
+`parameter values notebook <https://github.com/pybamm-team/PyBaMM/blob/develop/examples/notebooks/parameterization/parameter-values.ipynb>`_.
 
 Adding a set of parameters values
 ---------------------------------

examples/notebooks/parameterization/parameterization.ipynb

@@ -12,8 +12,8 @@
     "\n",
     "- The API documentation of [Parameters](https://pybamm.readthedocs.io/en/latest/source/parameters/index.html) can be found at [pybamm.readthedocs.io](https://pybamm.readthedocs.io/)\n",
     "- [Setting parameter values](https://github.com/pybamm-team/PyBaMM/blob/develop/examples/notebooks/Getting%20Started/Tutorial%204%20-%20Setting%20parameter%20values.ipynb) can be found at `pybamm/examples/notebooks/Getting Started/Tutorial 4 - Setting parameter values.ipynb`. This explains the basics of how to set the parameters of a model (in less detail than here).\n",
-    "- [parameter-management.ipynb](https://github.com/pybamm-team/PyBaMM/blob/develop/examples/notebooks/parameter-management.ipynb) can be found at `pybamm/examples/notebooks/parameter-management.ipynb`. This explains how to add and edit parameters in the `pybamm/input` directories\n",
-    "- [parameter-values.ipynb](https://github.com/pybamm-team/PyBaMM/blob/develop/examples/notebooks/parameter-values.ipynb) can be found at `pybamm/examples/notebooks/parameter-values.ipynb`. This explains the basics of the `ParameterValues` class.\n"
+    "- [parameter-management.ipynb](https://github.com/pybamm-team/PyBaMM/blob/develop/examples/notebooks/parameterization/parameter-management.ipynb) can be found at `pybamm/examples/notebooks/parameterization/parameter-management.ipynb`. This explains how to add and edit parameters in the `pybamm/input` directories\n",
+    "- [parameter-values.ipynb](https://github.com/pybamm-team/PyBaMM/blob/develop/examples/notebooks/parameterization/parameter-values.ipynb) can be found at `pybamm/examples/notebooks/parameterization/parameter-values.ipynb`. This explains the basics of the `ParameterValues` class.\n"
    ]
   },
   {

examples/scripts/DFN_half_cell.py

@@ -9,58 +9,48 @@
 
 # load model
 options = {"working electrode": "positive"}
-model = pybamm.lithium_ion.BasicDFNHalfCell(options=options)
+model1 = pybamm.lithium_ion.DFN(options=options)
+model2 = pybamm.lithium_ion.BasicDFNHalfCell(options=options)
 
-# create geometry
-geometry = model.default_geometry
+sols = []
+for model in [model1, model2]:
+    # create geometry
+    geometry = model.default_geometry
 
-# load parameter values
-chemistry = pybamm.parameter_sets.Chen2020
-param = pybamm.ParameterValues(chemistry=chemistry)
+    # load parameter values
+    chemistry = pybamm.parameter_sets.Xu2019
+    param = pybamm.ParameterValues(chemistry=chemistry)
 
-# add lithium counter electrode parameter values
-param.update(
-    {
-        "Lithium counter electrode exchange-current density [A.m-2]": 12.6,
-        "Lithium counter electrode conductivity [S.m-1]": 1.0776e7,
-        "Lithium counter electrode thickness [m]": 250e-6,
-    },
-    check_already_exists=False,
-)
-
-param["Initial concentration in negative electrode [mol.m-3]"] = 1000
-param["Current function [A]"] = 2.5
-
-# process model and geometry
-param.process_model(model)
-param.process_geometry(geometry)
+    # process model and geometry
+    param.process_model(model)
+    param.process_geometry(geometry)
 
-# set mesh
-var_pts = model.default_var_pts
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
+    # set mesh
+    var_pts = model.default_var_pts
+    mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
 
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
+    # discretise model
+    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+    disc.process_model(model)
 
-# solve model
-t_eval = np.linspace(0, 7200, 1000)
-solver = pybamm.CasadiSolver(mode="safe", atol=1e-6, rtol=1e-3)
-solution = solver.solve(model, t_eval)
+    # solve model
+    t_eval = np.linspace(0, 7200, 1000)
+    solver = pybamm.CasadiSolver(mode="safe", atol=1e-6, rtol=1e-3)
+    solution = solver.solve(model, t_eval)
+    sols.append(solution)
 
 # plot
 plot = pybamm.QuickPlot(
-    solution,
+    sols,
     [
-        "Working particle concentration [mol.m-3]",
         "Electrolyte concentration [mol.m-3]",
         "Current [A]",
-        "Working electrode potential [V]",
+        "Positive electrode potential [V]",
         "Electrolyte potential [V]",
-        "Total electrolyte concentration",
-        "Total lithium in working electrode [mol]",
-        "Working electrode open circuit potential [V]",
-        ["Terminal voltage [V]", "Voltage drop in the cell [V]"],
+        "Total lithium in electrolyte [mol]",
+        "Positive electrode open circuit potential [V]",
+        ["Terminal voltage [V]"],
+        "Negative electrode potential drop [V]",
     ],
     time_unit="seconds",
     spatial_unit="um",

pybamm/expression_tree/concatenations.py

@@ -409,6 +409,8 @@ def intersect(s1, s2):
 
 def simplified_concatenation(*children):
     """Perform simplifications on a concatenation."""
+    # remove children that are None
+    children = list(filter(lambda x: x is not None, children))
     # Create Concatenation to easily read domains
     concat = Concatenation(*children)
     # Simplify concatenation of broadcasts all with the same child to a single

pybamm/expression_tree/symbol.py

@@ -35,8 +35,6 @@ def domain_size(domain):
         "negative particle size": 29,
         "positive particle size": 31,
     }
-    if isinstance(domain, str):
-        domain = [domain]
     if domain in [[], None]:
         size = 1
     elif all(dom in fixed_domain_sizes for dom in domain):
@@ -443,8 +441,6 @@ def relabel_tree(self, symbol, counter):
             name = "+"
         elif name == "**":
             name = "^"
-        elif name == "epsilon_s":
-            name = "ɛ"
 
         new_node = anytree.Node(str(counter), label=name)
         counter += 1

pybamm/expression_tree/unary_operators.py

@@ -598,7 +598,7 @@ def __init__(self, child, integration_variable):
             domain = child.domain
             auxiliary_domains = {
                 "secondary": child.auxiliary_domains["secondary"],
-                "tertiary": child.auxiliary_domains["tertiary"]
+                "tertiary": child.auxiliary_domains["tertiary"],
             }
         if any(isinstance(var, pybamm.SpatialVariable) for var in integration_variable):
             name += " {}".format(child.domain)
@@ -1268,17 +1268,19 @@ def x_average(symbol):
     if isinstance(symbol, (pybamm.PrimaryBroadcast, pybamm.FullBroadcast)):
         return symbol.reduce_one_dimension()
     # If symbol is a concatenation of Broadcasts, its average value is its child
-    elif (
-        isinstance(symbol, pybamm.Concatenation)
-        and all(isinstance(child, pybamm.Broadcast) for child in symbol.children)
-        and symbol.domain == ["negative electrode", "separator", "positive electrode"]
+    elif isinstance(symbol, pybamm.Concatenation) and all(
+        isinstance(child, pybamm.Broadcast) for child in symbol.children
     ):
-        a, b, c = [orp.orphans[0] for orp in symbol.orphans]
         geo = pybamm.geometric_parameters
         l_n = geo.l_n
         l_s = geo.l_s
         l_p = geo.l_p
-        out = (l_n * a + l_s * b + l_p * c) / (l_n + l_s + l_p)
+        if symbol.domain == ["negative electrode", "separator", "positive electrode"]:
+            a, b, c = [orp.orphans[0] for orp in symbol.orphans]
+            out = (l_n * a + l_s * b + l_p * c) / (l_n + l_s + l_p)
+        elif symbol.domain == ["separator", "positive electrode"]:
+            b, c = [orp.orphans[0] for orp in symbol.orphans]
+            out = (l_s * b + l_p * c) / (l_s + l_p)
         # To respect domains we may need to broadcast the child back out
         child = symbol.children[0]
         # If symbol being returned doesn't have empty domain, return it
@@ -1300,16 +1302,14 @@ def x_average(symbol):
         # Even if domain is "negative electrode", "separator", or
         # "positive electrode", and we know l, we still compute it as Integral(1, x)
         # as this will be easier to identify for simplifications later on
-        if (
-            symbol.domain == ["negative particle"] or
-            symbol.domain == ["negative particle size"]
-        ):
+        if symbol.domain == ["negative particle"] or symbol.domain == [
+            "negative particle size"
+        ]:
             x = pybamm.standard_spatial_vars.x_n
             l = geo.l_n
-        elif (
-            symbol.domain == ["positive particle"] or
-            symbol.domain == ["positive particle size"]
-        ):
+        elif symbol.domain == ["positive particle"] or symbol.domain == [
+            "positive particle size"
+        ]:
             x = pybamm.standard_spatial_vars.x_p
             l = geo.l_p
         else:
@@ -1486,17 +1486,15 @@ def size_average(symbol):
 
     # If symbol is a primary broadcast to "particle size", take the orphan
     elif isinstance(symbol, pybamm.PrimaryBroadcast) and symbol.domain in [
-        ["negative particle size"], ["positive particle size"]
+        ["negative particle size"],
+        ["positive particle size"],
     ]:
         return symbol.orphans[0]
     # If symbol is a secondary broadcast to "particle size" from "particle",
     # take the orphan
-    elif (
-        isinstance(symbol, pybamm.SecondaryBroadcast) and
-        symbol.domains["secondary"] in [
-            ["negative particle size"], ["positive particle size"]
-        ]
-    ):
+    elif isinstance(symbol, pybamm.SecondaryBroadcast) and symbol.domains[
+        "secondary"
+    ] in [["negative particle size"], ["positive particle size"]]:
         return symbol.orphans[0]
     # Otherwise, perform the integration in R
     else: