From 0dec1ad26d829b30237820895f8ee09ef4e910b5 Mon Sep 17 00:00:00 2001 From: carlosej Date: Sun, 14 Dec 2014 21:51:01 -0500 Subject: [PATCH 1/3] Create 2D Variable Diffusion Coefficent --- .../2D Variable Diffusion Coefficent | 155 ++++++++++++++++++ 1 file changed, 155 insertions(+) create mode 100644 Final projects/2D Variable Diffusion Coefficent diff --git a/Final projects/2D Variable Diffusion Coefficent b/Final projects/2D Variable Diffusion Coefficent new file mode 100644 index 0000000..53abf9e --- /dev/null +++ b/Final projects/2D Variable Diffusion Coefficent @@ -0,0 +1,155 @@ + +Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 Louis Camacho. +2D Diffusion with Variable Diffusion Coefficent +Background +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! +Alright lets jump right in. First a quick review of the diffusion in 2D: \begin{aligned} \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}}}) \end{aligned} +The Problem +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. +In [109]: +import numpy as np +import matplotlib.pyplot as plt +%matplotlib inline +from matplotlib import rcParams +rcParams['font.family']='serif' +rcParams['font.size']=17 +In [110]: +L = 1.0e-2 +R = 1.0e-2 + +nx = 21 +ny = 21 +nt = 500 + +dx = L/(nx-1) +dy = R/(ny-1) + +x = numpy.linspace(0,L,nx) +y = numpy.linspace(0,R,ny) + + +Ti = numpy.ones((ny, nx))*10 +Ti[0,:]= 100 +Ti[:,0] = 100 +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: +In [111]: +vi=np.linspace(.0001,.0003,3) +In [112]: +def fc(T, nt, vi, dt, dx, dy): + + j = (np.shape(T)[1])/2 + i = (np.shape(T)[0])/2 + + for n in range(nt): + + + Tn = T.copy() + T[1:-1,1:-1] = Tn[1:-1,1:-1] + vi *\ + (dt/dx**2 * (Tn[1:-1,2:] - 2*Tn[1:-1,1:-1] + Tn[1:-1,:-2]) +\ + dt/dy**2 * (Tn[2:,1:-1] - 2*Tn[1:-1,1:-1] + Tn[:-2,1:-1])) + + # Don't forget our Neumann BCs + T[-1,:] = T[-2,:] + T[:,-1] = T[:,-2] + + + if T[j, i] >= 85: + print ("Center of plate reached 85C at time {0:.2f}s.".format(dt*n)) + break + + if T[j, i]<85: + print ("Center has not reached 85C yet, it is only {0:.2f}C.".format(T[j, i])) + + +In [113]: +sigma = 0.25 +dt = sigma * min(dx, dy)**2 / vi +T=Ti.copy() +T = fc(T, nt, vi, dt, dx, dy) +--------------------------------------------------------------------------- +ValueError Traceback (most recent call last) + in () + 2 dt = sigma * min(dx, dy)**2 / vi + 3 T=Ti.copy() +----> 4 T = fc(T, nt, vi, dt, dx, dy) + + in fc(T, nt, vi, dt, dx, dy) + 8 + 9 Tn = T.copy() +---> 10 T[1:-1,1:-1] = Tn[1:-1,1:-1] + vi * (dt/dx**2 * (Tn[1:-1,2:] - 2*Tn[1:-1,1:-1] + Tn[1:-1,:-2]) + dt/dy**2 * (Tn[2:,1:-1] - 2*Tn[1:-1,1:-1] + Tn[:-2,1:-1])) + 11 + 12 # Don't forget our Neumann BCs + +ValueError: operands could not be broadcast together with shapes (3,) (19,19) +The Solution +NOOOO! What happened!? +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. +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: +In [114]: +def fc(T, nt, vi, dt, dx, dy): + + j = (np.shape(T)[1])/2 + i = (np.shape(T)[0])/2 + + + Ti=T.copy() #Copy blank T outside + + for vi in (v): #Our v loop + T=Ti.copy() #Copy blank T inside + 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! + + for n in range(nt): + + + Tn = T.copy() + T[1:-1,1:-1] = Tn[1:-1,1:-1] + vi *\ + (dt/dx**2 * (Tn[1:-1,2:] - 2*Tn[1:-1,1:-1] + Tn[1:-1,:-2]) +\ + dt/dy**2 * (Tn[2:,1:-1] - 2*Tn[1:-1,1:-1] + Tn[:-2,1:-1])) + + # Don't forget our Neumann BCs + T[-1,:] = T[-2,:] + T[:,-1] = T[:,-2] + + + if T[j, i] >= 85: + print ("Center of plate reached 85C at time {0:.2f}s.".format(dt*n)) + break + + if T[j, i]<85: + print ("Center has not reached 85C yet, it is only {0:.2f}C.".format(T[j, i])) + + +In [115]: +sigma = 0.25 +T=Ti.copy() +T = fc(T, nt, vi, dt, dx, dy) +Center of plate reached 85C at time 0.31s. +Center of plate reached 85C at time 0.16s. +Center of plate reached 85C at time 0.10s. + +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! +Some Pics! +In [9]: +from IPython.display import Image +Image(url='http://mathbench.umd.edu/modules/cell-processes_diffusion/Finalgraphics/graph-flux-vs-gradient-rollover.gif') +Out[9]: + +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: \begin{aligned} {\mathbf{flux}} = {\mathbf{-D}} (\frac{{\mathbf{Du}}}{{\mathbf{Dx}}}) \end{aligned} +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. If we wanted to we could have this go on until the chip fails, but we'll save that for another day. +Citations & Recongition (Shout Outs!): +(http://mathbench.umd.edu/modules/cell-processes_diffusion/page10.htm) University of Maryland, Mathbench +(Numerical MOOC Module 4 Lesson 4) George Washington University, Professor Lorena Barba & Gilbert Forsyth +(http://en.wikipedia.org/wiki/Fick%27s_laws_of_diffusion) Wiki Article, Wikipedia +In [10]: +from IPython.core.display import HTML +css_file = '../../styles/numericalmoocstyle.css' +HTML(open(css_file, "r").read()) +--------------------------------------------------------------------------- +IOError Traceback (most recent call last) + in () + 1 from IPython.core.display import HTML + 2 css_file = '../../styles/numericalmoocstyle.css' +----> 3 HTML(open(css_file, "r").read()) + +IOError: [Errno 2] No such file or directory: '../../styles/numericalmoocstyle.css' +In []: From 0f34f780c149b732af230a35c35f4e7b09dd53bc Mon Sep 17 00:00:00 2001 From: carlosej Date: Mon, 15 Dec 2014 13:34:48 -0500 Subject: [PATCH 2/3] Revert "Create 2D Variable Diffusion Coefficent" This reverts commit 0dec1ad26d829b30237820895f8ee09ef4e910b5. --- .../2D Variable Diffusion Coefficent | 155 ------------------ 1 file changed, 155 deletions(-) delete mode 100644 Final projects/2D Variable Diffusion Coefficent diff --git a/Final projects/2D Variable Diffusion Coefficent b/Final projects/2D Variable Diffusion Coefficent deleted file mode 100644 index 53abf9e..0000000 --- a/Final projects/2D Variable Diffusion Coefficent +++ /dev/null @@ -1,155 +0,0 @@ - -Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 Louis Camacho. -2D Diffusion with Variable Diffusion Coefficent -Background -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! -Alright lets jump right in. First a quick review of the diffusion in 2D: \begin{aligned} \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}}}) \end{aligned} -The Problem -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. -In [109]: -import numpy as np -import matplotlib.pyplot as plt -%matplotlib inline -from matplotlib import rcParams -rcParams['font.family']='serif' -rcParams['font.size']=17 -In [110]: -L = 1.0e-2 -R = 1.0e-2 - -nx = 21 -ny = 21 -nt = 500 - -dx = L/(nx-1) -dy = R/(ny-1) - -x = numpy.linspace(0,L,nx) -y = numpy.linspace(0,R,ny) - - -Ti = numpy.ones((ny, nx))*10 -Ti[0,:]= 100 -Ti[:,0] = 100 -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: -In [111]: -vi=np.linspace(.0001,.0003,3) -In [112]: -def fc(T, nt, vi, dt, dx, dy): - - j = (np.shape(T)[1])/2 - i = (np.shape(T)[0])/2 - - for n in range(nt): - - - Tn = T.copy() - T[1:-1,1:-1] = Tn[1:-1,1:-1] + vi *\ - (dt/dx**2 * (Tn[1:-1,2:] - 2*Tn[1:-1,1:-1] + Tn[1:-1,:-2]) +\ - dt/dy**2 * (Tn[2:,1:-1] - 2*Tn[1:-1,1:-1] + Tn[:-2,1:-1])) - - # Don't forget our Neumann BCs - T[-1,:] = T[-2,:] - T[:,-1] = T[:,-2] - - - if T[j, i] >= 85: - print ("Center of plate reached 85C at time {0:.2f}s.".format(dt*n)) - break - - if T[j, i]<85: - print ("Center has not reached 85C yet, it is only {0:.2f}C.".format(T[j, i])) - - -In [113]: -sigma = 0.25 -dt = sigma * min(dx, dy)**2 / vi -T=Ti.copy() -T = fc(T, nt, vi, dt, dx, dy) ---------------------------------------------------------------------------- -ValueError Traceback (most recent call last) - in () - 2 dt = sigma * min(dx, dy)**2 / vi - 3 T=Ti.copy() -----> 4 T = fc(T, nt, vi, dt, dx, dy) - - in fc(T, nt, vi, dt, dx, dy) - 8 - 9 Tn = T.copy() ----> 10 T[1:-1,1:-1] = Tn[1:-1,1:-1] + vi * (dt/dx**2 * (Tn[1:-1,2:] - 2*Tn[1:-1,1:-1] + Tn[1:-1,:-2]) + dt/dy**2 * (Tn[2:,1:-1] - 2*Tn[1:-1,1:-1] + Tn[:-2,1:-1])) - 11 - 12 # Don't forget our Neumann BCs - -ValueError: operands could not be broadcast together with shapes (3,) (19,19) -The Solution -NOOOO! What happened!? -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. -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: -In [114]: -def fc(T, nt, vi, dt, dx, dy): - - j = (np.shape(T)[1])/2 - i = (np.shape(T)[0])/2 - - - Ti=T.copy() #Copy blank T outside - - for vi in (v): #Our v loop - T=Ti.copy() #Copy blank T inside - 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! - - for n in range(nt): - - - Tn = T.copy() - T[1:-1,1:-1] = Tn[1:-1,1:-1] + vi *\ - (dt/dx**2 * (Tn[1:-1,2:] - 2*Tn[1:-1,1:-1] + Tn[1:-1,:-2]) +\ - dt/dy**2 * (Tn[2:,1:-1] - 2*Tn[1:-1,1:-1] + Tn[:-2,1:-1])) - - # Don't forget our Neumann BCs - T[-1,:] = T[-2,:] - T[:,-1] = T[:,-2] - - - if T[j, i] >= 85: - print ("Center of plate reached 85C at time {0:.2f}s.".format(dt*n)) - break - - if T[j, i]<85: - print ("Center has not reached 85C yet, it is only {0:.2f}C.".format(T[j, i])) - - -In [115]: -sigma = 0.25 -T=Ti.copy() -T = fc(T, nt, vi, dt, dx, dy) -Center of plate reached 85C at time 0.31s. -Center of plate reached 85C at time 0.16s. -Center of plate reached 85C at time 0.10s. - -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! -Some Pics! -In [9]: -from IPython.display import Image -Image(url='http://mathbench.umd.edu/modules/cell-processes_diffusion/Finalgraphics/graph-flux-vs-gradient-rollover.gif') -Out[9]: - -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: \begin{aligned} {\mathbf{flux}} = {\mathbf{-D}} (\frac{{\mathbf{Du}}}{{\mathbf{Dx}}}) \end{aligned} -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. If we wanted to we could have this go on until the chip fails, but we'll save that for another day. -Citations & Recongition (Shout Outs!): -(http://mathbench.umd.edu/modules/cell-processes_diffusion/page10.htm) University of Maryland, Mathbench -(Numerical MOOC Module 4 Lesson 4) George Washington University, Professor Lorena Barba & Gilbert Forsyth -(http://en.wikipedia.org/wiki/Fick%27s_laws_of_diffusion) Wiki Article, Wikipedia -In [10]: -from IPython.core.display import HTML -css_file = '../../styles/numericalmoocstyle.css' -HTML(open(css_file, "r").read()) ---------------------------------------------------------------------------- -IOError Traceback (most recent call last) - in () - 1 from IPython.core.display import HTML - 2 css_file = '../../styles/numericalmoocstyle.css' -----> 3 HTML(open(css_file, "r").read()) - -IOError: [Errno 2] No such file or directory: '../../styles/numericalmoocstyle.css' -In []: From 1492680e34f39a9033d107fc2355bd9277b3efe6 Mon Sep 17 00:00:00 2001 From: carlosej Date: Mon, 15 Dec 2014 13:36:39 -0500 Subject: [PATCH 3/3] Adding Final Project --- .../2D Variable Diffusion Coefficent.ipynb | 422 ++++++++++++++++++ 1 file changed, 422 insertions(+) create mode 100644 Final projects/2D Variable Diffusion Coefficent.ipynb 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