thermal_explosion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from bokeh.io import  output_notebook, push_notebook, show\n",
    "from bokeh.plotting import figure\n",
    "from bokeh.layouts import column, row\n",
    "from bokeh.palettes import Category10, Category20\n",
    "from bokeh.models import PrintfTickFormatter\n",
    "\n",
    "from ipywidgets import FloatProgress, IntProgress\n",
    "from IPython.display import display\n",
    "\n",
    "from mylib.thermal_explosion_model import thermal_explosion_model, fuel_thermal_explosion_model\n",
    "\n",
    "import mylib.integration as integration\n",
    "\n",
    "output_notebook(hide_banner=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Themal explosion "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model with fuel consumption neglected\n",
    "\n",
    "\\begin{equation}\n",
    "\\partial_\\tau \\theta - {\\displaystyle 1 \\over \\displaystyle \\lambda} \\partial_{zz} \\theta = \\exp(\\theta)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integration with Radau5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def thermal_explosion_radau():\n",
    "\n",
    "    lamb = 0.86\n",
    "    xmin = 0.\n",
    "    xmax = 2.\n",
    "    # nb of points including boundary conditions\n",
    "    nxib = 1001\n",
    "    nx = nxib-2\n",
    "    tini = 0.0\n",
    "    tend = 1.0\n",
    "\n",
    "    tem = thermal_explosion_model(lamb=lamb, xmin=xmin, xmax=xmax, nx=nx)\n",
    "    fcn_radau = tem.fcn_radau\n",
    "\n",
    "    dx = (xmax-xmin)/(nxib-1)\n",
    "    x = np.linspace(xmin+dx, xmax-dx, nx)\n",
    "    yini = np.zeros(nx)\n",
    "    fig_sol = figure(x_range=(xmin, xmax), plot_height=500, plot_width=950, \n",
    "                     title=\"Solution (click on legend to hide correponding graph)\")\n",
    "    fig_sol.x(x, yini, color=\"Grey\", legend=\"sol at t = 0\")\n",
    "    fig_sol.legend.click_policy=\"hide\"\n",
    "    \n",
    "    nt_plot = 6\n",
    "    t_plot = np.linspace(tini, tend, nt_plot)\n",
    "    for it in range(nt_plot-1):\n",
    "        sol = integration.radau5(t_plot[it], t_plot[it+1], yini, fcn_radau, njac=1, atol=1.e-6, rtol=1.e-6)\n",
    "        yini = sol.y\n",
    "        fig_sol.x(x, sol.y, legend=f\"sol at t={t_plot[it+1]:.2f}\", color=Category10[10][it])\n",
    "        \n",
    "    show(fig_sol, notebook_handle=True)\n",
    "     \n",
    "thermal_explosion_radau()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integration with Rock4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def thermal_explosion_rock():\n",
    "\n",
    "    lamb = 0.86\n",
    "    xmin = 0.\n",
    "    xmax = 2.\n",
    "    # nb of points including boundary conditions\n",
    "    nxib = 1001\n",
    "    nx = nxib-2\n",
    "    tini = 0.0\n",
    "    tend = 1.0\n",
    "\n",
    "    progress_bar = FloatProgress(min=tini, max=tend, value=0, description='Progress:',)\n",
    "    display(progress_bar)\n",
    "\n",
    "    tem = thermal_explosion_model(lamb=lamb, xmin=xmin, xmax=xmax, nx=nx)\n",
    "    fcn_rock = tem.fcn_rock\n",
    "    \n",
    "    dx = (xmax-xmin)/(nxib-1)\n",
    "    x = np.linspace(xmin+dx, xmax-dx, nx)\n",
    "    yini = np.zeros(nx)\n",
    "    fig_sol = figure(x_range=(xmin, xmax), plot_height=500, plot_width=950, \n",
    "                     title=\"Solution (click on legend to hide correponding graph)\")\n",
    "    fig_sol.x(x, yini, color=\"Grey\", legend=\"sol at t = 0\")\n",
    "    fig_sol.legend.click_policy=\"hide\"\n",
    "\n",
    "    nt_plot = 6\n",
    "    t_plot = np.linspace(tini, tend, nt_plot)\n",
    "    for it in range(nt_plot-1):\n",
    "        sol = integration.rock4(t_plot[it], t_plot[it+1], yini, fcn_rock, tol=1.e-6)\n",
    "        yini = sol.y\n",
    "        fig_sol.x(x, sol.y, legend=f\"sol at t={t_plot[it+1]:.2f}\", color=Category10[10][it])\n",
    "        progress_bar.value = t_plot[it+1]\n",
    "        \n",
    "    show(fig_sol, notebook_handle=True)\n",
    "     \n",
    "thermal_explosion_rock()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model with fuel consumption\n",
    "\n",
    "$$\n",
    "\\left\\{\\begin{aligned}\n",
    "\\partial_\\tau \\theta - {\\displaystyle 1 \\over \\displaystyle \\lambda} \\partial_{zz} \\theta & = \\exp(\\theta)\\, Y\\\\\n",
    "\\partial_\\tau Y - {\\displaystyle 1 \\over \\displaystyle \\lambda} \\partial_{zz} Y & = - \\varepsilon \\exp(\\theta)\\, Y\n",
    "\\end{aligned}\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fuel_thermal_explosion_radau():\n",
    "\n",
    "    lamb = 0.9\n",
    "    eps = 0.002\n",
    "    xmin = 0.\n",
    "    xmax = 2.\n",
    "    # nb of points including boundary conditions\n",
    "    nxib = 1001\n",
    "    nx = nxib-2\n",
    "    tini = 0.0\n",
    "    tend = 10.0\n",
    "    neq = 2\n",
    "\n",
    "    ftem = fuel_thermal_explosion_model(lamb=lamb, eps=eps, xmin=xmin, xmax=xmax, nx=nx)\n",
    "    fcn_radau = ftem.fcn_radau\n",
    "\n",
    "    yini = np.zeros((nx, neq))\n",
    "    for inx in range(nx):\n",
    "        yini[inx, 1] = 0.5+inx/(nx)\n",
    "        #yini[inx, 1] = 1.0\n",
    "    yini = yini.flatten()\n",
    "    \n",
    "    dx = (xmax-xmin)/(nxib-1)\n",
    "    x = np.linspace(xmin, xmax, nx)\n",
    "    fig_sol_01 = figure(x_range=(xmin, xmax), plot_height=450, plot_width=950, \n",
    "                         title=\"Solution (click on legend to hide correponding graph)\")\n",
    "    fig_sol_02 = figure(x_range=(xmin, xmax), plot_height=450, plot_width=950, \n",
    "                         title=\"Solution (click on legend to hide correponding graph)\")\n",
    "    fig_sol_01.x(x, yini.reshape(nx,neq)[:,0], color=\"Grey\", legend=\"sol at t = 0\")\n",
    "    fig_sol_02.x(x, yini.reshape(nx,neq)[:,1], color=\"Grey\", legend=\"sol at t = 0\")\n",
    "    fig_sol_01.legend.click_policy=\"hide\"\n",
    "    fig_sol_02.legend.click_policy=\"hide\"\n",
    "    fig_sol_02.legend.location = \"top_left\"\n",
    "\n",
    "    # radau integration\n",
    "    nt_plot = 6\n",
    "    t_plot = np.linspace(tini, tend, nt_plot)\n",
    "    for it in range(nt_plot-1):\n",
    "        sol = integration.radau5(t_plot[it], t_plot[it+1], yini, fcn_radau, njac=2, atol=1.e-6, rtol=1.e-6)\n",
    "        yini = sol.y\n",
    "        theta = sol.y.reshape(nx,neq)[:,0]\n",
    "        Y = sol.y.reshape(nx,neq)[:,1]\n",
    "        fig_sol_01.x(x, theta, legend=f\"sol at t={t_plot[it+1]:.2f}\", color=Category10[10][it])\n",
    "        fig_sol_02.x(x, Y, legend=f\"sol at t={t_plot[it+1]:.2f}\", color=Category10[10][it])\n",
    "        \n",
    "    show(column(fig_sol_01, fig_sol_02))\n",
    "\n",
    "fuel_thermal_explosion_radau()"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}