Modeling Thermal Diffusivity II 20a.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling Thermal Diffusivity II\n",
    "## Versions\n",
    "- 19a \n",
    "  - Totally redefine approach. The main class is ThermalMesh with the temperature \n",
    "    stored in a 2D `numpy` array.\n",
    "  - `savePGM` to write an image file\n",
    "  - `setBorderAdiabatic` to set no heat flow out of the side\n",
    "  - Implemented _Newton's Law of Cooling_\n",
    "  - `iterate` using equation below\n",
    "  \n",
    "- 19b\n",
    "  - Add heat input to certain nodes\n",
    "  - Add constant Temperature nodes\n",
    "  \n",
    "- 19c\n",
    "  - Clean up code and documentation\n",
    "  - Change image & animation code\n",
    " \n",
    "- 20a\n",
    "  - Add explicity _template_ to class. A template describes what a node is\n",
    "    or does. A template must be set before calculation can be done Examples are:\n",
    "    - .: a normal temperature evolving node\n",
    "    - Q: constant heat input to node\n",
    "    - T: constant termperature node\n",
    "    - P: a periodic boundary node\n",
    "    - A: an adiabatic border node\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thermal Modeling\n",
    "\n",
    "## Definitions\n",
    "_Heat_ is the _transfer_ of thermal energy. Note that it is not the temperature of \n",
    "an object, or a property of an object, but a _change_ or movement of thermal energy. \n",
    "The main mechanisms of heat transfer are _conduction_ , _convection_ , and _radiation_ .\n",
    "\n",
    "\n",
    "## Heat Transfer Physics\n",
    "Heat is thermal energy associated with the temperature-dependent \n",
    "motion of particles. \n",
    "The macroscopic energy equation for infinitesimal volume used in \n",
    "heat transfer analysis is\n",
    "$$\n",
    "\\nabla \\cdot \\mathbf{q} = -c_p \\rho \\frac{\\partial T}{\\partial t} + \n",
    "\\sum_i{\\dot s_i}\n",
    "$$\n",
    "where $\\mathbf{q}$ is the _heat flux vector_ , $T$ is the temperature, $-c \\rho_p \\frac{\\partial T}{\\partial t}$ \n",
    "is the temporal change of internal energy ($\\rho$ is density, $c_p$ is  \n",
    "specific heat capacity at constant pressure, $T$ is temperature and $t$ is time.)\n",
    "\n",
    "## Fourier's Model of Heat Conduction\n",
    "Fourier's law of thermal conduction is\n",
    "$$\n",
    "\\mathbf{q} = -k \\nabla T\n",
    "$$\n",
    "where $k$ is the material's conductivity in $W/({m \\cdot K})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Heat Equation\n",
    "Working from the equation in _Heat Transfer Physics_ , you can derive the _Heat Equation_\n",
    "$$\n",
    "\\frac{\\partial T}{\\partial t} = \\frac{k c_p}{\\rho} \\nabla^2 T = \\alpha \\nabla^2 T\n",
    "$$\n",
    "where \n",
    "$$\\alpha = \\frac{k}{c_p \\rho}$$\n",
    "is the _thermal diffusivity_ of the medium. The table below shows some values \n",
    "for steel and aluminum.\n",
    "\n",
    "| Material | Thermal<br />Conductivity | Heat<br />Capacity | Density | $\\alpha$ |\n",
    "| --- | ---| --- | --- | --- |\n",
    "| - | ($W / (m \\cdot K$)) | ($J/(kg \\cdot K$) | ($kg/m^3$) | ($m^2/s$) |\n",
    "| Aluminum | 240 | 900 | 2700 |  80 |\n",
    "| Steel    |  35 | 480 | 7900 | 2.1 |\n",
    "\n",
    "### General Heat Loss\n",
    "If there are sources or sinks of thermal energy, the heat equation becomes\n",
    "$$\n",
    "\\rho c_p \\frac{\\partial T}{\\partial t} - k \\nabla^2 T = \\dot{q_V}\n",
    "$$\n",
    "where $\\dot{q_V}$ is the volumetric heat source (a negative $\\dot{q_V}$ is a _heat sink_.)\n",
    "\n",
    "\n",
    "You can model the general heat loss by assuming _Newton's Law of Cooling_\n",
    "$$\n",
    "\\frac{d T}{d t} = \\beta (T - T_0)\n",
    "$$\n",
    "where $\\beta$ is a constant and $T_0$ is the ambient termperature. The heat equation \n",
    "then becomes \n",
    "$$\n",
    "\\frac{\\partial T}{\\partial t} - \\alpha \\nabla^2 T = \\dot{q_V} -\\beta (T - T_0)\n",
    "$$\n",
    "Note that if $T > T_0$ there is cooling, and if $T < T_0$ there is warming."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applying the Heat Equation\n",
    "We want to apply the heat equation to simulate the time evolution of \n",
    "the temperature of a plate, so we can write the equation as\n",
    "$$\n",
    "\\frac{\\partial T}{\\partial t} = \\alpha \\nabla^2 T -\\beta (T - T_0) + \\dot q(x,y)\n",
    "$$\n",
    "where $\\dot q(x,y) $ is where heat is flowing into ($\\dot q > 0$) or out of \n",
    "($\\dot q < 0$) the plate. The $\\dot q$ term is used to model a heat gun pointed at \n",
    "a certain area of the plate, or maybe and ice pack against the plate.\n",
    "\n",
    "Another possibility is that you can have the plate in contact with a constant \n",
    "_temperature_ source, so the software ought to be able to handle that situation also.\n",
    "\n",
    "\n",
    "Now make the times discrete, so \n",
    "$$\n",
    "\\begin{align}\n",
    "T^{i+1} &= T^i + \\frac{\\partial T}{\\partial t} \\vert_{t^i} \\Delta t + \n",
    "           \\dot q(x,y) \\Delta t \\\\\n",
    "        & = T^i + \\left ( \\alpha \\nabla^2 T^i -\\beta (T^i - T_0) \\right ) \\Delta t \n",
    "           + \\dot q(x, y) \\Delta t\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Next we need to make the spatial coordinates, $x$ and $y$, discrete. \n",
    "We will simplify the grid by requiring it to be square so $\\Delta x = \\Delta y$.\n",
    "I will use the \n",
    "indices $j$ and $k$ for $x$ and $y$ respectively. \n",
    "The $\\nabla^2$ operator is basically the sum of the second \n",
    "derivatives in Cartesian coordinates\n",
    "$$\n",
    "\\nabla^2 = \\frac{\\partial^2 T}{\\partial x^2} + \\frac{\\partial^2 T}{\\partial y^2}\n",
    "$$\n",
    "When discretized, the second derivative becomes, in the $x$ direction\n",
    "$$\n",
    "\\frac{\\partial^2 T}{\\partial x^2}|_{t_i} \\approx \n",
    "  \\frac{T_{j+1,k}^i - 2 T_{j,k}^i + T_{j-1,k}^i}{{\\Delta x}^2}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the final result is\n",
    "$$\n",
    "T_{j,k}^{i+1} = T_{j,k}^i + \\frac{\\alpha \\Delta t}{{\\Delta x}^2}\n",
    "\\left({T_{j+1,k}^i + T_{j-1,k}^i + T_{j,k+1}^i + T_{j,k-1}^i - 4 T_{j,k}^i}\\right) \n",
    "- \\beta \\Delta t\\left({T_{j,k}^i - T_0}\\right) \n",
    "+ \\Delta t \\; \\dot q_{j,k}\n",
    "$$\n",
    "\n",
    "An important point to note is that _ALL_ of the right side terms are at time $t_i$ \n",
    "and the only left side term is at $t_{i+1}, so the spatial data at one time is \n",
    "used to calculate the spatial data at the next time. This is a huge simplifying \n",
    "factor in this model of thermal diffusion.\n",
    "\n",
    "All we have to do now is program this!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Implementation Hints\n",
    "From _Python Heat Transfer_ I found [here](http://www.cs.jhu.edu/~joanne/cs112/assign/p6.pdf).\n",
    "\n",
    "Use a 2D `numpy` array for the spatial temperature on a plate. Hints:\n",
    "1. _To implement __adiabatic__ edges_ (no heat flows \n",
    "off of the plate at the borders) add one row or column to the edge and give the values \n",
    "the _same_ values as the last row or column. If there is no temperature difference at the \n",
    "the edge, heat won't flow that direction. \n",
    "\n",
    "2. _To increase efficiency_ you can use the `np.roll` to rotate the array in a direction. \n",
    "then you can do the sum $T_{i,j+1,k} + T_{i,j-1,k} + T_{i,j,k+1} + T_{i,j,k-1}$ by rolling \n",
    "and adding. \n",
    "Since we are going to have extra rows and columns on the edges, we do not have \n",
    "to worry about the temerature at the bottom affecting the temperature at the top,\n",
    "And similarly for the other three sides. \n",
    "The missing elements in the corners don't matter, so define a rectangular array.\n",
    "\n",
    "BTW, if you want periodic boundary conditions, leave out the borders and the solution \n",
    "will smoothly roll across the boundaries.\n",
    "\n",
    "<img style=\"float: right;\" src=\"Modeling_Numpy_Hint.png\" alt=\"Figure\" width=\"400\"/>\n",
    "\n",
    "In the figure, \"B\" are border rows and columns used to implement (or not) an \n",
    "adiabatic edge to the plate. \n",
    "The contents of the boxes labeled \"C\" are arbitrary; they are not used in \n",
    "the calculation. The python coordinates of the corner pixels are shown.\n",
    "\n",
    "3. To add _Heat Sources and Sinks_, create both a template and \n",
    "a `heatSource` array the size of your \n",
    " array that adds (or subtracts) heat from particular cells. \n",
    " The non-zero values are the amount of heat added or subtracted from a\n",
    " given node.\n",
    " We are assuming \n",
    " the heat sources or sinks are constant with time. At each time step, the \n",
    " temperature \n",
    " of those cells is changed by the value $\\dot q_{j,k} \\Delta t$.\n",
    "\n",
    "4. To implement constant temperature nodes, modify the template and create \n",
    " an array the size of the nodes array where the non-zero values are the \n",
    " fixed temperatures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#23456789013456789013456789013456789013456789013456789013456789072\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from copy import copy, deepcopy\n",
    "from PIL import Image\n",
    "from ThermalMesh import ThermalMesh\n",
    "import pathlib\n",
    "# from ThermalMesh9Pt import ThermalMesh\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation Set-Up\n",
    "You will see in the next cell, that is takes some careful work to set up \n",
    "a simulation. You have to:\n",
    "1. Define the physics parameters.\n",
    "  1. The _thermal diffusivity_ is the main one.\n",
    "1. Define the size and resolution of the simulation.\n",
    "  1. The spatial spacing, `DX`.\n",
    "  1. The mesh size, `X_SIZE` and `Y_SIZE`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DeltaT = 0.000500 s\n",
      "The Mesh Fourier Number, MFN, is 0.0250\n",
      "nIterations = 40000\n",
      "With DX = 0.02 m and a 40x30 array the physical size is 0.8 x 0.6 m.\n",
      "template:\n",
      "   012345678911234567892123456789312345678941\n",
      " 0 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\n",
      " 1 A........................................A\n",
      " 2 A........................................A\n",
      " 3 A........................................A\n",
      " 4 A........................................A\n",
      " 5 A........................................A\n",
      " 6 A........................................A\n",
      " 7 A........................................A\n",
      " 8 A........................................A\n",
      " 9 A........................................A\n",
      "10 A........................................A\n",
      "11 A........................................A\n",
      "12 A........................................A\n",
      "13 A........................................A\n",
      "14 A........................................A\n",
      "15 A........................................A\n",
      "16 A........................................A\n",
      "17 A........................................A\n",
      "18 A........................................A\n",
      "19 A........................................A\n",
      "20 A........................................A\n",
      "21 A........................................A\n",
      "22 A........................................A\n",
      "23 A........................................A\n",
      "24 A........................................A\n",
      "25 A........................................A\n",
      "26 A...................................Q....A\n",
      "27 A........................................A\n",
      "28 A........................................A\n",
      "29 A........................................A\n",
      "30 A....T...................................A\n",
      "31 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\n"
     ]
    }
   ],
   "source": [
    "testPrint = False\n",
    "\n",
    "# Physical parameters\n",
    "# The thermal diffusivity\n",
    "alpha_Al = 80.0\n",
    "alpha_Steel = 2.1\n",
    "alpha_test = 0.02\n",
    "# Newton's cooling coefficient\n",
    "BETA = 0.02\n",
    "\n",
    "# Simulation size\n",
    "DX = 0.02 # mesh size in meters\n",
    "timeSim = 20.0 # Simulation for 20 seconds\n",
    "# Note, that this particular simulation turns off the constant \n",
    "# temperature and heat nodes and lets the plate cool for \n",
    "# `timeSim`\n",
    "\n",
    "# A very important parameter in simluating this is the\n",
    "# *Mesh Fourier Number* = MFN ≡ α∆t/(∆x)2.\n",
    "# For this algorithm to be stable MNF < 0.25, so I will set it to 0.25/10\n",
    "DeltaT = (.25/10) * DX**2 / alpha_test\n",
    "print(f\"DeltaT = {DeltaT:.6f} s\")\n",
    "MFN = alpha_test * DeltaT / DX**2\n",
    "print(f\"The Mesh Fourier Number, MFN, is {MFN:.4f}\")\n",
    "nIterations = int(timeSim / DeltaT + 0.5)\n",
    "print(f\"nIterations = {nIterations}\")\n",
    "\n",
    "X_SIZE = 40\n",
    "Y_SIZE = 30\n",
    "print(f\"With DX = {DX} m and a {X_SIZE}x{Y_SIZE} array the \\\n",
    "physical size is {X_SIZE*DX} x {Y_SIZE*DX} m.\")\n",
    "\n",
    "# Location and values for one constant heat node and \n",
    "# one contant temperature node.\n",
    "X_Q = 35 # 0 based index relative to array, not borders\n",
    "Y_Q = 25\n",
    "Q_VALUE = 2532.1\n",
    "X_T= 4\n",
    "Y_T = 29\n",
    "T_VALUE = 100.0\n",
    "\n",
    "# Ambient and initial temperature\n",
    "TAmb = 25.0\n",
    "initT = 27.0\n",
    "\n",
    "# Iteration parameters\n",
    "nImages = 100 # Number (actual the same number are computed for cooling)\n",
    "nImageSkip = nIterations // nImages\n",
    "# Name formats for the image files\n",
    "pgmCQSName = r\"Images/MTDII%04d.pgm\"\n",
    "pngCQSName = r\"Images/MTDII%04d.png\"\n",
    "\n",
    "# Build the template\n",
    "tempRow = list((X_SIZE + 2) * '.')\n",
    "template = []\n",
    "for i in range(Y_SIZE+2):\n",
    "    template.append(copy(tempRow))\n",
    "\n",
    "# Add adiabatic boundaries\n",
    "aRow = []\n",
    "for i in range(X_SIZE+2):\n",
    "    aRow.append('A')\n",
    "template[0] = copy(aRow)\n",
    "template[-1] = copy(aRow)\n",
    "# Set columns\n",
    "for i in range(Y_SIZE + 2):\n",
    "    template[i][0] = 'A'\n",
    "    template[i][-1] = 'A'\n",
    "\n",
    "# Add a constant heat point\n",
    "template[Y_Q+1][X_Q+1] = 'Q'\n",
    "\n",
    "# Add a constant temperature point\n",
    "template[Y_T+1][X_T+1] = 'T'\n",
    "\n",
    "# Compress template into strings\n",
    "for i in range(Y_SIZE + 2):\n",
    "    template[i] = \"\".join(template[i])\n",
    "\n",
    "# Pretty print the template array\n",
    "print(f\"template:\")\n",
    "print(\"   012345678911234567892123456789312345678941\")\n",
    "for i,row in enumerate(template):\n",
    "    print(f\"{i:2d}\",row)\n",
    "\n",
    "\n",
    "# Make contant Q array\n",
    "CQSources = np.zeros((Y_SIZE+2, X_SIZE+2), dtype=float)\n",
    "CQSources[Y_Q+1, X_Q+1] = Q_VALUE\n",
    "if testPrint:\n",
    "    print(\"CQSources\")\n",
    "    for i in range(Y_SIZE+2):\n",
    "        print(CQSources[i,:])\n",
    "\n",
    "# Make contant temperature array\n",
    "ConstT = np.zeros((Y_SIZE+2, X_SIZE+2), dtype=float)\n",
    "ConstT[Y_T+1, X_T+1] = T_VALUE\n",
    "if testPrint:\n",
    "    print(\"ConstT\")\n",
    "    for i in range(Y_SIZE+2):\n",
    "        print(ConstT[i,:])\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now make `myMesh` and Iterate\n",
    "Also, after the the simulation is done, let it cool for the same \n",
    "amount of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "setTemplate:  True True True True\n",
      "Mesh(40,30): \n",
      "MFN = 0.025, beta = 0.02, TAmb = 25.0\n",
      "[[27. 27. 27. ... 27. 27. 27.]\n",
      " [27. 27. 27. ... 27. 27. 27.]\n",
      " [27. 27. 27. ... 27. 27. 27.]\n",
      " ...\n",
      " [27. 27. 27. ... 27. 27. 27.]\n",
      " [27. 27. 27. ... 27. 27. 27.]\n",
      " [27. 27. 27. ... 27. 27. 27.]]\n",
      "S T A R T   I T E R A T I O N\n",
      "After iteration:\n",
      " Mesh(40,30): \n",
      "MFN = 0.025, beta = 0.02, TAmb = 25.0\n",
      "[[26.99998 26.99998 26.99998 ... 26.99998 26.99998 26.99998]\n",
      " [26.99998 26.99998 26.99998 ... 26.99998 26.99998 26.99998]\n",
      " [26.99998 26.99998 26.99998 ... 26.99998 26.99998 26.99998]\n",
      " ...\n",
      " [26.99998 26.99998 26.99998 ... 26.99998 26.99998 26.99998]\n",
      " [26.99998 26.99998 26.99998 ... 26.99998 26.99998 26.99998]\n",
      " [26.99998 26.99998 26.99998 ... 26.99998 26.99998 26.99998]]\n",
      "(30, 40)\n",
      "Writing Images/MTDII0000.pgm, Writing Images/MTDII0001.pgm, Writing Images/MTDII0002.pgm, \n",
      "Writing Images/MTDII0003.pgm, Writing Images/MTDII0004.pgm, Writing Images/MTDII0005.pgm, \n",
      "Writing Images/MTDII0006.pgm, Writing Images/MTDII0007.pgm, Writing Images/MTDII0008.pgm, \n",
      "Writing Images/MTDII0009.pgm, Writing Images/MTDII0010.pgm, Writing Images/MTDII0011.pgm, \n",
      "Writing Images/MTDII0012.pgm, Writing Images/MTDII0013.pgm, Writing Images/MTDII0014.pgm, \n",
      "Writing Images/MTDII0015.pgm, Writing Images/MTDII0016.pgm, Writing Images/MTDII0017.pgm, \n",
      "Writing Images/MTDII0018.pgm, Writing Images/MTDII0019.pgm, Writing Images/MTDII0020.pgm, \n",
      "Writing Images/MTDII0021.pgm, Writing Images/MTDII0022.pgm, Writing Images/MTDII0023.pgm, \n",
      "Writing Images/MTDII0024.pgm, Writing Images/MTDII0025.pgm, Writing Images/MTDII0026.pgm, \n",
      "Writing Images/MTDII0027.pgm, Writing Images/MTDII0028.pgm, Writing Images/MTDII0029.pgm, \n",
      "Writing Images/MTDII0030.pgm, Writing Images/MTDII0031.pgm, Writing Images/MTDII0032.pgm, \n",
      "Writing Images/MTDII0033.pgm, Writing Images/MTDII0034.pgm, Writing Images/MTDII0035.pgm, \n",
      "Writing Images/MTDII0036.pgm, Writing Images/MTDII0037.pgm, Writing Images/MTDII0038.pgm, \n",
      "Writing Images/MTDII0039.pgm, Writing Images/MTDII0040.pgm, Writing Images/MTDII0041.pgm, \n",
      "Writing Images/MTDII0042.pgm, Writing Images/MTDII0043.pgm, Writing Images/MTDII0044.pgm, \n",
      "Writing Images/MTDII0045.pgm, Writing Images/MTDII0046.pgm, Writing Images/MTDII0047.pgm, \n",
      "Writing Images/MTDII0048.pgm, Writing Images/MTDII0049.pgm, Writing Images/MTDII0050.pgm, \n",
      "Writing Images/MTDII0051.pgm, Writing Images/MTDII0052.pgm, Writing Images/MTDII0053.pgm, \n",
      "Writing Images/MTDII0054.pgm, Writing Images/MTDII0055.pgm, Writing Images/MTDII0056.pgm, \n",
      "Writing Images/MTDII0057.pgm, Writing Images/MTDII0058.pgm, Writing Images/MTDII0059.pgm, \n",
      "Writing Images/MTDII0060.pgm, Writing Images/MTDII0061.pgm, Writing Images/MTDII0062.pgm, \n",
      "Writing Images/MTDII0063.pgm, Writing Images/MTDII0064.pgm, Writing Images/MTDII0065.pgm, \n",
      "Writing Images/MTDII0066.pgm, Writing Images/MTDII0067.pgm, Writing Images/MTDII0068.pgm, \n",
      "Writing Images/MTDII0069.pgm, Writing Images/MTDII0070.pgm, Writing Images/MTDII0071.pgm, \n",
      "Writing Images/MTDII0072.pgm, Writing Images/MTDII0073.pgm, Writing Images/MTDII0074.pgm, \n",
      "Writing Images/MTDII0075.pgm, Writing Images/MTDII0076.pgm, Writing Images/MTDII0077.pgm, \n",
      "Writing Images/MTDII0078.pgm, Writing Images/MTDII0079.pgm, Writing Images/MTDII0080.pgm, \n",
      "Writing Images/MTDII0081.pgm, Writing Images/MTDII0082.pgm, Writing Images/MTDII0083.pgm, \n",
      "Writing Images/MTDII0084.pgm, Writing Images/MTDII0085.pgm, Writing Images/MTDII0086.pgm, \n",
      "Writing Images/MTDII0087.pgm, Writing Images/MTDII0088.pgm, Writing Images/MTDII0089.pgm, \n",
      "Writing Images/MTDII0090.pgm, Writing Images/MTDII0091.pgm, Writing Images/MTDII0092.pgm, \n",
      "Writing Images/MTDII0093.pgm, Writing Images/MTDII0094.pgm, Writing Images/MTDII0095.pgm, \n",
      "Writing Images/MTDII0096.pgm, Writing Images/MTDII0097.pgm, Writing Images/MTDII0098.pgm, \n",
      "Writing Images/MTDII0099.pgm, Writing Images/MTDII0100.pgm, \n",
      "\n",
      "### Now Start COOLING ###\n",
      "Writing Images/MTDII0101.pgm, \n",
      "Writing Images/MTDII0102.pgm, Writing Images/MTDII0103.pgm, Writing Images/MTDII0104.pgm, \n",
      "Writing Images/MTDII0105.pgm, Writing Images/MTDII0106.pgm, Writing Images/MTDII0107.pgm, \n",
      "Writing Images/MTDII0108.pgm, Writing Images/MTDII0109.pgm, Writing Images/MTDII0110.pgm, \n",
      "Writing Images/MTDII0111.pgm, Writing Images/MTDII0112.pgm, Writing Images/MTDII0113.pgm, \n",
      "Writing Images/MTDII0114.pgm, Writing Images/MTDII0115.pgm, Writing Images/MTDII0116.pgm, \n",
      "Writing Images/MTDII0117.pgm, Writing Images/MTDII0118.pgm, Writing Images/MTDII0119.pgm, \n",
      "Writing Images/MTDII0120.pgm, Writing Images/MTDII0121.pgm, Writing Images/MTDII0122.pgm, \n",
      "Writing Images/MTDII0123.pgm, Writing Images/MTDII0124.pgm, Writing Images/MTDII0125.pgm, \n",
      "Writing Images/MTDII0126.pgm, Writing Images/MTDII0127.pgm, Writing Images/MTDII0128.pgm, \n",
      "Writing Images/MTDII0129.pgm, Writing Images/MTDII0130.pgm, Writing Images/MTDII0131.pgm, \n",
      "Writing Images/MTDII0132.pgm, Writing Images/MTDII0133.pgm, Writing Images/MTDII0134.pgm, \n",
      "Writing Images/MTDII0135.pgm, Writing Images/MTDII0136.pgm, Writing Images/MTDII0137.pgm, \n",
      "Writing Images/MTDII0138.pgm, Writing Images/MTDII0139.pgm, Writing Images/MTDII0140.pgm, \n",
      "Writing Images/MTDII0141.pgm, Writing Images/MTDII0142.pgm, Writing Images/MTDII0143.pgm, \n",
      "Writing Images/MTDII0144.pgm, Writing Images/MTDII0145.pgm, Writing Images/MTDII0146.pgm, \n",
      "Writing Images/MTDII0147.pgm, Writing Images/MTDII0148.pgm, Writing Images/MTDII0149.pgm, \n",
      "Writing Images/MTDII0150.pgm, Writing Images/MTDII0151.pgm, Writing Images/MTDII0152.pgm, \n",
      "Writing Images/MTDII0153.pgm, Writing Images/MTDII0154.pgm, Writing Images/MTDII0155.pgm, \n",
      "Writing Images/MTDII0156.pgm, Writing Images/MTDII0157.pgm, Writing Images/MTDII0158.pgm, \n",
      "Writing Images/MTDII0159.pgm, Writing Images/MTDII0160.pgm, Writing Images/MTDII0161.pgm, \n",
      "Writing Images/MTDII0162.pgm, Writing Images/MTDII0163.pgm, Writing Images/MTDII0164.pgm, \n",
      "Writing Images/MTDII0165.pgm, Writing Images/MTDII0166.pgm, Writing Images/MTDII0167.pgm, \n",
      "Writing Images/MTDII0168.pgm, Writing Images/MTDII0169.pgm, Writing Images/MTDII0170.pgm, \n",
      "Writing Images/MTDII0171.pgm, Writing Images/MTDII0172.pgm, Writing Images/MTDII0173.pgm, \n",
      "Writing Images/MTDII0174.pgm, Writing Images/MTDII0175.pgm, Writing Images/MTDII0176.pgm, \n",
      "Writing Images/MTDII0177.pgm, Writing Images/MTDII0178.pgm, Writing Images/MTDII0179.pgm, \n",
      "Writing Images/MTDII0180.pgm, Writing Images/MTDII0181.pgm, Writing Images/MTDII0182.pgm, \n",
      "Writing Images/MTDII0183.pgm, Writing Images/MTDII0184.pgm, Writing Images/MTDII0185.pgm, \n",
      "Writing Images/MTDII0186.pgm, Writing Images/MTDII0187.pgm, Writing Images/MTDII0188.pgm, \n",
      "Writing Images/MTDII0189.pgm, Writing Images/MTDII0190.pgm, Writing Images/MTDII0191.pgm, \n",
      "Writing Images/MTDII0192.pgm, Writing Images/MTDII0193.pgm, Writing Images/MTDII0194.pgm, \n",
      "Writing Images/MTDII0195.pgm, Writing Images/MTDII0196.pgm, Writing Images/MTDII0197.pgm, \n",
      "Writing Images/MTDII0198.pgm, Writing Images/MTDII0199.pgm, Writing Images/MTDII0200.pgm, \n",
      "Mesh(40,30): \n",
      "MFN = 0.025, beta = 0.02, TAmb = 25.0\n",
      "[[66.16737112 66.16737112 66.16738182 ... 66.17082991 66.1708406\n",
      "  66.1708406 ]\n",
      " [66.16737112 66.16737112 66.16738182 ... 66.17082991 66.1708406\n",
      "  66.1708406 ]\n",
      " [66.16737118 66.16737118 66.16738188 ... 66.17082997 66.17084067\n",
      "  66.17084067]\n",
      " ...\n",
      " [66.16738247 66.16738247 66.16739317 ... 66.17084126 66.17085196\n",
      "  66.17085196]\n",
      " [66.16738253 66.16738253 66.16739323 ... 66.17084132 66.17085202\n",
      "  66.17085202]\n",
      " [66.16738253 66.16738253 66.16739323 ... 66.17084132 66.17085202\n",
      "  66.17085202]]\n",
      "#### DONE ####\n"
     ]
    }
   ],
   "source": [
    "myMesh = ThermalMesh( nX=X_SIZE, nY=Y_SIZE, dx=DX, dt=DeltaT, MFN=MFN)\n",
    "myMesh.setTemplate(template)\n",
    "myMesh.setInitT(initT)\n",
    "myMesh.setBeta(BETA)\n",
    "myMesh.setTAmb(TAmb)\n",
    "myMesh.setCQSources(CQSources)\n",
    "myMesh.setConstT(ConstT)\n",
    "print(myMesh)\n",
    "\n",
    "print(\"S T A R T   I T E R A T I O N\")\n",
    "myMesh.iterate()\n",
    "print(\"After iteration:\\n\", myMesh)\n",
    "\n",
    "# Now iterate until done\n",
    "iIteration = 0\n",
    "iImage = 0\n",
    "time = 0.0\n",
    "\n",
    "# Arrays to save so I can make images\n",
    "iterArrays = [myMesh.nodes[1:1+myMesh.nY,1:1+myMesh.nX]]\n",
    "print(iterArrays[0].shape)\n",
    "\n",
    "while time < timeSim:\n",
    "    if iIteration % nImageSkip == 0:\n",
    "        print(\"Writing\", pgmCQSName%(iImage), end=', ')\n",
    "        iterArrays.append(myMesh.nodes[1:1+myMesh.nY,1:1+myMesh.nX])\n",
    "        myMesh.savePGM(pgmCQSName%(iImage))\n",
    "        myMesh.savePNG(pngCQSName%(iImage))\n",
    "        if iImage % 3 == 2:\n",
    "            print()\n",
    "        iImage += 1\n",
    "        \n",
    "    myMesh.iterate()\n",
    "    iIteration += 1\n",
    "    time += DeltaT\n",
    "\n",
    "# Now let it cool\n",
    "print(\"\\n\\n### Now Start COOLING ###\")\n",
    "del myMesh.CQSources\n",
    "del myMesh.constT\n",
    "time = 0.0\n",
    "while time < timeSim:\n",
    "    if iIteration % nImageSkip == 0:\n",
    "        print(\"Writing\", pgmCQSName%(iImage), end=', ')\n",
    "        iterArrays.append(myMesh.nodes[1:1+myMesh.nY,1:1+myMesh.nX])\n",
    "        myMesh.savePGM(pgmCQSName%(iImage))\n",
    "        myMesh.savePNG(pngCQSName%(iImage))\n",
    "        if iImage % 3 == 2:\n",
    "            print()\n",
    "        iImage += 1\n",
    "        \n",
    "    myMesh.iterate()\n",
    "    iIteration += 1\n",
    "    time += DeltaT\n",
    "\n",
    "print(myMesh)\n",
    "print(\"#### DONE ####\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DD OO NN EE\n"
     ]
    }
   ],
   "source": [
    "print(\"DD OO NN EE\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save as Images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " # # D O N E # #\n"
     ]
    }
   ],
   "source": [
    "normEach = False\n",
    "paletteFile = \"IronPalette.csv\"\n",
    "imgFormat = r\"Images/MTDII%04d.gif\"\n",
    "animFile = \"MTDII_Animated.gif\"\n",
    "durationGif = 25\n",
    "\n",
    "# Get color map\n",
    "palette = np.genfromtxt(paletteFile, dtype='uint8', delimiter=\",\")\n",
    "                        \n",
    "# Keep originals protected\n",
    "itArr = deepcopy(iterArrays)\n",
    "gifList = []\n",
    "\n",
    "# Norm the images\n",
    "if normEach:\n",
    "    for i in range(len(itArr)):\n",
    "        img = copy(itArr[i])\n",
    "        min0 = np.min(img)\n",
    "        max0 = np.max(img)\n",
    "        print(min0, ' -> ', max0)\n",
    "        if min0 == max0:\n",
    "            img[:] = 0\n",
    "        else:\n",
    "            img = 255 * (img - min0)/(max0 - min0)\n",
    "        Image0 = Image.fromarray(img, mode='L')\n",
    "        Image0.putpalette(palette)\n",
    "        gifList.append(Image0)\n",
    "        \n",
    "        gifFile = imgFormat % (i)\n",
    "        # print(f\"{gifFile} from {np.min(img)} -> {np.max(img)}\")\n",
    "        Image0.save(gifFile)\n",
    "else:\n",
    "    # All images normed to the same absolute min & max\n",
    "    \n",
    "    min0 = 1e30\n",
    "    max0 = -1e30\n",
    "    for i in range(len(itArr)):\n",
    "        min0 = min(min0, np.min(itArr[i]))\n",
    "        max0 = max(max0, np.max(itArr[i]))\n",
    "    for i in range(len(itArr)):\n",
    "        img = copy(itArr[i])\n",
    "        img = 255 * (img - min0)/(max0 - min0)\n",
    "        img = np.array(img, dtype='uint8')\n",
    "        Image0 = Image.fromarray(img, mode='L')\n",
    "        Image0.putpalette(palette)\n",
    "        gifList.append(Image0)\n",
    "        \n",
    "        gifFile = imgFormat % (i)\n",
    "        # print(f\"{gifFile} from {np.min(img)} -> {np.max(img)}\")\n",
    "        Image0.save(gifFile)\n",
    "\n",
    "gifList[0].save(animFile, save_all=True, append_images=gifList[1:], \\\n",
    "               duration=durationGif, loop=0)\n",
    "\n",
    "print(\" # # D O N E # #\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final GIF\n",
    "<img src=\"MTDII_Animated.gif\" alt=\"Figure\" width=\"400\"/>\n",
    "<img src=\"http://i.stack.imgur.com/SBv4T.gif\" alt=\"this slowpoke moves\"  width=250/>"
   ]
  },
  {
   "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
}