docs/source/examples/notebooks/creating_models/4-comparing-full-and-reduced-order-models.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparing full and reduced-order models"
   ]
  },
  {
   "attachments": {},
   "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",
    "$$\n",
    "  \\frac{\\textrm{d} c}{\\textrm{d} t} = -3\\frac{j}{RF}\n",
    "$$\n",
    "with the initial condition:\n",
    "$$\n",
    "\\left.c\\right\\vert_{t=0} = c_0,\n",
    "$$\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}$                            |"
   ]
  },
  {
   "attachments": {},
   "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": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[33mWARNING: You are using pip version 22.0.4; however, version 22.3.1 is available.\n",
      "You should consider upgrading via the '/home/mrobins/git/PyBaMM/env/bin/python -m pip install --upgrade pip' command.\u001b[0m\u001b[33m\n",
      "\u001b[0mNote: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
    "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\")"
   ]
  },
  {
   "attachments": {},
   "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": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "models = [full_model, reduced_model]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then define our parameters, as seen previously, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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]\")"
   ]
  },
  {
   "attachments": {},
   "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": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = pybamm.Variable(\"Concentration [mol.m-3]\", domain=\"negative particle\")\n",
    "c_av = pybamm.Variable(\"Average concentration [mol.m-3]\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define our model equations, initial and boundary conditions (where appropriate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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 = {\n",
    "    c: {\"left\": (lbc, \"Neumann\"), \"right\": (rbc, \"Neumann\")}\n",
    "}"
   ]
  },
  {
   "attachments": {},
   "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 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](../expression_tree/broadcasts.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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",
    "}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we provide our parameter values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "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",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then define and process our geometry, and process both of the models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# geometry\n",
    "r = pybamm.SpatialVariable(\n",
    "    \"r\", domain=[\"negative particle\"], coord_sys=\"spherical polar\"\n",
    ")\n",
    "geometry = {\"negative particle\": {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)"
   ]
  },
  {
   "attachments": {},
   "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": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mesh\n",
    "submesh_types = {\"negative particle\": 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)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both models are now discretised and ready to be solved."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving the model"
   ]
  },
  {
   "attachments": {},
   "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": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "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].\n",
      "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].\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6fa9ebee57924ba5b79a2c51313fba25",
       "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",
    "t = np.linspace(0, 3600, 600)\n",
    "solutions = [None] * len(models)  # create list to hold solutions\n",
    "for i, model in enumerate(models):\n",
    "    solver = pybamm.ScipySolver()\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 reduced 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\"].nodes  # radial position\n",
    "\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",
    "\n",
    "import ipywidgets as widgets\n",
    "\n",
    "widgets.interact(plot, t=widgets.FloatSlider(min=0, max=3600, step=1, value=0));"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [next notebook](./5-half-cell-model.ipynb) we will show how to set up and solve a model which contains multiple domains."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "The relevant papers for this notebook are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[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.\n",
      "[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.\n",
      "[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.\n",
      "[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.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pybamm.print_citations()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env",
   "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.9.15"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "vscode": {
   "interpreter": {
    "hash": "19e5ebaa8d5a3277b4deed2928f02ad0cad6c3ab0b2beced644d557f155bce64"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}