Tip

An interactive online version of this notebook is available, which can be accessed via

Alternatively, you may download this notebook and run it offline.

Attention

You are viewing this notebook on the latest version of the documentation, where these notebooks may not be compatible with the stable release of PyBaMM since they can contain features that are not yet released. We recommend viewing these notebooks from the stable version of the documentation. To install the latest version of PyBaMM that is compatible with the latest notebooks, build PyBaMM from source.

Comparing full and reduced-order models#

In the previous notebook we saw how to solve the problem of diffusion on a sphere, motivated by the problem in the negative particle in battery modelling. In this notebook we consider a reduced-order ODE model for the particle behaviour, suitable in the limit of fast diffusion. We also show how to compare the results of the full and reduced-order models.

In the limit of fast diffusion in the particles the concentration is uniform in $r$. This result in the following ODE model for the (uniform) concentration in the particle

\[\frac{\textrm{d} c}{\textrm{d} t} = -3\frac{j}{RF}\]

with the initial condition:

\[\left.c\right\vert_{t=0} = c_0,\]

where $c$$ is the concentration, $r$ the radial coordinate, $t$ time, $R$ the particle radius, $D$ the diffusion coefficient, $j$ the interfacial current density, $F$ Faraday’s constant, and $c_0$ the initial concentration.

As in the previous example we use the following parameters:

Symbol	Units	Value
$R$	m	$10 \times 10^{-6}$
$D$	m${^2}$ s$^{-1}$	$3.9 \times 10^{-14}$
$j$	A m$^{-2}$	$1.4$
$F$	C mol$^{-1}$	$96485$
$c_0$	mol m$^{-3}$	$2.5 \times 10^{4}$

Setting up the models#

As in the single particle diffusion example, we begin by importing the pybamm library into this notebook, along with any other packages we require. In this notebook we want to compare the results of the full and reduced-order models, so we create two empty pybamm.BaseModel objects. We can pass in a name when we create the model, for easy reference.

[11]:

%pip install "pybamm[plot,cite]" -q    # install PyBaMM if it is not installed
import pybamm
import numpy as np
import matplotlib.pyplot as plt

full_model = pybamm.BaseModel(name="full model")
reduced_model = pybamm.BaseModel(name="reduced model")

WARNING: You are using pip version 22.0.4; however, version 22.3.1 is available.
You should consider upgrading via the '/home/mrobins/git/PyBaMM/env/bin/python -m pip install --upgrade pip' command.
Note: you may need to restart the kernel to use updated packages.

It can be useful to add the models to a list so that we can perform the same operations on each model easily

[12]:

models = [full_model, reduced_model]

We then define our parameters, as seen previously,

[13]:

R = pybamm.Parameter("Particle radius [m]")
D = pybamm.Parameter("Diffusion coefficient [m2.s-1]")
j = pybamm.Parameter("Interfacial current density [A.m-2]")
F = pybamm.Parameter("Faraday constant [C.mol-1]")
c0 = pybamm.Parameter("Initial concentration [mol.m-3]")

The reduced order model solves and ODE for the (uniform) concentration in the particle. In the parameter regime where this is valid, we expect that the solution of the ODE model should agree with the average concentration in the PDE mode. In anticipation of this, we create two variables: the concentration (which we will use in the PDE model), and the average concentration (which we will use in the ODE model). This will make it straightforward to compare the results in a consistent way. Note that the average concentration doesn’t have a domain since it is a scalar quantity.

[14]:

c = pybamm.Variable("Concentration [mol.m-3]", domain="negative particle")
c_av = pybamm.Variable("Average concentration [mol.m-3]")

Now we define our model equations, initial and boundary conditions (where appropriate)

[15]:

# governing equations for full model
N = -D * pybamm.grad(c)  # flux
dcdt = -pybamm.div(N)
full_model.rhs = {c: dcdt}

# governing equations for reduced model
dc_avdt = -3 * j / R / F
reduced_model.rhs = {c_av: dc_avdt}

# initial conditions (these are the same for both models)
full_model.initial_conditions = {c: c0}
reduced_model.initial_conditions = {c_av: c0}

# boundary conditions (only required for full model)
lbc = pybamm.Scalar(0)
rbc = -j / F / D
full_model.boundary_conditions = {
    c: {"left": (lbc, "Neumann"), "right": (rbc, "Neumann")}
}

We can now populate the variables dictionary of both models with any variables of interest. We can compute the average concentration in the full model using the operator pybamm.r_average. We may also wish to compare the concentration profile predicted by the full model with the uniform concentration profile predicted by the reduced model. We can use the operator pybamm.PrimaryBroadcast to broadcast the scalar valued uniform concentration across the particle domain so that it can be visualised as a function of $r$.

Note: the “Primary” refers to the fact the we are broadcasting in only one dimension. For some models, such as the DFN, variables may depend on a “pseudo-dimension” (e.g. the position in $x$ across the cell), but spatial operators only act in the “primary dimension” (e.g. the position in $r$ within the particle). If you are unfamiliar with battery models, do not worry, the details of this are not important for this example. For more information see the broadcasts notebook.

