demo/ign_sens_excel.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "adiabatic kinetics simulation with sensitivity analysis\n",
    "\"\"\"\n",
    "import sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import cantera as ct\n",
    "import matplotlib.pyplot as plt\n",
    "import timeit\n",
    "\n",
    "def ign_uq(g, factor, temp, pres, fuel, phi, T_equi, simtype='UV', bootplot=False):\n",
    "\tgas.set_multiplier(1.0) # reset all multipliers\n",
    "\tfor i in range(gas.n_reactions):\n",
    "\t\tgas.set_multiplier(factor[i],i)\n",
    "\t\t# print(gas.reaction_equation(i)+' index_reaction:',i,'multi_factor:',factor[i])\n",
    "\t\n",
    "\tgas.set_equivalence_ratio( phi, fuel, 'O2:1.0, N2:3.76' )\n",
    "\tgas.TP = temp, pres*ct.one_atm\n",
    "\n",
    "\t# here it is constant volume with IdealGasReactor\n",
    "\tif simtype == 'UV':\n",
    "\t\tr = ct.IdealGasReactor(gas)\n",
    "\telse:\n",
    "\t\tr = ct.IdealGasConstPressureReactor(gas)\n",
    "\n",
    "\tsim = ct.ReactorNet([r])\n",
    "\t\n",
    "\t# set the tolerances for the solution and for the sensitivity coefficients\n",
    "\tsim.rtol = 1.0e-6\n",
    "\tsim.atol = 1.0e-15\n",
    "\n",
    "\tt_end = 10;\n",
    "\ttime = []\n",
    "\ttemp = []\n",
    "\tstates = ct.SolutionArray(gas, extra=['t'])\n",
    "\n",
    "\t# stateArray = []\n",
    "\twhile sim.time < t_end and r.T < T_equi - 0.1:\n",
    "\t\tsim.step()\n",
    "\t\ttime.append(sim.time)\n",
    "\t\ttemp.append(r.T)\n",
    "\t\tstates.append(r.thermo.state, t=sim.time)\n",
    "\t\t# stateArray.append(r.thermo.X)\n",
    "\n",
    "\ttime = np.array(time)\n",
    "\ttemp = np.array(temp)\n",
    "\tdiff_temp = np.diff(temp)/np.diff(time)\n",
    "\n",
    "\tif 0:\n",
    "\t\t# print(r.T, T_equi)\n",
    "\t\tplt.figure()\n",
    "\t\tplt.plot(states.t, states.T, '-ok')\n",
    "\n",
    "\t\t# plt.plot( 1.0, states.T[ign_pos], 'rs',markersize=6  )\n",
    "\t\tplt.xlabel('Time [ms]')\n",
    "\t\tplt.ylabel('T [K]')\n",
    "\t\t# plt.xlim( 0.99, 1.01 )\n",
    "\t\t#plt.ylim( -0.001,0.001 )\n",
    "\t\tplt.tight_layout()\n",
    "\t\tplt.show()\n",
    "\n",
    "\tign_pos = np.argmax( diff_temp )\n",
    "\n",
    "\tign = time[ign_pos]\n",
    "\t\n",
    "\tgas.set_multiplier(1.0) # reset all multipliers\n",
    "\n",
    "\treturn ign\n",
    "\n",
    "def ign_sens(g, temp, pres, fuel, phi, ign0, simtype='UV', bootplot=False):\n",
    "\t'''\n",
    "\tign_sens(g, temp, pres, fuel, phi, ign0, simtype='UV', bootplot=False)\n",
    "\t'''\n",
    "\tgas.set_multiplier(1.0) # reset all multipliers\n",
    "\n",
    "\tgas.set_equivalence_ratio( phi, fuel, 'O2:1.0, N2:3.76' )\n",
    "\tgas.TP = temp, pres*ct.one_atm\n",
    "\n",
    "\t# here it is constant volume with IdealGasReactor\n",
    "\tif simtype == 'UV':\n",
    "\t\tr = ct.IdealGasReactor(gas)\n",
    "\telse:\n",
    "\t\tr = ct.IdealGasConstPressureReactor(gas)\n",
    "\n",
    "\tsim = ct.ReactorNet([r])\n",
    "\tfor i in range(gas.n_reactions):\n",
    "\t\tr.add_sensitivity_reaction(i)\n",
    "\n",
    "\t# set the tolerances for the solution and for the sensitivity coefficients\n",
    "\tsim.rtol = 1.0e-6\n",
    "\tsim.atol = 1.0e-15\n",
    "\tsim.rtol_sensitivity = 1.0e-6\n",
    "\tsim.atol_sensitivity = 1.0e-6\n",
    "\tsens_T = []\n",
    "\tstates = ct.SolutionArray(gas, extra=['t'])\n",
    "\n",
    "\tt_end = 10;\n",
    "\ttime = []\n",
    "\ttemp = []\n",
    "\t# stateArray = []\n",
    "\twhile sim.time < ign0*1.5:\n",
    "\t\tsim.step()\n",
    "\t\tsens_T.append( sim.sensitivities()[2,:] )\n",
    "\t\tstates.append(r.thermo.state, t=sim.time)\n",
    "\t\ttime.append( sim.time )\n",
    "\t\ttemp.append( r.T )\n",
    "\n",
    "\ttime = np.array(time)\n",
    "\ttemp = np.array(temp)\n",
    "\tdiff_temp = np.diff(temp)/np.diff(time)\n",
    "\tsens_T = np.array( sens_T )\n",
    "\n",
    "\tign_pos = np.argmax( diff_temp )\n",
    "\tign = time[ign_pos]\n",
    "\t\n",
    "\tif bootplot:\n",
    "\t\tplt.figure()\n",
    "\t\tplt.plot(states.t/ign0, states.T, '-ok')\n",
    "\n",
    "\t\tplt.plot( 1.0, states.T[ign_pos], 'rs',markersize=6  )\n",
    "\t\tplt.xlabel('Time [ms]')\n",
    "\t\tplt.ylabel('T [K]')\n",
    "\t\tplt.xlim( 0.99, 1.01 )\n",
    "\t\t#plt.ylim( -0.001,0.001 )\n",
    "\t\tplt.tight_layout()\n",
    "\t\tplt.show()\n",
    "\n",
    "\tsens_T_ign = sens_T[ign_pos, :]\n",
    "\treturn sens_T_ign"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CH4_gri30.xml_0.5_20_1000.0_0.05_UV\n",
      "Ignition Delay is: 34.8857 ms\n",
      "Start Adjoint\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH/RJREFUeJzt3XuQXOWd3vHvo5E0uq0MaMZY6Dayg9c7wldNFG1c3tgr\nymi93hWVVDbaErdgIMywXvClvMZUYW+qYFnwJaZiCQioQKiDrMVgyJbZ5RKXSSqR5RE3XUCWDMhI\nCJAlOxKFEZb45Y/zDjStmenRTJ/u063nU3VKZ95z6V83Z+bhPeftcxQRmJmZFc24RhdgZmY2GAeU\nmZkVkgPKzMwKyQFlZmaF5IAyM7NCckCZmVkh5RZQkuZI+rGkbZK2Srq8YvmXJIWkjrK2KyXtlLRd\n0lll7QslbU7LbpSkvOo2M7NiyLMHdQT4UkR0A4uByyR1QxZewKeBXw6snJYtBxYAS4GVktrS4lXA\nxcDpaVqaY91mZlYAuQVUROyNiMfS/CHgaWBWWvwd4CtA+beElwHrIuJwRDwH7AQWSZoJTI+IDZF9\nq3gNcHZedZuZWTGMr8eLSOoCPgr8VNIyYE9EPFlxpm4WsKHs592p7XdpvrJ9sNe5BLgEYOrUqQs/\n8IEP1OgdmJlZpU2bNv0qIjrz2n/uASVpGvAD4Aqy035fIzu9V3MRcQtwC0BPT0/09/fn8TJmZgZI\n2pXn/nMdxSdpAlk4lSLiHuB9wHzgSUnPA7OBxyS9B9gDzCnbfHZq25PmK9vNzKyF5TmKT8BtwNMR\n8W2AiNgcEe+OiK6I6CI7XfexiHgJuB9YLqld0nyywRAbI2IvcFDS4rTP84D78qrbzMyKIc9TfB8H\nzgU2S3oitX0tIn402MoRsVXSemAb2anAyyLiaFrcB9wOTAYeSJOZmbUwterjNnwNyswsX5I2RURP\nXvv3nSTMzKyQHFBmBVcqlWhvb0dSXae2tjYkMX78eCTR0dFBR0fHqNq6urro6+ujq6uLcePG0dXV\nRalUeuv9DdZuRkS05LRw4cIwa1ZkX2Jv6WnKlCnR29sbU6ZMeUf7hAkTYsaMGSEp5s2bF729vTFv\n3rwhf167dm2sXbv2mLYBgy1rRNtgdVca7n0UcR2gP3L8O97wIMlrckBZs2p0cNRzGjdu3Jj30dbW\nFhMmTHhHW3t7e1x99dVx9dVXR3t7e9X1x48fn3tb5TRp0qS45ppr4vHHH4/HH388rr322pg0adIx\n61x77bXx5JNPxlNPPRV/93d/N+g61113XWzZsiW2bNkS11133aDrXH/99bFt27bYtm1bXH/99YOu\nc8MNN8T27dtj+/bt8c1vfnPQdb71rW/Fjh07YseOHQE8Gzn+HfcgCbMC8X2Qrcm8ERHtee3c16DM\nCmLWrEHv4NXSxo3L70/QwLW0IpPEvffeyz333DNkrZK4++67ufvuu4ddZ/369axfv37YddatW8dd\nd901bE2lUqnqdcA777yTNWvWAEwcdsWxyrN71sjJp/is2VCAU271nIa6BlWrad68eTFv3ryGv89q\nNQ4YqtYirwMcjhz/jjc8SPKaHFDWbBr9x7JyGrg+1NbWFkDMmDEjZsyYMaq24QYJlF+InzFjRkyc\nOPG46pwwYcIx20yZMuWtgQqDDcKoXL8ebZXTQI0DBqu16OuQ8zWohgdJXpMDyppJb2/viP8gn3ba\naY0uN1fVRr95FF9x1iHnUXweJGFWANWulfT29rJy5co6VWM2Mr6ThJk5nOyE5IAyM7NCckCZFdyS\nJUsaXYJZQzigzAru4YcfbnQJZg3hgDIzs0JyQJmZWSE5oMzMrJAcUGYNdiLeg89sJBxQZg324osv\nNroEs0JyQJmZWSE5oMwKrLu7u9ElmDWMA8qswLZu3droEswaxgFlZmaFlFtASZoj6ceStknaKuny\n1H6DpGckPSXpXkknlW1zpaSdkrZLOqusfaGkzWnZjSr6YzLNzGzM8uxBHQG+FBHdwGLgMkndwEPA\nGRHxIeDnwJUAadlyYAGwFFgpqS3taxVwMXB6mpbmWLeZmRVAbgEVEXsj4rE0fwh4GpgVEQ9GxJG0\n2gZgdppfBqyLiMMR8RywE1gkaSYwPSI2pAdkrQHOzqtuMzMrhrpcg5LUBXwU+GnFoguBB9L8LOCF\nsmW7U9usNF/ZPtjrXCKpX1L/vn37xl64mZk1TO4BJWka8APgiog4WNZ+FdlpwFKtXisibomInojo\n6ezsrNVuzcysAcbnuXNJE8jCqRQR95S1XwB8FlgSbz9zfg8wp2zz2altD2+fBixvNzOzFpbnKD4B\ntwFPR8S3y9qXAl8B/jwiXivb5H5guaR2SfPJBkNsjIi9wEFJi9M+zwPuy6tuMzMrhjx7UB8HzgU2\nS3oitX0NuBFoBx5Ko8U3RMSlEbFV0npgG9mpv8si4mjarg+4HZhMds3qAczMrKXp7TNsraWnpyf6\n+/sbXYZZVcN9ra9Vfz+tNUjaFBE9ee3fd5IwM7NCckCZmVkhOaDMzKyQHFBmZlZIDigzMyskB5SZ\nmRWSA8qsgUqlmt3py6zlOKDMGuicc85pdAlmheWAMjOzQnJAmZlZITmgzApq8uTJjS7BrKEcUGYF\n9dprr1VfyayFOaDMzKyQHFBmZlZIDigzMyskB5SZmRWSA8rMzArJAWVmZoXkgDIzs0JyQJmZWSE5\noMzMrJAcUGZmVki5BZSkOZJ+LGmbpK2SLk/tp0h6SNKO9O/JZdtcKWmnpO2SziprXyhpc1p2oyTl\nVbeZmRVDnj2oI8CXIqIbWAxcJqkb+CrwSEScDjySfiYtWw4sAJYCKyW1pX2tAi4GTk/T0hzrNjOz\nAsgtoCJib0Q8luYPAU8Ds4BlwB1ptTuAs9P8MmBdRByOiOeAncAiSTOB6RGxISICWFO2jZmZtai6\nXIOS1AV8FPgpcGpE7E2LXgJOTfOzgBfKNtud2mal+cp2MzNrYbkHlKRpwA+AKyLiYPmy1COKGr7W\nJZL6JfXv27evVrs1M7MGyDWgJE0gC6dSRNyTml9Op+1I/76S2vcAc8o2n53a9qT5yvZjRMQtEdET\nET2dnZ21eyNmZlZ3eY7iE3Ab8HREfLts0f3A+Wn+fOC+svblktolzScbDLExnQ48KGlx2ud5ZduY\nmVmLGp/jvj8OnAtslvREavsacB2wXtLngF3AXwBExFZJ64FtZCMAL4uIo2m7PuB2YDLwQJrMzKyF\nKbsM1Hp6enqiv7+/0WWYDWu4r/S16u+mtQ5JmyKiJ6/9+04SZg20l2yUUOW0d7iNzE4QDiizBnrP\ncbabnUgcUGZmVkgOKDMzKyQHlJmZFZIDyszMCskBZdZALx1nu9mJJM8v6ppZFTOHWeZvQdmJzj0o\nMzMrJAeUmZkVkgPKzMwKyQFlZmaF5IAyM7NCckCZmVkhOaDMzKyQHFBmZlZIDigzMyskB5SZmRWS\nA8rMzArJAWXWIKVSqdElmBWaA8qsQS666KJGl2BWaA4oswZ5/fXXG12CWaE5oMwKaOrUqY0uwazh\ncgsoSaslvSJpS1nbRyRtkPSEpH5Ji8qWXSlpp6Ttks4qa18oaXNadqMk5VWzWVHcfPPNjS7BrOHy\n7EHdDiytaLse+NuI+AhwdfoZSd3AcmBB2malpLa0zSrgYuD0NFXu06zlrFixotElmDVcbgEVEY8C\nByqbgelp/l3Ai2l+GbAuIg5HxHPATmCRpJnA9IjYEBEBrAHOzqtmMzMrjno/8v0K4J8lfZMsHP91\nap8FbChbb3dq+12ar2wflKRLgEsA5s6dW7uqzcys7uo9SKIX+EJEzAG+ANxWy51HxC0R0RMRPZ2d\nnbXctZmZ1Vm9A+p84J40/w/AwCCJPcCcsvVmp7Y9ab6y3czMWly9A+pF4N+k+T8GdqT5+4Hlktol\nzScbDLExIvYCByUtTqP3zgPuq3PNZmbWALldg5J0F/BJoEPSbuDrZKPxvitpPPA66XpRRGyVtB7Y\nBhwBLouIo2lXfWQjAicDD6TJzMxanLLBca2np6cn+vv7G12G2ZCG+0pfq/5eWmuRtCkievLav+8k\nYWZmheSAMjOzQnJAmZlZITmgzMyskBxQZmZWSMMOM5d0/wj2cSAiLqhNOWZmZplq34P6A2C4x34K\n+F7tyjEzM8tUC6irIuInw60g6W9rWI+ZmRlQ/RrU9qEWSOoFiIj1Na3IzMyM6gF1r6SFlY2p13Rx\nPiWZmZlVD6h/D/yDpD8EUOYm4I/I7rNnZmaWi2EDKiI2kT3Bdq2kpcDdQCewNCIO1qE+s5ZUKpUa\nXYJZ4Q0bUJJOIXuK7fnAWrIn3P4nYGpaZmajcPnllze6BLPCqzaKbxMwcFvlQ8C/AjaSDS8P4L35\nlWbWuvbv39/oEswKb9iAioj59SrEzDIzZsxodAlmhVDtFN97qu1gJOuY2ch997vfbXQJZoVQbRTf\nj0awj5GsY2YjtGLFikaXYFYI1a5BfVjScKP1BHg0n5mZ1Vy1a1Bt9SrEzMysnB+3YWZmhVRtkES1\nU4BmZma5qNaD2liXKszMzCpUCyiNdseSVkt6RdKWivbPS3pG0lZJ15e1Xylpp6Ttks4qa18oaXNa\ndqOkUddkZmbNo9opvE5JXxxqYUR8e5htbwf+K7BmoEHSp4BlwIcj4rCkd6f2bmA5sAA4DXhY0vsj\n4iiwiuzO6T8lG9K+FHigSt1mheX78JmNTLUeVBswDfi9IaYhRcSjwIGK5l7guog4nNZ5JbUvA9ZF\nxOGIeA7YCSySNBOYHhEbIiLIwu7skb45syK69NJLG12CWVOo1oPaGxH/uYav937gE5KuAV4HvhwR\nPwNmARvK1tud2n6X5ivbzZrWq6++2ugSzJpCtYCq9fWe8cApwGLgXwLrJdXshrOSLgEuAZg7d26t\ndmtmZg1Q7RTfkhq/3m7gnshsBN4EOoA9wJyy9Wantj1pvrJ9UBFxS0T0RERPZ2dnjUs3y19vb2+j\nSzArjGoPLKy8hjRWPwQ+BSDp/cBE4FfA/cBySe2S5gOnAxsjYi9wUNLiNHrvPOC+GtdkVhgrV65s\ndAlmhZHbF3El3UX2WPgOSbuBrwOrgdVp6PkbwPlp8MNWSeuBbcAR4LI0gg+gj2xE4GSy0XsewWdm\ndgJQlg+tp6enJ/r7+xtdhtkxhvsqX6v+PlprkrQpInry2r/vxWdWR319fY0uwaxpOKDM6mjVqlWN\nLsGsaTigzMyskBxQZgUxderURpdgVigOKLOCuPnmmxtdglmhOKDMCmLFihWNLsGsUBxQZnXiu5ib\nHR8HlFmd+C7mZsfHAWVWJ76LudnxcUCZFYBvEmt2LAeUWR1Uu4OEbxJrdiwHlFkd3HTTTY0uwazp\nOKDM6sA3gTU7fg4os5xVG17u609mg/PjNsxy1tHRwf79+4dc3qq/g9b6/LgNsyY3XDiZ2dAcUGYN\n5BvEmg3NAWWWo2rXn3yDWLOh+RqUWY58/clama9BmTUxX38yGz0HlFlOPLzcbGx8is8sJ5MnT+b1\n118fcnmr/u7ZicOn+MyaUKlUGjacJNWxGrPmlFtASVot6RVJWwZZ9iVJIamjrO1KSTslbZd0Vln7\nQkmb07Ib5d9sawIXXXTRsMv9bCiz6vLsQd0OLK1slDQH+DTwy7K2bmA5sCBts1JSW1q8CrgYOD1N\nx+zTrEiq9Z7Ady83G4ncAioiHgUODLLoO8BXgPIT8MuAdRFxOCKeA3YCiyTNBKZHxIbITtivAc7O\nq2azWrjggguGXe7BEWYjU9drUJKWAXsi4smKRbOAF8p+3p3aZqX5yvah9n+JpH5J/fv27atR1WYj\n19fXx5EjR4Zdx70ns5GpW0BJmgJ8Dbg6r9eIiFsioiciejo7O/N6GbMhrVq1atjlvrWR2ciNr+Nr\nvQ+YDzyZxjnMBh6TtAjYA8wpW3d2atuT5ivbzQqn2lNzwbc2MjsedetBRcTmiHh3RHRFRBfZ6bqP\nRcRLwP3AckntkuaTDYbYGBF7gYOSFqfRe+cB99WrZrPjUa33tGTJElasWFGnasyaX57DzO8C/i/w\n+5J2S/rcUOtGxFZgPbAN+Cfgsog4mhb3AbeSDZz4BfBAXjWbjdaCBQuqrvPwww/XoRKz1pHbKb6I\n+Msqy7sqfr4GuGaQ9fqBM2panFkNlUoltm3bNuw6EydOrFM1Zq3Dd5IwG6Nqw8oBVq9enX8hZi3G\nAWU2BmeeeWbVYeW+9mQ2Og4os1E688wzeeSRR4Zd57TTTvO1J7NRckCZjcJIwglgzx5/K8JstBxQ\nZsepVCqNKJx8SyOzsXFAmR2HUqnEOeecU3W9trY239LIbIwcUGYj1NfXN6JwArjjjjtyrsas9Tmg\nzEagr6+v6p0iBvT29nrUnlkN1PNefGZNaaQDIiALJ5/aM6sNB5TZMBYsWFD1LhEDlixZ4nAyqyEH\nlNkQZs2axYsvvjiidbu7u/19J7Mac0CZDeLkk0/mN7/5zYjW7e7uZuvWrTlXZHbi8SAJszKlUglJ\nIw6n0047zeFklhP3oMyS4xkMAXDSSSf5ThFmOXIPyk54pVKJ9vb24wqn7u5ufv3rX+dYlZk5oOyE\nNvDl2zfeeGPE2/T29vq0nlkd+BSfnZBKpRIXXnjhcQUT+HtOZvXkgLITymiDSRJ33nmn7xBhVkcO\nKDshjDaYwMPIzRrF16CspfX19SHpuK8zAUyaNIm1a9c6nMwaxAFlLadUKjFt2jQkjfgGr5V6e3v5\n7W9/61N6Zg3kU3zWMsZyGm/ApEmTuPXWWx1MZgXgHpQ1tVKpREdHx6hP45Vzr8msWHILKEmrJb0i\naUtZ2w2SnpH0lKR7JZ1UtuxKSTslbZd0Vln7Qkmb07IbJSmvmq05VIbS/v37x7S/3t5eIsLDx80K\nJs8e1O3A0oq2h4AzIuJDwM+BKwEkdQPLgQVpm5WS2tI2q4CLgdPTVLlPOwGUSiW6urpqFkqSHExm\nBZdbQEXEo8CBirYHI+JI+nEDMDvNLwPWRcThiHgO2AkskjQTmB4RGyIigDXA2XnVbMVR3ksaCKVd\nu3aNeb8DI/PefPNNB5NZwTXyGtSFwANpfhbwQtmy3altVpqvbB+UpEsk9Uvq37dvX43LtbzV+tRd\nuWnTprF27VpfYzJrIg0JKElXAUeAUi33GxG3RERPRPR0dnbWcteWg/LTduPGjat5KMHb15cOHTrk\nYDJrMnUfZi7pAuCzwJJ02g5gDzCnbLXZqW0Pb58GLG+3JlEqlbjqqqv45S9/ySmnnALA/v37kcTb\n//l5x/xYSeLSSy/1KTyzJlfXHpSkpcBXgD+PiNfKFt0PLJfULmk+2WCIjRGxFzgoaXEavXcecF89\na7bjM1ivaNeuXUQE+/fvf6uHVMtAGjBjxgxfXzJrIbn1oCTdBXwS6JC0G/g62ai9duChNFp8Q0Rc\nGhFbJa0HtpGd+rssIo6mXfWRjQicTHbN6gGsoRrRKxrMvHnzuOaaa3zqzqxFKe8/Io3S09MT/f39\njS6jJQwE0q5du44JoXoZN24cb775pkPJrEAkbYqInrz27ztJnMAGTseNGzeOjo4OOjo6GDduHF1d\nXfT19b11qu7cc899a4h3PcNp4JRdRHD06FEigueff97hZHaCcA+qxZWfjps7dy6f+cxn+NGPftTQ\n3tBg3EMyaz7uQdmgyns/5T2e8t5Qee8nIti1axerVq1qSG8IYOrUqcyYMQOAtrbsRiHz5s1j7dq1\n7iGZ2THcgyqYoXo81QYkFJF7RWatLe8elAMqRyMJmwMHDjRd8JSbOnUqkyZNYv/+/bS1tXH06FEH\nktkJwgE1SrUIqMGGU5cHyoEDB5q+lzNSA+/FIWRmAxxQo/Te97433nzzzWPCZSQ9mYH5Q4cOjen5\nQs1gIHgGrg0N9hnNnTvXYWRmx3BAjZKk1nxjozAQQvPmzXPwmFnN5B1QfuR7ExsqeCp7jA4hM2tG\nDqgCcvCYmTmgcjeSsHHwmJkdywE1jAkTJjB9+vRBB1AMzA836MJhY2Y2ei0dUCMZoeaejJlZMbVs\nQE2cOJHVq1c7YMzMmlTL3ovvgx/8oMPJzKyJtWxAmZlZc3NAmZlZITmgzMyskBxQZmZWSA4oMzMr\nJAeUmZkVkgPKzMwKKbeAkrRa0iuStpS1nSLpIUk70r8nly27UtJOSdslnVXWvlDS5rTsRknKq2Yz\nMyuOPHtQtwNLK9q+CjwSEacDj6SfkdQNLAcWpG1WSmpL26wCLgZOT1PlPs3MrAXlFlAR8ShwoKJ5\nGXBHmr8DOLusfV1EHI6I54CdwCJJM4HpEbEhsicrrinbxszMWli9r0GdGhF70/xLwKlpfhbwQtl6\nu1PbrDRf2T4oSZdI6pfUv2/fvtpVbWZmddewQRKpR1TTx7JHxC0R0RMRPZ2dnbXctZmZ1Vm9A+rl\ndNqO9O8rqX0PMKdsvdmpbU+ar2w3M7MWV++Auh84P82fD9xX1r5cUruk+WSDITam04EHJS1Oo/fO\nK9vGzMxaWG7Pg5J0F/BJoEPSbuDrwHXAekmfA3YBfwEQEVslrQe2AUeAyyLiaNpVH9mIwMnAA2ky\nM7MWp+xSUOvp6emJ/v7+RpdhZtayJG2KiJ689u87SZiZWSE5oMzMrJAcUGZmVkgOKDMzKyQHlJmZ\nFVLLjuKTdAjY3ug6RqgD+FWjixihZqoVmqte15qPZqoVmqve34+I38tr57l9D6oAtuc5/LGWJPW7\n1nw0U72uNR/NVCs0V72Scv0uj0/xmZlZITmgzMyskFo5oG5pdAHHwbXmp5nqda35aKZaobnqzbXW\nlh0kYWZmza2Ve1BmZtbEHFBmZlZIhQ0oSUslbZe0U9JXB1l+sqR7JT0laaOkM8qWXS5pi6Stkq4o\naz9F0kOSdqR/Ty5bdmV6re2SzipArTdIeiZtc6+kk1J7l6TfSnoiTTcVoNZvSNpTVtNnypaN+nPN\nsd7vl9X6vKQnUvuoP1tJqyW9ImnLEMsl6cb0Pp6S9LFq7zGv4zXHevM6ZvOoNZdjNqdaa3681qje\nQbet6XEbEYWbgDbgF8B7gYnAk0B3xTo3AF9P8x8AHknzZwBbgClk3/N6GPgXadn1wFfT/FeBv0/z\n3ek12oH56bXbGlzrp4Hxaf7vy2rtArYU7HP9BvDlQV5v1J9rnvVWbP8t4OoafLZ/BHxsqO2Bz5A9\ny0zAYuCn1d5jHsdrzvXW/JjNsda8jtma15rH8TrWeofbtpbHbVF7UIuAnRHxbES8AawDllWs0w38\nT4CIeAboknQq8AdkH+JrEXEE+Anwb9M2y4A70vwdwNll7esi4nBEPAfsTDU0rNaIeDC1AWwge9z9\nWOX1uQ5lLJ9r7vVKEtlDM+86jpoGFRGPAgeGWWUZsCYyG4CTJM2s8h7zOF5zqzenYzavz3a4fY36\ns82z1loerzWod7hta3bcFjWgZgEvlP28O7WVe5L0B0fSImAe2S/EFuATkmZImkL2fwBz0janRvYY\neYCXgFOP4/XqXWu5C3nnk4Tnpy79TyR9YoR15l3r59MpgNVlXfqxfK551wvwCeDliNhR1jbaz7aa\nod7LcO8xj+M1z3rL1eqYzbPWPI7ZvGqF+h6v1eodTs2O26IG1EhcR5bmTwCfBx4HjkbE02SnFx4E\n/gl4AjhauXFkfc56jbEfda2SrgKOAKXUtBeYGxEfAb4I/HdJ0xtc6yqyUxMfSfV9q4b15FHvgL/k\nnf83mvdnO2p1Pl7HpAHH7Gg08pgdraY5XgeM9bgt6r349vDO/9udndreEhEHgf8Ib3V9nwOeTctu\nA25Ly64lS2qAlyXNjIi9qZv6ykhfrwG1IukC4LPAkvQfmog4DBxO85sk/QJ4PzCSe2LlUmtEvFxW\n838D/nGkr9eIetPP48l6XgvL9jWWz3a072XCMO8xj+M1z3rzOGZzqTXHY7bmtaYa6328Vqt3OLU7\nboe7QNWoiSw4nyW7kDZwwXBBxTonARPT/MVk50kHlr07/TsXeAY4Kf18A++8eHd9ml/AOy/ePcvI\nB0nkVetSYBvQWbGvzoHayP4PcA9wSoNrnVm2zhfIzjOP6XPNs96yz/cntfps0zZdDH2x+U9558Xm\njdXeYx7Ha8711vyYzbHWXI7ZPGrN63gdS73DbVvL4/a43kg9J7JrBj8nG+lxVWq7FLg0zf9hWr4d\nuAc4uWzb/5V+UZ4k+z+5gfYZwCPADrJRXaeULbsqvdZ24E8KUOtOsvO1T6TpptT+74Ctqe0x4M8K\nUOudwGbgKeB+3vnLP+rPNa9607LbB/ZR1jbqz5bs1Mte4HdkPbXPVdQp4HvpfWwGeoZ7j3kerznW\nm9cxm0etuRyzedSax/Fao3qP2bbWx61vdWRmZoXUzIMkzMyshTmgzMyskBxQZmZWSA4oMzMrJAeU\nmZkVkgPKbATSLZMG7hr9kt55J+z/k8PrXSBpn6Rba7Cv/5DuIP2P1dc2K46i3knCrFAiYj/ZbXGQ\n9A3g1Yj4Zs4v+/2I+Kux7iQivi/pZeDLNajJrG7cgzIbI0mvpn8/mW7aeZ+kZyVdJ2mFsudUbZb0\nvrRep6QfSPpZmj4+gte4QNIPlT1f53lJfyXpi5Iel7RB0ilpvb+WtC3dBHVdvu/cLF/uQZnV1ofJ\nHvVxgOxWLrdGxCJJl5PdzPYK4LvAdyLif0uaC/xz2qaaM4CPApPI7trwNxHxUUnfAc4D/gvZrWXm\nR8RhpQcGmjUrB5RZbf0s0qMG0s07H0ztm4FPpfkzge7s3rYATJc0LSJerbLvH0fEIeCQpP8H/I+y\nfX8ozT8FlCT9EPjhmN+NWQM5oMxq63DZ/JtlP7/J279v44DFEfF6Dvv+U7Innf4ZcJWkD8bbDxE0\nayq+BmVWfw+Sne4DQNJHarFTSeOAORHxY+BvgHcB02qxb7NGcA/KrP7+GviepKfIfgcfJbuD9Fi1\nAWslvYvsLtQ3RsRvarBfs4bw3czNCig9+K+nFsPM0/4+CXw5Ij5bi/2Z1YNP8ZkV02+BP6nVF3WB\nlcCvx1yVWR25B2VmZoXkHpSZmRWSA8rMzArJAWVmZoXkgDIzs0JyQJmZWSH9f4rQbe5UN5VZAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106281c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish Adjoint 27.15\n",
      "Start Brute Force\n",
      "Finish Brute Force 35.23 [s]\n",
      "14 4.45180845443\n",
      "22 11.6854922906\n"
     ]
    }
   ],
   "source": [
    "pres = 20\n",
    "temp = 1000.0\n",
    "phi = .5\n",
    "simtype = 'UV'\n",
    "\n",
    "fuel = 'CH4'\n",
    "mech = 'gri30.xml'\n",
    "dk = 5.e-2\n",
    "\n",
    "# mech = 'h2_li_19.xml'\n",
    "# fuel = 'nc7h16'\n",
    "# mech = 'nc7sk88.cti'\n",
    "# gas = ct.Solution('mech/dme_sk39.cti')\n",
    "# gas = ct.Solution('mech/ic8sk143.cti')\n",
    "# mech = \"c4_49.xml\"\n",
    "\n",
    "string_list = [fuel, mech, str(phi), str(pres), str(temp), str(dk), simtype]\n",
    "string_list = '_'.join(string_list)\n",
    "print(string_list)\n",
    "\n",
    "\n",
    "mechfile = 'mech/'+mech\n",
    "gas = ct.Solution(mechfile)\n",
    "gas.set_equivalence_ratio( phi, fuel, 'O2:1.0, N2:3.76' )\n",
    "gas.TP = temp, pres*ct.one_atm\n",
    "\n",
    "# print('species_names:')\n",
    "# print(gas.species_names)\n",
    "# print(gas.species_index('OH'))\n",
    "\n",
    "# Get equilibrium temperature for ignition break\n",
    "gas.equilibrate(simtype)\n",
    "T_equi = gas.T\n",
    "\n",
    "m = gas.n_reactions\n",
    "\n",
    "# Create a dataframe to store sensitivity-analysis data\n",
    "ds = pd.DataFrame(data=[], index=gas.reaction_equations(range(m)))\n",
    "pd.options.display.float_format = '{:,.2e}'.format\n",
    "\n",
    "# Create an empty column to store the sensitivities data. \n",
    "# baseCase for brute force method\n",
    "ds[\"index\"] = \"\"\n",
    "ds[\"bruteforce\"] = \"\"\n",
    "ds[\"adjoint\"] = \"\"\n",
    "ds[\"ratio\"] = \"\"\n",
    "\n",
    "# base case\n",
    "factor = np.ones( gas.n_reactions )\n",
    "ign0 = ign_uq(gas, factor, temp, pres, fuel, phi, T_equi, simtype, False)\n",
    "print(\"Ignition Delay is: {:.4f} ms\".format(ign0*1000))\n",
    "\n",
    "\n",
    "start_time = timeit.default_timer()\n",
    "print('Start Adjoint')\n",
    "\n",
    "sens_T_ign = ign_sens(gas, temp, pres, fuel, phi, ign0, simtype, True)\n",
    "\n",
    "elapsed = timeit.default_timer() - start_time\n",
    "print('Finish Adjoint {:.2f}'.format(elapsed))\n",
    "\n",
    "\n",
    "factor = np.ones( ( gas.n_reactions, 1 ) )\n",
    "\n",
    "start_time = timeit.default_timer()\n",
    "print('Start Brute Force')\n",
    "\n",
    "for i in range(m):\n",
    "    factor[i] = 1+dk\n",
    "    ign = ign_uq(gas, factor, temp, pres, fuel, phi, T_equi, simtype, False)\n",
    "    factor[i] = 1.0\n",
    "    ds[\"bruteforce\"][i] = (ign-ign0)/(ign0*dk)\n",
    "\n",
    "elapsed = timeit.default_timer() - start_time\n",
    "print('Finish Brute Force {:.2f} [s]'.format(elapsed))\n",
    "\n",
    "for i in range(m):\n",
    "    ds['index'][i] = i\n",
    "    ds[\"adjoint\"][i] = -sens_T_ign[i]\n",
    "    ds[\"ratio\"][i] = ds[\"bruteforce\"][i] /(ds[\"adjoint\"][i] + 1.E-30)\n",
    "\n",
    "ds['adjointabs'] = ds['adjoint'].abs()\n",
    "ds = ds.sort_values( by=[\"adjointabs\"] ,ascending=False)\n",
    "# print( ds.to_string(index=True) )\n",
    "\n",
    "\n",
    "index = ds['adjoint'].abs() > 0.05*ds['adjoint'].abs().max()\n",
    "rxn1 = sum(index)\n",
    "error1 = (ds[index]['ratio'].max() - ds[index]['ratio'].min())/ds[index]['ratio'].max()*100\n",
    "print( sum(index) ,error1 )\n",
    "\n",
    "index = ds['adjoint'].abs() > 0.01*ds['adjoint'].abs().max()\n",
    "rxn2 = sum(index)\n",
    "error2 = (ds[index]['ratio'].max() - ds[index]['ratio'].min())/ds[index]['ratio'].max()*100\n",
    "print( sum(index) ,error2 )\n",
    "\n",
    "ubf = ds['bruteforce']/np.linalg.norm(ds['bruteforce'])\n",
    "uad = ds['adjoint']/np.linalg.norm(ds['adjoint'])\n",
    "results_df = pd.DataFrame(data = {\"maxerror_0.1\":[error1, rxn1],\"maxerror_0.01\":[error2, rxn2], \"ign0|cos\":[ign0, np.dot(ubf.values, uad.values)]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                             index bruteforce   adjoint     ratio adjointabs\n",
      "CH3 + HO2 <=> CH3O + OH        118  -6.26e-01 -2.05e+03  3.05e-04   2.05e+03\n",
      "2 CH3 (+M) <=> C2H6 (+M)       157   5.18e-01  1.64e+03  3.17e-04   1.64e+03\n",
      "CH2O + O2 <=> HCO + HO2         31  -3.48e-01 -1.12e+03  3.11e-04   1.12e+03\n",
      "CH3 + O2 <=> CH2O + OH         155  -2.45e-01 -7.81e+02  3.13e-04   7.81e+02\n",
      "2 HO2 <=> H2O2 + O2            115   2.34e-01  7.38e+02  3.17e-04   7.38e+02\n",
      "CH3 + O2 <=> CH3O + O          154  -1.95e-01 -6.17e+02  3.16e-04   6.17e+02\n",
      "CH2O + HO2 <=> H2O2 + HCO      120  -1.35e-01 -4.32e+02  3.13e-04   4.32e+02\n",
      "2 OH (+M) <=> H2O2 (+M)         84  -1.08e-01 -3.53e+02  3.06e-04   3.53e+02\n",
      "CH2O + CH3 <=> CH4 + HCO       160  -9.79e-02 -3.11e+02  3.15e-04   3.11e+02\n",
      "CH2O + OH <=> H2O + HCO        100  -7.26e-02 -2.32e+02  3.13e-04   2.32e+02\n",
      "CH4 + OH <=> CH3 + H2O          97   6.44e-02  2.07e+02  3.11e-04   2.07e+02\n",
      "2 HO2 <=> H2O2 + O2            114   6.49e-02  2.04e+02  3.19e-04   2.04e+02\n",
      "CH4 + H <=> CH3 + H2            52   4.30e-02  1.35e+02  3.19e-04   1.35e+02\n",
      "H + O2 <=> O + OH               37  -3.88e-02 -1.22e+02  3.17e-04   1.22e+02\n",
      "CH3 + H2O2 <=> CH4 + HO2       156  -3.05e-02 -1.01e+02  3.00e-04   1.01e+02\n",
      "CH3 + CH3OH <=> CH3O + CH4     162   2.93e-02  9.18e+01  3.19e-04   9.18e+01\n",
      "CH3 + HO2 <=> CH4 + O2         117   2.73e-02  8.02e+01  3.40e-04   8.02e+01\n",
      "CH3O + O2 <=> CH2O + HO2       169  -1.97e-02 -6.34e+01  3.11e-04   6.34e+01\n",
      "H + N2 + O2 <=> HO2 + N2        35  -1.15e-02 -3.75e+01  3.05e-04   3.75e+01\n",
      "HO2 + OH <=> H2O + O2           86   1.11e-02  3.48e+01  3.19e-04   3.48e+01\n",
      "CH2O + H (+M) <=> CH3O (+M)     56  -8.99e-03 -2.91e+01  3.08e-04   2.91e+01\n",
      "C2H6 + OH <=> C2H5 + H2O       112  -6.80e-03 -2.15e+01  3.16e-04   2.15e+01\n",
      "2 CH3 <=> C2H5 + H             158  -4.21e-03 -1.33e+01  3.18e-04   1.33e+01\n",
      "CH3 + OH (+M) <=> CH3OH (+M)    94   3.42e-03  1.06e+01  3.22e-04   1.06e+01\n",
      "CH2O + H <=> H2 + HCO           57   3.16e-03  1.00e+01  3.16e-04   1.00e+01\n",
      "CH3 + H (+M) <=> CH4 (+M)       51   2.77e-03  8.59e+00  3.23e-04   8.59e+00\n",
      "C2H6 + CH3 <=> C2H5 + CH4      164  -2.66e-03 -8.34e+00  3.19e-04   8.34e+00\n",
      "H + 2 O2 <=> HO2 + O2           33  -2.60e-03 -8.09e+00  3.21e-04   8.09e+00\n",
      "CH2O + O <=> HCO + OH           14  -2.41e-03 -7.62e+00  3.16e-04   7.62e+00\n",
      "H + HO2 <=> H2 + O2             44   2.41e-03  7.56e+00  3.19e-04   7.56e+00\n",
      "...                            ...        ...       ...       ...        ...\n",
      "CH2 + NO <=> HCN + OH          249  -1.57e-06 -3.44e-09  4.55e+02   3.44e-09\n",
      "H + NH <=> H2 + N              190   3.48e-06 -3.16e-09 -1.10e+03   3.16e-09\n",
      "CN + O2 <=> NCO + O            219   1.39e-06 -2.74e-09 -5.08e+02   2.74e-09\n",
      "CH2(S) + NO <=> H + HNCO       251  -1.32e-06  2.30e-09 -5.76e+02   2.30e-09\n",
      "NH3 + O <=> NH2 + OH           278   2.44e-06  2.27e-09  1.08e+03   2.27e-09\n",
      "HNCO + OH <=> CO2 + NH2        267   3.62e-06 -1.93e-09 -1.87e+03   1.93e-09\n",
      "N + OH <=> H + NO              179   1.32e-06 -1.92e-09 -6.89e+02   1.92e-09\n",
      "CH + H2 (+M) <=> CH3 (+M)      288   1.17e-06  1.81e-09  6.45e+02   1.81e-09\n",
      "CH + CH3 <=> C2H3 + H          128  -1.71e-05  1.76e-09 -9.70e+03   1.76e-09\n",
      "CN + OH <=> H + NCO            217   3.57e-07  1.15e-09  3.12e+02   1.15e-09\n",
      "HCN + OH <=> H + HOCN          233   2.28e-06 -8.20e-10 -2.78e+03   8.20e-10\n",
      "HCN + OH <=> CO + NH2          235   2.32e-07 -6.52e-10 -3.56e+02   6.52e-10\n",
      "H + HCNO <=> H + HNCO          269   2.84e-06 -2.78e-10 -1.02e+04   2.78e-10\n",
      "HCN + OH <=> H + HNCO          234   4.07e-06 -2.15e-10 -1.89e+04   2.15e-10\n",
      "H + HOCN <=> H + HNCO          272   2.21e-06 -1.22e-10 -1.81e+04   1.22e-10\n",
      "CN + O <=> CO + N              216   3.73e-06  6.96e-11  5.36e+04   6.96e-11\n",
      "HCN + M <=> CN + H + M         229   6.37e-07  3.97e-11  1.60e+04   3.97e-11\n",
      "C + CH2 <=> C2H + H            122   2.74e-06  1.48e-11  1.85e+05   1.48e-11\n",
      "CH + NO <=> H + NCO            246   2.06e-06  1.09e-11  1.89e+05   1.09e-11\n",
      "NCO + NO <=> CO2 + N2          228  -1.56e-06  7.26e-12 -2.15e+05   7.26e-12\n",
      "H + HCNO <=> CO + NH2          271   2.65e-06  7.20e-12  3.68e+05   7.20e-12\n",
      "H + HCNO <=> HCN + OH          270   5.33e-06  2.81e-12  1.90e+06   2.81e-12\n",
      "C + NO <=> CN + O              243  -1.37e-06 -1.12e-12  1.23e+06   1.12e-12\n",
      "CH + NO <=> HCO + N            247   6.39e-06 -9.61e-13 -6.66e+06   9.61e-13\n",
      "HNCO + O <=> CO + HNO          262   3.70e-06  5.32e-14  6.96e+07   5.32e-14\n",
      "CN + NO2 <=> NCO + NO          280   2.64e-06 -2.81e-14 -9.39e+07   2.81e-14\n",
      "NCO + NO2 <=> CO2 + N2O        281   2.23e-06 -3.87e-15 -5.78e+08   3.87e-15\n",
      "C + NO <=> CO + N              244   7.55e-07  3.28e-15  2.30e+08   3.28e-15\n",
      "AR + H + O2 <=> AR + HO2        36   0.00e+00 -0.00e+00  0.00e+00   0.00e+00\n",
      "AR + CH2(S) <=> AR + CH2       142   0.00e+00 -0.00e+00  0.00e+00   0.00e+00\n",
      "\n",
      "[325 rows x 5 columns]\n"
     ]
    }
   ],
   "source": [
    "print(ds)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}