PDF de programación - Simulaciones de Física Estadística usando GPGPU: tres herramientas

Simulaciones de Física Estadística
usando GPGPU: tres herramientas

Alejandro. B. Kolton

CONICET, Centro Atómico Bariloche

WHPC13

Gracias:

● Ezequiel Ferrero [Grenoble]

● S. Bustingorry [Centro Atómico Bariloche].
● GPGPU@CAB [Centro Atómico Bariloche].
● Curso ICNPG 2012-2013 [Instituto Balseiro]
● NVIDIA Academic Partnership
● SNCAD

* Estamos organizando la 3er EAGPGPU en Bariloche (2014)

Mensaje de un inexperto

entusiasta en HPC ...

● Principiantes: cuando diga rápidamente en la charla que algo
es fácil, realmente lo es !

● Expertos: pueden ayudar muchísimo mas conociendo un
poquito de la problemática y perspectiva particular de cada
científico (esfuerzo mutuo de comunicación).

● Objetivo: intentar convencerlos, con ejemplos simples
concretos, de que es relativamente fácil empezar a desarrollar
códigos de alto rendimiento para innovar la investigación
científica, y de que definitivamente debemos fomentar mas la
sinergia interdisciplinaria en el país.

Física Estadística Computacional con GPGPU
Variedad de Modelos dinámicos

Con aplicación a la materia condensada

con desorden y fluctuaciones térmicas

● Vórtices en superconductores.
● Paredes de dominio Magnéticas.
● Paredes de dominio ferroeléctricas.
● Lineas de contacto de mojado.
● Modelos de terremotos/avalanchas.

●

“Hopping” de electrones localizados..

● Etc.

● Dificultad:

La física interesante “emerge” a
gran escala: tamaño/número de
partículas del sistema grande.

Necesidad de resolver numéricamente y de
realizar cómputos de alto rendimiento en
una variedad de modelos

Ej: modelar la dinámica de un
vórtice en un superconductor
con corriente

Una perspectiva

● Ciencia Básica: Preparar rápidamente un código de
simulación numérica como si se tratara de un “setup
experimental” para testear ideas...

Input
Ideas:
modeIos y
teorías

● Motivación
● Intuición
● Abstracción
● Creatividad.

Output:
¿ Entendemos mas ?
¿ Predicciones ?
¿ Comunicación ?

¿Implementación?
-->
¿ Que ?
¿ Como ?
¿ Con que ?
¿ Cuando ?

Optimizar a alto nivel!!

Computer Simulation of Liquids
Allen-Tildesley (1987)

Una perspectiva

Desafío:
Escribir rápidamente un código de simulación:
Correcto, Robusto, y de Buen rendimiento: Que permita obtener resultados
(estadísticos) exactos y suficientemente precisos para testear las hipótesis...
pero no mas !
Alto nivel y en un lenguaje “expresivo”: Pensar mas en “que” quiero calcular,
que en el “como” hacerlo... y que sea fácil exponerlo e implementarlo.

● GPGPU:

➢ Buenos speed-ups para muchos problemas con poco esfuerzo!.
➢ Relativamente baratas para el poder de cálculo que brindan.
➢ Toda una variedad creciente de librerías activamente mantenidas, muchas

de ellas gratis, útiles para aplicaciones científicas (expresivas y productivas).


Pero antes de largarse a programar ...

● BUSCAR, DISCUTIR, PROBAR (95% del laburo):

Que herramientas de GPGPU se adecuan mejor a estas
necesidades particulares ?

Un tour sobre tres herramientas

C++ template library

Random Numbers
Random123:
a Library of Counter-Based
Random Number Generators

GPU y CPU multicore

CUFFT

GPU y CPU multicore

y su uso en dos ejemplos concretos

de investigación científica

Que es Thrust ?

http://thrust.github.io/