[16]:

# full model
full_model.variables = {
    "Concentration [mol.m-3]": c,
    "Surface concentration [mol.m-3]": pybamm.surf(c),
    "Average concentration [mol.m-3]": pybamm.r_average(c),
}

# reduced model
reduced_model.variables = {
    "Concentration [mol.m-3]": pybamm.PrimaryBroadcast(c_av, "negative particle"),
    "Surface concentration [mol.m-3]": c_av,  # in this model the surface concentration is just equal to the scalar average concentration
    "Average concentration [mol.m-3]": c_av,
}

Using the model#

As before, we provide our parameter values

[17]:

param = pybamm.ParameterValues(
    {
        "Particle radius [m]": 10e-6,
        "Diffusion coefficient [m2.s-1]": 3.9e-14,
        "Interfacial current density [A.m-2]": 1.4,
        "Faraday constant [C.mol-1]": 96485,
        "Initial concentration [mol.m-3]": 2.5e4,
    }
)

We then define and process our geometry, and process both of the models

[18]:

# geometry
r = pybamm.SpatialVariable(
    "r", domain=["negative particle"], coord_sys="spherical polar"
)
geometry = {"negative particle": {r: {"min": pybamm.Scalar(0), "max": R}}}
param.process_geometry(geometry)

# models
for model in models:
    param.process_model(model)

We can now set up our mesh, choose a spatial method, and discretise our models. Note that, even though the reduced-order model is an ODE model, we discretise using the mesh for the particle so that our PrimaryBroadcast operator is discretised in the correct way.

[19]:

# mesh
submesh_types = {"negative particle": pybamm.Uniform1DSubMesh}
var_pts = {r: 20}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

# discretisation
spatial_methods = {"negative particle": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)

# process models
for model in models:
    disc.process_model(model)

Both models are now discretised and ready to be solved.

Solving the model#

As before, we choose a solver and times at which we want the solution returned. We then solve both models, post-process the results, and create a slider plot to compare the results.

[20]:

# loop over models to solve
t = np.linspace(0, 3600, 600)
solutions = [None] * len(models)  # create list to hold solutions
for i, model in enumerate(models):
    solver = pybamm.ScipySolver()
    solutions[i] = solver.solve(model, t)

# post-process the solution of the full model
c_full = solutions[0]["Concentration [mol.m-3]"]
c_av_full = solutions[0]["Average concentration [mol.m-3]"]


# post-process the solution of the reduced model
c_reduced = solutions[1]["Concentration [mol.m-3]"]
c_av_reduced = solutions[1]["Average concentration [mol.m-3]"]

# plot
r = mesh["negative particle"].nodes  # radial position


def plot(t):
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))

    # Plot concetration as a function of r
    ax1.plot(r * 1e6, c_full(t=t, r=r), label="Full Model")
    ax1.plot(r * 1e6, c_reduced(t=t, r=r), label="Reduced Model")
    ax1.set_xlabel("Particle radius [microns]")
    ax1.set_ylabel("Concentration [mol.m-3]")
    ax1.legend()

    # Plot average concentration over time
    t_hour = np.linspace(0, 3600, 600)  # plot over full hour
    c_min = c_av_reduced(t=3600) * 0.98  # minimum axes limit
    c_max = param["Initial concentration [mol.m-3]"] * 1.02  # maximum axes limit

    ax2.plot(t_hour, c_av_full(t=t_hour), label="Full Model")
    ax2.plot(t_hour, c_av_reduced(t=t_hour), label="Reduced Model")
    ax2.plot([t, t], [c_min, c_max], "k--")  # plot line to track time
    ax2.set_xlabel("Time [s]")
    ax2.set_ylabel("Average concentration [mol.m-3]")
    ax2.legend()

    plt.tight_layout()
    plt.show()


import ipywidgets as widgets

widgets.interact(plot, t=widgets.FloatSlider(min=0, max=3600, step=1, value=0));

2022-12-12 12:41:59.589 - [WARNING] processed_variable.get_spatial_scale(518): No length scale set for negative particle. Using default of 1 [m].
2022-12-12 12:41:59.609 - [WARNING] processed_variable.get_spatial_scale(518): No length scale set for negative particle. Using default of 1 [m].

From the results we observe that the reduced-order model does a good job of predicting the average concentration, but, since it is only an ODE model, cannot predicted the spatial distribution.

In the next notebook we will show how to set up and solve a model which contains multiple domains.

References#

The relevant papers for this notebook are:

[21]:

pybamm.print_citations()

[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.
[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.
[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.
[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.

Symbol	Units	Value
\(R\)	m	\(10 \times 10^{-6}\)
\(D\)	m\({^2}\) s\(^{-1}\)	\(3.9 \times 10^{-14}\)
\(j\)	A m\(^{-2}\)	\(1.4\)
\(F\)	C mol\(^{-1}\)	\(96485\)
\(c_0\)	mol m\(^{-3}\)	\(2.5 \times 10^{4}\)