diff --git a/Final projects/2D Variable Diffusion Coefficent.ipynb b/Final projects/2D Variable Diffusion Coefficent.ipynb new file mode 100644 index 0000000..928c492 --- /dev/null +++ b/Final projects/2D Variable Diffusion Coefficent.ipynb @@ -0,0 +1,422 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:3531a2310566805b91dc170ac7a1eef960dda1a13007aca64ca09115d7b2b0a8" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 Louis Camacho.*" + ] + }, + { + "cell_type": "heading", + "level": 1, + "metadata": {}, + "source": [ + "2D Diffusion with Variable Diffusion Coefficent" + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Background" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "I hope everyone is ready for a new lesson! First a little background on what we are going to do. Since we already know what diffusion is, we know that after a certain amount of time and space the concentration of a substance will diminish into the surrounding area until there is equilibrium. Assumming we know what our initial conditions are, we can calculate the time elapsed.The diffusion coefficient is unique in every situation. It is measured at a particular temperature in a particualr medium for a particular material. Up until this point we have held the diffusion coefficent constant making our lives easy. Now its time to instill fear into our hearts and make everything harder, but not to worry when this over we will have a realisitic view on the world around us and have conquered our fears! \n" + ] + }, + { + "cell_type": "heading", + "level": 5, + "metadata": {}, + "source": [ + "Alright lets jump right in. First a quick review of the diffusion in 2D:\n", + "\n", + "\n", + "\\begin{aligned}\n", + "\\frac{\\partial{\\mathbf{u}}}{\\partial t} & = {\\mathbf{v}} (\\frac{\\partial^2{\\mathbf{u}}}{\\partial{\\mathbf{x^2}}}+\\,\\frac{\\partial^2{\\mathbf{u}}}{\\partial{\\mathbf{y^2}}})\n", + "\\end{aligned}" + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "The Problem" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What if our silicon chip was defective, like if it contained several impurities like dust and other contaminants, this will definitely affect our diffusivity since it is no longer just silicone. What are looking to do is have the **v** vary through our experiment while in progress. To do that lets first setup our 2D equation as we have before." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "from matplotlib import rcParams\n", + "rcParams['font.family']='serif'\n", + "rcParams['font.size']=17" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 109 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "L = 1.0e-2\n", + "R = 1.0e-2\n", + "\n", + "nx = 21\n", + "ny = 21\n", + "nt = 500\n", + "\n", + "dx = L/(nx-1)\n", + "dy = R/(ny-1)\n", + "\n", + "x = numpy.linspace(0,L,nx)\n", + "y = numpy.linspace(0,R,ny)\n", + "\n", + "\n", + "Ti = numpy.ones((ny, nx))*10\n", + "Ti[0,:]= 100\n", + "Ti[:,0] = 100" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 110 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "While we can use our code for foward time central space. We will have to play with it a little bit. The first time we have to do is what are we going to do with \"V\"? Well, before that we have ask ourselves what is going on inside the code. We could just slap on a numpy array with certain values and have the code run those values but look what happens:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "vi=np.linspace(.0001,.0003,3) " + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 111 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def fc(T, nt, vi, dt, dx, dy):\n", + "\n", + " j = (np.shape(T)[1])/2\n", + " i = (np.shape(T)[0])/2\n", + " \n", + " for n in range(nt):\n", + " \n", + " \n", + " Tn = T.copy()\n", + " T[1:-1,1:-1] = Tn[1:-1,1:-1] + vi *\\\n", + " (dt/dx**2 * (Tn[1:-1,2:] - 2*Tn[1:-1,1:-1] + Tn[1:-1,:-2]) +\\\n", + " dt/dy**2 * (Tn[2:,1:-1] - 2*Tn[1:-1,1:-1] + Tn[:-2,1:-1]))\n", + " \n", + " # Don't forget our Neumann BCs\n", + " T[-1,:] = T[-2,:]\n", + " T[:,-1] = T[:,-2]\n", + " \n", + " \n", + " if T[j, i] >= 85:\n", + " print (\"Center of plate reached 85C at time {0:.2f}s.\".format(dt*n))\n", + " break\n", + " \n", + " if T[j, i]<85:\n", + " print (\"Center has not reached 85C yet, it is only {0:.2f}C.\".format(T[j, i]))\n", + " \n", + " " + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 112 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "sigma = 0.25\n", + "dt = sigma * min(dx, dy)**2 / vi\n", + "T=Ti.copy()\n", + "T = fc(T, nt, vi, dt, dx, dy)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "operands could not be broadcast together with shapes (3,) (19,19) ", + "output_type": "pyerr", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 2\u001b[0m \u001b[0mdt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msigma\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mmin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdy\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;36m2\u001b[0m \u001b[1;33m/\u001b[0m \u001b[0mvi\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 3\u001b[0m \u001b[0mT\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTi\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m \u001b[0mT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnt\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvi\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdt\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdy\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[1;32m\u001b[0m in \u001b[0;36mfc\u001b[1;34m(T, nt, vi, dt, dx, dy)\u001b[0m\n\u001b[0;32m 8\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 9\u001b[0m \u001b[0mTn\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mT\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 10\u001b[1;33m \u001b[0mT\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mTn\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mvi\u001b[0m \u001b[1;33m*\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mdt\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mdx\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;36m2\u001b[0m \u001b[1;33m*\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mTn\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m-\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mTn\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mTn\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mdt\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mdy\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;36m2\u001b[0m \u001b[1;33m*\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mTn\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m-\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mTn\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mTn\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 11\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 12\u001b[0m \u001b[1;31m# Don't forget our Neumann BCs\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;31mValueError\u001b[0m: operands could not be broadcast together with shapes (3,) (19,19) " + ] + } + ], + "prompt_number": 113 + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "The Solution" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**NOOOO!** **What happened!?** \n", + "\n", + "Well, what happened is that Python is using all values for **v** or *vi* (as we used it) and putting them in to the function and not refreshing **T** like we want. In other words, the first **v** is being used and taking up all the spaces, so that the rest of the **v** values are left out with no spaces from them.\n", + "\n", + "So, how do we fix you say? First, we can use a for loop for **v** values before the **n** loop so that all **v**'s are used. With that done we have to have a copy of that array both outside and inside the **v** loop. That way when the first value of **v** are used **T** can be cleared and ready for the next **v** value until the there are none left. Let's take a look:\n", + " \n", + " " + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def fc(T, nt, vi, dt, dx, dy):\n", + "\n", + " j = (np.shape(T)[1])/2\n", + " i = (np.shape(T)[0])/2\n", + " \n", + " \n", + " Ti=T.copy() #Copy blank T outside\n", + " \n", + " for vi in (v): #Our v loop\n", + " T=Ti.copy() #Copy blank T inside\n", + " dt = sigma * min(dx, dy)**2 / vi #Don't forget to move dt into the v loop, so that dt can also be updated with the new T!\n", + " \n", + " for n in range(nt):\n", + " \n", + " \n", + " Tn = T.copy()\n", + " T[1:-1,1:-1] = Tn[1:-1,1:-1] + vi *\\\n", + " (dt/dx**2 * (Tn[1:-1,2:] - 2*Tn[1:-1,1:-1] + Tn[1:-1,:-2]) +\\\n", + " dt/dy**2 * (Tn[2:,1:-1] - 2*Tn[1:-1,1:-1] + Tn[:-2,1:-1]))\n", + " \n", + " # Don't forget our Neumann BCs\n", + " T[-1,:] = T[-2,:]\n", + " T[:,-1] = T[:,-2]\n", + " \n", + " \n", + " if T[j, i] >= 85:\n", + " print (\"Center of plate reached 85C at time {0:.2f}s.\".format(dt*n))\n", + " break\n", + " \n", + " if T[j, i]<85:\n", + " print (\"Center has not reached 85C yet, it is only {0:.2f}C.\".format(T[j, i]))\n", + " \n", + " " + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 114 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "sigma = 0.25\n", + "T=Ti.copy()\n", + "T = fc(T, nt, vi, dt, dx, dy)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Center of plate reached 85C at time 0.31s.\n", + "Center of plate reached 85C at time 0.16s." + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "Center of plate reached 85C at time 0.10s." + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n" + ] + } + ], + "prompt_number": 115 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Yes!** We get our expected answers! The physics is correct the higher our diffusion coefficent the less time it takes to spread out. The fact that our code is telling us that our varying **v** has a profound effect on the solution is incredible. Basically, it is modeling something like, for example *turbulence*, albeit very miniscule, can effect any experiment. We now have the basis of any real world scenario ready to go!" + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Some Pics!" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from IPython.display import Image\n", + "Image(url='http://mathbench.umd.edu/modules/cell-processes_diffusion/Finalgraphics/graph-flux-vs-gradient-rollover.gif')" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "html": [ + "" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 9, + "text": [ + "" + ] + } + ], + "prompt_number": 9 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As this sophisticated graph shows that a higher D, **v** for us, the higher the slope and that less time is required for complete diffusion. Remember that our flux is dependant on our coefficent and steepness of the gradient as shown below: \n", + "\\begin{aligned}\n", + "{\\mathbf{flux}} = {\\mathbf{-D}} (\\frac{{\\mathbf{Du}}}{{\\mathbf{Dx}}})\n", + "\\end{aligned}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So the amount of flux (rate of heat) moving through the chip is dependant on the coefficent which is not constant and the gradient. Causing the time to shift slow to fast and vice versa. Another thing to note is this formula: " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\\begin{aligned}\n", + "{\\mathbf{Time}} = (\\frac{{\\mathbf{x^2}}}{{\\mathbf{2D}}})\n", + "\\end{aligned}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we increase our distance evenly in small intervals our time will increase dramatically. Just by looking at it we see a quadratic equation that means our solution when graphed we be a parabola. That will show us that as distance increases time will explode upward. Now add to the fact that now we have a varying coefficent our time will be will increase or decrease further depending on the value.\n", + "As we see there are a lot of factors, but it is not impossible to deal with them. We have our tools let's use them!" + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "References:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "(http://mathbench.umd.edu/modules/cell-processes_diffusion/page10.htm) University of Maryland, Mathbench\n", + "\n", + "\n", + "(Numerical MOOC \"Spreading Out\" Module 4 Lesson 4) George Washington University, Professor Lorena Barba & Gilbert Forsyth\n", + "\n", + "(http://en.wikipedia.org/wiki/Fick%27s_laws_of_diffusion) Wiki Article, Wikipedia" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from IPython.core.display import HTML\n", + "css_file = '../../styles/numericalmoocstyle.css'\n", + "HTML(open(css_file, \"r\").read())" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "ename": "IOError", + "evalue": "[Errno 2] No such file or directory: '../../styles/numericalmoocstyle.css'", + "output_type": "pyerr", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mIOError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mIPython\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcore\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdisplay\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mHTML\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 2\u001b[0m \u001b[0mcss_file\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'../../styles/numericalmoocstyle.css'\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0mHTML\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mopen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcss_file\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"r\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mread\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[1;31mIOError\u001b[0m: [Errno 2] No such file or directory: '../../styles/numericalmoocstyle.css'" + ] + } + ], + "prompt_number": 10 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file