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.

Tutorial 10 - Creating a model

In Tutorial 9 we showed how to change the mesh using on of the built-in battery models in PyBaMM. In this tutorial we show how to create a simple model from scratch in PyBaMM.

As simple example, we consider the problem of linear diffusion on a unit sphere with a flux at the boundary that depends on the concentration. We solve

\[\frac{\partial c}{\partial t} = \nabla \cdot (\nabla c),\]

with the following boundary and initial conditions:

\[\left.\frac{\partial c}{\partial r}\right\vert_{r=0} = 0, \quad \left.\frac{\partial c}{\partial r}\right\vert_{r=1} = -j, \quad \left.c\right\vert_{t=0} = c_0,\]

where

\[j = \left.j_0(1-c)^{1/2}c^{1/2}\right\vert_{r=1}\]

Here \(c_0\) and \(j_0\) are parameters we can control. In this example we will assume that everything is non-dimensional and focus on how to set up and solve the model rather than any specific physical interpretation.

%pip install "pybamm[plot,cite]" -q    # install PyBaMM if it is not installed
import pybamm

Note: you may need to restart the kernel to use updated packages.

Setting up the model

First we load an empty model. We use the BaseModel class that sets up all the basic framework on which our model will be built.

model = pybamm.BaseModel()

We then define our variables and parameters using the Variable and Parameter classes, respectively. Since we are solving a PDE we need to tell PyBaMM the domain each variable belongs to so that it can be discretised in space in the correct way. This is done by passing the keyword argument domain, and in this example we arbitrarily choose the domain “negative particle”.

c = pybamm.Variable("Concentration", domain="negative particle")
c0 = pybamm.Parameter("Initial concentration")
j0 = pybamm.Parameter("Flux parameter")

We then state out governing equations. In PyBaMM we distinguish between Ordinary Differential Equations of the form \(dy/dt = \text{rhs}\) and Algebraic Equations of the form \(f(y) = 0\). The model equations are stored in dictionaries where the key is the variable and the value is the rhs for ODEs and the residual (\(f(y)\)) for algebraic equations.

Sometime it is useful to define intermediate quantities in order to express the governing equations more easily. In this example we define the flux, then define the rhs to be minus the divergence of the flux. The equation is then added to the dictionary model.rhs

N = -pybamm.grad(c)  # define the flux
dcdt = -pybamm.div(N)  # define the rhs equation

model.rhs = {c: dcdt}  # add the equation to rhs dictionary with the variable as the key

Next we add the necessary boundary and initial conditions to the model. These are also stored in dictionaries called model.boundary_conditions and model.initial_conditions, respectively.

# boundary conditions
c_surf = pybamm.surf(c)  # concentration at the surface of the sphere
j = j0 * (1 - c_surf) ** (1 / 2) * c_surf ** (1 / 2)  # prescribed boundary flux
model.boundary_conditions = {c: {"left": (0, "Neumann"), "right": (-j, "Neumann")}}

# initial conditions
model.initial_conditions = {c: c0}

We can add any variables of interest to the dictionary model.variables. These can simply be the variables we solve for (in this case \(c\)) or any other user-defined quantities.

model.variables = {
    "Concentration": c,
    "Surface concentration": c_surf,
    "Flux": N,
    "Boundary flux": j,
}

Setting up the geometry and mesh

In order to solve the model we need to define the geometry and choose how we are going to discretise the equations in space. We first define our radial coordinate using pybamm.SpatialVariable. When we define our spatial variable we pass in a name, the domain on which the variable lives, and the coordinate system we want to use.

r = pybamm.SpatialVariable(
    "r", domain=["negative particle"], coord_sys="spherical polar"
)

We can then define our geometry using a dictionary. The key is the name of the domain, and the value is another dictionary which gives the coordinate to use and the limits. In this case we solve on a unit sphere, so we pass out SpatialVariable, r, and the limit 0 and 1.

geometry = {"negative particle": {r: {"min": 0, "max": 1}}}

Finally we choose how we are going to discretise in space. We choose to use the Finite Volume method on a uniform mesh with 20 volumes.

spatial_methods = {"negative particle": pybamm.FiniteVolume()}
submesh_types = {"negative particle": pybamm.Uniform1DSubMesh}
var_pts = {r: 20}
# create a mesh of our geometry, using a uniform grid with 20 volumes
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

Solving the model

Now we are ready to solve the model. First we need to provide values for the parameters in our model. We do this by passing a dictionary of parameter names and values to the pybamm.ParameterValues class.

parameter_values = pybamm.ParameterValues(
    {
        "Initial concentration": 0.9,
        "Flux parameter": 0.8,
    }
)

Next we choose a solver. Since this is a system of ODEs we can use the ScipySolver which uses a Runge-Kutta scheme by default.

solver = pybamm.ScipySolver()

We can then create a simulation by passing information about the model, geometry, parameters, discretisation and solver to the pybamm.Simulation class.

sim = pybamm.Simulation(
    model,
    geometry=geometry,
    parameter_values=parameter_values,
    submesh_types=submesh_types,
    var_pts=var_pts,
    spatial_methods=spatial_methods,
    solver=solver,
)

Finally we can solve the model

sim.solve([0, 1])  # solve up to a time of 1

<pybamm.solvers.solution.Solution at 0x7f41afa86b20>

The easiest way to quickly plot the results is to call sim.plot to create a slider plot.

Note that at present the plot method is set up to plot dimensional results from battery simulations, so the labels include units which can be ignored (the model assumes a default length scale of 1m and default time scale of 1s). Alternatively we could extract the solution data as seen in Tutorial 6 and create the plots manually. You can find out more about customising plots in this notebook.

# pass in a list of the variables we want to plot
sim.plot(["Concentration", "Surface concentration", "Flux", "Boundary flux"])

2022-07-11 13:21:23.957 - [WARNING] processed_variable.get_spatial_scale(521): No length scale set for negative particle. Using default of 1 [m].
2022-07-11 13:21:23.975 - [WARNING] processed_variable.get_spatial_scale(521): No length scale set for negative particle. Using default of 1 [m].

<pybamm.plotting.quick_plot.QuickPlot at 0x7f41afa9a9d0>

Here we have seen how to create a basic model from scratch in PyBaMM. In the next tutorial we will see how to split this model up into separate submodels.

References

The relevant papers for this notebook are:

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.