{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Shock tube modeling with constant volume\n", "\n", "This example is available as an ipynb (Jupyter Notebook) file in the main GitHub repository at https://github.com/pr-omethe-us/PyKED/blob/main/docs/shock-tube-example.ipynb\n", "\n", "The ChemKED file used in this example can be found in the `tests` directory of the PyKED\n", "repository at https://github.com/pr-omethe-us/PyKED/blob/main/tests/testfile_st_p5.yaml.\n", "The data in this file comes from [Stranic et al.][Stranic2012], describing\n", "shock-tube ignition delays for _tert_-butanol.\n", "\n", "[Stranic2012]: https://doi.org/10.1016/j.combustflame.2011.08.014\n", "\n", "The file starts with the meta information about the ChemKED file itself:\n", "\n", "```yaml\n", "file-authors:\n", " - name: Morgan Mayer\n", " ORCID: 0000-0001-7137-5721\n", "file-version: 0\n", "chemked-version: 0.0.1\n", "```\n", "\n", "This is followed by the information about the reference for the data, and the information about the experimental apparatus:\n", "\n", "```yaml\n", "reference:\n", " doi: 10.1016/j.combustflame.2011.08.014\n", " authors:\n", " - name: Ivo Stranic\n", " - name: Deanna P. Chase\n", " - name: Joseph T. Harmon\n", " - name: Sheng Yang\n", " - name: David F. Davidson\n", " - name: Ronald K. Hanson\n", " journal: Combustion and Flame\n", " year: 2012\n", " volume: 159\n", " pages: 516-527\n", "experiment-type: ignition delay\n", "apparatus:\n", " kind: shock tube\n", " institution: High Temperature Gasdynamics Laboratory, Stanford University\n", " facility: stainless steel shock tube\n", "```\n", "\n", "This ChemKED file specifies multiple data points with some common\n", "conditions, including a common mixture composition and common definition of\n", "ignition delay. Therefore, a `common-properties` section is specified, followed by the\n", "`datapoints` list:\n", "\n", "```yaml\n", "common-properties:\n", " composition: &comp\n", " kind: mole fraction\n", " species:\n", " - species-name: t-butanol\n", " InChI: 1S/C4H10O/c1-4(2,3)5/h5H,1-3H3\n", " amount:\n", " - 0.003333333\n", " - species-name: O2\n", " InChI: 1S/O2/c1-2\n", " amount:\n", " - 0.04\n", " - species-name: Ar\n", " InChI: 1S/Ar\n", " amount:\n", " - 0.956666667\n", " ignition-type: &ign\n", " target: OH*\n", " type: 1/2 max\n", "datapoints:\n", " - temperature:\n", " - 1459 kelvin\n", " ignition-delay:\n", " - 347 us\n", " pressure:\n", " - 1.60 atm\n", " composition: *comp\n", " ignition-type: *ign\n", " equivalence-ratio: 0.5\n", " - temperature:\n", " - 1389 kelvin\n", " ignition-delay:\n", " - 756 us\n", " pressure:\n", " - 1.67 atm\n", " composition: *comp\n", " ignition-type: *ign\n", " equivalence-ratio: 0.5\n", " - temperature:\n", " - 1497 kelvin\n", " ignition-delay:\n", " - 212 us\n", " pressure:\n", " - 1.55 atm\n", " composition: *comp\n", " ignition-type: *ign\n", " equivalence-ratio: 0.5\n", " - temperature:\n", " - 1562 kelvin\n", " ignition-delay:\n", " - 105 us\n", " pressure:\n", " - 1.50 atm\n", " composition: *comp\n", " ignition-type: *ign\n", " equivalence-ratio: 0.5\n", "```\n", "\n", "In this example, we will run constant-volume simulations at each\n", "pressure and temperature condition in the `datapoints` list. Once again,\n", "the ChemKED file specifies all the information required for the simulations except\n", "for the chemical kinetic model, and Cantera can be used to simulate autoignition. For these simulations, we will be using the butanol mechanism developed by [Sarathy et al.](https://doi.org/10.1016/j.combustflame.2011.12.017) and available from the [LLNL website](https://combustion.llnl.gov/mechanisms/alcohols/butanol-isomers). In Python, additional functionality can be imported into a script or session by the `import`\n", "keyword. Cantera, NumPy, and PyKED must be imported into the session so that we can work with the\n", "code. In the case of Cantera and NumPy, we will use many functions from these libraries, so we\n", "assign them abbreviations (`ct` and `np`, respectively) for convenience. From PyKED, we\n", "will only be using the `ChemKED` class, so this is all that is imported. In this example, we also import\n", "Python's built-in `multiprocessing` library so that we can run the simulations on multiple cores. We\n", "import the `Pool` class, which offers facilities to manage a jobs on a pool of processors:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from multiprocessing import Pool\n", "\n", "import cantera as ct\n", "import numpy as np\n", "\n", "from pyked import ChemKED" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, we define a function that will be mapped onto each job. This function takes the initial\n", "temperature, pressure, and mole fractions as input and returns the simulated ignition delay time,\n", "defined as the time the temperature increases by 400 K over the initial temperature, a\n", "simplified definition for this example. In general, the user could process the mole fraction of\n", "OH\\* (provided that the kinetic model includes accurate chemistry for OH\\*) to match the\n", "definition of ignition delay in the experiments. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Suppress warnings from loading the mechanism file\n", "ct.suppress_thermo_warnings()\n", "\n", "\n", "def run_simulation(T, P, X):\n", " gas = ct.Solution(\"LLNL_sarathy_butanol.cti\")\n", " gas.TPX = T, P, X\n", " reac = ct.IdealGasReactor(gas)\n", " netw = ct.ReactorNet([reac])\n", " while reac.T < T + 400:\n", " netw.step()\n", "\n", " return netw.time" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, we load the ChemKED file and generate a list of initial conditions that will be mapped onto\n", "the `run_simulation()` function. We first define a convenience function to collect the input\n", "from a single datapoint and return a Python tuple with the conditions (tuples are lightweight groups\n", "of data in Python). In this function, we define a conversion between the species names in the ChemKED file and the species in the mechanism using a dictionary. Then, we use the built-in `map` function to apply the\n", "`collect_input` function to each of the elements in the `ck.datapoints` list:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from urllib.request import urlopen\n", "\n", "import yaml\n", "\n", "st_link = \"https://raw.githubusercontent.com/pr-omethe-us/PyKED/main/tests/testfile_st_p5.yaml\"\n", "with urlopen(st_link) as response:\n", " testfile_st = yaml.safe_load(response.read())\n", "ck = ChemKED(dict_input=testfile_st)\n", "\n", "\n", "def collect_input(dp):\n", " T_initial = dp.temperature.to(\"K\").magnitude\n", " P_initial = dp.pressure.to(\"Pa\").magnitude\n", " species_conversion = {\"t-butanol\": \"tc4h9oh\", \"O2\": \"o2\", \"Ar\": \"ar\"}\n", " X_initial = dp.get_cantera_mole_fraction(species_conversion)\n", " return (T_initial, P_initial, X_initial)\n", "\n", "\n", "initial_conditions = list(map(collect_input, ck.datapoints))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we create the processor `Pool` (with 4 processes) and send the jobs out to run:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The ignition delay for T_initial=1459 K, P_initial=162120.0 Pa is: 0.0006703341602991564 seconds\n", "The ignition delay for T_initial=1389 K, P_initial=169212.75 Pa is: 0.0012624985834179291 seconds\n", "The ignition delay for T_initial=1497 K, P_initial=157053.75 Pa is: 0.0005355282509854273 seconds\n", "The ignition delay for T_initial=1562 K, P_initial=151987.5 Pa is: 0.00043835671071621276 seconds\n" ] } ], "source": [ "with Pool(processes=4) as pool:\n", " ignition_delays = pool.starmap(run_simulation, initial_conditions)\n", "\n", "for (T, P, X), tau in zip(initial_conditions, ignition_delays):\n", " print(f\"The ignition delay for T_initial={T} K, P_initial={P} Pa is: {tau} seconds\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The simulated ignition delay results are returned in the `ignition_delays` list. The results\n", "are printed to the screen using a Python formatted string (f-string). In addition, the results can be plotted against the experimental data on an Arrhenius plot" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " this.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('