Panini_Model_Jupyter.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DE Coursework\n",
    "# Shahmeer Hassan & Tsz Yin Yung\n",
    "\n",
    "# This code takes a square chicken panini and heats it in a sandwich grill\n",
    "\n",
    "# The loading bar in this code only works in Jupyter\n",
    "# Use Panini_Model.py if NOT running in Jupyter\n",
    "\n",
    "\n",
    "#import libraries\n",
    "import numpy as np\n",
    "import sys # Used for a clean exit\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "\n",
    "# This imports the library to create a laoding bar\n",
    "from tqdm.auto import tqdm\n",
    "\n",
    "# This is to open an image at the end\n",
    "from PIL import Image\n",
    "\n",
    "# This is to open an mp4 animation at the end\n",
    "from os import startfile\n",
    "\n",
    "\n",
    "# This function asks the user for a mesh size\n",
    "def get_h():\n",
    "    # Ask for a step size\n",
    "    h = input(\"Enter the mesh size in millimetres (numerical value only): \")\n",
    "    h = float(h) # Convert to a float number\n",
    "    h = h/1000 # convert to metres\n",
    "    return h\n",
    "\n",
    "\n",
    "def get_k(h, alpha_bread, alpha_chicken):\n",
    "    \n",
    "    # Find largest possible time step for both bread and chicken\n",
    "    k_maxbread = 0.25*(h**2)/alpha_bread\n",
    "    k_maxchicken = 0.25*(h**2)/alpha_chicken\n",
    "    \n",
    "    # Select lowest of the two\n",
    "    k_max = min([k_maxbread, k_maxchicken])\n",
    "    \n",
    "    print()\n",
    "    print(\"The maximum time step that will give a stable model is \"+str(k_max)+\" seconds\") \n",
    "    \n",
    "    # Ask user what they want as a time step\n",
    "    k = input(\"Please enter a time step in seconds (numerical value only): \")\n",
    "    k = float(k) # convert to a number\n",
    "    \n",
    "    return k\n",
    "    \n",
    "\n",
    "# This function recieves a boolean variable, execute\n",
    "    # and asks the user if they want to continue with an unstable system\n",
    "def get_input(execute):\n",
    "    \n",
    "    # If the Courant condition is satisfied\n",
    "    if execute == True:\n",
    "        print(\"Proceeding with calculation...\")\n",
    "        return execute\n",
    "    \n",
    "    # If the Courant condition is not satisfied\n",
    "    elif execute == False:\n",
    "        ask_user = input(\"Do you wish to continue? (y/n):\") #Ask the user\n",
    "        if ask_user == 'y' or ask_user=='Y' or ask_user=='yes':\n",
    "            #Set the boolean variable to True\n",
    "            execute = True\n",
    "            print(\"Proceeding with calculation...\")\n",
    "            return execute\n",
    "            \n",
    "        elif ask_user == 'n' or ask_user=='N' or ask_user=='no':\n",
    "            #End the program\n",
    "            sys.exit(\"The user terminated the program\")\n",
    "            \n",
    "        else:\n",
    "            #The user mistyped. Ask again.\n",
    "            print(\"Invalid response\")\n",
    "            execute = get_input(execute)\n",
    "            \n",
    "\n",
    "# This function opens an image\n",
    "def open_image(filename):\n",
    "    image = Image.open(filename)\n",
    "    image.show()\n",
    "\n",
    "\n",
    "#~~~~~~~~~~~~~~~~~~~~ Define the inputs ~~~~~~~~~~~~~~~~~~~~#\n",
    "\n",
    "# Thickness of bread\n",
    "bread_thickness = 0.003 # [m]\n",
    "\n",
    "#Thickness of filling\n",
    "filling_thickness = 0.01 # [m]\n",
    "\n",
    "# Dimensions of square sandwich\n",
    "XL = 0.1 # side length of sandwich [m]\n",
    "YL = (2*bread_thickness) + filling_thickness # height of sandwich [m]\n",
    "\n",
    "#set initial temperature of sandwich\n",
    "initial_temperature = 0 # degC\n",
    "\n",
    "# Time forecast\n",
    "ttl_time = input(\"How long do you want to put the panini in the press for in seconds? \")\n",
    "ttl_time = float(ttl_time)\n",
    "\n",
    "# Define thermal properties\n",
    "# Define heat transfer coefficient between sandwich and air\n",
    "h_conv = 10 # [W/m^2 K]\n",
    "\n",
    "# Define conductivity\n",
    "k_bread = 0.068 # [W/m K]\n",
    "k_chicken = 1.6 # [W/m K]\n",
    "\n",
    "#Define density\n",
    "rho_bread = 200 # [kg/m^3]\n",
    "rho_chicken = 1068 # [kg/m^3]\n",
    "\n",
    "#Define specfic heat capacity\n",
    "cp_bread = 2800 # [J/kg K]\n",
    "cp_chicken = 3220 # [J/kg K]\n",
    "\n",
    "# Define Thermal Diffusivity\n",
    "alpha_bread = k_bread/(rho_bread*cp_bread)\n",
    "alpha_chicken = k_chicken/(rho_chicken * cp_chicken)\n",
    "\n",
    "# Weighted thermal diffusivity for use in chicken-filling boundaries\n",
    "alpha_breadbc = (3*alpha_bread + alpha_chicken)/4\n",
    "alpha_chickenbc = (3*alpha_chicken + alpha_bread)/4\n",
    "\n",
    "#Define ambient temperature\n",
    "T_amb = 21 # degC\n",
    "\n",
    "#Define hotplate temperature\n",
    "T_plate = 150 # degC\n",
    "\n",
    "#Speed of video relative to real time\n",
    "# E.g. speed_factor = 5 delivers video x5 faster than realtime\n",
    "speed_factor = 2\n",
    "\n",
    "#step in space\n",
    "h = get_h() # units: [m]\n",
    "\n",
    "#step in time\n",
    "k = get_k(h, alpha_bread, alpha_chicken) # units: [s]\n",
    "\n",
    "#~~~~~~~~~~~~Initialise the Temperatures 3D matrix~~~~~~~~~~~~~~~~~~~~#\n",
    "\n",
    "#First initialise the 2D Matrix at t=0\n",
    "xdim = int(XL/h)+1\n",
    "ydim = int(YL/h)+1\n",
    "T0 = np.zeros((ydim,xdim))\n",
    "\n",
    "#Then intialise the 3D temperatures matrix\n",
    "tdim = int(ttl_time/k)+1\n",
    "T = np.zeros((tdim,ydim,xdim))\n",
    "\n",
    "# Apply this initial temperature of sandwich to the initial temperature matrix\n",
    "for i in range(0, len(T0)):\n",
    "    for j in range(0,len(T0[0])):\n",
    "        T0[i,j]=initial_temperature\n",
    "\n",
    "# Add the temp matrix of t=0 to the temperatures 3D matrix\n",
    "T[0]=T0\n",
    "\n",
    "#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#\n",
    "\n",
    "# Define indices for sandwich\n",
    "breaddim = int(bread_thickness/h)\n",
    "filldim = ydim - 2 * breaddim\n",
    "\n",
    "#~~~~~~~~~~~~~~~~~Check for Convergence~~~~~~~~~~~~~~~~~~~~~~~~~~~#\n",
    "\n",
    "#Calculate Courant values for bread and chicken\n",
    "courant_bread=alpha_bread*k/(h**2)\n",
    "courant_chicken=alpha_chicken*k/(h**2)\n",
    "\n",
    "# Calculate weighted Courant values for later use\n",
    "courant_3b1c = (3*courant_bread + courant_chicken)/4\n",
    "courant_1b3c = (3*courant_chicken + courant_bread)/4\n",
    "\n",
    "print() # Blank line\n",
    "\n",
    "if courant_bread <= 0.25 and courant_chicken <= 0.25:\n",
    "    #Courant satisfied, stable\n",
    "    print('Stable, Courant condition satisfied')\n",
    "    execute = True\n",
    "    \n",
    "elif courant_bread <= 0.25 and courant_chicken > 0.25:\n",
    "    #Courant not satisfied, unstable\n",
    "    print('Potentially unstable, Courant condition not satisfied for panini filling')\n",
    "    execute = False\n",
    "    \n",
    "elif courant_bread > 0.25 and courant_chicken <= 0.25:\n",
    "    print('Potentially unstable, Courant condition not satisfied for bread')\n",
    "    execute = False\n",
    "    \n",
    "else:\n",
    "    print('Potentially unstable, Courant condition not satisfied for bread or panini filling')\n",
    "    execute = False\n",
    "\n",
    "# Tell user what the courant values are\n",
    "print(\"Bread: alpha dt / (dx)^2 =\",courant_bread)\n",
    "print(\"Filling: alpha dt / (dx)^2 =\", courant_chicken)\n",
    "\n",
    "#Ask the user whether to continue to run the program\n",
    "execute = get_input(execute)\n",
    "\n",
    "#~~~~~~~~~~~~~~~~~~~~ Solve problem ~~~~~~~~~~~~~~~~~~~~#\n",
    "\n",
    "for t in tqdm(range(1,len(T))): #for every instance in time, starting at t=0+k\n",
    "\n",
    "    #Define boundary conditions\n",
    "\n",
    "    #Top and bottom\n",
    "    for j in range(0,len(T[0,0])): # for each column in the matrix\n",
    "\n",
    "        # Top of food: Dirichlet\n",
    "        T[t,0,j] = T_plate    #top of food = T_plate\n",
    "\n",
    "        # Bottom of food: Dirichlet\n",
    "        T[t,ydim-1,j]= T_plate #bottom of food = T_plate\n",
    "        \n",
    "        \n",
    "    # Boundary conditions on sides - convection to environment\n",
    "    for i in range(1,ydim-1):\n",
    "        \n",
    "        #For bread\n",
    "        if i <= breaddim-1 or i >= filldim + breaddim:\n",
    "            \n",
    "            # Left side\n",
    "            T[t,i,0] =(h/(k_bread+h_conv*h))*(h_conv*T_amb+T[t-1,i,1]*k_bread/h)\n",
    "\n",
    "            #Right side\n",
    "            T[t,i,xdim-1] =(h/(k_bread+h_conv*h))*(h_conv*T_amb+T[t-1,i,xdim-2]*k_bread/h)\n",
    "        \n",
    "        #For filling\n",
    "        elif i>=breaddim and i < breaddim+filldim:\n",
    "            \n",
    "            # Left side\n",
    "            T[t,i,0] =(h/(k_chicken+h_conv*h))*(h_conv*T_amb+T[t-1,i,1]*k_chicken/h)\n",
    "\n",
    "            #Right side\n",
    "            T[t,i,xdim-1] =(h/(k_chicken+h_conv*h))*(h_conv*T_amb+T[t-1,i,xdim-2]*k_chicken/h)\n",
    "    \n",
    "\n",
    "    #Centre of sandwich\n",
    "    for i in range(1,ydim-1): # For every row in the matrix\n",
    "        for j in range(1,xdim-1): # For every column in the matrix\n",
    "            \n",
    "            #if inside bread region and inside bread boundaries\n",
    "            if i < breaddim-1 or i > breaddim + filldim:\n",
    "                T[t,i,j]=T[t-1,i,j]+courant_bread*(T[t-1,i+1,j]+T[t-1,i-1,j]+T[t-1,i,j+1]+T[t-1,i,j-1]-4*T[t-1,i,j])\n",
    "                #equivalent statement would be:\n",
    "                '''T[t,i,j]= (k*alpha_bread/(h**2)) * (T[t-1,i+1,j] - 2*T[t-1,i,j] + T[t-1,i-1,j] + \\\n",
    "                 T[t-1,i,j+1] - 2*T[t-1,i,j] + T[t-1,i,j-1]) + T[t-1,i,j]'''\n",
    "                \n",
    "            # if inside filling region and inside filling boundaries\n",
    "            if i > breaddim and i < breaddim+filldim-1:\n",
    "                T[t,i,j]=T[t-1,i,j]+courant_chicken*(T[t-1,i+1,j]+T[t-1,i-1,j]+T[t-1,i,j+1]+T[t-1,i,j-1]-4*T[t-1,i,j])\n",
    "                \n",
    "            # For undefined bread boundaries\n",
    "            if i==breaddim-1 or i==breaddim+filldim:\n",
    "                T[t,i,j]= (k*alpha_breadbc/(h**2)) * (T[t-1,i+1,j] - 2*T[t-1,i,j] + T[t-1,i-1,j] + \\\n",
    "                 T[t-1,i,j+1] - 2*T[t-1,i,j] + T[t-1,i,j-1]) + T[t-1,i,j]\n",
    "            \n",
    "            # For undefined filling boundaries\n",
    "            if i==breaddim or i==breaddim+filldim-1:\n",
    "                T[t,i,j]= (k*alpha_chickenbc/(h**2)) * (T[t-1,i+1,j] - 2*T[t-1,i,j] + T[t-1,i-1,j] + \\\n",
    "                 T[t-1,i,j+1] - 2*T[t-1,i,j] + T[t-1,i,j-1]) + T[t-1,i,j]\n",
    "    \n",
    "    #Update loading bar (only works in Jupyter)\n",
    "    print(\" \", end='\\r')\n",
    "\n",
    "print(\"Calculation completed successfully!\")\n",
    "print()\n",
    "\n",
    "#~~~~~~~~~~~~~~~~~~~~~~~~Generate Contour Plot~~~~~~~~~~~~~~~~~~~~~~~~~~~#\n",
    "\n",
    "plt.style.use('ggplot')\n",
    "\n",
    "print(\"Generating colour plot...\")\n",
    "\n",
    "H=1000*h\n",
    "# This step is required to reduce the effects of floating point error\n",
    "\n",
    "#create x and y arrays for plotting in millimetres\n",
    "x=np.arange(0, XL*1000+H, H)\n",
    "y=np.arange(0, YL*1000+H, H)\n",
    "# The multiplication by 1000 caused problems due to floating point error\n",
    "# (YL+h)*1000 multiplies floating point error more than YL*1000*H\n",
    "\n",
    "print(\"Final Temperature Distribtuion:\")\n",
    "\n",
    "# Plot x, y, final temperature matrix with the inferno colour scheme\n",
    "plt.contourf(x,y,T[-1,:,:], cmap='inferno')\n",
    "\n",
    "# Plot lines to indicate where bread and filling meet\n",
    "plt.plot([0, XL*1000],[bread_thickness*1000,bread_thickness*1000], color='gray', linestyle='dotted')\n",
    "plt.plot([0,XL*1000],[(bread_thickness+filling_thickness)*1000,(bread_thickness+filling_thickness)*1000], color='gray', linestyle='dotted')\n",
    "\n",
    "# Generate a colour bar\n",
    "plt.colorbar()\n",
    "\n",
    "plt.title(\"Temperature Distribution in Panini after \"+str(ttl_time)+\" seconds\", pad=15)\n",
    "\n",
    "# Define axes labels\n",
    "plt.xlabel(\"Length of panini [mm]\")\n",
    "plt.ylabel(\"Height of panini [mm]\")\n",
    "\n",
    "# Save then show figure\n",
    "plt.savefig('panini.png', dpi=300, bbox_inches='tight')\n",
    "plt.show()\n",
    "print() # Blank line\n",
    "\n",
    "#~~~~~~~~~~~~~~~~~~~~~~~~Generate Animation~~~~~~~~~~~~~~~~~~~~~~~~~~~#\n",
    "\n",
    "print(\"Generating animation...\")\n",
    "\n",
    "# Set up the figure and axes by calling plt.subplots\n",
    "# Also defining the figure size\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "\n",
    "#Set the axes range so that they do not auto scale during the animation\n",
    "ax.set(xlim=(0,(XL)*1000),ylim=(0,(YL)*1000))\n",
    "\n",
    "#Label the axes\n",
    "ax.set_xlabel(\"Length of panini [mm]\")\n",
    "ax.set_ylabel(\"Height of panini [mm]\")\n",
    "\n",
    "#Initialise the heatmap\n",
    "heatmap = ax.pcolormesh(x,y,T[0,:-1,:-1], vmin=initial_temperature, vmax=T_plate+10, cmap='inferno')\n",
    "# use a colourmesh, defining length, height and intial temperature matrix\n",
    "# Set the minimum & maximum temperature value on the colourscale\n",
    "# Use the colourscheme 'inferno' (because it looks coolest)\n",
    "\n",
    "#Generate the colourbar for the heatmap\n",
    "fig.colorbar(heatmap)\n",
    "\n",
    "\n",
    "# define a function that will be called later\n",
    "def animate(i):\n",
    "    heatmap.set_array(T[i,:-1,:-1].flatten())\n",
    "    #.flatten() will collapse a matrix into a 1D array\n",
    "    ax.set_title(\"Temperature Distribution in Panini at Time = \"+str(round(i*k,1))+\" s\")\n",
    "    # Apply a dynamic title, which will show the time to one decimal place\n",
    "    \n",
    "    # This plots lines indicating where bread and filling meet\n",
    "    ax.plot([0, XL*1000],[bread_thickness*1000,bread_thickness*1000], color='gray', linestyle='dotted')\n",
    "    ax.plot([0,XL*1000],[(bread_thickness+filling_thickness)*1000,(bread_thickness+filling_thickness)*1000], color='gray', linestyle='dotted')\n",
    "\n",
    "anim = FuncAnimation(fig, animate, interval=(k*1000/speed_factor), frames=len(T))\n",
    "# Use the predefined FuncAnimation\n",
    "# input the 'fig' defined in line 244\n",
    "# and input the animate func defined in line 263\n",
    "# this interval gives x5 faster video (argument is in milliseconds)\n",
    "print(\"Animation generated successfully!\", end='\\n')\n",
    "\n",
    "# This code will generate an mp4 file\n",
    "print(\"Generating mp4...\")\n",
    "anim.save('panini.mp4', dpi=150) # save as mp4\n",
    "print(\"Successfully saved as panini.mp4\", end='\\n')\n",
    "\n",
    "print() # Blank line\n",
    "\n",
    "# This code will generate a gif\n",
    "#print(\"Generating gif...\")\n",
    "#anim.save('panini.gif', writer='imagemagick', dpi=150) # save as gif\n",
    "#print(\"Successfully saved as panini.gif\")\n",
    "#print()\n",
    "\n",
    "#~~~~~~~~~~~~~~~~Deliver summary of info to user~~~~~~~~~~~~~~~~~#\n",
    "\n",
    "print(\"The panini was cooked for \"+str(ttl_time)+\" seconds at \"+str(T_plate)+\" degC\")\n",
    "\n",
    "# Calculate minimum temperature in final panini\n",
    "min_temp = np.min(T[-1])\n",
    "min_temp = round(min_temp, 2) # round the value\n",
    "\n",
    "# Calculate maximum temperature in the final panini\n",
    "flat_temps = T[-1,:,:].flatten()\n",
    "max_temp = initial_temperature\n",
    "for i in range(0,len(flat_temps)):\n",
    "    if flat_temps[i] > max_temp and flat_temps[i] != T_plate:\n",
    "        max_temp = flat_temps[i]\n",
    "\n",
    "print(\"The minimum temperature in the panini is \"+str(min_temp)+\" degC\")\n",
    "print(\"The maximum temperature in the panini is \"+str(round(max_temp,2))+\" degC\")\n",
    "\n",
    "# Deliver warnings to user:\n",
    "hotwarning = False\n",
    "coldwarning = True\n",
    "\n",
    "if max_temp > 50:\n",
    "    hotwarning = True\n",
    "\n",
    "if min_temp > 20:\n",
    "    coldwarning = False\n",
    "    \n",
    "if hotwarning == True and coldwarning == True:\n",
    "    print(\"Warning! The panini is really badly cooked. Risk of a terrible food experience.\")\n",
    "\n",
    "if hotwarning == False and coldwarning == True:\n",
    "    print(\"Warning! The panini is cold\")\n",
    "\n",
    "if hotwarning == True and coldwarning == False:\n",
    "    print(\"Warning! The panini is hot!\")\n",
    "\n",
    "#~~~~~~~~~~~~~~~~~~~~Celebrate~~~~~~~~~~~~~~~~~~#\n",
    "    \n",
    "print()\n",
    "\n",
    "open_image(\"panini.png\")\n",
    "\n",
    "startfile(\"panini.mp4\")\n",
    "\n",
    "open_image('Dont-Open-Surprise-Inside.jpg')\n",
    "\n",
    "print(\"Here's your panini. Enjoy!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}