Issue 791 complete models

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added some basic models (BasicSPM and BasicDFN) in order to clearly demonstrate the PyBaMM model structure for battery models ([#795](https://github.com/pybamm-team/PyBaMM/pull/795))
 -   Added the harmonic mean to the Finite Volume method, which is now used when computing fluxes ([#783](https://github.com/pybamm-team/PyBaMM/pull/783))
 -   Refactored `Solution` to make it a dictionary that contains all of the solution variables. This automatically creates `ProcessedVariable` objects when required, so that the solution can be obtained much more easily. ([#781](https://github.com/pybamm-team/PyBaMM/pull/781))
 -   Added notebook to explain broadcasts ([#776](https://github.com/pybamm-team/PyBaMM/pull/776))

docs/source/models/lithium_ion/dfn.rst

@@ -3,3 +3,6 @@ Doyle-Fuller-Newman (DFN)
 
 .. autoclass:: pybamm.lithium_ion.DFN
     :members:
+
+.. autoclass:: pybamm.lithium_ion.BasicDFN
+    :members:

docs/source/models/lithium_ion/spm.rst

@@ -3,3 +3,6 @@ Single Particle Model (SPM)
 
 .. autoclass:: pybamm.lithium_ion.SPM
     :members:
+
+.. autoclass:: pybamm.lithium_ion.BasicSPM
+    :members:

examples/notebooks/README.md

@@ -48,7 +48,9 @@ The following notebooks are specific to different stages of the PyBaMM pipeline,
 
 ### Models
 
-The following models are implemented and can easily be used or [compared](./models/lead-acid.ipynb). We always welcome [new models](https://pybamm.readthedocs.io/en/latest/tutorials/add-model.html)!
+Several battery models are implemented and can easily be used or [compared](./models/lead-acid.ipynb). The notebooks below show the solution of each individual model. We always welcome [new models](https://pybamm.readthedocs.io/en/latest/tutorials/add-model.html)!
+
+Once you are comfortable with the expression tree structure, a good starting point to understand the models in PyBaMM is to take a look at the [basic SPM](https://github.com/pybamm-team/PyBaMM/blob/master/pybamm/models/full_battery_models/lithium_ion/basic_spm.py) and [basic DFN](https://github.com/pybamm-team/PyBaMM/blob/master/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py), since these define the entire model (variables, equations, initial and boundary conditions, events) in a single class and so are easier to understand. However, we recommend that you subsequently use the full models as they offer much greater flexibility for coupling different physical effects and visualising a greater range of variables.
 
 #### Lithium-ion models
 

pybamm/expression_tree/functions.py

@@ -159,7 +159,7 @@ def evaluate(self, t=None, y=None, u=None, known_evals=None):
                 evaluated_children = [None] * len(self.children)
                 for i, child in enumerate(self.children):
                     evaluated_children[i], known_evals = child.evaluate(
-                        t, y, known_evals=known_evals
+                        t, y, u, known_evals=known_evals
                     )
                 known_evals[self.id] = self._function_evaluate(evaluated_children)
             return known_evals[self.id], known_evals

pybamm/models/full_battery_models/lithium_ion/__init__.py

@@ -5,3 +5,5 @@
 from .spm import SPM
 from .spme import SPMe
 from .dfn import DFN
+from .basic_dfn import BasicDFN
+from .basic_spm import BasicSPM

pybamm/models/full_battery_models/lithium_ion/basic_dfn.py

@@ -0,0 +1,271 @@
+#
+# Basic Doyle-Fuller-Newman (DFN) Model
+#
+import pybamm
+from .base_lithium_ion_model import BaseModel
+
+
+class BasicDFN(BaseModel):
+    """Doyle-Fuller-Newman (DFN) model of a lithium-ion battery, from [2]_.
+
+    This class differs from the :class:`pybamm.lithium_ion.DFN` model class in that it
+    shows the whole model in a single class. This comes at the cost of flexibility in
+    comparing different physical effects, and in general the main DFN class should be
+    used instead.
+
+    Parameters
+    ----------
+    name : str, optional
+        The name of the model.
+
+    References
+    ----------
+    .. [2] SG Marquis, V Sulzer, R Timms, CP Please and SJ Chapman. “An asymptotic
+           derivation of a single particle model with electrolyte”. In: arXiv preprint
+           arXiv:1905.12553 (2019).
+
+
+    **Extends:** :class:`pybamm.lithium_ion.BaseModel`
+    """
+
+    def __init__(self, name="Doyle-Fuller-Newman model"):
+        super().__init__({}, name)
+        # `param` is a class containing all the relevant parameters and functions for
+        # this model. These are purely symbolic at this stage, and will be set by the
+        # `ParameterValues` class when the model is processed.
+        param = self.param
+
+        ######################
+        # Variables
+        ######################
+        # Variables that depend on time only are created without a domain
+        Q = pybamm.Variable("Discharge capacity [A.h]")
+        # Variables that vary spatially are created with a domain
+        c_e_n = pybamm.Variable(
+            "Negative electrolyte concentration", domain="negative electrode",
+        )
+        c_e_s = pybamm.Variable(
+            "Separator electrolyte concentration", domain="separator",
+        )
+        c_e_p = pybamm.Variable(
+            "Positive electrolyte concentration", domain="positive electrode",
+        )
+        # Concatenations combine several variables into a single variable, to simplify
+        # implementing equations that hold over several domains
+        c_e = pybamm.Concatenation(c_e_n, c_e_s, c_e_p)
+
+        # Electrolyte potential
+        phi_e_n = pybamm.Variable(
+            "Negative electrolyte potential", domain="negative electrode",
+        )
+        phi_e_s = pybamm.Variable(
+            "Separator electrolyte potential", domain="separator",
+        )
+        phi_e_p = pybamm.Variable(
+            "Positive electrolyte potential", domain="positive electrode",
+        )
+        phi_e = pybamm.Concatenation(phi_e_n, phi_e_s, phi_e_p)
+
+        # Electrode potential
+        phi_s_n = pybamm.Variable(
+            "Negative electrode potential", domain="negative electrode",
+        )
+        phi_s_p = pybamm.Variable(
+            "Positive electrode potential", domain="positive electrode",
+        )
+        # Particle concentrations are variables on the particle domain, but also vary in
+        # the x-direction (electrode domain) and so must be provided with auxiliary
+        # domains
+        c_s_n = pybamm.Variable(
+            "Negative particle concentration",
+            domain="negative particle",
+            auxiliary_domains={"secondary": "negative electrode"},
+        )
+        c_s_p = pybamm.Variable(
+            "Positive particle concentration",
+            domain="positive particle",
+            auxiliary_domains={"secondary": "positive electrode"},
+        )
+
+        # Constant temperature
+        T = param.T_init
+
+        ######################
+        # Other set-up
+        ######################
+
+        # Current density
+        i_cell = param.current_with_time
+
+        # Porosity
+        # Primary broadcasts are used to broadcast scalar quantities across a domain
+        # into a vector of the right shape, for multiplying with other vectors
+        eps_n = pybamm.PrimaryBroadcast(param.epsilon_n, "negative electrode")
+        eps_s = pybamm.PrimaryBroadcast(param.epsilon_s, "separator")
+        eps_p = pybamm.PrimaryBroadcast(param.epsilon_p, "positive electrode")
+        eps = pybamm.Concatenation(eps_n, eps_s, eps_p)
+
+        # Tortuosity
+        tor = pybamm.Concatenation(
+            eps_n ** param.b_e_n, eps_s ** param.b_e_s, eps_p ** param.b_e_p
+        )
+
+        # Interfacial reactions
+        # Surf takes the surface value of a variable, i.e. its boundary value on the
+        # right side. This is also accessible via `boundary_value(x, "right")`, with
+        # "left" providing the boundary value of the left side
+        c_s_surf_n = pybamm.surf(c_s_n)
+        j0_n = (
+            param.m_n(T)
+            / param.C_r_n
+            * c_e_n ** (1 / 2)
+            * c_s_surf_n ** (1 / 2)
+            * (1 - c_s_surf_n) ** (1 / 2)
+        )
+        j_n = (
+            2
+            * j0_n
+            * pybamm.sinh(
+                param.ne_n / 2 * (phi_s_n - phi_e_n - param.U_n(c_s_surf_n, T))
+            )
+        )
+        c_s_surf_p = pybamm.surf(c_s_p)
+        j0_p = (
+            param.gamma_p
+            * param.m_p(T)
+            / param.C_r_p
+            * c_e_p ** (1 / 2)
+            * c_s_surf_p ** (1 / 2)
+            * (1 - c_s_surf_p) ** (1 / 2)
+        )
+        j_s = pybamm.PrimaryBroadcast(0, "separator")
+        j_p = (
+            2
+            * j0_p
+            * pybamm.sinh(
+                param.ne_p / 2 * (phi_s_p - phi_e_p - param.U_p(c_s_surf_p, T))
+            )
+        )
+        j = pybamm.Concatenation(j_n, j_s, j_p)
+
+        ######################
+        # State of Charge
+        ######################
+        I = param.dimensional_current_with_time
+        # The `rhs` dictionary contains differential equations, with the key being the
+        # variable in the d/dt
+        self.rhs[Q] = I * param.timescale / 3600
+        # Initial conditions must be provided for the ODEs
+        self.initial_conditions[Q] = pybamm.Scalar(0)
+
+        ######################
+        # Particles
+        ######################
+
+        # The div and grad operators will be converted to the appropriate matrix
+        # multiplication at the discretisation stage
+        N_s_n = -param.D_n(c_s_n, T) * pybamm.grad(c_s_n)
+        N_s_p = -param.D_p(c_s_p, T) * pybamm.grad(c_s_p)
+        self.rhs[c_s_n] = -(1 / param.C_n) * pybamm.div(N_s_n)
+        self.rhs[c_s_p] = -(1 / param.C_p) * pybamm.div(N_s_p)
+        # Boundary conditions must be provided for equations with spatial derivatives
+        self.boundary_conditions[c_s_n] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (-param.C_n * j_n / param.a_n, "Neumann"),
+        }
+        self.boundary_conditions[c_s_p] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"),
+        }
+        self.initial_conditions[c_s_n] = param.c_n_init
+        self.initial_conditions[c_s_p] = param.c_p_init
+        # Events specify points at which a solution should terminate
+        self.events.update(
+            {
+                "Minimum negative particle surface concentration": (
+                    pybamm.min(c_s_surf_n) - 0.01
+                ),
+                "Maximum negative particle surface concentration": (1 - 0.01)
+                - pybamm.max(c_s_surf_n),
+                "Minimum positive particle surface concentration": (
+                    pybamm.min(c_s_surf_p) - 0.01
+                ),
+                "Maximum positive particle surface concentration": (1 - 0.01)
+                - pybamm.max(c_s_surf_p),
+            }
+        )
+        ######################
+        # Current in the solid
+        ######################
+        i_s_n = -param.sigma_n * (1 - eps_n) ** param.b_s_n * pybamm.grad(phi_s_n)
+        sigma_eff_p = param.sigma_p * (1 - eps_p) ** param.b_s_p
+        i_s_p = -sigma_eff_p * pybamm.grad(phi_s_p)
+        # The `algebraic` dictionary contains differential equations, with the key being
+        # the main scalar variable of interest in the equation
+        self.algebraic[phi_s_n] = pybamm.div(i_s_n) + j_n
+        self.algebraic[phi_s_p] = pybamm.div(i_s_p) + j_p
+        self.boundary_conditions[phi_s_n] = {
+            "left": (pybamm.Scalar(0), "Dirichlet"),
+            "right": (pybamm.Scalar(0), "Neumann"),
+        }
+        self.boundary_conditions[phi_s_p] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (i_cell / pybamm.boundary_value(-sigma_eff_p, "right"), "Neumann"),
+        }
+        # Initial conditions must also be provided for algebraic equations, as an
+        # initial guess for a root-finding algorithm which calculates consistent initial
+        # conditions
+        self.initial_conditions[phi_s_n] = pybamm.Scalar(0)
+        self.initial_conditions[phi_s_p] = param.U_p(
+            param.c_p_init, param.T_init
+        ) - param.U_n(param.c_n_init, param.T_init)
+
+        ######################
+        # Current in the electrolyte
+        ######################
+        i_e = (param.kappa_e(c_e, T) * tor * param.gamma_e / param.C_e) * (
+            param.chi(c_e) * pybamm.grad(c_e) / c_e - pybamm.grad(phi_e)
+        )
+        self.algebraic[phi_e] = pybamm.div(i_e) - j
+        self.boundary_conditions[phi_e] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (pybamm.Scalar(0), "Neumann"),
+        }
+        self.initial_conditions[phi_e] = -param.U_n(param.c_n_init, param.T_init)
+
+        ######################
+        # Electrolyte concentration
+        ######################
+        N_e = -tor * param.D_e(c_e, T) * pybamm.grad(c_e)
+        self.rhs[c_e] = (1 / eps) * (
+            -pybamm.div(N_e) / param.C_e + (1 - param.t_plus) * j / param.gamma_e
+        )
+        self.boundary_conditions[c_e] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (pybamm.Scalar(0), "Neumann"),
+        }
+        self.initial_conditions[c_e] = param.c_e_init
+        self.events["Zero electrolyte concentration cut-off"] = pybamm.min(c_e) - 0.002
+
+        ######################
+        # (Some) variables
+        ######################
+        voltage = pybamm.boundary_value(phi_s_p, "right")
+        # The `variables` dictionary contains all variables that might be useful for
+        # visualising the solution of the model
+        self.variables = {
+            "Negative particle surface concentration": c_s_surf_n,
+            "Electrolyte concentration": c_e,
+            "Positive particle surface concentration": c_s_surf_p,
+            "Current [A]": I,
+            "Negative electrode potential": phi_s_n,
+            "Electrolyte potential": phi_e,
+            "Positive electrode potential": phi_s_p,
+            "Terminal voltage": voltage,
+        }
+        self.events["Minimum voltage"] = voltage - param.voltage_low_cut
+        self.events["Maximum voltage"] = voltage - param.voltage_high_cut
+
+    @property
+    def default_geometry(self):
+        return pybamm.Geometry("1D macro", "1+1D micro")