{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Two years ago, I've had the pleasure to take a Harvard class entitled [Science & Cooking: From Haute Cuisine to Soft Matter Science (part 1)](https://www.edx.org/course/science-cooking-from-haute-cuisine-to-soft-matter) on edx, a MOOC platform.\n",
"\n",
"This has lead me to do some fun kitchen science, in particular trying to [cook an egg using math formulas](https://flothesof.github.io/cooking-egg-math.html) (spoiler: I was not successful).\n",
"\n",
"One of the tools introduced by the class is the [Cook my Meat App](http://up.csail.mit.edu/science-of-cooking/). The app allows you to specify a cooking protocol for meat modelled as a 1D medium. Ever since seing this app, I thought it would be nice to replicate it with my own tools and learn from it. This post is the attempt to do just this. \n",
"\n",
"I'll try to go from the theory (the heat equation in 1D) to the implementation using the Crank-Nicolson time stepping method, in Python. In the end, we'll try to compare our results with these from the app (and they'll hopefully be similar). \n",
"\n",
"I'll also apply these simulations to some other cooking problems."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The 1D heat equation and the Crank-Nicolson method "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The heat equation using the Crank-Nicolson discretization\n",
"\n",
"The problem tackled by the Cook my Meat app is that of cooking a steak. The steak has thickness $L$ and can be heated on one side or both sides by a heat source, while being in contact with the outside air on the other side. It also starts at a given initial temperature.\n",
"\n",
"Our problem then consists of solving for the temperature distribution as a function of time. The temperature evolution is governed by the 1D heat equation with a constant diffusion constant, written as \n",
"\n",
"$$\n",
"\\frac{\\partial T}{\\partial t} = D \\frac{\\partial^2 T}{\\partial x^2} \n",
"$$\n",
"\n",
"The physical interpretation of this equation is that heat flows through space along the $x$ coordinate at a rate that is proportional to the diffusion constant $D$ and to the local second derivative with respect to space of temperature. \n",
"\n",
"The Crank-Nicolson method consists, when discretizing the continuous temperature field to a $N + 1$ nodes grid (with notation $T_i^n$, where $i$ represents the space dimension and $n$ the time dimension), to use the following equations\n",
"\n",
"$$\n",
"\\frac{T_i^{n+1} - T_i^{n}}{\\Delta t} = \\dfrac{1}{2} \\left [ D \\frac{T_{i+1}^{n} - 2 T_{i}^{n} + T_{i-1}^{n}}{\\Delta x^2} + D \\frac{T_{i+1}^{n+1} - 2 T_{i}^{n+1} + T_{i-1}^{n+1}}{\\Delta x^2}\\right ] \\, \\forall i \\in [1, N-1]\n",
"$$\n",
"\n",
"Rearranging terms to have time $n+1$ on the left and time $n$ on the right and using the notation $\\sigma = \\frac{D \\Delta t}{2 \\Delta x ^2}$, this equation becomes\n",
"\n",
"$$\n",
"- \\sigma T_{i-1}^{n+1} + (1 + 2 \\sigma) T_{i}^{n+1} - \\sigma T_{i+1}^{n+1} = - \\sigma T_{i-1}^{n} + (1 - 2 \\sigma) T_{i}^{n} - \\sigma T_{i+1}^{n}\n",
"$$\n",
"\n",
"This is the general term in the medium that will be used to update temperatures across time and space."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Boundary conditions\n",
"\n",
"However, our problem still lacks boundary conditions. The two edges can have either a fixed temperature or a zero heat flux condition, depending on whether there is heat transfer (through a pan for instance, fixed temperature) or not (contact with air, zero heat flux).\n",
"\n",
"To write the boundary conditions properly, we will use a little trick, which I read about [here](http://georg.io/2013/12/Crank_Nicolson), based on the above equation. Let's tackle the left boundary first. Writing the equation for $i=0$, we obtain:\n",
"\n",
"$$\n",
"- \\sigma T_{-1}^{n+1} + (1 + 2 \\sigma) T_{0}^{n+1} - \\sigma T_{1}^{n+1} = - \\sigma T_{-1}^{n} + (1 - 2 \\sigma) T_{0}^{n} - \\sigma T_{1}^{n}\n",
"$$\n",
"\n",
"The peculiar thing here is the appearance of a *virtual* node of index -1. \n",
"\n",
"Since we wish to have a zero heat flux condition at the left boundary, we would write $\\frac{\\partial T}{\\partial x} |_{x=0} = 0$ in the continuous system.\n",
"\n",
"In the discrete system, we can use the virtual node and a finite difference to write\n",
"\n",
"\n",
"$$\n",
"\\frac{T_{0}^{n+1} - T_{-1}^{n+1}}{\\Delta x} = \\frac{T_{0}^{n} - T_{-1}^{n}}{\\Delta x} = 0\n",
"$$\n",
"\n",
"This allows us to obtain\n",
"\n",
"\n",
"$$\n",
"(1 + \\sigma) T_{0}^{n+1} -\\sigma T_{1}^{n+1} = (1 - \\sigma) T_{0}^{n} - \\sigma T_{1}^{n} \n",
"$$\n",
"\n",
"In the case of a fixed temperature on the left, denoted $T_l$, we can just write $T_{0}^{n+1} = T_l$ and put it on the first line of the system.\n",
"\n",
"\n",
"Doing the same thing at the right boundary yields the following formula in the case of zero heat flux right boundary\n",
"\n",
"\n",
"$$\n",
"-\\sigma T_{N}^{n+1} (1 + \\sigma) T_{N+1}^{n+1} = -\\sigma T_{N}^{n} (1 - \\sigma) T_{N+1}^{n}\n",
"$$\n",
"\n",
"and $T_{N}^{n+1} = T_r$ on the last line of the linear system.\n",
"\n",
"## Assembling the whole system\n",
"\n",
"Finally, the whole system becomes \n",
"\n",
"$$\n",
"A \\underline{T}^{n+1} = B \\underline{T}^{n} + C\n",
"$$\n",
"\n",
"To see what the $A$, $B$ and $C$ matrices look like in an exemple, we can write them out in the case of a fixed temperature on the left and a zero heat flux boundary condition on the right \n",
"\n",
"$$\n",
"\\begin{bmatrix}\n",
"1 & 0 & 0 & 0 & 0 & \\cdots & 0 & 0 & 0 & 0\\\\\\\\\n",
"-\\sigma & 1+2\\sigma & -\\sigma & 0 & 0 & \\cdots & 0 & 0 & 0 & 0 \\\\\\\\\n",
"0 & -\\sigma & 1+2\\sigma & -\\sigma & \\cdots & 0 & 0 & 0 & 0 & 0 \\\\\\\\\n",
"0 & 0 & \\ddots & \\ddots & \\ddots & \\ddots & 0 & 0 & 0 & 0 \\\\\\\\\n",
"0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\sigma & 1+2\\sigma & -\\sigma \\\\\\\\\n",
"0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\sigma & 1+\\sigma\n",
"\\end{bmatrix}\n",
"\\begin{bmatrix}\n",
"T_0^{n+1} \\\\\\\\\n",
"T_1^{n+1} \\\\\\\\\n",
"T_2^{n+1} \\\\\\\\\n",
"\\vdots \\\\\\\\\n",
"T_{N}^{n+1} \\\\\\\\\n",
"T_{N+1}^{n+1}\n",
"\\end{bmatrix} =\n",
"\\begin{bmatrix}\n",
"0 & 0 & 0 & 0 & 0 & \\cdots & 0 & 0 & 0 & 0\\\\\\\\\n",
"\\sigma & 1-2\\sigma & \\sigma & 0 & 0 & \\cdots & 0 & 0 & 0 & 0 \\\\\\\\\n",
"0 & \\sigma & 1-2\\sigma & \\sigma & \\cdots & 0 & 0 & 0 & 0 & 0 \\\\\\\\\n",
"0 & 0 & \\ddots & \\ddots & \\ddots & \\ddots & 0 & 0 & 0 & 0 \\\\\\\\\n",
"0 & 0 & 0 & 0 & 0 & 0 & 0 & \\sigma & 1-2\\sigma & \\sigma \\\\\\\\\n",
"0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\sigma & 1-\\sigma\n",
"\\end{bmatrix}\n",
"\\begin{bmatrix}\n",
"T_0^{n} \\\\\\\\\n",
"T_1^{n} \\\\\\\\\n",
"T_2^{n} \\\\\\\\\n",
"\\vdots \\\\\\\\\n",
"T_{N}^{n} \\\\\\\\\n",
"T_{N+1}^{n}\n",
"\\end{bmatrix} +\n",
"\\begin{bmatrix}\n",
"T_l \\\\\\\\\n",
"0 \\\\\\\\\n",
"0 \\\\\\\\\n",
"\\vdots \\\\\\\\\n",
"0 \\\\\\\\\n",
"0\n",
"\\end{bmatrix} \n",
"$$\n",
"\n",
"Assuming the matrix $B$ is invertible, the solution at timestep $n+1$ is obtained by \n",
"\n",
"$$\n",
" \\underline{T}^{n+1} = A^{-1} (B \\underline{T}^n + C)\n",
"$$\n",
"\n",
"Since material parameters are constant, we can compute $B^{-1} A$ only once and then apply it repeatedly, which should be computationally efficient."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Implementation using Python "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The units we will use for the implementation of our simulations are as follows: millimeters, seconds, ° Kelvin.\n",
"\n",
"In this unit, the heat diffusion constant of water (to which steak and other food is assimilated) is $D = 0.14$ (in mm^2/sec)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"application/javascript": [
"\n",
"(function(root) {\n",
" function now() {\n",
" return new Date();\n",
" }\n",
"\n",
" var force = true;\n",
"\n",
" if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n",
" root._bokeh_onload_callbacks = [];\n",
" root._bokeh_is_loading = undefined;\n",
" }\n",
"\n",
" var JS_MIME_TYPE = 'application/javascript';\n",
" var HTML_MIME_TYPE = 'text/html';\n",
" var EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n",
" var CLASS_NAME = 'output_bokeh rendered_html';\n",
"\n",
" /**\n",
" * Render data to the DOM node\n",
" */\n",
" function render(props, node) {\n",
" var script = document.createElement(\"script\");\n",
" node.appendChild(script);\n",
" }\n",
"\n",
" /**\n",
" * Handle when an output is cleared or removed\n",
" */\n",
" function handleClearOutput(event, handle) {\n",
" var cell = handle.cell;\n",
"\n",
" var id = cell.output_area._bokeh_element_id;\n",
" var server_id = cell.output_area._bokeh_server_id;\n",
" // Clean up Bokeh references\n",
" if (id != null && id in Bokeh.index) {\n",
" Bokeh.index[id].model.document.clear();\n",
" delete Bokeh.index[id];\n",
" }\n",
"\n",
" if (server_id !== undefined) {\n",
" // Clean up Bokeh references\n",
" var cmd = \"from bokeh.io.state import curstate; print(curstate().uuid_to_server['\" + server_id + \"'].get_sessions()[0].document.roots[0]._id)\";\n",
" cell.notebook.kernel.execute(cmd, {\n",
" iopub: {\n",
" output: function(msg) {\n",
" var id = msg.content.text.trim();\n",
" if (id in Bokeh.index) {\n",
" Bokeh.index[id].model.document.clear();\n",
" delete Bokeh.index[id];\n",
" }\n",
" }\n",
" }\n",
" });\n",
" // Destroy server and session\n",
" var cmd = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n",
" cell.notebook.kernel.execute(cmd);\n",
" }\n",
" }\n",
"\n",
" /**\n",
" * Handle when a new output is added\n",
" */\n",
" function handleAddOutput(event, handle) {\n",
" var output_area = handle.output_area;\n",
" var output = handle.output;\n",
"\n",
" // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n",
" if ((output.output_type != \"display_data\") || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n",
" return\n",
" }\n",
"\n",
" var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n",
"\n",
" if (output.metadata[EXEC_MIME_TYPE][\"id\"] !== undefined) {\n",
" toinsert[toinsert.length - 1].firstChild.textContent = output.data[JS_MIME_TYPE];\n",
" // store reference to embed id on output_area\n",
" output_area._bokeh_element_id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n",
" }\n",
" if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n",
" var bk_div = document.createElement(\"div\");\n",
" bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n",
" var script_attrs = bk_div.children[0].attributes;\n",
" for (var i = 0; i < script_attrs.length; i++) {\n",
" toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\n",
" }\n",
" // store reference to server id on output_area\n",
" output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n",
" }\n",
" }\n",
"\n",
" function register_renderer(events, OutputArea) {\n",
"\n",
" function append_mime(data, metadata, element) {\n",
" // create a DOM node to render to\n",
" var toinsert = this.create_output_subarea(\n",
" metadata,\n",
" CLASS_NAME,\n",
" EXEC_MIME_TYPE\n",
" );\n",
" this.keyboard_manager.register_events(toinsert);\n",
" // Render to node\n",
" var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n",
" render(props, toinsert[toinsert.length - 1]);\n",
" element.append(toinsert);\n",
" return toinsert\n",
" }\n",
"\n",
" /* Handle when an output is cleared or removed */\n",
" events.on('clear_output.CodeCell', handleClearOutput);\n",
" events.on('delete.Cell', handleClearOutput);\n",
"\n",
" /* Handle when a new output is added */\n",
" events.on('output_added.OutputArea', handleAddOutput);\n",
"\n",
" /**\n",
" * Register the mime type and append_mime function with output_area\n",
" */\n",
" OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n",
" /* Is output safe? */\n",
" safe: true,\n",
" /* Index of renderer in `output_area.display_order` */\n",
" index: 0\n",
" });\n",
" }\n",
"\n",
" // register the mime type if in Jupyter Notebook environment and previously unregistered\n",
" if (root.Jupyter !== undefined) {\n",
" var events = require('base/js/events');\n",
" var OutputArea = require('notebook/js/outputarea').OutputArea;\n",
"\n",
" if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n",
" register_renderer(events, OutputArea);\n",
" }\n",
" }\n",
"\n",
" \n",
" if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n",
" root._bokeh_timeout = Date.now() + 5000;\n",
" root._bokeh_failed_load = false;\n",
" }\n",
"\n",
" var NB_LOAD_WARNING = {'data': {'text/html':\n",
" \"
\\n\"+\n",
" \"
\\n\"+\n",
" \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n",
" \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n",
" \"
\\n\"+\n",
" \"
\\n\"+\n",
" \"
re-rerun `output_notebook()` to attempt to load from CDN again, or
\"}};\n",
"\n",
" function display_loaded() {\n",
" var el = document.getElementById(null);\n",
" if (el != null) {\n",
" el.textContent = \"BokehJS is loading...\";\n",
" }\n",
" if (root.Bokeh !== undefined) {\n",
" if (el != null) {\n",
" el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n",
" }\n",
" } else if (Date.now() < root._bokeh_timeout) {\n",
" setTimeout(display_loaded, 100)\n",
" }\n",
" }\n",
"\n",
"\n",
" function run_callbacks() {\n",
" try {\n",
" root._bokeh_onload_callbacks.forEach(function(callback) {\n",
" if (callback != null)\n",
" callback();\n",
" });\n",
" } finally {\n",
" delete root._bokeh_onload_callbacks\n",
" }\n",
" console.debug(\"Bokeh: all callbacks have finished\");\n",
" }\n",
"\n",
" function load_libs(css_urls, js_urls, callback) {\n",
" if (css_urls == null) css_urls = [];\n",
" if (js_urls == null) js_urls = [];\n",
"\n",
" root._bokeh_onload_callbacks.push(callback);\n",
" if (root._bokeh_is_loading > 0) {\n",
" console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n",
" return null;\n",
" }\n",
" if (js_urls == null || js_urls.length === 0) {\n",
" run_callbacks();\n",
" return null;\n",
" }\n",
" console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n",
" root._bokeh_is_loading = css_urls.length + js_urls.length;\n",
"\n",
" function on_load() {\n",
" root._bokeh_is_loading--;\n",
" if (root._bokeh_is_loading === 0) {\n",
" console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n",
" run_callbacks()\n",
" }\n",
" }\n",
"\n",
" function on_error() {\n",
" console.error(\"failed to load \" + url);\n",
" }\n",
"\n",
" for (var i = 0; i < css_urls.length; i++) {\n",
" var url = css_urls[i];\n",
" const element = document.createElement(\"link\");\n",
" element.onload = on_load;\n",
" element.onerror = on_error;\n",
" element.rel = \"stylesheet\";\n",
" element.type = \"text/css\";\n",
" element.href = url;\n",
" console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n",
" document.body.appendChild(element);\n",
" }\n",
"\n",
" for (var i = 0; i < js_urls.length; i++) {\n",
" var url = js_urls[i];\n",
" var element = document.createElement('script');\n",
" element.onload = on_load;\n",
" element.onerror = on_error;\n",
" element.async = false;\n",
" element.src = url;\n",
" console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n",
" document.head.appendChild(element);\n",
" }\n",
" };\n",
"\n",
" function inject_raw_css(css) {\n",
" const element = document.createElement(\"style\");\n",
" element.appendChild(document.createTextNode(css));\n",
" document.body.appendChild(element);\n",
" }\n",
"\n",
" \n",
" var js_urls = [];\n",
" var css_urls = [];\n",
" \n",
"\n",
" var inline_js = [\n",
" function(Bokeh) {\n",
" /* BEGIN bokeh.min.js */\n",
" /*!\n",
" * Copyright (c) 2012 - 2019, Anaconda, Inc., and Bokeh Contributors\n",
" * All rights reserved.\n",
" * \n",
" * Redistribution and use in source and binary forms, with or without modification,\n",
" * are permitted provided that the following conditions are met:\n",
" * \n",
" * Redistributions of source code must retain the above copyright notice,\n",
" * this list of conditions and the following disclaimer.\n",
" * \n",
" * Redistributions in binary form must reproduce the above copyright notice,\n",
" * this list of conditions and the following disclaimer in the documentation\n",
" * and/or other materials provided with the distribution.\n",
" * \n",
" * Neither the name of Anaconda nor the names of any contributors\n",
" * may be used to endorse or promote products derived from this software\n",
" * without specific prior written permission.\n",
" * \n",
" * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\"\n",
" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE\n",
" * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE\n",
" * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE\n",
" * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR\n",
" * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF\n",
" * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS\n",
" * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN\n",
" * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)\n",
" * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF\n",
" * THE POSSIBILITY OF SUCH DAMAGE.\n",
" */\n",
" (function(root, factory) {\n",
" root[\"Bokeh\"] = factory();\n",
" })(this, function() {\n",
" var define;\n",
" var parent_require = typeof require === \"function\" && require\n",
" return (function(modules, entry, aliases, externals) {\n",
" if (aliases === undefined) aliases = {};\n",
" if (externals === undefined) externals = {};\n",
"\n",
" var cache = {};\n",
"\n",
" var normalize = function(name) {\n",
" if (typeof name === \"number\")\n",
" return name;\n",
"\n",
" if (name === \"bokehjs\")\n",
" return entry;\n",
"\n",
" var prefix = \"@bokehjs/\"\n",
" if (name.slice(0, prefix.length) === prefix)\n",
" name = name.slice(prefix.length)\n",
"\n",
" var alias = aliases[name]\n",
" if (alias != null)\n",
" return alias;\n",
"\n",
" var trailing = name.length > 0 && name[name.lenght-1] === \"/\";\n",
" var index = aliases[name + (trailing ? \"\" : \"/\") + \"index\"];\n",
" if (index != null)\n",
" return index;\n",
"\n",
" return name;\n",
" }\n",
"\n",
" var require = function(name) {\n",
" var mod = cache[name];\n",
" if (!mod) {\n",
" var id = normalize(name);\n",
"\n",
" mod = cache[id];\n",
" if (!mod) {\n",
" if (!modules[id]) {\n",
" if (parent_require && externals[id]) {\n",
" try {\n",
" mod = {exports: parent_require(id)};\n",
" cache[id] = cache[name] = mod;\n",
" return mod.exports;\n",
" } catch (e) {}\n",
" }\n",
"\n",
" var err = new Error(\"Cannot find module '\" + name + \"'\");\n",
" err.code = 'MODULE_NOT_FOUND';\n",
" throw err;\n",
" }\n",
"\n",
" mod = {exports: {}};\n",
" cache[id] = cache[name] = mod;\n",
" modules[id].call(mod.exports, require, mod, mod.exports);\n",
" } else\n",
" cache[name] = mod;\n",
" }\n",
"\n",
" return mod.exports;\n",
" }\n",
"\n",
" var main = require(entry);\n",
" main.require = require;\n",
"\n",
" main.register_plugin = function(plugin_modules, plugin_entry, plugin_aliases, plugin_externals) {\n",
" if (plugin_aliases === undefined) plugin_aliases = {};\n",
" if (plugin_externals === undefined) plugin_externals = {};\n",
"\n",
" for (var name in plugin_modules) {\n",
" modules[name] = plugin_modules[name];\n",
" }\n",
"\n",
" for (var name in plugin_aliases) {\n",
" aliases[name] = plugin_aliases[name];\n",
" }\n",
"\n",
" for (var name in plugin_externals) {\n",
" externals[name] = plugin_externals[name];\n",
" }\n",
"\n",
" var plugin = require(plugin_entry);\n",
"\n",
" for (var name in plugin) {\n",
" main[name] = plugin[name];\n",
" }\n",
"\n",
" return plugin;\n",
" }\n",
"\n",
" return main;\n",
" })\n",
" ([\n",
" function _(n,o,r){n(1),function(n){for(var o in n)r.hasOwnProperty(o)||(r[o]=n[o])}(n(102))},\n",
" function _(n,c,f){n(2),n(11),n(14),n(21),n(49),n(52),n(87),n(94),n(100)},\n",
" function _(e,n,a){e(3)()||Object.defineProperty(Object,\"assign\",{value:e(4),configurable:!0,enumerable:!1,writable:!0})},\n",
" function _(r,t,o){t.exports=function(){var r,t=Object.assign;return\"function\"==typeof t&&(t(r={foo:\"raz\"},{bar:\"dwa\"},{trzy:\"trzy\"}),r.foo+r.bar+r.trzy===\"razdwatrzy\")}},\n",
" function _(t,r,n){var o=t(5),c=t(10),a=Math.max;r.exports=function(t,r){var n,f,h,i=a(arguments.length,2);for(t=Object(c(t)),h=function(o){try{t[o]=r[o]}catch(t){n||(n=t)}},f=1;f= 0\");if(!isFinite(r))throw new RangeError(\"Count must be < ∞\");for(n=\"\";r;)r%2&&(n+=t),r>1&&(t+=t),r>>=1;return n}},\n",
" function _(t,i,n){var r=t(18),a=Math.abs,o=Math.floor;i.exports=function(t){return isNaN(t)?0:0!==(t=Number(t))&&isFinite(t)?r(t)*o(a(t)):t}},\n",
" function _(n,t,i){t.exports=n(19)()?Math.sign:n(20)},\n",
" function _(n,t,o){t.exports=function(){var n=Math.sign;return\"function\"==typeof n&&(1===n(10)&&-1===n(-20))}},\n",
" function _(n,r,t){r.exports=function(n){return n=Number(n),isNaN(n)||0===n?n:n>0?1:-1}},\n",
" function _(e,r,a){e(22)()||Object.defineProperty(Array,\"from\",{value:e(23),configurable:!0,enumerable:!1,writable:!0})},\n",
" function _(n,o,r){o.exports=function(){var n,o,r=Array.from;return\"function\"==typeof r&&(o=r(n=[\"raz\",\"dwa\"]),Boolean(o&&o!==n&&\"dwa\"===o[1]))}},\n",
" function _(e,l,r){var n=e(24).iterator,t=e(44),a=e(45),i=e(46),u=e(47),o=e(10),f=e(8),c=e(48),v=Array.isArray,h=Function.prototype.call,y={configurable:!0,enumerable:!0,writable:!0,value:null},s=Object.defineProperty;l.exports=function(e){var l,r,A,g,p,w,b,d,x,j,O=arguments[1],m=arguments[2];if(e=Object(o(e)),f(O)&&u(O),this&&this!==Array&&a(this))l=this;else{if(!O){if(t(e))return 1!==(p=e.length)?Array.apply(null,e):((g=new Array(1))[0]=e[0],g);if(v(e)){for(g=new Array(p=e.length),r=0;r
\n",
""
],
"text/plain": [
":QuadMesh [time,x] (temperature)"
]
},
"execution_count": 4,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1393"
}
},
"output_type": "execute_result"
}
],
"source": [
"solver_params = dict(x=np.linspace(0, 30, num=nx), \n",
" T_initial=np.ones_like(np.linspace(0, 30, num=nx)) * (23 + KELVIN) + 10 * np.exp(-0.01 * (np.linspace(0, 30, num=nx) - 10)**2), \n",
" dt=.0001, \n",
" D=0.14)\n",
"\n",
"boundary_conditions = dict(Tl=23 + KELVIN, \n",
" Tr=23 + KELVIN, \n",
" Tair=23 + KELVIN)\n",
"\n",
"time_grid, snapshots = make_snapshots(solver_params, boundary_conditions, tmax, n_snapshots)\n",
"\n",
"hv.QuadMesh((time_grid, solver_params['x'], np.array(snapshots).T - KELVIN), \n",
" kdims=['time', 'x'], \n",
" vdims='temperature')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here again, this behaviour is what we expect: a bump in temperature slowly diffuses to the surrounding medium making it more homogeneous as time progresses."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Using cook my meat app inputs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that the basic machinery is in place, we can turn to replicating the data from the Cook My Meat app. It relies on a simple text format to describe cooking procedures.\n",
"\n",
"The default recipe proposed by the app is as follows:\n",
"\n",
"```\n",
"3cm Steak starts at 23°C\n",
"150°C and 23°C for 4:00\n",
"23°C and 150°C for 4:00\n",
"23°C and 23°C for 5:00\n",
"```\n",
"It consists of a header and an initial temperature followed by boundary temperatures and the duration of cooking with these in place. \n",
"\n",
"We can build a simple parser using regular expressions that extracts the information from such an input. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'thickness_mm': 30.0,\n",
" 'initial_temperature_C': 296.15,\n",
" 'steps': [(423.15, 296.15, 240),\n",
" (296.15, 423.15, 240),\n",
" (296.15, 296.15, 300)]}"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import re\n",
"\n",
"four_minutes_a_side = \"\"\"3cm Steak starts at 23°C\n",
"150°C and 23°C for 4:00\n",
"23°C and 150°C for 4:00\n",
"23°C and 23°C for 5:00\"\"\"\n",
"\n",
"def parse_input(recipe):\n",
" \"\"\"Parses the Cook my meat text format and returns a dictionary with data properly formatted.\"\"\"\n",
" header, *lines = recipe.split('\\n')\n",
" thickness_cm, T_initial = list(map(float, re.findall(\"([\\d.]+)cm Steak starts at ([-\\d.]+)°[C]\", header)[0]))\n",
" T_initial += KELVIN\n",
" steps = []\n",
" for line in lines:\n",
" Tleft, Tright, mins, secs = re.findall(\"([-\\d.]+)°[C] and ([-\\d.]+)°[C] for ([\\d]+):(\\d{2})\", line)[0]\n",
" steps.append((float(Tleft) + KELVIN, float(Tright) + KELVIN, 60 * int(mins) + int(secs)))\n",
" parsed_params = dict(thickness_mm=thickness_cm * 10, initial_temperature_C=T_initial, steps=steps)\n",
" return parsed_params\n",
"\n",
"parse_input(four_minutes_a_side)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we can program a simulation that runs each substep of the cooking while using the final temperature as input for the next simulation step."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":QuadMesh [time,x] (temperature)"
]
},
"execution_count": 6,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1514"
}
},
"output_type": "execute_result"
}
],
"source": [
"def griddify(parsed_params, fine_time_grid):\n",
" \"\"\"Transforms the input temperature data to a regular time grid.\"\"\"\n",
" data = np.array(parsed_params['steps'])\n",
" temps, t = data[:, :2], data[:, 2]\n",
" t = np.cumsum(t)\n",
" index = 0\n",
" boundary_conditions = []\n",
" for tt in fine_time_grid:\n",
" if tt > t[index]:\n",
" index += 1\n",
" if index >= data.shape[0]:\n",
" index -= 1\n",
" boundary_conditions.append(temps[index])\n",
" return boundary_conditions\n",
"\n",
"def cook_recipe(recipe, n_snapshots=100, M=100, nx=100):\n",
" \"\"\"Takes a Cook my meat recipe as input, simulates it and returns temperature snapshots.\"\"\"\n",
" parsed_params = parse_input(recipe)\n",
" total_time = sum([step_time for Tl, Tr, step_time in parsed_params['steps']])\n",
" coarse_time_grid = np.linspace(0, total_time, num=n_snapshots)\n",
" coarse_time_grid_indices = np.arange(n_snapshots) * M\n",
" fine_time_grid = np.linspace(0, total_time, num=n_snapshots * M)\n",
" dt = fine_time_grid[1] - fine_time_grid[0]\n",
" boundary_conditions = griddify(parsed_params, fine_time_grid)\n",
" x = np.linspace(0, parsed_params['thickness_mm'], num=nx)\n",
" Tair = parsed_params['initial_temperature_C']\n",
" T = np.ones_like(x) * Tair\n",
" solver = CrankNicolsonSolver(x, T, dt, D=0.14)\n",
" snapshots = []\n",
" for ind, (t, (Tl, Tr)) in enumerate(zip(fine_time_grid, boundary_conditions)):\n",
" solver.set_boundary_conditions(Tl, Tr, Tair)\n",
" solver.step()\n",
" if ind in coarse_time_grid_indices:\n",
" T = solver.T.copy()\n",
" snapshots.append(T)\n",
" return x, coarse_time_grid, snapshots\n",
"\n",
"x, time_grid, snapshots = cook_recipe(four_minutes_a_side, n_snapshots=50, M=500, nx=50)\n",
"\n",
"hv.QuadMesh((time_grid, x, np.array(snapshots).T - KELVIN), \n",
" kdims=['time', 'x'], \n",
" vdims='temperature').redim.range(temperature=(23, 100))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The above plot is qualitatively similar to the original app.\n",
"\n",
"We can compare the final temperature from our solver with the final output of the app (obtained with some Javascript hacking):"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"x_app = np.linspace(0, 30, num=30)\n",
"final_temps_app = np.array([86.18586243519505,\n",
" 85.86449361301695,\n",
" 85.25487377408177,\n",
" 84.37978835482498,\n",
" 83.27157213371525,\n",
" 81.97050758592101,\n",
" 80.52286390721143,\n",
" 78.97869610729124,\n",
" 77.38953387284053,\n",
" 75.80609063774105,\n",
" 74.27611499608037,\n",
" 72.84249054685091,\n",
" 71.54166834143024,\n",
" 70.40249052953442,\n",
" 69.44543686657468,\n",
" 68.68229960704663,\n",
" 68.11626877794897,\n",
" 67.74239024927502,\n",
" 67.54834421381581,\n",
" 67.51548195134583,\n",
" 67.62005391346293,\n",
" 67.83456168692639,\n",
" 68.12916948444723,\n",
" 68.473116548245,\n",
" 68.83607928479354,\n",
" 69.1894401963763,\n",
" 69.50742898563637,\n",
" 69.76810900571606,\n",
" 69.954189125051,\n",
" 70.05364687780765])"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":Overlay\n",
" .Curve.Cook_my_meat :Curve [x] (y)\n",
" .Curve.Ours :Curve [x] (y)"
]
},
"execution_count": 8,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1634"
}
},
"output_type": "execute_result"
}
],
"source": [
"hv.Curve((x_app, final_temps_app[::-1]), label='Cook my meat') * hv.Curve((x, snapshots[-1] - KELVIN), label='ours')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There seems to be a little offset between them.\n",
"\n",
"I can think of several reasons for which the resultats are not in *very* good agreement:\n",
"\n",
"- time and space grids are not the same in the two versions\n",
"- I used a method to project our system to a grid of times instead of applying a constant dt=1 throughout the simulation\n",
"\n",
"However, the obtained temperatures still look quite close and the spatial shape is also quite resembling. So I'd say this validates our replication attempt."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Trying the different cooking methods "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's now try the other default recipes featured in the app using our solver."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":QuadMesh [time,x] (temperature)"
]
},
"execution_count": 9,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "1854"
}
},
"output_type": "execute_result"
}
],
"source": [
"every_15_secs = \"\"\"3cm Steak starts at 23°C\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\n",
"150°C and 23°C for 0:15\n",
"23°C and 150°C for 0:15\"\"\"\n",
"\n",
"x, time_grid, snapshots = cook_recipe(every_15_secs)\n",
"\n",
"hv.QuadMesh((time_grid, x, np.array(snapshots).T - KELVIN), \n",
" kdims=['time', 'x'], \n",
" vdims='temperature').redim.range(temperature=(23, 100))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
""
],
"text/plain": [
":QuadMesh [time,x] (temperature)"
]
},
"execution_count": 11,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "2096"
}
},
"output_type": "execute_result"
}
],
"source": [
"sous_vide_liquid_nitrogen = \"\"\"3cm Steak starts at 23°C\n",
"53°C and 53°C for 60:00\n",
"-200°C and -200°C for 0:30\n",
"200°C and 200°C for 2:00\"\"\"\n",
"\n",
"x, time_grid, snapshots = cook_recipe(sous_vide_liquid_nitrogen)\n",
"\n",
"hv.QuadMesh((time_grid, x, np.array(snapshots).T - KELVIN), \n",
" kdims=['time', 'x'], \n",
" vdims='temperature').redim.range(temperature=(23, 100))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All of these cooking methods look quite similar to these produced by the Cook my meat app. Let's then turn to some analysis using our tool."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Post-processing to get the meat state "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As in the app, let's first add some post-processing that allows to determine the meat quality instead of the temperature.\n",
"\n",
"This essentially consists of monitoring the maximum temperature seen at a given location as a function of time.\n",
"\n",
"This in turn determines the state of the proteins according to this table:\n",
"\n",
"| Cooking state \t| Protein and temperature of change \t|\n",
"|---------------\t|---------------------------------------\t|\n",
"| Raw \t| No protein denaturization below 40°C \t|\n",
"| Rare \t| Myosin denatures 40-55°C \t|\n",
"| Medium Rare \t| Glycogen denatures 55-60°C \t|\n",
"| Medium \t| Myoglobin denatures 60-70°C \t|\n",
"| Well \t| Actin denatures 70-120°C \t|\n",
"| Browned \t| Browning reactions 120-180°C \t|\n",
"| Charred \t| Charring at >180°C \t|\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can implement this using the matplotlib color bar logic which already has a function for mapping from values to integers `matplotlib.colors.BoundaryNorm`."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-1 0 1 2 3 4 5 6\n"
]
},
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
""
],
"text/plain": [
":Layout\n",
" .QuadMesh.Flip_every_15_seconds :QuadMesh [time,x] (protein_state)\n",
" .QuadMesh.Sear_then_cook_low :QuadMesh [time,x] (protein_state)\n",
" .QuadMesh.Sous_vide_plus_liquid_nitrogen :QuadMesh [time,x] (protein_state)"
]
},
"execution_count": 13,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "2511"
}
},
"output_type": "execute_result"
}
],
"source": [
"layouts = []\n",
"for recipe, label in zip([every_15_secs, sear_then_cook_low, sous_vide_liquid_nitrogen],\n",
" ['flip every 15 seconds', 'sear then cook low', 'sous vide + liquid nitrogen']):\n",
" x, time_grid, snapshots = cook_recipe(recipe)\n",
" qm = hv.QuadMesh((time_grid, x, np.array(snapshots_to_protein_state(snapshots)).T),\n",
" kdims=['time', 'x'], \n",
" vdims='protein_state', label=label).options(cmap='reds_r')\n",
" layouts.append(qm)\n",
"hv.Layout(layouts).cols(1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Comparison of cooking procedures "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"While discussing this post with colleagues, I found out that the solver we've developed is a nice tool to make hypotheses and understand the physics of cooking. In particular, I found two experiments interesting. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Why the liquid nitrogen matters "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we have another look at this recipe\n",
"\n",
"```\n",
"3cm Steak starts at 23°C\n",
"53°C and 53°C for 60:00\n",
"-200°C and -200°C for 0:30\n",
"200°C and 200°C for 2:00\n",
"```\n",
"\n",
"we could wonder if the liquid nitrogen is really needed. Does this really make a difference if we leave it out?\n",
"\n",
"Let's see:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
""
],
"text/plain": [
":Layout\n",
" .QuadMesh.Sous_vide_with_liquid_nitrogen :QuadMesh [time,x] (protein_state)\n",
" .QuadMesh.Modified_sous_vide_without_liquid_nitrogen :QuadMesh [time,x] (protein_state)"
]
},
"execution_count": 15,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "3426"
}
},
"output_type": "execute_result"
}
],
"source": [
"(qm1[3750, :] + qm2[3750 - 30, :]).cols(1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What we do see from the above comparison is that the cooling acts as a sort of \"heat buffer\" that prevents the searing at the end from adding too much heat to the steak center and thus overcooking much of the inner part. So in a sense, the liquid nitrogen method is quite optimal, the core being medium rare cooked and at the same having a \"crust\" which is slightly more cooked."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Several hours low-temperature roasts "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A good example of a slow-roast recipe can be found here: https://www.thespruceeats.com/prime-rib-roast-slow-method-995284.\n",
"\n",
"\n",
"What is interesting is that the recipe has quite a lot of steps:\n",
"\n",
"- searing the roast\n",
"- cooking at low temperature (200°F/93°C) until it reaches 128°F/53°C\n",
"- letting it sit until it reaches 120°C/48°C\n",
"\n",
"I tried to transcribe the recipe using the Cook my meat format:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":Layout\n",
" .QuadMesh.Slow_roast_left_parenthesis_temperature_right_parenthesis :QuadMesh [time,x] (temperature)\n",
" .QuadMesh.Slow_roast_left_parenthesis_protein_state_right_parenthesis :QuadMesh [time,x] (protein_state)"
]
},
"execution_count": 16,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "3816"
}
},
"output_type": "execute_result"
}
],
"source": [
"slow_roast = \"\"\"12cm Steak starts at 23°C\n",
"200°C and 23°C for 3:00\n",
"23°C and 200°C for 3:00\n",
"93°C and 93°C for 100:00\n",
"23°C and 23°C for 25:00\"\"\"\n",
"\n",
"x, time_grid, snapshots = cook_recipe(slow_roast)\n",
"qmt = hv.QuadMesh((time_grid, x, np.array(snapshots).T - KELVIN), \n",
" kdims=['time', 'x'], \n",
" vdims='temperature', label='slow roast (temperature)').redim.range(temperature=(30, 70))\n",
"\n",
"qmp = hv.QuadMesh((time_grid, x, np.array(snapshots_to_protein_state(snapshots)).T),\n",
" kdims=['time', 'x'], \n",
" vdims='protein_state', label='slow roast (protein state)').options(cmap='reds_r')\n",
"\n",
"hv.Layout([qmt, qmp]).cols(1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As can be seen, the resting part is really important too since it allows the heat to change the proteins from rare to medium-rate. \n",
"Now the interesting part is to play with the recipe. Let's say we cook the steak 10°C higher. What happens?"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":Layout\n",
" .QuadMesh.Slow_roast_left_parenthesis_temperature_right_parenthesis :QuadMesh [time,x] (temperature)\n",
" .QuadMesh.Slow_roast_left_parenthesis_protein_state_right_parenthesis :QuadMesh [time,x] (protein_state)"
]
},
"execution_count": 17,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "4206"
}
},
"output_type": "execute_result"
}
],
"source": [
"slow_roast_10deg_higher = \"\"\"12cm Steak starts at 23°C\n",
"200°C and 23°C for 3:00\n",
"23°C and 200°C for 3:00\n",
"103°C and 103°C for 100:00\n",
"23°C and 23°C for 25:00\"\"\"\n",
"\n",
"x, time_grid, snapshots = cook_recipe(slow_roast_10deg_higher)\n",
"qmt = hv.QuadMesh((time_grid, x, np.array(snapshots).T - KELVIN), \n",
" kdims=['time', 'x'], \n",
" vdims='temperature', label='slow roast (temperature)').redim.range(temperature=(30, 70))\n",
"\n",
"qmp = hv.QuadMesh((time_grid, x, np.array(snapshots_to_protein_state(snapshots)).T),\n",
" kdims=['time', 'x'], \n",
" vdims='protein_state', label='slow roast (protein state)').options(cmap='reds_r')\n",
"\n",
"hv.Layout([qmt, qmp]).cols(1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Too bad, the roast is now medium instead of medium rare! \n",
"\n",
"An interesting application of our simulation tool is that we can use it to trace the right cooking time as a function of steak thickness. We will do this by a simple grid search."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "ba3c2e3805f3450cbec84b526d0cab43",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HBox(children=(FloatProgress(value=0.0, max=17.0), HTML(value='')))"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"from tqdm.notebook import tqdm\n",
"\n",
"thicknesses = np.arange(8.5, 17, 0.5)\n",
"durations = np.linspace(30, 215, num=30)\n",
"results = []\n",
"for thickness in tqdm(thicknesses):\n",
" for duration in durations:\n",
" recipe = \"\"\"{}cm Steak starts at 23°C\n",
"200°C and 23°C for 3:00\n",
"23°C and 200°C for 3:00\n",
"103°C and 103°C for {:.0f}:00\n",
"23°C and 23°C for 25:00\"\"\".format(thickness, duration)\n",
" x, time_grid, snapshots = cook_recipe(recipe)\n",
" protein_state = snapshots_to_protein_state(snapshots)[-1]\n",
" val = np.sum(protein_state == 1) / protein_state.size\n",
" results.append((thickness, duration, val * 100))"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.holoviews_exec.v0+json": "",
"text/html": [
"
\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" \n",
"
\n",
""
],
"text/plain": [
":HeatMap [thickness,duration] (percent_medium_rare)"
]
},
"execution_count": 19,
"metadata": {
"application/vnd.holoviews_exec.v0+json": {
"id": "4480"
}
},
"output_type": "execute_result"
}
],
"source": [
"hv.HeatMap(results, kdims=['thickness', 'duration'], vdims='percent_medium_rare').opts(tools=['hover'], width=500)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The heatmap is quite intuitive: there is a sweet spot in terms of cooking time where the roast becomes just right and if you pass it, it's overcooked. Below the sweet spot you get a little of the meat cooked the way you want it, but not the most of it. Also, an interesting observation from this plot is that the thicker your piece of meat, the wider the margin of error you have in terms of duration for cooking your roast right."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Conclusion "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this blog post, we went from implementing the heat equation with several boundary conditions to checking its results to doing some kitchen science. \n",
"I hope you've enjoyed reading about this as much as I had pleasure developing it. Now back to cooking eggs :)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*This post was entirely written using the Jupyter Notebook. Its content is BSD-licensed. You can see a static view or download this notebook with the help of nbviewer at [20200109_CookMyMeatClone.ipynb](http://nbviewer.ipython.org/urls/raw.github.com/flothesof/posts/master/20200109_CookMyMeatClone.ipynb).*"
]
}
],
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}