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.

Creating a submodel

In Tutorial 3 we showed how to create a simple model from scratch in PyBaMM. In this tutorial we will solve a very similar problem, but using separate submodels for the linear diffusion problem and the model for the surface flux. In this simple example the surface flux is just some known function of the concentration, so we could just explicitly define it in the model for diffusion. However, we write it as a separate model to show how submodels interact.

We solved the problem of linear diffusion on a unit sphere with a flux at the boundary that depends on the concentration

\[\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. Again 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

Setting up the model

Again we start with an empty BaseModel class

model = pybamm.BaseModel()

Next we set up the submodel for diffusion in the particle using a pybamm.BaseSubModel. Each submodel has methods that define the variables, equations, initial and boundary conditions corresponding to the physics the model describes.

First get_fundamental_variables defines any variables that can be defined independently of any other submodels that may be included in our model. Here we can define the concentration, surface concentration and flux, and add them to the dictionary of variables.

Next we can use get_coupled_variables to define any variables that do depend on variables defined in another submodel. In this simple example we don’t have any variables to define here. However, if we had included a temperature dependent diffusivity, for example, then we would have needed to define the flux in get_coupled_variables since it would now depend on the temperature which would be defined in a separate submodel.

Once we have defined all the variables we need we can write down our equations. Any equations that include time derivatives will turn into ordinary differential equations after discretisation. We set the right hand sides of such equations in the set_rhs method. In this example we add the right hand side of the diffusion equation to the rhs dictionary. Equations that don’t contain time derivatives give algebraic constraints in our model. These equations are set in the set_algebraic method. In this example we don’t have any algebraic equations, so we can skip this method.

Finally we set the boundary and initial conditions using the methods set_boundary_conditions and set_initial_conditions, respectively.

class Particle(pybamm.BaseSubModel):
    def __init__(self, param, domain, options=None):
        super().__init__(param, domain, options=options)

    def get_fundamental_variables(self):
        # create concentration variable
        c = pybamm.Variable("Concentration", domain="negative particle")

        # define concentration at the surface of the sphere
        c_surf = pybamm.surf(c)

        # define flux
        N = -pybamm.grad(c)

        # create dictionary of model variables
        variables = {
            "Concentration": c,
            "Surface concentration": c_surf,
            "Flux": N,
        }

        return variables

    def get_coupled_variables(self, variables):
        return variables

    def set_rhs(self, variables):
        # extract the variables we need
        c = variables["Concentration"]
        N = variables["Flux"]

        # define the rhs of the PDE
        dcdt = -pybamm.div(N)

        # add it to the submodel dictionary
        self.rhs = {c: dcdt}

    def set_algebraic(self, variables):
        pass

    def set_boundary_conditions(self, variables):
        # extract the variables we need
        c = variables["Concentration"]
        j = variables["Boundary flux"]

        # add the boundary conditions to the submodel dictionary
        self.boundary_conditions = {
            c: {"left": (0, "Neumann"), "right": (-j, "Neumann")}
        }

    def set_initial_conditions(self, variables):
        # extract the variable we need
        c = variables["Concentration"]

        # define the initial concentration parameter
        c0 = pybamm.Parameter("Initial concentration")

        # add the initial conditions to the submodel dictionary
        self.initial_conditions = {c: c0}

class BoundaryFlux(pybamm.BaseSubModel):
    def __init__(self, param, domain, options=None):
        super().__init__(param, domain, options=options)

    def get_coupled_variables(self, variables):
        # extract the variable we need
        c_surf = variables["Surface concentration"]

        # define the flux parameter
        j0 = pybamm.Parameter("Flux parameter")
        j = j0 * (1 - c_surf) ** (1 / 2) * c_surf ** (1 / 2)  # prescribed boundary flux

        # update dictionary of model variables
        variables.update({"Boundary flux": j})

        return variables

We can now set the submodels in a model by assigning a dictionary to model.submodels. The dictionary key is the name we want to give to the submodel and the value is an instance of the submodel class we want to use.

When we instantiate a submodel we are required to pass in param, a class of parameter symbols we are going to call, and domain, the domain on which the submodel lives. In this example we will simply set param to None and hard-code the definition of our parameters into the submodel. When writing lots of submodels it is more efficient to define all the parameters in a shared class, and pass this to each submodel. For the domain we will choose “Negative”.

model.submodels = {
    "Particle": Particle(None, "Negative"),
    "Boundary flux": BoundaryFlux(None, "Negative"),
}

At this stage we have just told the model which submodels it is constructed from, but the variables and equations have not yet been created. For example if we look at the rhs dictionary it is empty.

model.rhs

{}

To populate the model variables, equations, boundary and initial conditions we need to “build” the model. To do this we call build_model

model.build_model()

This loops through all of the submodels, first creating the “fundamental variables”, followed by the “coupled variables” and finally the equations (rhs and algebraic) and the boundary and initial conditions. Now we see that model.rhs contains our diffusion equation.

model.rhs

{Variable(0x48a91fa41dea72d3, Concentration, children=[], domains={'primary': ['negative particle']}): Divergence(0x4befa5e7cf20d0c5, div, children=['grad(Concentration)'], domains={'primary': ['negative particle']})}

Using the model

We can now use our model as in the previous tutorial.

We first set up our geometry and mesh

r = pybamm.SpatialVariable(
    "r", domain=["negative particle"], coord_sys="spherical polar"
)
geometry = {"negative particle": {r: {"min": 0, "max": 1}}}
spatial_methods = {"negative particle": pybamm.FiniteVolume()}
submesh_types = {"negative particle": pybamm.Uniform1DSubMesh}
var_pts = {r: 20}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

and then set up our simulation, remembering to set values for our parameters

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

solver = pybamm.ScipySolver()

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

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

and plot the results

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

2022-07-11 13:20:04.194 - [WARNING] processed_variable.get_spatial_scale(521): No length scale set for negative particle. Using default of 1 [m].
2022-07-11 13:20:04.209 - [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 0x7f9dc1b2ec10>

In this notebook we saw how to split a model up into submodels. Although this was a simple example it let us understand how to construct submodels and see how they interact via coupled variables.

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.