//#includes...
int main(void)
{
// generate 32M random numbers serially
thrust::host_vector<int> h_vec(32 << 20);
std::generate(h_vec.begin(), h_vec.end(), rand);

// transfer data to the device
thrust::device_vector<int> d_vec = h_vec;

// sort data on the device (846M keys per second on GeForce GTX 480)
thrust::sort(d_vec.begin(), d_vec.end());

// transfer data back to host
thrust::copy(d_vec.begin(), d_vec.end(), h_vec.begin());

return 0;
}



nvcc -O2 -arch=sm_21 programa.cu -o ejecutable



Que tipo de algoritmos tiene Thrust ?

Fundamentalmente
“patrones paralelos”



Reducción en árbol

3

1

7

0

4

1

6

3

+

4

+

7

+

5

+

9

+

11

+

14

+

25

Calcula la suma

Vale para int, float, double,o cualquier tipo de dato con la operación “+” bien definida...

Reducción

3

1

7

0

“>”

3

“>”

7

“>”

7

4

“>”

1

6

“>”

3

4

6

“>”

6

A “>” B = (A>B)?(A):(B)

“>”

7

Calcula el máximo

Vale para int, float, double,o cualquier tipo de dato con la operación “>” bien definida...

Reducción tenística

“play”

“play”

“play”

Mutua Madrid Open
Spain - 07.05.2012-13.05.2012
Draw: 64
Surface: Clay
Prize Money: €3,090,150
63 partidos en 6 jornadas

“play”

Calcula el campeón

ATP MADRID:
64 players
63 partidos en
6 jornadas

Name: Roland Garros
Category: Grand Slam
Place: Paris, France
Date: 27.05.2012-
10.06.2012
Draw Size: S-128 D-64

log_2(N) operaciones:
Grand-Grand-Slam: 2^20 = 1048576 jugadores → solamente 20 “jornadas” !!!

Reducción en Thrust

Cualquier operacion binaria asoc. y conmutativa
(+,-,*, max, min, etc, o user-defined)

Generic algorithms

… definida sobre cualquier
Estructura de datos
(int, float, etc, o user-defined),
Viva en el HOST o en el DEVICE (GPU)

Thrust

Thrust Content from GTC 2012: Nathan Bell



Principales desarroladores: Jared Hoberock y Nathan Bell

Thrust

Thrust Content from GTC 2012: Nathan Bell



Thrust

Thrust Content from GTC 2012: Nathan Bell



Thrust

MISMO código
Corre en CPU
Multicore!!

Thrust Content from GTC 2012: Nathan Bell

Delegamos la implementación
de bajo nivel o “mapeo al hardware”
a la librería



Thrust

CUDA Toolkit 5.0 Performance Report



Thrust

Thrust Content from GTC 2012: Nathan Bell

La mayoría son
“Parallel primitives”

Thrust se puede usar
para alivianar la escritura
de gran parte del código:

Alocacion de memoria,
copias, y composición de
patrones paralelas.



Limitaciones de Thrust

Review de M. Harris en:
https://developer.nvidia.com/content/expressive-algorithmic-programming-thrust

GTC 2012 Nathan Bell

o
g
o

l

o
t
u
p
m
o
c
-
o
c

z
a
f
r
e
t
n

I

i

í

s
F



Usualmente en Cuda C/C++



Un ejemplo concreto de

aplicación científica

• Modelo para la dinámica de una interfaz en un
➢ Matemáticamente: PDE estocástica en función del

medio desordenado:

tiempo para un campo escalar unidimensional en un
potencial desordenado.

➢ Con los dedos: una cuerdita elástica que se mueve

en un paisaje desordenado.

Interfaces en medios desordenados

Magnetic domain wall in a Pt/Co/Pt film in the creep regime.
Lemerle, Jamet, Ferre, et al (LPS Orsay).

Un modelo mínimo para capturar la
dinámica “universal” de las paredes

Desorden congelado

Desorden dinámico (térmico)



Un modelo mínimo para capturar la
dinámica “universal” de las paredes

u

u

u

