{
"metadata": {
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "raw",
"metadata": {},
"source": [
"Texto y c\u00f3digo sujeto bajo Creative Commons Attribution license, CC-BY-SA. (c) Original por Lorena A. Barba, 2013, traducido por F.J. Navarro-Brull para CAChemE.org "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[@LorenaABarba](https://twitter.com/LorenaABarba) - \n",
"[@CAChemEorg](https://twitter.com/cachemeorg)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"12 pasos para Navier-Stokes\n",
"=====\n",
"***"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Este notebook complementa el primer m\u00f3dulo interectivo online de [CFD con Python](https://bitbucket.org/franktoffel/cfd-python-class-es) y clases dadas por la profesora Lorena A. Barba, llamadas **12 Pasos para Navier-Stokes.**. En particular, estos notebook fueron escritos con la ayuda del estudiante de grado Gilbert Forsyth."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Operaciones matriciales (arrays) con NumPy\n",
"----------------\n",
"\n",
"Para los programas de c\u00e1lculo m\u00e1s intensivos, el uso de las funciones incorporadas en la librer\u00eda Numpy puede proporcionar varios aumentos en la velocidad de ejecuci\u00f3n. Como un ejemplo sencillo, considera la siguiente ecuaci\u00f3n:\n",
"\n",
"$$u^{n+1}_i = u^n_i-u^n_{i-1}$$\n",
"\n",
"Ahora, dado un vector $u^n = [0, 1, 2, 3, 4, 5]\\ \\ $, podemos calcular los valores de $u^{n+1}$ por iteraci\u00f3n sobre los valores de $u^n$ con un bucle `for`. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import numpy as np\n",
"\n",
"u = np.array((0, 1, 2, 3, 4, 5))\n",
"\n",
"for i in range(1,len(u)):\n",
" print u[i]-u[i-1]"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"1\n",
"1\n",
"1\n",
"1\n",
"1\n"
]
}
],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Este es el resultado esperado y el tiempo de ejecuci\u00f3n fue casi instant\u00e1neo. Si realizamos la misma operaci\u00f3n como una operaci\u00f3n vectorizada mediante el uso de arrays, en lugar de calcular $u^n_i-u^n_{i-1}\\ $ cinco veces por separado, se puede cortar (*slice*) el array $u$ y calcular cada operaci\u00f3n con un s\u00f3lo comando:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"u[1:]-u[0:-1]\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 7,
"text": [
"array([1, 2, 3, 4, 5])"
]
}
],
"prompt_number": 7
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Lo que dice este comando es lo siguiente, restar los cinco primeros elementos de $u$ a los cinco \u00faltimos del mismo.\n",
"\n",
"**Nota**: Recuerda que en Python se indican los primeros elementos del vector con el \u00edndice 0 y que los \u00edndices negativos dan \"la vuelta\" al vector. A diferencia de MATLAB, el \u00faltimo indice de una selecci\u00f3n no se incluye en la misma siendo en notaci\u00f3n matem\u00e1tica $[i,j)$.\n",
"\n",
"###Aumento de velocidad\n",
"\n",
"Para un array de seis elementos, los beneficios de las operaciones vectoriales (mediante arrays) son bastante escasos. No habr\u00e1 ninguna diferencia apreciable en el tiempo de ejecuci\u00f3n debido a las pocas operaciones que tienen lugar. Pero si volvemos a la convecci\u00f3n lineal 2D, podemos ver algunos aumentos significativos de velocidad."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"nx = 81\n",
"ny = 81\n",
"nt = 100\n",
"c = 1\n",
"dx = 2.0/(nx-1)\n",
"dy = 2.0/(ny-1)\n",
"sigma = .2\n",
"dt = sigma*dx\n",
"\n",
"x = np.linspace(0,2,nx)\n",
"y = np.linspace(0,2,ny)\n",
"\n",
"u = np.ones((ny,nx)) ## Crea un vector 1xn de unos\n",
"un = np.ones((ny,nx)) ##\n",
"\n",
"###Asigna las condiciones iniciales\n",
"\n",
"u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Con nuestras condiciones iniciales est\u00e1 todo listo, primero intentaremos ejecutar nuestro c\u00f3digo con el bucle anidado original, mediante la funci\u00f3n IPython \"magica\" `%%timeit`, que nos ayudar\u00e1 a evaluar el rendimiento de nuestro c\u00f3digo.\n",
"\n",
"** Nota **: La funci\u00f3n m\u00e1gica `%%timeit` ejecutar\u00e1 el c\u00f3digo varias veces y luego dar\u00e1 un tiempo promedio de ejecuci\u00f3n como resultado. Si tienes que representar alguna gr\u00e1fica en la celda donde se ejecuta `%%timeit`, se dibujar\u00e1n las figuras repetidas ocasiones por lo que puede ser un poco lioso. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%timeit\n",
"u = np.ones((ny,nx))\n",
"u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2\n",
"\n",
"for n in range(nt+1): ##loop across number of time steps\n",
" un[:] = u[:]\n",
" for i in range(1, len(u)):\n",
" for j in range(1, len(u)):\n",
" u[i,j] = un[i, j] - (c*dt/dx*(un[i,j] - un[i-1,j]))-(c*dt/dy*(un[i,j]-un[i,j-1]))\n",
" u[0,:] = 1\n",
" u[-1,:] = 1\n",
" u[:,0] = 1\n",
" u[:,-1] = 1"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"1 loops, best of 3: 6.37 s per loop\n"
]
}
],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Con el c\u00f3digo Python \"puro\" anterior, el mejor tiempo de ejecuci\u00f3n alcanzado fue de 6,24 segundos. Ten en cuenta que con estos tres bucles anidados, las sentencias dentro del bucle **j** se est\u00e1n evaluando m\u00e1s de 650.000 veces. Vamos a compararlo con el rendimiento del mismo c\u00f3digo vecotorizado (mediante arrays):"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%timeit\n",
"u = np.ones((ny,nx))\n",
"u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2\n",
"\n",
"for n in range(nt+1): ##loop across number of time steps\n",
" un[:] = u[:]\n",
" u[1:,1:]=un[1:,1:]-(c*dt/dx*(un[1:,1:]-un[0:-1,1:]))-(c*dt/dy*(un[1:,1:]-un[1:,0:-1]))\n",
" u[0,:] = 1\n",
" u[-1,:] = 1\n",
" u[:,0] = 1\n",
" u[:,-1] = 1"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"100 loops, best of 3: 8.24 ms per loop\n"
]
}
],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Como puedes ver, el aumento de velocidad es sustancial. El mismo c\u00e1lculo se reduce de 6,24 segundos a 9,59 milisegundos. 6 segundos no es una gran cantidad de tiempo de espera, pero las mejoras de velocidad aumentar\u00e1n exponencialmente con el tama\u00f1o y la complejidad del problema que se est\u00e1 evaluando. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = open(\"../styles/custom.css\", \"r\").read()\n",
" return HTML(styles)\n",
"css_styling()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
""
],
"output_type": "pyout",
"prompt_number": 1,
"text": [
""
]
}
],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"> (La celda de arriba establece el formato de este notebook)"
]
}
],
"metadata": {}
}
]
}