{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PDE session 5\n",
    "## Parabolic PDE, stability and convergence of explicit methods"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAACQCAYAAABZL1J4AAAAAXNSR0IArs4c6QAAAERlWElmTU0AKgAAAAgAAYdpAAQAAAABAAAAGgAAAAAAA6ABAAMAAAABAAEAAKACAAQAAAABAAACEKADAAQAAAABAAAAkAAAAACAZ6rjAAALQ0lEQVR4Ae3dXYhN3R8H8OX9Jc0o5F0uuVNcUV4iF26UCxempChX4oorebuRW0lJ0SQpochbFJNIorgSJYPczY0UM4nn/6z995+YkT2rf+ucOWs+u6bm7L3Xs9b6/Fbn+Tqzz96j/vl3CzYCBAgQIECAQILA6IRznUqAAAECBAgQqAQECAuBAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYQIJLJNCBAgAABAgQECGuAAAECBAgQSBYYm9xiBDW4f/9+iD82AgQIEBieAqtXrw7xx9Z4AQHiL+YxPBw6dOgvZzhEgAABAs0WECCaUwEBYgjuq1atknCH4OQUAgQINEog/gOvq6urUd3p5w8CAsQfUAbuiun24MGDA3d7TYAAAQJNEojvyQJEk/B/dusiyub6650AAQIECLSkgADRkmUzaAIECBAg0FwBAaK5/nonQIAAAQItKSBAtGTZDJoAAQIECDRXQIBorr/eCRAgQIBASwoIEC1ZNoMmQIAAAQLNFRAgmuuvdwIECBBoAYFjx46FzZs3hxs3boTv37+3wIjzD1GAyG+sBwIECBBocYHdu3eHtWvXhiNHjoR58+aFvXv3hpcvX7b4rP6/4buR1BD84h3P3EhqCFBOIUCAQIMEGv2cogkTJoSdO3dWP69evQpnz56tAkUME9u2bQsdHR2hvb29QbMfHt2M+uffbXgMZfiNIoYGz8IYfnUxIgIECPxJYPTo0eHHjx+/HRq4L+X1r+fGOxLXhZb4v9Oenp5w6tSp0NbWFh4/fhyOHj0aLl++HMaPHx9u3rwZTp48GebMmfPbGFv1hU8g/lK5uGBigDhw4MBfznKIAAECBJohEN+f7927F+J7dSO3+AnEmTNnQmdnZ5g/f371CcSWLVtCDBDHjx8P+/fvD2PHjg1fvnwJK1asqIJDd3d3uHr1ati1a1cxAcInEDWrbtSoUdWiqDnNYQIECBBosEAj35/7+vqq0BD/dPHu3buwdevWKjgsXry4f9bXr18Ps2bNCkuXLq327du3Lzx69Cg8ePCgChNPnjxpeNjpH1yGXwSIGtRGLtCaoThMgAABAr8INPL9OX4L49mzZ1VoWL9+fRgzZswvI/nvr9++fQvjxo3r3x+f5Bx/Dh8+3L+vpF8EiJpqNnKB1gzFYQIECBD4RWA4vz/39vaGqVOnhvipRPz2Ronb6BInZU4ECBAgQKCZAg8fPqz+/L18+fL+YdRdhNl/Yov8IkC0SKEMkwABAgSGt8CePXuq+0TEUV67di0sWrQoTJo0qRr0x48fw5s3b4b3BBJHJ0AkgjmdAAECBAj8SeDixYth9uzZ1bUS06ZNC/GaiM+fP4f379+H06dPV9dP/Kldq+5zDURN5Ybz39hqhu4wAQIEihYYbu/Pr1+/Dnfv3g0LFy4MGzZsCG/fvq1ufT1jxoywadOm6qudJRVEgKip5nBboDXDdZgAAQIjRsD7c3NL7U8YzfXXOwECBAgQaEkBAaIly2bQBAgQIECguQICRI2/h2jVADlMgACBJgl4f24S/M9uXQPRXH+9EyBAgACBlhTwCURLls2gCRAgQIBAcwUEiOb6650AAQIECLSkgMd5D6Fs8clr8Znu8YYg8fnw8UYh8RHfkydPHkJrpxAgQIBAboH4qOzbt2+HEydO5O7Kf/+ngABRsxQ+ffoUduzYEc6dOxdmzpxZnf306dPQ0dERrly5UtPaYQIECBDIKdDV1RXOnz8f2traqhs35ezLf/t3AQHid49Bry5cuFA99/1/4SGesGzZsjB37tzw/PnzsGTJkkFt7CBAgACBxgjE9+CVK1eGeFOpjRs3NqZTvVQCroGoWQgxJMTnuQ/c4r4XL14M3O01AQIECDRQoL29vQoPDexSVz8FBIiapdDT0xPifcwHbtOnTw/xmI0AAQIECIxEAQGipuoTJ04Mvb29g87q6+sLEyZMGLTfDgIECBAgMBIEBIiaKsdrH+K3MAZu3d3d/RdVDjzmNQECBAgQKF1AgKip8Jo1a8KtW7cGnXXnzp3qwp1BB+wgQIAAAQIjQECAqCnyunXrwqVLl8KHDx/6z4zPe49X/P76zYz+g34hQIAAAQIjQMDXOGuKHK9z6OzsDNu3bw8LFiwIX79+rX7Onj1b09JhAgQIEGikQPyHna1xAh6mlWAdv3URL6qcMmVKQiunEiBAgACB8gQEiPJqakYECBAgQCC7gGsgshPrgAABAgQIlCcgQJRXUzMiQIAAAQLZBQSI7MQ6IECAAAEC5QkIEOXV1IwIECBAgEB2AQEiO7EOCBAgQIBAeQICRHk1NSMCBAgQIJBdQIDITqwDAgQIECBQnoAAUV5NzYgAAQIECGQXECCyE+uAAAECBAiUJyBAlFdTMyJAgAABAtkFBIjsxDogQIAAAQLlCQgQ5dXUjAgQIECAQHYBASI7sQ4IECBAgEB5AgJEeTU1IwIECBAgkF1AgMhOrAMCBAgQIFCegABRXk3NiAABAgQIZBcQILIT64AAAQIECJQnIECUV1MzIkCAAAEC2QUEiOzEOiBAgAABAuUJCBDl1dSMCBAgQIBAdgEBIjuxDggQIECAQHkCAkR5NTUjAgQIECCQXUCAyE6sAwIECBAgUJ6AAFFeTc2IAAECBAhkFxAgshPrgAABAgQIlCcgQJRXUzMiQIAAAQLZBQSI7MQ6IECAAAEC5QkIEOXV1IwIECBAgEB2AQEiO7EOCBAgQIBAeQICRHk1NSMCBAgQIJBdQIDITqwDAgQIECBQnoAAUV5NzYgAAQIECGQXECCyE+uAAAECBAiUJyBAlFdTMyJAgAABAtkFBIjsxDogQIAAAQLlCQgQ5dXUjAgQIECAQHYBASI7sQ4IECBAgEB5AgJEeTU1IwIECBAgkF1AgMhOrAMCBAgQIFCegABRXk3NiAABAgQIZBcQILIT64AAAQIECJQnIECUV1MzIkCAAAEC2QUEiOzEOiBAgAABAuUJCBDl1dSMCBAgQIBAdoH/AKrYibbQSPPVAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this exercise we are looking at the temperature distribution in a 1D bar. In 1D, the heat conduction problem is represented by:\n",
    "\n",
    "$$ \\frac{\\partial T}{\\partial t} = \\alpha \\frac{\\partial^2 T}{\\partial x^2}  \\tag{1}$$\n",
    "\n",
    "The underlying problem is illustrated below where we are interested by the temperature distribution in a bar of length 1m within the timeframe [0,1s]. The parameter $\\alpha$ should take the value 2.3e-1 $m^2/s$\n",
    "\n",
    "![bar](attachment:image.png)\n",
    "\n",
    "To complete our problem we need to define initial conditions in addition to boundary conditions.\n",
    "\n",
    "As boundary conditions we will enforce that:\n",
    "\n",
    "1. $\\partial T / \\partial x = q(t) $ for $x = 0$  and for all instance of time $t$\n",
    "2. $\\partial T / \\partial x = 0$ for $x = 1 \\text{ m}$  and for all instance of time $t$\n",
    "\n",
    "In addition we will use as initial condition:\n",
    "\n",
    "3. $T = 25^\\circ\\text{ deg} $ for $ 0 \\ge x \\ge1 \\text{ m}$ at time $t = 0$\n",
    "\n",
    "The flux on the left-hand side of the bar is given as function of time $t$ and has the following expression:\n",
    "\n",
    "$$ q(t) = -\\left(q_0 + \\frac{q_m}{t_s\\sqrt{2 \\pi}}\\exp{\\left(-\\frac{1}{2} \\left(\\frac{t-t_m}{t_s} \\right)^2 \\right)} \\right) \\tag{2}$$\n",
    "\n",
    "where $q_0 = 25$, $q_m = 200$, $t_m = 0.4$ and $t_s = 0.1$\n",
    "\n",
    "To solve this problem, we will use the Dufort-Frankel scheme. As discussed in the lecture, this scheme has three levels in time and we thus need some special treatment when we start the computations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this exercise, you have to:\n",
    "\n",
    "0. Discretize the PDE using the Dufort-Frankel scheme.\n",
    "1. Import the necessary python modules\n",
    "2. Implement a function which computes the flux on the left-hand side of the bar. As input we need the time only. The parameters governing equation (2) should be defined inside this function.\n",
    "3. Implement a function which solve the PDE using the Dufort-Frankel scheme. For the first time increment you should use the FTCS scheme. As input we want the initial conditions as well as the number of points along the x axis and the number of time steps. The Fourier number $\\gamma$ (defined as $\\gamma = \\alpha \\Delta t / \\Delta x^2$) will be computed inside this function, and we should then supply the timestep, the grid size and the thermal diffusivity parameter $\\alpha$ in addition. This function should return the complete history of the temperature $T$ (i.e. a matrix).\n",
    "4. Run the problem with 100 grid points along $x$ and a timesteps $\\Delta t$ of 0.0001 s. Plot the distribution of the temperature along $x$ for a time of approximately 0.1s, 0.5s and 1s. You should also plot the evolution of the temperature over time in the start, the middle and the end of the bar. Use subplots to visualize theses plots side-by-side.\n",
    "5. Run the problem with 100 grid points along $x$ and a timesteps $\\Delta t$ of 0.01 s. Plot the distribution of the temperature along $x$ for a time of approximately 0.1s, 0.5s and 1s. You should also plot the evolution of the temperature over time in the start, the middle and the end of the bar. Use subplots to visualize theses plots side-by-side.\n",
    "6. What can you comment on your results when comparing the plots of question 4 and 5?\n",
    "7. Define a new function where you solve the PDE problem using the Dufort-Frankel scheme but where the first timestep is computed using the BTCS scheme.\n",
    "8. Re-run the same problem as question 5 using the new function from question 7.\n",
    "9. What can you comment on your results when comparing the plots of question 5 and 8?\n",
    "10. Re-run question 5 but this time the parameter $q_0$ of the flux term should be equal to 0. What can you say about your results?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}