u

u

u'>u

VS

VS



suaviza

distorsiona



empuja

Implementación numérica del modelo

interface~polimero

Implementación mas simple: diferencias finitas, Euler explicito



Implementación numérica del modelo

Implementación mas simple: explicit finite-difference

“Casi Embarasosamente Paralelo”

1) “Parallel For”(X=0,...,L-1)

Thread “X” calcula la fuerza en X, Ftot(X) → Array “Ftot” (device memory)

2) “Parallel For”(X=0,...,L-1)

Thread “X” calcula u(X,n+1) → Nuevo array “u” (device memory)

3) “parallel reductions” : calculo cantidades de interés (o visualización ...)

4) Vuelvo a (1), en un “For(n=0,1,...,Trun)” de evolución temporal



Implementación mas naive en GPU

// varios #include<...>
using namespace thrust;
int main()
{

// vectores que viven en la device (GPU global memory)
// declaracion y alocacion:
device_vector<REAL> u(L); // mi pared (discretizada)
device_vector<REAL> Ftot(L); // la fuerza sobre mi pared

// iteradores sobre GPU global memory: definen rangos
// necesarios para aplicar algoritmos paralelos genericos
device_vector<REAL>::iterator u_begin = u.begin();
device_vector<REAL>::iterator Ftot_begin = Ftot.begin();

// genero una evolucion temporal de la pared
for(long n=0;n<Trun;n++)
{

....
// (1) calculo paralelo Ftot(X,n) en la GPU
....
// (2) calculo paralelo u(X,n+1) en la GPU
.....
// (3) calculo paralelo propiedades en GPU
.....
// copio de mGPU -> mCPU, lo mínimo necesario!

}
….

}



Implementación en GPU

Paso de Euler Paralelo

int main()
{

Functor

....

#define Dt
struct paso_de_Euler
{
__device__
REAL operator()(REAL u_old, REAL Ftot)
{

// u_old == u(X,n)
// Ftot == Ftot(X,n)
return (u_old + Ftot*Dt); // → u(X,n+1)

}
}

// vectores
device_vector<REAL> u(L); // mi pared (discretizada)
device_vector<REAL> Ftot(L); // la fuerza sobre mi pared

// iteradores
device_vector<REAL>::iterator u_begin = u.begin();
device_vector<REAL>::iterator Ftot_begin = Ftot.begin();

// genero una evolución temporal de la pared
for(long n=0;n<Trun;n++)
{

....
// (1) calculo Ftot(X,n) en el device (GPU)
....
// (2) avanzo un paso de Euler:
transform(
u_begin,u_end, Ftot_begin,
u_begin,
paso_de_Euler()
);
.....
// calculo propiedades o visualizo (GPU)
.....
// copio de memGPU -> memCPU, lo mínimo necesario!

Inputs (y su rango)
Output (in-place)
functor

“Parecido” a un kernel de cuda pero:

● La implementación de mas bajo nivel
esta delegada (ej: en threads/blocks
de Cuda).

● La interfaz es uniforme y genérica: el
input/output y las operaciones pueden

ser user-defined (classes).

}
….

}

Implementación en GPU

Calculo de las fuerzas

struct fuerza
{
long t;
fuerza(long _t):t(_t){}; // t == "estado interno"

__device__
REAL operator()(tuple<int,REAL,REAL,REAL> tt)
{

int tid = get<0>(tt); // tid==X

REAL um1 = get<1>(tt); // u(X-1)
REAL u = get<2>(tt); // u(X)
REAL up1 = get<3>(tt);

// u(X+1)

?

// ... desorden = RN(tid,[u])+RN(tid,[u]+1)(u-[u]);
// ... ruido = RN(tid,t);

REAL fuerza =
(up1+um1-2*u) + desorden + ruido;

return fuerza; // → Ftot(X)

}
}



}
….
}

int main()
{
// vectores
device_vector<REAL> u(L); // mi pared (discretizada)
devic