diff --git a/alexabaumer/.ipynb_checkpoints/Baumer_finalProject-checkpoint.ipynb b/alexabaumer/.ipynb_checkpoints/Baumer_finalProject-checkpoint.ipynb new file mode 100644 index 0000000..a1512f4 --- /dev/null +++ b/alexabaumer/.ipynb_checkpoints/Baumer_finalProject-checkpoint.ipynb @@ -0,0 +1,119 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "We are going to be looking at crime hotspots- how they develop and disspate. \n", + "\n", + "We are going to be looking at pedestrian flow, specifically how people exit crowded spaces. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There is relatively little understood about how crime hotspots develop, grow, and disspate in communities across the country. In order to develop more robust and effective crime prevention strategies, mathematicians have attempted to model crime as a reaction-diffusion model. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# The Short et al Model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$ {dB}/{dt} = \\eta D (del square) B - \\omega B + \\kappa \\rho A $$\n", + "\n", + "$$ {d\\rho}/{dt} = D (grad) * [(grad) \\rho - 2 \\rho (grad) lnA] - \\rho A + \\gamma $$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Variable declarations:\n", + " \n", + "B : the risk of victimization \n", + "D : the diffusion coefficient \n", + "\n", + "$$\\eta$$ : control rate, between 0 and 1\n", + "$$\\omega$$ : rate at which the risk decays towards a fixed, environmental value\n", + "$$\\kappa$$ : attractive force pulling offenders back to previous successful locations, inferred from empirical evidence\n", + "$$\\rho(x)$$ : the density of criminal agents \n", + "A : \n", + "$$\\gamma$$ : \n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Simplifying to specific crime types " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Short's model operates independently of crime type. However, for our purposes, we will simplify this model to focus on home burglary. Our crime operatives remain mobile while our targets (homes) are stationary." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "References! \n", + "\n", + "[1] Short et al 2010 PNAS " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python [Root]", + "language": "python", + "name": "Python [Root]" + }, + "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.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/alexabaumer/Lahar_simple.ipynb b/alexabaumer/Lahar_simple.ipynb new file mode 100644 index 0000000..0156e0b --- /dev/null +++ b/alexabaumer/Lahar_simple.ipynb @@ -0,0 +1,524 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Lahars: rivers of concrete" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\"Lahar\" is an Indonesian term that describes a hot of cold mixture of water and rock fragments that flows down the slopes of a volcano and typically enters a river valley [1]. One of the major volcanic hazards, they typically occur around stratovolcanoes- such as the Cascade Range in the Western United States. Lahar events can happen with or without a volcanic eruption- smaller events are often termed \"debris flows.\"\n", + "\n", + "A quick search on YouTube reveals many lahar videos, like this one: https://www.youtube.com/watch?v=TpwiFtVRBTs\n", + "\n", + "Lahars can best be equated to concrete, while it moves- it is a wet, slurry mixture, when it rests, it solidifies into a rocky deposit. Mitigating volcanic risks is a high priority goal for geolgists around the world; if we are able to better predict the behavior of future eruptions, we can warn at risk populations living in hazard zones as well as think about \n", + "engineering structures to protect those populations [2]. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For this particular lahar model, we will use what we learned in modules 2 and 3 about convection problems. \n", + "\n", + "When we begin to think about modeling the height of lahars in a numerical space, we first will make some assumptions about the characteristics of our lahar. The first is that our lahar will be of constant width. This assumption limits the applications of this model to lahars that do not experience bulking and debulking. We also are able to reduce the probelm to one dimension with this assumption- simplifying our calculations considerably.The second assumption we will make is that the rheological behavior of our lahar is similar to water. Since we are modeling the lahar while in motion, this is not too bad of an assumption to make. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To begin our modeling process, we will import our favorite libraries" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import numpy \n", + "from matplotlib import pyplot\n", + "from matplotlib import rcParams\n", + "rcParams['font.family'] = 'serif'\n", + "rcParams['font.size'] = 16\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Derivation of lahar governing equation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The majority of this derivation comes from the textbook, Mathematical Modeling of Earth's Dynamical Systems by Lee Kump and Rudy Slingerland [3]. \n", + "\n", + "We will start by balancing the lahar mass with conservation. We set up a control volume where the volumetric flow rate is the volumetric flux, q, multiplied by the width. The mass rate out of the section is obtained from a Taylor series. And assuming width and bulk density do not vary in space or time, we get..." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\\frac{\\partial h}{\\partial t} + \\frac{\\partial q}{\\partial x} = 0$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This equation looks familiar to us, as it takes the form of our \"traffic\" model used in Module 2 and 3. We now need to find a way to relate q to h, x, and t. A force balance on a segment of the lahar should do the trick. Starting with the downslope gravity force... A is our cross-sectional area of the flow and we know that as long as the bed slope is small, sin(alpha) equals the bed slope, S. The force opposing gravity is the shear stress multiplied by the area over which it acts, where P is the wetted perimeter..." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$ F_g = \\sigma A dx g sin \\alpha$$\n", + "$$ F_s = \\tau_0 P dx $$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we equate the two in our force balance, we reach an expression for shear stress- where A/P is R- the hydraulic radius. As long as the flow width is 20 times the flow depth, the hydraulic radis is not different from the flow depth, h ; We aso know that a turbulent flow of water exerts a shear stress on its boundaries. This is proportional to the square of the depth averaged velocity multiplied by the drag coefficient... " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$ \\tau_0 = \\sigma g R S = C_f \\sigma v^2 $$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "When we equate these, we can find an expression for velocity of the flow; we will call the first term \"Beta\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\\nu = \\sqrt{\\frac{g S}{C_f}}\\sqrt{h} $$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The expression of velocity is a generalization of the Chezy equation, predicting the cross-sectional average velocity of steady, uniform flows. Substituting our equation for velocity into our definition of q, we find the following expression: " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$ q = \\beta h^k $$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And now we have enough information to return to our governing equations and begin discretizing!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discretization of lahar governing equations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For this problem, we will be using the MacCormack scheme. We first learned about this scheme in Module 3, Lesson 2. As a reminder, it is a two step method utilizing a predictor step and a corrector step. It achieves second-order accuracy in both space and time [4]. Our discretized equation looks as follows:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$ h_i^* = h_i^n - \\frac{\\Delta t}{\\Delta x}(q_{i+1}^n - q_i^n) $$ \n", + "$$ h_i^{n+1} = \\frac{1}{2} (h_i^n + h_i^* - \\frac{\\Delta t}{\\Delta x}(q_i^* - q_{i-1}^*)$$" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def computeQ(beta, h, k):\n", + " \"\"\"Parameters:\n", + " beta: coeeficient relating gravity, slope, and drag coefficient\n", + " h: height of lahar\n", + " k: coefficient, given in text\n", + " \"\"\" \n", + " return beta*h**k" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def lahar_maccormack(h, nt, dt, dx, beta, k):\n", + " \n", + " \"\"\"Computes the height of the lahar with the MacCormack scheme \n", + " \n", + " Parameters:\n", + " h: height of the lahar\n", + " nt: number of time steps\n", + " dt: time-step size\n", + " dx: mesh spacing\n", + " alpha: \"\"\"\n", + " \n", + " #initialize results \n", + " h_n = numpy.zeros((nt,len(h)))\n", + " h_star = numpy.empty_like(h)\n", + " \n", + " #copy initial array into new array\n", + " h_n[:,:] = h.copy()\n", + " h_star = h.copy()\n", + " \n", + " #compute h_n at time step t+1 and then set h = h_n\n", + " \n", + " for t in range(1,nt):\n", + " q = computeQ(beta, h, k)\n", + " h_star[:-1] = h[:-1] - dt/dx * (q[1:] - q[:-1])\n", + " q_star = computeQ(beta, h_star, k)\n", + " h_n[t,1:] = 0.5 * (h[1:] + h_star[1:] - dt/dx * (q_star[1:] - q_star[:-1]))\n", + " h = h_n[t].copy()\n", + " \n", + " vel = beta * (h**(0.5))\n", + "\n", + " return h_n, vel" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Defining our lahar field " + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#Lahar and river valley characteristics \n", + "\n", + "g = 9.8 #gravity, m/s^2\n", + "S = 0.1 #bed slope, can find this by looking at topographical maps\n", + "C_f = 10.0 #drag coefficient of river beds\n", + "h_initial = 2.0 #initial height of bed, m\n", + "h_max = 30.0 #initial height of source of lahar, m\n", + "k = 1.5 #Vignaux and Weir 1990 observed k values for lahars to range from 1.24 to 1.47\n", + "\n", + "beta = ((g * S)/C_f)**0.5" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def lahar_initial(nx, h_max, h_initial):\n", + " \"\"\"Returns the initial height of the lahar by describing its lateral extent\n", + " \n", + " nx: number of grid points in x\n", + " h_max: maximum height of the lahar at its source\n", + " h_initial: the height of the bed slope\n", + " \n", + " h: height of the lahar as it propogates\n", + " \"\"\"\n", + " h = h_initial*numpy.ones(nx)\n", + " h[:int((nx-1)*1./10.)] = h_max\n", + " \n", + " return h" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 36, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAENCAYAAAD+CUlOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE4hJREFUeJzt3X+0ZWV93/H3h/JLYAZQGYEqSQYjGCM2JlpYCF5r1YbB\nVi0kWV1RNKuVstoGjFOVIdUbfxCgDE1CjQvtyrI/omnX6MLwK4ZIbrFkoiBGk+pSCShFDZBk2pkk\naMD59o+zB+4e98xzh7n7njMz79das+7ez/mxv+fhcD937+fZe6eqkCRpdw6adgGSpNlnWEiSmgwL\nSVKTYSFJajIsJElNhoUkqengMd88yWbgESDAY1X1iiTHAlcA9wLPBjZU1cNj1iFJ2jujhgVwS1W9\ne6e2y4Fbq2pTknOBjcAbRq5DkrQXMuZJeUk2AZ8BjgDurKqbk9wPnFFV3+z2Mu6pqqeNVoQkaa+N\nvWdxRVXdleQg4PYk24DjgG3d41uBY5IcVFXbR65FkvQkjTrAXVV3dT+3A58GXgY8BKzqnrIa2GJQ\nSNJsG23PIskpwJlV9Rtd0w8DHwduAs4ANgFndutDr/eiVZL0JFRVlvs9xzwMtRU4J8kJwNHA/VX1\nkSS3AFd0YbIWWL+rN/Aih8tnfn6e+fn5aZexX7Avl5f9ubySZc8JYMSwqKpvA+cNtG8BLhxru5Kk\n5edJeZKkJsPiADE3NzftEvYb9uXysj/3DaOeZ7E3ktSs1iZJsyrJKAPc7llIkpoMC0lSk2EhSWoy\nLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNC\nktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJ\nTYaFJKnJsJAkNY0eFkkOT/KFJFd168cmuS7J25N8KMlxY9cgSdo7K7Fn8V7g7kXrlwO3VtWVwCeA\njStQgyRpL4waFkl+FvhfwNcXNa8DNnfLd3TrkqQZNlpYJHkucGpVXb/TQ2uAbd3yVuCYJI6dSNIM\nO3jE934t8J0kbwdeAhyS5GLgQWAVk6BYDWypqu1DbzA/P//48tzcHHNzcyOWK0n7noWFBRYWFkbf\nTqpq/I0k7wKOrKq3Jfl14Laq2pTkXOD8qrpg4DW1ErVJ0v4kCVWV5X7fMfcsAEjyOuAs4NAkPw1s\nAK5McgqwFlg/dg2SpL2zInsWT4Z7FpK058bas3BgWZLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJ\nsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwL\nSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAk\nNRkWkqSmg8d64yQBbgD+EDgMOBl4E3AEcAVwL/BsYENVPTxWHZKkvTdaWHT+oKouB0hyPfBPgbOA\nW6tqU5JzgY3AG0auQ5K0F1JV428kORjYDFwIXA+cUVXfTHIscE9VPW3gNbUStUnS/iQJVZXlft+x\n9yxI8krgLcCNVXV3kjXAtu7hrcAxSQ6qqu07v/blF10zdnm79by1JzL/5lfz1KOPnGodkjRto4dF\nVf0u8LtJ/nOSi4AHgVVMgmI1sGUoKABuu/MrY5e3W7fd+RVOePrRXPqmn5xqHZI0bWMOcD8X+KGq\nurlrug/4IeAm4AxgE3Bmtz7sW3c9sbzqxMm/FfaNb//lim9TkpZqYWGBhYWF0bcz2phFkrXAVcDd\nwKHAqcDPA48ymQ11P7AWeMfQbKgk9Xuf+dIotbXc+Ok/5lc++ikA3vzas7jusp+dSh2StKf2uTGL\nqroXOG8XD1+4lPd4+Yufu3wF7YE/feCJ7CocZJckT8obMDlFZMIJWZJkWAxaHBaSJMOiyXM9JMmw\nGLR4v8KskCTDYlB/zMK0kCTDYsDiIQtnQ0mSYSFJWgLDYoBTZyWpz7AY4JiFJPUZFgOcDSVJfYaF\nJKnJsBjgYShJ6ltSWHQ3MNqxfFqSS8crafqcOitJfUvdszh9x0JVfRE4apxyZoOzoSSpb7eXKE9y\nMXAJcHSSN3bN24HPjVyXJGmG7DYsqupXgV9N8i+q6kMrVNPUOWYhSX1LOgy1c1Akedk45cwGp85K\nUt+S7pSX5CeBi5iMVQQ4CTh5xLqmyj0LSepb6m1VLwMuBv6cSVhcMFpFkqSZs9SwuLuqHh/UTvJf\nR6pnJjh1VpL6WrOh3tktHt8FxD3d+lnAPxyzsGly6qwk9bX2LE4DbgC+sVP788cpZzY4ZiFJfa2w\nuKSqHkhyVFX91Y7GJLeMXJckaYa0zrN4oFv8tSQfXvxQkkMXPb5fceqsJPUtdYD7OcB7gfuAtcB3\nAJL8dlVdO1JtU+NhKEnqW+q1oW6pqrOr6oKqOgu4tapeATxjxNqmxtlQktS31LA4caf1HSfkbV3G\nWiRJM2qpYfFYkpuS/EqSm4G/TfITwEtHrG1qnDorSX1LGrOoqouTnAM8j8khqJu6h9aNVtkUOWYh\nSX1LHeCmqm4GbgZI8tNV9d9Hq2rKnA0lSX2tM7g/UlX/LMl98PhIb4DVwH4bFpKkvtaexXz385rF\nU2STXDRaRTPAw1CS1LfbAe6q+mr389okByV5epJU1QdWprzpcOqsJPUt9X4WrwQ+CHwR+K0kq6rq\nusZr1jI5ke9zwLOAv6iq9yQ5FrgCuBd4NrChqh7ei8+w7JwNJUl9S506+2rgVOCOqvoI33/exZCn\nAh+tqo1VdQnwM0l+DLicyYyqK4FPABufRN2SpBW01LB4oKq+wxOD3N9tvaCq7qqqGxY1BfhrJtNt\nN3dtdzCD028ds5CkvqWGxXOSvAP4kST/GnjmnmwkyWuAT3ZjIGuAbd1DW4Fjkiy1jhXh1FlJ6mtN\nnX1md2XZS4BLgacDxwPvW+oGkswBc92hKIAHgVVMgmI1sKWqtg+9dn5+/vHlubk55ubmlrrZvdLb\ns3CAW9IMW1hYYGFhYfTtZHeHWZL8N+AjOzcDP1VVzftwJ1kHvKSqLk1yAvADwBuA26pqU5JzgfOH\n3msy6Wo6v6hvueNPOOfiyUzhV53xI/zOtRdPpQ5J2lNJqKq0n7lnWrOhzgT+brd8MpPbqgY4qfXG\nSV4I/BZwZ5LfB44A3g9sAK5McgqTy52vf3KljyfL3s2StG9rhcX6qvoYQJL/UFVv6ZZf03rjqrqb\nyeGmIRfuUZUrzKmzktTXOinvY4tXF7VfP1pFM8DZUJLUt9uwSHL6LtpfPE45s8eskKT2Yairk/xB\nt3xWkqu65dOBs8cra7ocspCkvlZYPMrkRDqAG3dq3285dVaS+lph8baqunPnxiQ/PlI9M6F3IUGP\nQ0lSc4D7+4Kia//cOOXMHrNCkpZ+uY8DSjzRQpJ6DIsBTp2VpD7DYoAXEpSkPsOiwdlQkmRYDHLM\nQpL6DIsBTp2VpD7DYoAXEpSkPsOiwT0LSTIsBjlmIUl9hsUAp85KUp9hMcALCUpSn2HR4JiFJBkW\ngxyykKQ+w2KAU2clqc+waPAwlCQZFoPcs5CkPsNigEMWktRnWAxw6qwk9RkWDY5ZSJJhMah/1dnp\n1SFJs8KwGOBtVSWpz7AY4IUEJanPsGhwx0KSDItBvavOOhtKkgyLIZ6UJ0l9hsUAhywkqc+waHA2\nlCSNGBZJnpHkQ0k+u6jt2CTXJXl799hxY21/bzh1VpL6xtyzOBP4xE5tlwO3VtWV3WMbR9z+k+aY\nhST1jRYWVfVxYNtOzeuAzd3yHd36zHHIQpL6VnrMYg1PBMhW4JgkMz1u4tRZSYKDV3h7DwKrmATF\namBLVW3f1ZPn5+cfX56bm2Nubm7k8iY8DCVpX7GwsMDCwsLo28mYA7hJXgpcXVUv6tZ/HbitqjYl\nORc4v6ou2MVra1qDy1+691s876d+CYBTf/B4vrzpl6ZShyTtqSRU1bIfTR9zNtTZwOuB45NsSHIY\ncBnwiiSXAa8F1o+1/b3htaEkqW+0w1BVdTtw+07N3wUuHGubY3DqrCR5Ut4gxywkqc+wGOCFBCWp\nz7AY4JiFJPUZFg0ehpIkw2JQ/x7cpoUkGRYDvJCgJPUZFgMcs5CkPsOiwf0KSTIsBvWmznoYSpIM\niyGelCdJfYbFAIcsJKnPsGjwMJQkGRaDnDorSX2GxYBeWEyxDkmaFYbFAIcsJKnPsGjwMJQkGRaD\nnDorSX2GxQAvJChJfYbFAGdDSVKfYSFJajIsBjh1VpL6DIsBXkhQkvoMiwHOhpKkPsNCktRkWAxw\n6qwk9RkWA5w6K0l9B0+7gFm0OCy+t7345kNbpliNJE2fYdHwf7f9Dc885x3TLkOSpsrDUAOOPPxQ\nDj/skGmXIUkzwz2LAU85/FA2XnIe1/zm7/HIdx+ddjmStGTfGul9M6sDuElqVmuTpFmVhKpa9tvy\neBhKktQ0lcNQSV4OvA54EKCq3j2NOiRJS7Pih6GSPAX4IvDcqnosySbg/VX1+zs9z8NQkrSH9qfD\nUGcAX6+qx7r1O4B1U6hDkrRE0wiLNcC2RetbuzZJ0oyaRlg8BKxetL66a5MkzahpDHBvBk5KckhV\nPQqcCbx/6Inz8/OPL8/NzTE3N7cS9UnSPmNhYYGFhYXRtzOV8yy62VDnM9mjeLSq3jPwHAe4JWkP\njTXA7Ul5krQf2Z9mQ0mS9jGGhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUk\nqcmwkCQ1GRaSpCbDQpLUZFgcIFbievcHCvtyedmf+wbD4gDh/5DLx75cXvbnvsGwkCQ1GRaSpKaZ\nvlPetGuQpH3RAXVbVUnS7PAwlCSpybCQJDUdPO0Cdpbk5cDrgAcBqurd061odiXZDDwCBHisql6R\n5FjgCuBe4NnAhqp6uHv+emA1cAxwa1Xd0LW/APhXwH3AGuCtVbV9pT/PSkvyDOC9wAuq6sVd27L1\nX5LDgKuBb3bvdWVVfW0FP+KK2UVfvgt46aKnva+qPtU9Zl/uRpK1TPrzc8CzgL+oqvdM9ftZVTPz\nD3gK8DXg4G59E/Cyadc1q/+Adw60fQA4r1s+F/gv3fKLgRu75b8DfBVY1a3/MXBct3w18KZpf7YV\n6r/XdX302TH6D3g7sL5b/lHg9ml/5hXuy+/7ftqXS+7PnwBevWj9fwM/Ns3v56wdhjoD+HpVPdat\n3wGsm2I9s+60JP82ybuSnNO1rQM2d8t3ADvaz93RXlXfA74MvLT7C+bw6v464QDq86r6OLBtp+bl\n7L91i17zJ0z+ex01wkeZul30ZZJsSPLWJG9L8pSu3b5sqKq7qtsz6AT4a6b4/Zy1w1Br6H/htjJJ\nUw27oqruSnIQcHuSbcBxPNGHW4Fju8fXAF9a9NqtXduf8/19vmb0ymfX4u/g3vbfzt/nbV3bXy1/\n2TPpfzD54++RJBcB1wL/HPtyjyR5DfDJqvpqkql9P2dtz+IhJsfcdljdtWlAVd3V/dwOfBp4GZP+\nWtU9ZTWwpXt8cfuOxx7CPt/Zgyxf/+3qNQeEqvpyVT3Srd7G5PsJ9uWSJZkD5qrqLV3T1L6fsxYW\nm4GTkhzSrZ8J3DTFemZWklOS/Nyiph8G7mHSX2d0bS/hif57vL3r31OZHKO8F/ib7i8WODD7fPEJ\nTMvZfzcues3zgT+qqv39L+HH+zLJVYvanwP8abdsXy5BknXAq6rqkiQnJDmdKX4/Z+6kvG421PlM\nEu7RqnrPlEuaSUlOYLJb/3ngaCaTAn5h0WyJ+4G1wDvqidkSbwWeymS2xC1VdWPXfhrw88DXu8fX\n14ExG+ps4A3Aq5gMHG4EjmCZ+i/J4cC/B/4MOBm4vKruWbEPuIIG+vIa4N8x6c+HmQygvnPH57cv\ndy/JC4H/CdzJJICPAN4P/DZwJVP4fs5cWEiSZs+sHYaSJM0gw0KS1GRYSJKaDAtJUpNhIUlqMiwk\nSU2zdrkPaVkleRFwFXAo8Ekm89WLyaVS/l+SzwJ/v3YxhzzJPwE+X1X3r1TN0iwyLLRfq6o7kywA\nR1Z3ufsk/wi4LcmLqruc9m68BtjC5CQo6YBlWOiAU1W/k+SdwLokv8bkngsPAf+RydU6nwXcDvwf\n4O8Bb0xyelVdleSDwAPAUcCfVdU13WUYPgDcDfwtcBrwb6rq7u4qnhuZXOrieCb3JXhfkrOANzG5\n9PQpTM7E/cuV6gNpTxkWOlDdz+SX933d+inAC4D1TC4F/fzul/0fAR+uqtu7591QT9xU5vNJrquq\nP0xyPXBoVV2W5DzgAibhsQH4WlVd3b3mjd37fBR4UVV9O8kFwC8CvzDyZ5aeNMNCB6ofAL6xY6Wq\nvpDkA8DHgO8y+SU/5MQk72VyOedVwNOYhAtMbjgDk2sh7bia52nAf1q0nQ8neTqTa/S8PkmAY4Ed\n93CRZpJhoQNOklcChwGfAi7t2n4Q+ExV/UZ3I6l5JuMV3wMO6m4icxTwtqo6uXvNP97prYcGyb/A\n5CJtO7b95qr6YJKHgOu6QfanAqcv3yeUlp9hof1akh8HzgYOSbIBOJLJ9/4fMLlC6knAvwQ+DPxi\nks8zGbO4rnuLW5ncsGc7cCHwpW7c4ivAicDPJfnNbhs/muSTwOuZ3HXshcAvA1cnuZTuaqDd+/4M\ncFWSe7vtXTtaJ0jLwKvOSpKaPClPktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpKb/\nD/UkVLF1+ukaAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "nx = 1001 #number of spatial steps\n", + "nt = 500 #number of time steps\n", + "L = 20000 #length of river bed\n", + "dx = L/(nx-1)\n", + "x = numpy.linspace(0,L,nx) #x grid\n", + "\n", + "h = lahar_initial(nx, h_max, h_initial)\n", + "\n", + "pyplot.plot(x, h, linewidth = 3, ls = '-', color='#003366')\n", + "pyplot.ylim(-5,50);\n", + "pyplot.xlabel('Distance')\n", + "pyplot.ylabel('Height')" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "u_max = 4.0 #maximum velocity of the lahar, m/s\n", + "sigma = 1.0 \n", + "dt = sigma * dx / u_max" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "h_n, vel = lahar_maccormack(h, nt, dt, dx, beta, k)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Animating our results" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "from matplotlib import animation\n", + "from IPython.display import HTML" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def animate(data):\n", + " x = numpy.linspace(0, L, nx)\n", + " y = data\n", + " line.set_data(x,y)\n", + " return line," + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "fig = pyplot.figure();\n", + "ax = pyplot.axes(xlim=(0,L),ylim=(0,50),xlabel=('Distance'),ylabel=('Height'));\n", + "line, = ax.plot([],[],color='#003366',lw=3)\n", + "\n", + "anim = animation.FuncAnimation(fig, animate, frames = h_n, interval = 50)\n", + "HTML(anim.to_html5_video())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Calculations of velocity at bottom, size of deposit" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Velocity of deposit at bottom (m/s): 0.856168\n" + ] + } + ], + "source": [ + "#Velocity at the bottom calculation\n", + "vel=numpy.mean(vel)\n", + "print('Velocity of deposit at bottom (m/s): {0:f}'.format(vel))" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Size of deposit at bottom (m^3): 79552961.974933\n" + ] + } + ], + "source": [ + "#Size of the deposit at the bottom\n", + "deposit = numpy.sum(h_n*dx)\n", + "print('Size of deposit at bottom (m^3): {0:f}'.format(deposit))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Conclusions & challenge task" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have successfully modeled a lahar as a one dimensional, nonlinear advection equation. We used the MacCormack scheme to solve for the height of the lahar as a function of space and time. We were able to find the velocity and size of the lahar deposit and the bottom of our slope. For the challenge task, consider setting up the initial condition as a time dependent function for lahar height. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Dig Deeper ..." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This model of a lahar did not account for bulking and debulking. However, there are numerical models that have investigated adding in these terms-- see [5]. How could we manipulate our governing equation to include bulking and debulking? Would we continue to use the MacCormack scheme?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## References" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[1] Program, Volcano Hazards. \"USGS: Volcano Hazards Program.\" USGS: Volcano Hazards Program. N.p. n.d. Web. 5 Dec. 2016.\n", + "\n", + "[2] Pierson, Thomas C., Nathan J. Wood, and Carolyn L. Driedger. \"Reducing risk from lahar hazards: concepts, case studies, and roles for scientists.\" Journal of Applied Volcanology 3.1 (2014): 1\n", + "\n", + "[3] Slingerland, Rudy, and Lee R. Kump. Mathematical Modeling of Earth's Dynamical Systems: A Primer. Princeton, NJ: Princeton UP, 2011. Print\n", + "\n", + "[4] Numerical MOOC- Module 3, Lesson 2\n", + "\n", + "[5] Fagents, Sara A., and Stephen M. Baloga. \"Toward a Model for the Bulking and Debulking of Lahars.\" Journal of Geophysical Research 111.B10 (2006): n. pag. Web. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "from IPython.core.display import HTML\n", + "css_file = 'numericalmoocstyle.css'\n", + "HTML(open(css_file,\"r\").read())" + ] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python [Root]", + "language": "python", + "name": "Python [Root]" + }, + "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.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/alexabaumer/numericalmoocstyle.css b/alexabaumer/numericalmoocstyle.css new file mode 100644 index 0000000..93bbd51 --- /dev/null +++ b/alexabaumer/numericalmoocstyle.css @@ -0,0 +1,152 @@ + + + + + + + +