Texto y código 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

12 pasos para Navier-Stokes


Este notebook complementa el primer módulo interectivo online de CFD con Python 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.

Operaciones matriciales (arrays) con NumPy

Para los programas de cálculo más intensivos, el uso de las funciones incorporadas en la librería Numpy puede proporcionar varios aumentos en la velocidad de ejecución. Como un ejemplo sencillo, considera la siguiente ecuación:

$$u^{n+1}_i = u^n_i-u^n_{i-1}$$

Ahora, dado un vector $u^n = [0, 1, 2, 3, 4, 5]\ \ $, podemos calcular los valores de $u^{n+1}$ por iteración sobre los valores de $u^n$ con un bucle for.

In [2]:
import numpy as np

u = np.array((0, 1, 2, 3, 4, 5))

for i in range(1,len(u)):
    print  u[i]-u[i-1]
1
1
1
1
1

Este es el resultado esperado y el tiempo de ejecución fue casi instantáneo. Si realizamos la misma operación como una operación 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ón con un sólo comando:

In [7]:
u[1:]-u[0:-1]
Out[7]:
array([1, 2, 3, 4, 5])

Lo que dice este comando es lo siguiente, restar los cinco primeros elementos de $u$ a los cinco últimos del mismo.

Nota: Recuerda que en Python se indican los primeros elementos del vector con el índice 0 y que los índices negativos dan "la vuelta" al vector. A diferencia de MATLAB, el último indice de una selección no se incluye en la misma siendo en notación matemática $[i,j)$.

Aumento de velocidad

Para un array de seis elementos, los beneficios de las operaciones vectoriales (mediante arrays) son bastante escasos. No habrá ninguna diferencia apreciable en el tiempo de ejecución debido a las pocas operaciones que tienen lugar. Pero si volvemos a la convección lineal 2D, podemos ver algunos aumentos significativos de velocidad.

In [1]:
nx = 81
ny = 81
nt = 100
c = 1
dx = 2.0/(nx-1)
dy = 2.0/(ny-1)
sigma = .2
dt = sigma*dx

x = np.linspace(0,2,nx)
y = np.linspace(0,2,ny)

u = np.ones((ny,nx)) ## Crea un vector 1xn de unos
un = np.ones((ny,nx)) ##

###Asigna las condiciones iniciales

u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2

Con nuestras condiciones iniciales está todo listo, primero intentaremos ejecutar nuestro código con el bucle anidado original, mediante la función IPython "magica" %%timeit, que nos ayudará a evaluar el rendimiento de nuestro código.

Nota : La función mágica %%timeit ejecutará el código varias veces y luego dará un tiempo promedio de ejecución como resultado. Si tienes que representar alguna gráfica en la celda donde se ejecuta %%timeit, se dibujarán las figuras repetidas ocasiones por lo que puede ser un poco lioso.

In [2]:
%%timeit
u = np.ones((ny,nx))
u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2

for n in range(nt+1): ##loop across number of time steps
    un[:] = u[:]
    for i in range(1, len(u)):
        for j in range(1, len(u)):
            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]))
            u[0,:] = 1
            u[-1,:] = 1
            u[:,0] = 1
            u[:,-1] = 1
1 loops, best of 3: 6.37 s per loop

Con el código Python "puro" anterior, el mejor tiempo de ejecución alcanzado fue de 6,24 segundos. Ten en cuenta que con estos tres bucles anidados, las sentencias dentro del bucle j se están evaluando más de 650.000 veces. Vamos a compararlo con el rendimiento del mismo código vecotorizado (mediante arrays):

In [3]:
%%timeit
u = np.ones((ny,nx))
u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2

for n in range(nt+1): ##loop across number of time steps
    un[:] = u[:]
    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]))
    u[0,:] = 1
    u[-1,:] = 1
    u[:,0] = 1
    u[:,-1] = 1
100 loops, best of 3: 8.24 ms per loop

Como puedes ver, el aumento de velocidad es sustancial. El mismo cálculo 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án exponencialmente con el tamaño y la complejidad del problema que se está evaluando.

In [1]:
from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
Out[1]:

(La celda de arriba establece el formato de este notebook)