Issue 569 create model notebooks #877
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,383 @@ | ||
| { | ||
| "cells": [ | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "# Comparing full and reduced-order models" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "In the [previous notebook](./3-negative-particle-problem.ipynb) 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. \n", | ||
| "\n", | ||
| "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 \n", | ||
| "\n", | ||
| "\\begin{equation*}\n", | ||
| " \\frac{\\textrm{d} c}{\\textrm{d} t} = -3\\frac{j}{RF}\n", | ||
| "\\end{equation*}\n", | ||
| "with the initial condition:\n", | ||
| "\\begin{equation*}\n", | ||
| "\\left.c\\right\\vert_{t=0} = c_0,\n", | ||
| "\\end{equation*}\n", | ||
| "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. \n", | ||
| "\n", | ||
| "As in the previous example we use the following parameters:\n", | ||
| "\n", | ||
| "| Symbol | Units | Value |\n", | ||
| "|:-------|:-------------------|:-----------------------------------------------|\n", | ||
| "| $R$ | m | $10 \\times 10^{-6}$ |\n", | ||
| "| $D$ | m${^2}$ s$^{-1}$ | $3.9 \\times 10^{-14}$ |\n", | ||
| "| $j$ | A m$^{-2}$ | $1.4$ |\n", | ||
| "| $F$ | C mol$^{-1}$ | $96485$ |\n", | ||
| "| $c_0$ | mol m$^{-3}$ | $2.5 \\times 10^{4}$ |" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "## Setting up the models\n", | ||
| "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. " | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 1, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "import pybamm\n", | ||
| "import numpy as np\n", | ||
| "import matplotlib.pyplot as plt\n", | ||
| "\n", | ||
| "full_model = pybamm.BaseModel(name=\"full model\")\n", | ||
| "reduced_model = pybamm.BaseModel(name=\"reduced model\")" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "It can be useful to add the models to a list so that we can perform the same operations on each model easily" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 2, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "models = [full_model, reduced_model]" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We then define our parameters, as seen previously, " | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 3, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "R = pybamm.Parameter(\"Particle radius [m]\")\n", | ||
| "D = pybamm.Parameter(\"Diffusion coefficient [m2.s-1]\")\n", | ||
| "j = pybamm.Parameter(\"Interfacial current density [A.m-2]\")\n", | ||
| "F = pybamm.Parameter(\"Faraday constant [C.mol-1]\")\n", | ||
| "c0 = pybamm.Parameter(\"Initial concentration [mol.m-3]\")" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "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." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 4, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "c = pybamm.Variable(\"Concentration [mol.m-3]\", domain=\"negative particle\")\n", | ||
| "c_av = pybamm.Variable(\"Average concentration [mol.m-3]\")" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "Now we define our model equations, initial and boundary conditions (where appropriate)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 5, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# governing equations for full model\n", | ||
| "N = -D * pybamm.grad(c) # flux\n", | ||
| "dcdt = -pybamm.div(N)\n", | ||
| "full_model.rhs = {c: dcdt} \n", | ||
| "\n", | ||
| "# governing equations for reduced model\n", | ||
| "dc_avdt = -3 * j / R / F\n", | ||
| "reduced_model.rhs = {c_av: dc_avdt} \n", | ||
| "\n", | ||
| "# initial conditions (these are the same for both models)\n", | ||
| "full_model.initial_conditions = {c: c0}\n", | ||
| "reduced_model.initial_conditions = {c_av: c0}\n", | ||
| "\n", | ||
| "# boundary conditions (only required for full model)\n", | ||
| "lbc = pybamm.Scalar(0)\n", | ||
| "rbc = -j / F / D\n", | ||
| "full_model.boundary_conditions = {c: {\"left\": (lbc, \"Dirichlet\"), \"right\": (rbc, \"Neumann\")}}" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "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$. \n", | ||
| "\n", | ||
| "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 ion $r$ within the particle). If you are unfamiliar with battery models, do not worry, the details of this are not important for this example. " | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 6, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# full model\n", | ||
| "full_model.variables = {\n", | ||
| " \"Concentration [mol.m-3]\": c,\n", | ||
| " \"Surface concentration [mol.m-3]\": pybamm.surf(c),\n", | ||
| " \"Average concentration [mol.m-3]\": pybamm.r_average(c),\n", | ||
| "}\n", | ||
| "\n", | ||
| "# reduced model\n", | ||
| "reduced_model.variables = {\n", | ||
| " \"Concentration [mol.m-3]\": pybamm.PrimaryBroadcast(c_av, \"negative particle\"),\n", | ||
| " \"Surface concentration [mol.m-3]\": c_av, # in this model the surface concentration is just equal to the scalar average concentration \n", | ||
| " \"Average concentration [mol.m-3]\": c_av,\n", | ||
| "}" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "## Using the model" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "As before, we provide out parameter values" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 7, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "param = pybamm.ParameterValues(\n", | ||
| " {\n", | ||
| " \"Particle radius [m]\": 10e-6,\n", | ||
| " \"Diffusion coefficient [m2.s-1]\": 3.9e-14,\n", | ||
| " \"Interfacial current density [A.m-2]\": 1.4,\n", | ||
| " \"Faraday constant [C.mol-1]\": 96485,\n", | ||
| " \"Initial concentration [mol.m-3]\": 2.5e4,\n", | ||
| " }\n", | ||
| ")" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We then define and process our geometry, and process both of the models" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 8, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# geometry\n", | ||
| "r = pybamm.SpatialVariable(\"r\", domain=[\"negative particle\"], coord_sys=\"spherical polar\")\n", | ||
| "geometry = {\"negative particle\": {\"primary\": {r: {\"min\": pybamm.Scalar(0), \"max\": R}}}}\n", | ||
| "param.process_geometry(geometry)\n", | ||
| "\n", | ||
| "# models\n", | ||
| "for model in models:\n", | ||
| " param.process_model(model)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "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." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 9, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# mesh\n", | ||
| "submesh_types = {\"negative particle\": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh)}\n", | ||
| "var_pts = {r: 20}\n", | ||
| "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n", | ||
| "\n", | ||
| "# discretisation\n", | ||
| "spatial_methods = {\"negative particle\": pybamm.FiniteVolume()}\n", | ||
| "disc = pybamm.Discretisation(mesh, spatial_methods)\n", | ||
| "\n", | ||
| "# process models\n", | ||
| "for model in models:\n", | ||
| " disc.process_model(model)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "Both models are now discretised and ready to be solved." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "### Solving the model" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "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." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 23, | ||
| "metadata": {}, | ||
| "outputs": [ | ||
| { | ||
| "data": { | ||
| "application/vnd.jupyter.widget-view+json": { | ||
| "model_id": "b6a9b1125ea8432cacf0081560824fb4", | ||
| "version_major": 2, | ||
| "version_minor": 0 | ||
| }, | ||
| "text/plain": [ | ||
| "interactive(children=(FloatSlider(value=0.0, description='t', max=3600.0, step=1.0), Output()), _dom_classes=(…" | ||
| ] | ||
| }, | ||
| "metadata": {}, | ||
| "output_type": "display_data" | ||
| } | ||
| ], | ||
| "source": [ | ||
| "# loop over models to solve\n", | ||
| "solver = pybamm.ScipySolver()\n", | ||
| "t = np.linspace(0, 3600, 600)\n", | ||
| "solutions = [None] * len(models) # create list to hold solutions\n", | ||
| "for i, model in enumerate(models):\n", | ||
| " solutions[i] = solver.solve(model, t)\n", | ||
| "\n", | ||
| "# post-process the solution of the full model\n", | ||
| "c_full = solutions[0][\"Concentration [mol.m-3]\"]\n", | ||
| "c_av_full = solutions[0][\"Average concentration [mol.m-3]\"]\n", | ||
| "\n", | ||
| "\n", | ||
| "# post-process the solution of the full model\n", | ||
| "c_reduced = solutions[1][\"Concentration [mol.m-3]\"]\n", | ||
| "c_av_reduced = solutions[1][\"Average concentration [mol.m-3]\"]\n", | ||
| "\n", | ||
| "# plot\n", | ||
| "r = mesh[\"negative particle\"][0].nodes # radial position\n", | ||
| "\n", | ||
| "def plot(t):\n", | ||
| " fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))\n", | ||
| " \n", | ||
| " # Plot concetration as a function of r\n", | ||
| " ax1.plot(r * 1e6, c_full(t=t,r=r), label=\"Full Model\")\n", | ||
| " ax1.plot(r * 1e6, c_reduced(t=t,r=r), label=\"Reduced Model\") \n", | ||
| " ax1.set_xlabel(\"Particle radius [microns]\")\n", | ||
| " ax1.set_ylabel(\"Concentration [mol.m-3]\")\n", | ||
| " ax1.legend()\n", | ||
| " \n", | ||
| " # Plot average concentration over time\n", | ||
| " t_hour = np.linspace(0, 3600, 600) # plot over full hour\n", | ||
| " c_min = c_av_reduced(t=3600) * 0.98 # minimum axes limit \n", | ||
| " c_max = param[\"Initial concentration [mol.m-3]\"] * 1.02 # maximum axes limit \n", | ||
| " \n", | ||
| " ax2.plot(t_hour, c_av_full(t=t_hour), label=\"Full Model\")\n", | ||
| " ax2.plot(t_hour, c_av_reduced(t=t_hour), label=\"Reduced Model\") \n", | ||
| " ax2.plot([t, t], [c_min, c_max], \"k--\") # plot line to track time\n", | ||
| " ax2.set_xlabel(\"Time [s]\")\n", | ||
| " ax2.set_ylabel(\"Average concentration [mol.m-3]\") \n", | ||
| " ax2.legend()\n", | ||
| "\n", | ||
| " plt.tight_layout()\n", | ||
| " plt.show()\n", | ||
| " \n", | ||
| "import ipywidgets as widgets\n", | ||
| "widgets.interact(plot, t=widgets.FloatSlider(min=0,max=3600,step=1,value=0));\n", | ||
| " " | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "From the results we pbserve 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 dis" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "In the [next notebook](/5-a-simple-SEI-model.ipynb) we consider a simple model for SEI growth, and show how to correctly pose the model in non-dimensional form and then create and solve it using pybamm." | ||
|
There was a problem hiding this comment. this link is also broken |
||
| ] | ||
| } | ||
| ], | ||
| "metadata": { | ||
| "kernelspec": { | ||
| "display_name": "Python 3", | ||
| "language": "python", | ||
| "name": "python3" | ||
| }, | ||
| "language_info": { | ||
| "codemirror_mode": { | ||
| "name": "ipython", | ||
| "version": 3 | ||
| }, | ||
| "file_extension": ".py", | ||
| "mimetype": "text/x-python", | ||
| "name": "python", | ||
| "nbconvert_exporter": "python", | ||
| "pygments_lexer": "ipython3", | ||
| "version": "3.6.9" | ||
| } | ||
| }, | ||
| "nbformat": 4, | ||
| "nbformat_minor": 2 | ||
| } | ||
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
link to the broadcasts notebook (in the expression tree folder)?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think there is a typo where it says "(e.g. the position ion$r$ ...".