Laboratorios: Simple Navier-Stokes

Presenter Notes

Plan General de los Labs

  • CPU.
    • Paralelismo intra-CPU (ILP).
    • Jerarquía de memoria.
  • CPU Multicore.
    • Paralelización.
    • Jerarquía de memoria.
    • Escalabilidad.
  • GPU
    • Paralelización
    • Jerarquía de memoria.

Presenter Notes

Resumen:

  • El problema.
  • El programa.
  • Tareas Lab1.
  • Bibliografía

Nicolás Wolovick, $Date: 2012-04-17 14:27:51 -0300 (Tue, 17 Apr 2012) $, $Revision: 3418 $

Presenter Notes

El Problema

Equación de flujo de fluidos: Navier-Stokes

Navier-Stokes

  • Notación vectorial super-compacta.
  • Calcular como evoluciona un flujo (densidad y velocidad) en el tiempo.
  • Solo hay solución analítica para casos muy simples.

Presenter Notes

Qué es la ecuación?

Navier-Stokes

  • ∇: diferenciación por componente.
  • ∇^2: diferenciación doble por componente.
  • u: campo de velocidades.
  • ν: viscosidad.
  • f: campo fuente.
  • ρ: densidad.
  • κ: taza de difusión.
  • S: fuente.

Presenter Notes

Más sobre la ecuación

Navier-Stokes

  • La primera es no-linea.
  • La segunda es lineal.

Presenter Notes

El Programa

Presenter Notes

Discretización, la grilla

Grilla

  • Grilla bidimensional (puede ser 3d también).
  • Dentro de cada celda se supone u y ρ constante.

En C.

1 static u[size], v[size], u_prev[size], v_p rev[size];
2 static dens[size], dens_p rev[size];
3 
4 #define IX(i,j) ((i)+(N+2)*(j))

Presenter Notes

Evolución del ρ

Movimientos de densidades

  • Add Forces: + S
  • Diffuse: + κ∇^2ρ
  • Move: - (u ∇)ρ

Presenter Notes

Sumar la fuente, +S

1 void add_source ( int N, float * x, float * s, float dt )
2 {
3     int i, size=(N+2)*(N+2);
4     for ( i=0 ; i<size ; i++ ) x[i] += dt*s[i];
5 }

Presenter Notes

Difusión

Difusión

x0[IX(i-1,j)]+x0[IX(i+1,j)]+x0[IX(i,j-1)]+x0[IX(i,j+1)]-4*x0[IX(i,j)]

No sirve, es inestable.

 1 void diffuse_bad ( int N, int b, float * x, float * x0, float diff, float dt )
 2 {
 3     int i, j;
 4     float a=dt*diff*N*N;
 5     for ( i=1 ; i<=N ; i++ ) {
 6         for ( j=1 ; j<=N ; j++ ) {
 7             x[IX(i,j)] = x0[IX(i,j)] + a*(x0[IX(i-1,j)]+x0[IX(i+1,j)]+
 8                          x0[IX(i,j-1)]+x0[IX(i,j+1)]-4*x0[IX(i,j)]);
 9         }
10     }
11     set_bnd ( N, b, x );
12 }

Presenter Notes

Difusión (cont).

 1 void diffuse ( int N, int b, float * x, float * x0, float diff, float dt )
 2 {
 3     int i, j, k;
 4     float a=dt*diff*N*N;
 5     for ( k=0 ; k<20 ; k++ ) {
 6         for ( i=1 ; i<=N ; i++ ) {
 7             for ( j=1 ; j<=N ; j++ ) {
 8                 x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+
 9                 x[IX(i,j-1)]+x[IX(i,j+1)]))/(1+4*a);
10             }
11         }
12     set_bnd ( N, b, x );
13     }
14 }

Intercambia velocidad y exactitud por estabilidad.

Presenter Notes

Movimientos o advect

Advect

  • Linear backtrace.
  • Técnica estable, que no explota.

Presenter Notes

Advection

 1 void advect ( int N, int b, float * d, float * d0, float * u, float * v, float dt )
 2 {
 3     int i, j, i0, j0, i1, j1;
 4     float x, y, s0, t0, s1, t1, dt0;
 5     dt0 = dt*N;
 6     for ( i=1 ; i<=N ; i++ ) {
 7         for ( j=1 ; j<=N ; j++ ) {
 8             x = i-dt0*u[IX(i,j)]; y = j -dt0*v[IX(i,j)];
 9             if (x<0.5) x=0.5; if (x>N+0.5) x=N+ 0.5; i0=(int)x; i1=i0+ 1;
10             if (y<0.5) y=0.5; if (y>N+0.5) y=N+ 0.5; j0=(int)y; j1=j0+1;
11             s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1;
12             d[IX(i,j)] = s0*(t0*d0[IX(i0,j0)]+t1*d 0[IX(i0,j1)])+
13                 s1*(t0*d0[IX(i1,j0)]+t1*d 0[IX(i1,j1)]);
14         }
15     }
16     set_bnd ( N, b, d );
17 }

Presenter Notes

Paso de evolución de Densidad

1 void dens_step ( int N, float * x, float * x0, float * u, float * v, float diff,
2                  float dt )
3 {
4     add_source ( N, x, x0, dt );
5     SWAP ( x0, x ); diffuse ( N, 0, x, x0, diff, dt );
6     SWAP ( x0, x ); advect ( N, 0, x, x0, u, v, dt );
7 }

Donde SWAP intercambia dos punteros.

1 #define SWAP(x0,x) {float *tmp=x0;x0=x;x=tmp;}

Presenter Notes

Paso de evolución de Velocidad

Complicado de explicar, pero usa las mismas rutinas.

 1 void vel_step ( int N, float * u, float * v, float * u0, float * v0, float visc,
 2                 float dt )
 3 {
 4     add_source ( N, u, u0, dt ); add_source ( N, v, v0, dt );
 5     SWAP ( u0, u ); diffuse ( N, 1, u, u0, visc, dt );
 6     SWAP ( v0, v ); diffuse ( N, 2, v, v0, visc, dt );
 7     project ( N, u, v, u0, v0 );
 8     SWAP ( u0, u ); SWAP ( v0, v );
 9     advect ( N, 1, u, u0, u0, v0, dt ); advect ( N, 2, v, v0, u0, v0, dt );
10     project ( N, u, v, u0, v0 );
11 }

Faltaría project.

Presenter Notes

Resultado

Ver demo!

Presenter Notes

Tareas

  • Probar compiladores y sus opciones.
    • Recordar el problema del aliasing.
    • IPO.
    • Fast-math.
    • Autovectorización.
  • Explorar el espacio de N vs. performance.
  • Utilizar gprof para hacer profiling y ver que funciones toman la mayor parte del tiempo.
  • Opcionalmente medir ciclos, FPUops, instrucciones, etc., instrumentando el código con PAPI.
    • perf mide performance counters de punta a punta.

Presenter Notes

Tareas (cont.)

  • Optimizar el código utilizando estas técnicas entre otras:
    • Alineación.
    • Accesos cache-aware.
    • Caché blocking.
    • Loop unrolling manual.
    • Instrucciones SIMD.
    • AoS (Array of Structures), SoA (Structure of Arrays).
    • memset.

(Michael Gerndt, Efficient Programming of Multicore Processors and Supercomputers)

Presenter Notes

Bibliografía

Presenter Notes

Bibliografía

Presenter Notes

La clase que viene

Arquitecturas paralelas.

Presenter Notes