La tarea de N cuerpos o c贸mo volar una galaxia sin salir de la cocina




No hace mucho le铆 la novela de ciencia ficci贸n "El problema de los tres cuerpos" de Liu Qixin . En 茅l, algunos extraterrestres ten铆an un problema: no sab铆an c贸mo, con suficiente precisi贸n para ellos, calcular la trayectoria de su planeta de origen. A diferencia de nosotros, viv铆an en un sistema de tres estrellas, y el "clima" del planeta depend铆a en gran medida de su disposici贸n mutua, desde el calor incinerador hasta el fr铆o helado. Y decid铆 verificar si podemos resolver tales problemas.

F铆sica del fen贸meno.


Para comprender el problema, es necesario tratar con la f铆sica del fen贸meno. En el marco de la teor铆a cl谩sica, la fuerza de atracci贸n de dos cuerpos est谩 determinada por la ley de Newton :

 v e c F ( v e c r 1 , v e c r 2 ) = - G m 1 m 2 f r a c v e c r 1 - v e c r 2 | v e c r 1 - v e c r 2 | 3 ,       



donde  vecr1, vecr2 - la posici贸n de los cuerpos en el espacio, m1,m2 - masas de cuerpos, G - constante gravitacional.
En el sistema de N Los cuerpos en cada uno de ellos se ver谩n afectados por la fuerza de atracci贸n del resto, que se expresa mediante la ecuaci贸n:

 vecFn=G sumk neqnmnmk frac vecrn vecrk| vecrn vecrk|3.$



Usando la segunda ley de Newton, escribimos la aceleraci贸n para cada part铆cula:

 vecan= vecFn/mn=G sumk neqnmk frac vecrn vecrk| vecrn vecrk|3.$



Recordando que la aceleraci贸n es la segunda derivada del tiempo de la coordenada, obtenemos una ecuaci贸n diferencial parcial de segundo orden, que debe resolverse para obtener la trayectoria de cada cuerpo:

 frac partial2 vecrn partialt2=fn=G sumk neqnmk frac vecrn vecrk| vecrn vecrk|3.$


Es importante se帽alar aqu铆 que la complejidad de calcular una funci贸n fn es igual a O(N2) y aumenta significativamente con el aumento del n煤mero de cuerpos que interact煤an.

Matem谩ticas


El primer y m谩s simple m茅todo para resolver ecuaciones diferenciales es el m茅todo de Euler , que est谩 dise帽ado para resolver ecuaciones de la forma:

 fracdydx=f(x,y).$


Tras la transici贸n a la regi贸n discreta, obtenemos:

yi=yi1+hf(xi1,yi1), quadi=1,2,3, puntos,m,


donde h Es el paso de integraci贸n, y m - el n煤mero de pasos de integraci贸n. Por lo tanto, si necesitamos calcular la posici贸n de los cuerpos a la vez T entonces deber铆amos hacer m=T/h pasos de integraci贸n. Aqu铆 el primer problema es visible, si T grande, entonces debemos tomar una gran cantidad de pasos de integraci贸n.

Para aplicar el m茅todo de Euler a nuestro problema, debe reducirse a un sistema de primer orden. Para hacer esto, presentamos una variable adicional: velocidad de part铆culas:

 frac partial vecvn partialt=fn, frac partial vecrn partialt= vecvn.$



El segundo problema para resolver sistemas de ecuaciones diferenciales es la precisi贸n de la soluci贸n y su control. La precisi贸n se puede mejorar de dos maneras: disminuyendo el paso de integraci贸n y eligiendo un m茅todo con un mayor orden de precisi贸n. Ambos m茅todos conducen a una mayor complejidad computacional, pero de diferentes maneras. Por ejemplo, puede usar el m茅todo cl谩sico Runge-Kutta de cuarto orden ; requiere cuatro c谩lculos de funciones fn en cada paso, pero tiene un orden de precisi贸n O(h4) (a modo de comparaci贸n, el m茅todo de Euler tiene un orden de precisi贸n O(h) y requiere un c谩lculo fn ) La precisi贸n de la soluci贸n tambi茅n se puede controlar de varias maneras: compare con la soluci贸n anal铆tica, resuelva utilizando diferentes m茅todos o con diferentes pasos y compare los resultados, controle los par谩metros de terceros y las limitaciones que la soluci贸n debe cumplir.
Adem谩s, cada uno de estos m茅todos tiene sus inconvenientes. Las decisiones anal铆ticas pueden estar ausentes o, incluso en la mayor铆a de los casos, completamente ausentes. Por ejemplo, para nuestra tarea N la soluci贸n anal铆tica de tel es solo para N=2 , pero incluso esto es suficiente para probar la precisi贸n de los m茅todos. Resolver un problema mediante dos m茅todos o con diferentes pasos aumenta el tiempo de c谩lculo, pero este enfoque se puede aplicar a casi cualquier tarea. No todas las tareas tienen limitaciones, pero las nuestras las tienen: en cada paso de integraci贸n podemos controlar la implementaci贸n de las leyes de conservaci贸n . Este enfoque tambi茅n aumenta la complejidad del c谩lculo, pero hay mucho para elegir; calcular la suma de los momentos moment谩neos o angulares de todas las part铆culas tiene una complejidad del orden O(N) , mientras que calcular la energ铆a total de un sistema tiene complejidad de orden O(N2)

Nota sobre el c谩lculo de la energ铆a total
En nuestro caso, la energ铆a total del sistema consta de dos partes: energ铆a cin茅tica y energ铆a potencial. La energ铆a cin茅tica consiste en la suma de las energ铆as cin茅ticas de todos los cuerpos. Para calcular la energ铆a potencial, necesitamos agregar las energ铆as potenciales de cada part铆cula en el campo gravitacional de las part铆culas restantes, por lo que debemos agregar O(N2) t茅rminos La dificultad es que todos los t茅rminos son de un orden muy diferente, e incluso con c谩lculos de doble precisi贸n no es posible calcular este valor con una precisi贸n suficiente para comparar en diferentes pasos. Para superar este problema, fue necesario aplicar la suma de acuerdo con el algoritmo de Kahan .

Trayectoria el铆ptica

Fig 1. Un ejemplo de una trayectoria el铆ptica.




Considere el caso simple de un sat茅lite que se mueve en una 贸rbita el铆ptica alrededor de la Tierra. A medida que el sat茅lite se acerca a la Tierra, su velocidad aumenta, y al alejarse de la Tierra disminuye, en consecuencia, surge la posibilidad de disminuir el paso de integraci贸n en el tiempo en la parte verde de la 贸rbita, y aumentarlo en rojo y azul sin cambiar la precisi贸n de la soluci贸n. Intentemos comparar con m谩s detalle.

Tabla 1. M茅todos investigados para resolver ecuaciones diferenciales
MarcadoOrdenDescripci贸n
adams1-5M茅todo Adams-Bashfort
Euler1M茅todo Euler
rk44 4M茅todo cl谩sico de Runge-Kutta
rkck5 5M茅todo Kash-Karp
rkdp5 5M茅todo Dorman-Prince
rkdverk6 6M茅todo Werner 1) p. 181
rkf7 7M茅todo Felberg 1) p. 180
rkgl6 6M茅todo impl铆cito de Gauss-Legendre
rklc4 4M茅todo impl铆cito de Lobatto
trapecio2M茅todo trapezoidal

Para seleccionar el mejor m茅todo para nuestra tarea, compararemos varios m茅todos conocidos. Para hacer esto, simulamos la colisi贸n de dos sistemas de cuerpos. N=512 y mida el cambio relativo en la energ铆a total , el momento y su momento al final de la simulaci贸n (tiempo m谩ximo de simulaci贸n) Tmax=2.5 ) En este caso, variaremos el paso y los par谩metros de los m茅todos de integraci贸n y mediremos el n煤mero de llamadas a funciones fn , respectivamente, aquellos m茅todos que con un menor n煤mero de llamadas conducir谩n a una menor p茅rdida, consideraremos m谩s aceptables.





a)b)
Figura 2. Cambio relativo en energ铆a a), momento b), al final de la simulaci贸n del sistema N=512 cuerpos por varios m茅todos dependiendo del n煤mero de c谩lculos de funciones fn doble precisi贸n

De los gr谩ficos de la Figura 2 se puede ver que la mejor relaci贸n entre la cantidad de c谩lculo de la funci贸n fn y el cambio relativo en la energ铆a del sistema de cuerpos en los m茅todos Adams y Dorman-Prince de quinto orden. Tambi茅n se ve que para todos los m茅todos con un aumento en el n煤mero de c谩lculos fn El cambio relativo en el impulso del sistema aumenta. Para un cambio relativo en la energ铆a, esto tambi茅n es notable, pero solo para algunos m茅todos que podr铆an alcanzar el umbral dE/E0<1012 . Este efecto ya no se puede atribuir al error de los m茅todos, sino al error de los c谩lculos, y un aumento adicional en la precisi贸n solo es posible junto con un aumento en la precisi贸n de los c谩lculos con coma flotante.




a)b)
Fig. 3. Cambio relativo en energ铆a a), momento b), al final de la simulaci贸n del sistema N=512 cuerpos por varios m茅todos dependiendo del n煤mero de c谩lculos de funciones fn con precisi贸n cu谩druple (__float128)

Las Figuras 3a y 3b muestran que el uso de c谩lculos con precisi贸n cu谩druple puede reducir las p茅rdidas de energ铆a relativas hasta 1023 , pero debe comprender que el tiempo de c谩lculo aumenta en dos 贸rdenes de magnitud en comparaci贸n con la precisi贸n doble. Como en el caso de los c谩lculos de doble precisi贸n, la mejor relaci贸n de precisi贸n y el n煤mero de c谩lculos de funciones fn poseen m茅todos de quinto orden Adams y Dorman-Prince.

Los m茅todos Dorman-Prince y Werner pertenecen a la clase de m茅todos anidados y pueden calcular simult谩neamente dos soluciones con un orden de precisi贸n alto y bajo (para el m茅todo Dorman-Prince, 贸rdenes 5 y 4, y para el m茅todo Werner, 贸rdenes 6 y 5). Si estas dos soluciones son muy diferentes, entonces podemos dividir el paso de integraci贸n actual en otras m谩s peque帽as. Eso nos permite cambiar din谩micamente el paso de integraci贸n y reducirlo solo en aquellas 谩reas donde se requiere.

Comparemos los m茅todos de quinto orden de Dorman-Prince, Werner y Adams con m谩s detalle, durante un intervalo de simulaci贸n m谩s largo de nuestro sistema ( Tmax=300 )


a) El cambio relativo en energ铆a, momento y su momento en el proceso de modelado


b) N煤mero de c谩lculos de funciones fn en el intervalo de simulaci贸n  DeltaT=0.3
Fig. 4. Dependencias del cambio relativo en los par谩metros f铆sicos del sistema y el n煤mero de c谩lculos de funciones. fn de vez en cuando modelando con el m茅todo de paso variable Dorman-Prince h=105 puntos109



a) El cambio relativo en energ铆a, momento y su momento en el proceso de modelado


b) N煤mero de c谩lculos de funciones fn en el intervalo de simulaci贸n  DeltaT=0.3
Fig. 5. Dependencias del cambio relativo en los par谩metros f铆sicos del sistema y el n煤mero de c谩lculos de funciones. fn versus modelado de tiempo por el m茅todo de Werner con un paso variable h=105 puntos109



Fig. 6. Cambio relativo en energ铆a, momento y su momento en el proceso de modelado por el m茅todo Adams-Bashfort de quinto orden con un paso h=106



Fig 7. Dependencias del n煤mero de c谩lculos de funciones. fn para m茅todos Adams, Dorman-Prince y Werner de quinto orden a partir del tiempo de simulaci贸n

Se puede ver que en la etapa inicial de la evoluci贸n de nuestro sistema ( T<25 ) los tres m茅todos muestran caracter铆sticas similares, pero en etapas posteriores se producen algunos eventos en el sistema, como resultado de los cuales los errores en los par谩metros principales del sistema (energ铆a total, momento y su momento) saltan bruscamente. Pero los m茅todos de Dorman-Prince y Werner hacen frente a estos cambios mucho mejor debido a la capacidad de reducir el paso de integraci贸n en secciones "complejas", como resultado de lo cual aumenta el n煤mero de c谩lculos de funciones fn como se ve en las Figuras 4b y 5b, pero el n煤mero total de c谩lculos fn Los m茅todos anidados son m谩s peque帽os que el m茅todo Adams-Bashfort, como se puede ver en la Figura 7.

Me pregunto qu茅 pas贸 con el sistema en estos momentos.

Video 1. Modelado de un sistema de 512 cuerpos. M茅todo Dorman-Prince. Paso din谩mico h=105 puntos109


El video demuestra eso hasta el momento T=25 el movimiento es relativamente tranquilo, y despu茅s de una colisi贸n de los centros de las "galaxias", lo que conduce a un cambio brusco en las trayectorias y un fuerte aumento en las velocidades de algunas part铆culas. Adem谩s, para mantener la precisi贸n de la soluci贸n, es necesario reducir el paso de integraci贸n. Los m茅todos anidados pueden hacer esto autom谩ticamente; los gr谩ficos muestran que en algunas partes de la evoluci贸n del sistema, el paso de integraci贸n se redujo en casi dos 贸rdenes de magnitud con 105 antes h=107 . Al usar el m茅todo Adams y tal paso en todo el intervalo de la evoluci贸n del sistema, no hubi茅ramos esperado una soluci贸n.

Resumen


Para la soluci贸n, es mejor utilizar m茅todos anidados que le permitan controlar din谩micamente el paso de integraci贸n y reducirlo solo en secciones "complejas" de la trayectoria.

No persiga los m茅todos de orden superior. Incluso cuando usan el tipo de datos 'doble', no alcanzan sus capacidades potenciales, el uso de los tipos de datos con mayor precisi贸n aumenta en gran medida el tiempo que lleva resolver el problema.

Implementaci贸n de la CPU


Ahora que la elecci贸n de un m茅todo para resolver ecuaciones est谩 definida, intentemos calcular el c谩lculo de la fuerza de interacci贸n para cada part铆cula. Obtenemos un doble ciclo para todas las part铆culas:

C贸digo de implementaci贸n 'simple'
for(size_t body1 = 0; body1 < count; ++body1) { const nbvertex_t v1(rx[ body1 ], ry[ body1 ], rz[ body1 ]); nbvertex_t total_force; for(size_t body2 = 0; body2 != count; ++body2) { if(body1 == body2) { continue; } const nbvertex_t v2(rx[ body2 ], ry[ body2 ], rz[ body2 ]); const nbvertex_t force(m_data->force(v1, v2, mass[body1], mass[body2])); total_force += force; } frx[body1] = vx[body1]; fry[body1] = vy[body1]; frz[body1] = vz[body1]; fvx[body1] = total_force.x / mass[body1]; fvy[body1] = total_force.y / mass[body1]; fvz[body1] = total_force.z / mass[body1]; } 


Las fuerzas gravitacionales para cada cuerpo se calculan de forma independiente, y para usar todos los n煤cleos del procesador, es suficiente escribir la directiva OpenMP antes del primer ciclo:

Un fragmento de c贸digo de la implementaci贸n de 'openmp'
 #pragma omp parallel for for(size_t body1 = 0; body1 < count; ++body1) 


Porque Dado que cada cuerpo interact煤a con cada uno, para reducir la cantidad de interacciones del procesador con la RAM y mejorar el uso de la memoria cach茅, tenemos la capacidad de cargar parte de los datos en la memoria cach茅 y usarla repetidamente:

C贸digo de implementaci贸n 'openmp + block'
 #pragma omp parallel for for(size_t n1 = 0; n1 < count; n1 += BLOCK_SIZE) { nbcoord_t x1[BLOCK_SIZE]; nbcoord_t y1[BLOCK_SIZE]; nbcoord_t z1[BLOCK_SIZE]; nbcoord_t m1[BLOCK_SIZE]; nbvertex_t total_force[BLOCK_SIZE]; for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { size_t local_n1 = b1 + n1; x1[b1] = rx[local_n1]; y1[b1] = ry[local_n1]; z1[b1] = rz[local_n1]; m1[b1] = mass[local_n1]; } for(size_t n2 = 0; n2 < count; n2 += BLOCK_SIZE) { nbcoord_t x2[BLOCK_SIZE]; nbcoord_t y2[BLOCK_SIZE]; nbcoord_t z2[BLOCK_SIZE]; nbcoord_t m2[BLOCK_SIZE]; for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2) { size_t local_n2 = b2 + n2; x2[b2] = rx[local_n2]; y2[b2] = ry[local_n2]; z2[b2] = rz[local_n2]; m2[b2] = mass[n2 + b2]; } for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { const nbvertex_t v1(x1[ b1 ], y1[ b1 ], z1[ b1 ]); for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2) { const nbvertex_t v2(x2[ b2 ], y2[ b2 ], z2[ b2 ]); const nbvertex_t force(m_data->force(v1, v2, m1[b1], m2[b2])); total_force[b1] += force; } } } for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { size_t local_n1 = b1 + n1; frx[local_n1] = vx[local_n1]; fry[local_n1] = vy[local_n1]; frz[local_n1] = vz[local_n1]; fvx[local_n1] = total_force[b1].x / m1[b1]; fvy[local_n1] = total_force[b1].y / m1[b1]; fvz[local_n1] = total_force[b1].z / m1[b1]; } } 


Una optimizaci贸n adicional consiste en extraer el contenido de la funci贸n de calcular la fuerza en el ciclo principal y eliminar la divisi贸n y multiplicaci贸n por la masa corporal m1 [b1]. Adem谩s del hecho de que hemos reducido un poco los c谩lculos, el compilador podr谩 aplicar instrucciones vectoriales del procesador SSE y AVX en un ciclo tan extendido.

C贸digo de implementaci贸n 'openmp + block + optimization'
 #pragma omp parallel for for(size_t n1 = 0; n1 < count; n1 += BLOCK_SIZE) { nbcoord_t x1[BLOCK_SIZE]; nbcoord_t y1[BLOCK_SIZE]; nbcoord_t z1[BLOCK_SIZE]; nbcoord_t total_force_x[BLOCK_SIZE]; nbcoord_t total_force_y[BLOCK_SIZE]; nbcoord_t total_force_z[BLOCK_SIZE]; for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { size_t local_n1 = b1 + n1; x1[b1] = rx[local_n1]; y1[b1] = ry[local_n1]; z1[b1] = rz[local_n1]; total_force_x[b1] = 0; total_force_y[b1] = 0; total_force_z[b1] = 0; } for(size_t n2 = 0; n2 < count; n2 += BLOCK_SIZE) { nbcoord_t x2[BLOCK_SIZE]; nbcoord_t y2[BLOCK_SIZE]; nbcoord_t z2[BLOCK_SIZE]; nbcoord_t m2[BLOCK_SIZE]; for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2) { size_t local_n2 = b2 + n2; x2[b2] = rx[local_n2]; y2[b2] = ry[local_n2]; z2[b2] = rz[local_n2]; m2[b2] = mass[n2 + b2]; } for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2) { nbcoord_t dx = x1[b1] - x2[b2]; nbcoord_t dy = y1[b1] - y2[b2]; nbcoord_t dz = z1[b1] - z2[b2]; nbcoord_t r2(dx * dx + dy * dy + dz * dz); if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = (m2[b2]) / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; total_force_x[b1] -= dx; total_force_y[b1] -= dy; total_force_z[b1] -= dz; } } } for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { size_t local_n1 = b1 + n1; frx[local_n1] = vx[local_n1]; fry[local_n1] = vy[local_n1]; frz[local_n1] = vz[local_n1]; fvx[local_n1] = total_force_x[b1]; fvy[local_n1] = total_force_y[b1]; fvz[local_n1] = total_force_z[b1]; } } 


Tabla 2. Dependencia del tiempo de c谩lculo (en segundos) de la funci贸n fn sobre la cantidad de cuerpos que interact煤an N para diversas implementaciones de CPU
N2048409681921638432768
simple0,04250.16510.65942,6510,52
Openmp0.00780,02600,10790.4171.655
openmp + block + optimizaci贸n0.00370,01280,04950,1940,774

Par谩metros del sistema:
  • sistema: Debian 9, Intel Core i7-5820K (6 n煤cleos)
  • compilador: gcc 6.3.0


Se ve claramente que la versi贸n con soporte de OpenMP se acelera seis veces, exactamente en t茅rminos de la cantidad de n煤cleos, y la versi贸n optimizada es m谩s r谩pida un poco m谩s de dos veces. Entonces, con la optimizaci贸n, no debe confiar solo en la concurrencia. Curiosamente, en los c谩lculos de flujo 煤nico, el procesador trabaj贸 a una frecuencia de 3.6 GHz, en la versi贸n paralela (openmp) baj贸 la frecuencia a 3.4 GHz, y en el paralelo y optimizado (openmp + block + optimizaci贸n) baj贸 a 3.3 GHz, pero esto no le imped铆a trabajar 13,6 veces m谩s r谩pido. Tambi茅n se ve que un aumento en el tiempo de c谩lculo con un aumento en el tama帽o del problema es cuadr谩tico, y un aumento adicional N hace que la tarea sea insoluble en un tiempo razonable.

Implementaci贸n de GPU


Pero existe el deseo de hacer c谩lculos a煤n m谩s r谩pidos. Hay varias direcciones disponibles para la aceleraci贸n: c谩lculo de GPU, aproximaci贸n de funciones fn . Primero, para la inform谩tica de GPU, eleg铆 la tecnolog铆a OpenCL. Para un trabajo m谩s conveniente, se utiliz贸 la biblioteca CLHPP . La principal ventaja de OpenCL es que el c贸digo se puede ejecutar en el procesador y en la GPU, lo que simplifica la escritura y la depuraci贸n, adem谩s de expandir la lista de hardware para ejecutar. La herramienta Oclgrind ayuda en la depuraci贸n, que en tiempo de ejecuci贸n muestra llamadas incorrectas a la API OpenCL y problemas de acceso a la memoria.

Opencl


Para comenzar con OpenCL, debe obtener una lista de las plataformas disponibles. Las plataformas m谩s comunes son AMD, Intel y NVidia.

C贸digo
 std::vector<cl::Platform> platforms; cl::Platform::get(&platforms); 


Luego, despu茅s de elegir una plataforma, debe seleccionar el dispositivo inform谩tico que representa esta plataforma:

C贸digo
 const cl::Platform& platform(platforms[platform_n]); std::vector<cl::Device> all_devices; platform.getDevices(CL_DEVICE_TYPE_ALL, &all_devices); 


Y al final de la fase preparatoria, es necesario crear un contexto y colas dentro de las cuales se asignar谩 la memoria y se realizar谩n los c谩lculos. Por ejemplo, un contexto que combina todos los dispositivos inform谩ticos de una plataforma seleccionada se crea de la siguiente manera:

Contexto y c贸digo de creaci贸n de cola
 cl::Context context(all_devices); std::vector<cl::CommandQueue> queues; for(cl::Device device: all_devices) queues.push_back(cl::CommandQueue(context, device)); 


Para descargar el c贸digo fuente a un dispositivo inform谩tico, debe compilarse, la clase cl :: Program est谩 destinada a esto.

C贸digo de compilaci贸n del n煤cleo
 std::vector< std::string > source_data; cl::Program::Sources sources; for(int i = 0; i != files.size(); ++i) { source_data.push_back(load_program(files[i]));//   sources.push_back(std::make_pair(source_data.back().data(), source_data.back().size())); } cl::Program prog(context, sources); devices.push_back(all_devices); prog.build(devices, options); 


Para describir los par谩metros de una funci贸n (kernel) que se ejecuta en un dispositivo inform谩tico, hay una plantilla cl: make_kernel.

Ejemplo de n煤cleo de c谩lculo de fuerza de interacci贸n
 typedef cl::make_kernel< cl_int, cl_int, //Block offset cl::Buffer, //mass cl::Buffer, //y cl::Buffer, //f cl_int, cl_int, //yoff,foff cl_int, cl_int //points_count,stride > ComputeBlock; 


Adem谩s, todo es simple: declaramos una variable con el tipo de n煤cleo, pasamos el programa compilado y el nombre del n煤cleo computacional, podemos iniciar el n煤cleo casi como una funci贸n normal.

C贸digo de lanzamiento del kernel
 ComputeBlock fcompute(prog, "ComputeBlockLocal"); cl::NDRange global_range(device_data_size); cl::NDRange local_range(block_size); cl::EnqueueArgs eargs(ctx.m_queue, global_range, local_range); fcompute(eargs, ...  ); //,   . 


El n煤cleo computacional para OpenCL en s铆 es muy similar a la opci贸n 'openmp + block + optimization' para la CPU, solo que a diferencia de la versi贸n de la CPU, el primer ciclo se controla usando OpenCL (el rango del ciclo est谩 determinado por la variable global_range del c贸digo de inicio del kernel, y el n煤mero de iteraci贸n actual est谩 disponible desde el kernel utilizando la funci贸n get_global_id (0)). Primero, parte de los datos del cuerpo se cargan desde la memoria local y luego se procesan. La memoria local es com煤n a todos los subprocesos del grupo, por lo que la descarga se produce una vez para el grupo y cada subproceso del grupo la procesa y, dado que la memoria local es mucho m谩s r谩pida que la global, los c谩lculos son mucho m谩s r谩pidos.

C贸digo del n煤cleo
 __kernel void ComputeBlockLocal(int offset_n1, int offset_n2, __global const nbcoord_t* mass, __global const nbcoord_t* y, __global nbcoord_t* f, int yoff, int foff, int points_count, int stride) { int n1 = get_global_id(0) + offset_n1; __global const nbcoord_t* rx = y + yoff; __global const nbcoord_t* ry = rx + stride; __global const nbcoord_t* rz = rx + 2 * stride; __global const nbcoord_t* vx = rx + 3 * stride; __global const nbcoord_t* vy = rx + 4 * stride; __global const nbcoord_t* vz = rx + 5 * stride; __global nbcoord_t* frx = f + foff; __global nbcoord_t* fry = frx + stride; __global nbcoord_t* frz = frx + 2 * stride; __global nbcoord_t* fvx = frx + 3 * stride; __global nbcoord_t* fvy = frx + 4 * stride; __global nbcoord_t* fvz = frx + 5 * stride; nbcoord_t x1 = rx[n1]; nbcoord_t y1 = ry[n1]; nbcoord_t z1 = rz[n1]; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; __local nbcoord_t x2[NBODY_DATA_BLOCK_SIZE]; __local nbcoord_t y2[NBODY_DATA_BLOCK_SIZE]; __local nbcoord_t z2[NBODY_DATA_BLOCK_SIZE]; __local nbcoord_t m2[NBODY_DATA_BLOCK_SIZE]; // NB! get_local_size(0) == NBODY_DATA_BLOCK_SIZE for(int b2 = 0; b2 < points_count; b2 += NBODY_DATA_BLOCK_SIZE) { int n2 = b2 + offset_n2 + get_local_id(0); // Copy data block to local memory x2[ get_local_id(0) ] = rx[n2]; y2[ get_local_id(0) ] = ry[n2]; z2[ get_local_id(0) ] = rz[n2]; m2[ get_local_id(0) ] = mass[n2]; // Synchronize local work-items copy operations barrier(CLK_LOCAL_MEM_FENCE); nbcoord_t local_res_x = 0.0; nbcoord_t local_res_y = 0.0; nbcoord_t local_res_z = 0.0; for(int local_n2 = 0; local_n2 != NBODY_DATA_BLOCK_SIZE; ++local_n2) { nbcoord_t dx = x1 - x2[local_n2]; nbcoord_t dy = y1 - y2[local_n2]; nbcoord_t dz = z1 - z2[local_n2]; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = (m2[local_n2]) / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; local_res_x -= dx; local_res_y -= dy; local_res_z -= dz; } // Synchronize local work-items computations barrier(CLK_LOCAL_MEM_FENCE); res_x += local_res_x; res_y += local_res_y; res_z += local_res_z; } frx[n1] = vx[n1]; fry[n1] = vy[n1]; frz[n1] = vz[n1]; fvx[n1] = res_x; fvy[n1] = res_y; fvz[n1] = res_z; } 


Cuda


La implementaci贸n de la plataforma NVidia CUDA es un poco m谩s simple que OpenCL, no necesitamos crear el contexto del dispositivo y administrar la cola de ejecuci贸n nosotros mismos (al menos hasta que queramos hacer una implementaci贸n multi-GPU). Como en el caso de OpenCL, necesitamos asignar memoria en la GPU, copiar nuestros datos en ella y luego podemos iniciar el n煤cleo inform谩tico.

Puede leer m谩s sobre trabajar con CUDA aqu铆 .

C贸digo de lanzamiento del kernel de CUDA
 dim3 grid(count / block_size); dim3 block(block_size); size_t shared_size(4 * sizeof(nbcoord_t) * block_size); kfcompute <<< grid, block, shared_size >>> (... ...); 


A diferencia de OpenCL, CUDA no especifica el rango completo de iteraciones (en la implementaci贸n de OpenCL es global_range), pero establece el tama帽o de cuadr铆cula y los tama帽os de bloque en la cuadr铆cula, en consecuencia, el c谩lculo del n煤mero de cuerpo actual cambia ligeramente, de lo contrario, el n煤cleo es muy similar a OpenCL, con la excepci贸n de otros nombres funciones de sincronizaci贸n y especificaci贸n para memoria compartida. Otra caracter铆stica distintiva 煤til de CUDA es que podemos especificar el tama帽o requerido de memoria compartida cuando se inicia el n煤cleo. Al igual que en la implementaci贸n de OpenCL, al comienzo de cada bloque de iteraci贸n copiamos parte de los datos en la memoria compartida y luego trabajamos con esta memoria desde todos los hilos del bloque.

C贸digo del n煤cleo CUDA
 __global__ void kfcompute(int offset_n2, const nbcoord_t* y, int yoff, nbcoord_t* f, int foff, const nbcoord_t* mass, int points_count, int stride) { int n1 = blockDim.x * blockIdx.x + threadIdx.x; const nbcoord_t* rx = y + yoff; const nbcoord_t* ry = rx + stride; const nbcoord_t* rz = rx + 2 * stride; nbcoord_t x1 = rx[n1]; nbcoord_t y1 = ry[n1]; nbcoord_t z1 = rz[n1]; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; extern __shared__ nbcoord_t shared_xyzm_buf[]; nbcoord_t* x2 = shared_xyzm_buf; nbcoord_t* y2 = x2 + blockDim.x; nbcoord_t* z2 = y2 + blockDim.x; nbcoord_t* m2 = z2 + blockDim.x; for(int b2 = 0; b2 < points_count; b2 += blockDim.x) { int n2 = b2 + offset_n2 + threadIdx.x; // Copy data block to local memory x2[ threadIdx.x ] = rx[n2]; y2[ threadIdx.x ] = ry[n2]; z2[ threadIdx.x ] = rz[n2]; m2[ threadIdx.x ] = mass[n2]; // Synchronize local work-items copy operations __syncthreads(); nbcoord_t local_res_x = 0.0; nbcoord_t local_res_y = 0.0; nbcoord_t local_res_z = 0.0; for(int n2 = 0; n2 != blockDim.x; ++n2) { nbcoord_t dx = x1 - x2[n2]; nbcoord_t dy = y1 - y2[n2]; nbcoord_t dz = z1 - z2[n2]; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = (m2[n2]) / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; local_res_x -= dx; local_res_y -= dy; local_res_z -= dz; } // Synchronize local work-items computations __syncthreads(); res_x += local_res_x; res_y += local_res_y; res_z += local_res_z; } n1 += foff; f[n1 + 3 * stride] = res_x; f[n1 + 4 * stride] = res_y; f[n1 + 5 * stride] = res_z; } 


Tabla 3. La dependencia de la funci贸n de tiempo de c谩lculo (en segundos) fn sobre la cantidad de cuerpos que interact煤an N para diversas implementaciones de GPU y mejores implementaciones de CPU
N40968192163843276865536131072
openmp + block + optimizaci贸n0,01280,04950,1940,774------
OpenCL + half NVidia K800.0040.0080,0260,1340.3221,18
CUDA + half NVidia K800.0040.0080,02450,1150.2911.13

Donde conseguir NVidia K80
Las implementaciones de OpenCL y CUDA se lanzaron en el servicio gratuito Colab de Google, que proporciona acceso a las computadoras NVidia K80. Puede leer m谩s sobre c贸mo trabajar con este servicio en el concentrador

En general, el resultado es decepcionante, solo 5-6 veces m谩s r谩pido que la implementaci贸n de la CPU. Incluso si realizamos c谩lculos en todo el K80, obtendremos una aceleraci贸n de hasta 12 veces, pero desde Dado que la complejidad de la tarea es cuadr谩tica, en un tiempo razonable podremos procesar no 32768 cuerpos interactuando, sino 131072, que es solo 4 veces m谩s.

Aproximaci贸n de funciones fn


Si observa de cerca la funci贸n, que se establece por la fuerza de atracci贸n de dos cuerpos, est谩 claro que disminuye cuadr谩ticamente con la distancia. Por lo tanto, podemos calcular con precisi贸n la fuerza de interacci贸n entre cuerpos cercanos, y aproximadamente entre cuerpos distantes. Un enfoque famoso
es un algoritmo de c贸digo de 谩rbol propuesto por D. Barnes y P. Hat. Se construye un 谩rbol de octree en el espacio simulado, que contiene en sus hojas las coordenadas y masas de los cuerpos simulados. Los nodos principales contienen el centro de masa, la masa total de los nodos secundarios y el radio de la esfera descrita alrededor de los cuerpos de los nodos secundarios. La ra铆z del 谩rbol contiene el centro de masa de todos los cuerpos, su masa total y el radio de la esfera que se describe a su alrededor. Al calcular la fuerza de interacci贸n, la distancia a la ra铆z del 谩rbol se considera primero si la relaci贸n de la distancia al nodo a su radio es mayor que alguna constante  lambdacrit , entonces la ra铆z se considera un cuerpo con coordenadas iguales a las coordenadas del centro de masa de los cuerpos incluidos en 茅l, y una masa igual a la suma de las masas de los nodos hijos, si la relaci贸n es menor o igual a  lambdacrit , el procedimiento se repite recursivamente para cada nodo secundario. Si se alcanza la hoja de un 谩rbol, la fuerza de interacci贸n se considera de la manera habitual. Por lo tanto, si un cuerpo est谩 muy alejado del grupo compacto de otros cuerpos, este grupo se le presenta como un cuerpo, y la fuerza de interacci贸n se calcula no hasta cada cuerpo, sino solo hasta un cuerpo. Debido a esto, la complejidad del algoritmo disminuye con O(N2) antes O(N log(N)) a costa de alguna p茅rdida de precisi贸n.

En mi enfoque, decid铆 no usar el 谩rbol de octree , sino el 谩rbol de kd, porque es m谩s f谩cil de usar y tiene menos gastos de almacenamiento en comparaci贸n con el octree.

Volvamos a la implementaci贸n en la CPU. Un nodo de 谩rbol kd se puede representar en forma de una clase que contiene punteros a los descendientes izquierdo y derecho e informaci贸n sobre las coordenadas y la masa:

Nodo de 谩rbol Kd
 class node { node* m_left; //!<   node* m_right; //!<   nbvertex_t m_mass_center; //!<     nbcoord_t m_mass; //!<   nbcoord_t m_radius_sqr; //!<    ,   lambda_crit nbvertex_t m_bmin; //!<     nbvertex_t m_bmax; //!<     size_t m_body_n; //!<  ,    }; 


Con este m茅todo de almacenamiento del 谩rbol, tenemos dos opciones posibles para atravesar el 谩rbol: usar la recursividad expl铆cita o usar la pila nosotros mismos. Me decid铆 por la segunda opci贸n.

C谩lculo de la fuerza de interacci贸n atravesando un 谩rbol
 nbvertex_t force_compute(const nbvertex_t& v1, const nbcoord_t mass1) { nbvertex_t total_force; node* stack_data[MAX_STACK_SIZE] = {}; node** stack = stack_data; node** stack_head = stack; *stack++ = m_root; while(stack != stack_head) { node* curr = *--stack; const nbcoord_t distance_sqr((v1 - curr->m_mass_center).norm()); if(distance_sqr > curr->m_radius_sqr) {//  ,       total_force += force(v1, curr->m_mass_center, mass1, curr->m_mass); } else {//   ,      if(curr->m_right != NULL) { *stack++ = curr->m_right; } if(curr->m_left != NULL) { *stack++ = curr->m_left; } } } return total_force; } 


Como en el caso de la implementaci贸n de CPU "exacta", se llama a la funci贸n de c谩lculo de fuerza para cada cuerpo. El ciclo en todos los cuerpos se puede paralelizar f谩cilmente con las directivas OpenMP.

Pero las iteraciones de bucle vecino en este caso se referir谩n a partes completamente diferentes del 谩rbol, lo que no permite el uso eficiente de la memoria cach茅 del procesador. Para superar este problema, todos los cuerpos deben atravesarse no en el orden original, sino en el orden en que se ubican en las hojas del 谩rbol kd, entonces se producir谩n iteraciones vecinas para los cuerpos que est谩n cerca del espacio, y atravesar谩n el 谩rbol a lo largo de casi los mismos caminos.

Transversal de la hoja del 谩rbol
 template<class Visitor> void traverse(Visitor visit) { node* stack_data[MAX_STACK_SIZE] = {}; node** stack = stack_data; node** stack_head = stack; *stack++ = m_root; while(stack != stack_head) { node* curr = *--stack; if(curr->m_radius_sqr > 0) {//  .    . if(curr->m_left != NULL) { *stack++ = curr->m_left; } if(curr->m_right != NULL) { *stack++ = curr->m_right; } } else {//   .  . visit(curr->m_body_n, curr->m_mass_center, curr->m_mass); } } } 


Esta implementaci贸n tiene otro problema: no hay una forma universal de paralelizar tal recorrido de 谩rbol. Pero podemos cambiar completamente la forma en que se almacena el 谩rbol en la memoria: podemos almacenar todos los nodos en una matriz lineal y abandonar por completo el almacenamiento de punteros a los descendientes, por analog铆a con la construcci贸n de un mont贸n binario . Al comenzar la numeraci贸n de nodos con 1 si el nodo est谩 en el n煤mero de celda i , entonces su hijo izquierdo est谩 en la celda 2i , ni帽o derecho en la celda 2i+1 y el padre en la celda i/2 . El nodo derecho correspondiente a la izquierda con el n煤mero i tendr谩 un n煤mero i+1 . La matriz en s铆 tendr谩 una longitud 2N , y todos los nodos hoja se ubicar谩n en el 煤ltimo N Adem谩s, las celdas, los nodos muy espaciados se ubicar谩n en celdas cercanas de la matriz. La funci贸n transversal del 谩rbol cambiar谩 un poco, pero el marco sigue siendo el mismo:

C谩lculo de la fuerza atravesando un 谩rbol en una matriz
 nbvertex_t force_compute(const nbvertex_t& v1, const nbcoord_t mass1) { nbvertex_t total_force; size_t stack_data[MAX_STACK_SIZE] = {}; size_t* stack = stack_data; size_t* stack_head = stack; *stack++ = NBODY_HEAP_ROOT_INDEX; while(stack != stack_head) { size_t curr = *--stack; const nbcoord_t distance_sqr((v1 - m_mass_center[curr]).norm()); if(distance_sqr > m_radius_sqr[curr]) { total_force += force(v1, m_mass_center[curr], mass1, m_mass[curr]); } else { size_t left(left_idx(curr)); size_t rght(rght_idx(curr)); if(rght < m_body_n.size()) { *stack++ = rght; } if(left < m_body_n.size()) { *stack++ = left; } } } return total_force; } 


Pero estas no son todas las posibilidades que nos abre el almacenamiento de los nodos en la matriz: podemos rechazar la pila al rastrear. Para hacer esto, en la rama de c贸digo en la que vamos a los hijos del nodo, agregamos la funci贸n para calcular el siguiente nodo ( nextup ), y en la rama en la que calculamos la fuerza de interacci贸n, agregamos el c谩lculo del siguiente nodo con el sub谩rbol actual omitido ( s k i p d o w n )

Para omitir el sub谩rbol actual, debemos bajar a la ra铆z (direcci贸n D ), mientras estamos en el descendiente derecho, tan pronto como lleguemos a la izquierda, vamos al sub谩rbol derecho correspondiente (direcci贸n R ), si llegamos a la ra铆z, se completa el recorrido del 谩rbol.

C贸digo de funci贸n de omisi贸n de sub谩rbol
 index_t skip_down(index_t idx) { // While index is 'right' -> go down while(is_right(idx)) { index_t parent = parent_idx(idx); // We at root again. Stop traverse. if(parent == NBODY_HEAP_ROOT_INDEX) { return NBODY_HEAP_ROOT_INDEX; } idx = parent; } return left2right(idx); } 



Figura 8. Saltarse un sub谩rbol s k i p d o w n .

Para ir al siguiente nodo, si es posible, vaya al ni帽o izquierdo (direcci贸n U ), y si no hay descendiente, vaya al siguiente nodo 'desde abajo' utilizando la funci贸n s k i p d o w n .

Ir al siguiente c贸digo de funci贸n de nodo
 index_t next_up(index_t idx, index_t tree_size) { index_t left = left_idx(idx); if(left < tree_size) { return left; } return skip_down(idx); } 



Figura 9. Transiciones al siguiente nodo n e x t u p .

Puede parecer que reemplazamos la recursi贸n con un bucle w h i l e en funci贸n s k i p d o w n , y dicho reemplazo no da nada, pero veamos c贸mo determinar si un nodo con un n煤mero yo descendiente derecho Esto se puede hacer simplemente comprobando su n煤mero impar (el elemento secundario correcto est谩 en la celda 2 i + 1 ), para esto es suficiente calcular i \ & 1 . Es decir dividimos el n煤mero i dos si el bit menos significativo se establece en uno. Pero esto se puede hacer sin un bucle, en muchos procesadores hay una instrucci贸n find first set que devuelve la posici贸n del primer bit set, us谩ndolo convertimos el bucle en cuatro instrucciones de procesador.

C贸digo de funci贸n de sub谩rbol omitido optimizado
 index_t skip_down(index_t idx) { idx = idx >> (__builtin_ffs(~idx) - 1); return left2right(idx); } 


Despu茅s de eso, podemos excluir la pila de la funci贸n transversal del 谩rbol y reemplazarla con un par skipdown+nextup , despu茅s de eso la funci贸n se ve a煤n m谩s simple:

C谩lculo de la fuerza atravesando un 谩rbol en una matriz sin usar una pila
 nbvertex_t force_compute(const nbvertex_t& v1, const nbcoord_t mass1) const { nbvertex_t total_force; size_t curr = NBODY_HEAP_ROOT_INDEX; size_t tree_size = m_mass_center.size(); do { const nbcoord_t distance_sqr((v1 - m_mass_center[curr]).norm()); if(distance_sqr > m_radius_sqr[curr]) { total_force += force(v1, m_mass_center[curr], mass1, m_mass[curr]); curr = skip_down(curr); } else { curr = next_up(curr, tree_size); } } while(curr != NBODY_HEAP_ROOT_INDEX); return total_force; } 


En total, obtuvimos seis combinaciones de recorrido de 谩rbol y c谩lculo de fuerza. Compare estos enfoques en t茅rminos de tiempo de c谩lculo y calidad. Tomamos como medida de calidad el cambio relativo en la energ铆a total del sistema despu茅s de 100 iteraciones. Como sistema modelo, tomamos dos "galaxias" interactivas que consisten en16384cuerpos cada uno.

Tabla 4. Combinaciones de m茅todo de recorrido de 谩rbol y c谩lculo de fuerza
Tipo de c谩lculo transversal del 谩rbol / fuerza脕rbol con pila'Pila' con una pila'Mont贸n' sin pila
Iteraciones por n煤mero de cuerpociclo + 谩rbolciclo + mont贸nciclo + heapstackless
Desv铆o de hojas谩rbol + 谩rbol anidadonestedtree + heapnestedtree + heapstackless





a) Dependencia del tiempo de c谩lculo de la funci贸n fn desde la relaci贸n de la distancia cr铆tica al nodo del 谩rbol a su radio ( crit )b) La dependencia del cambio relativo en energ铆a en el sistema en la relaci贸n de la distancia cr铆tica al nodo del 谩rbol a su radio ( crit )
Figura 10. Resultados del c谩lculo de funciones fn de varias maneras en la CPU (n煤mero de cuerpos N=32768 )

Se puede ver que la implementaci贸n de '谩rbol anidado + 谩rbol' est谩 irremediablemente atrasado en velocidad, porque Carece de concurrencia. Las implementaciones con la ubicaci贸n de los nodos de 谩rbol en la matriz y la indexaci贸n como en un mont贸n binario est谩n a la cabeza. El cambio relativo en la energ铆a es insignificante para todas las variantes concrit>1 .Tambi茅n en la fig. 10a muestra que para (crit<10 ) todas las variantes del c谩lculo de funciones fn adelanta significativamente en velocidad la versi贸n m谩s optimizada del c谩lculo exacto ('openmp + block + optimization'), con un aumento adicional crit Las implementaciones con un 谩rbol pierden la versi贸n exacta.

Recorrido de 谩rbol de GPU


Trat茅 de evitar el 谩rbol en la GPU usando la tecnolog铆a OpenCL y CUDA. La opci贸n de almacenar nodos en forma de 谩rbol se descart贸 de inmediato, y solo quedaban las opciones para almacenar el 谩rbol en una matriz con indexaci贸n como en un mont贸n binario. En general, la implementaci贸n del n煤cleo inform谩tico no es muy diferente de la versi贸n de la CPU.

Kernel OpenCL para calcular la fuerza al atravesar un 谩rbol (transversal en el orden de los cuerpos de numeraci贸n)
 __kernel void ComputeTreeBH(int offset_n1, int points_count, int tree_size, __global const nbcoord_t* y, __global nbcoord_t* f, __global const nbcoord_t* tree_cmx, __global const nbcoord_t* tree_cmy, __global const nbcoord_t* tree_cmz, __global const nbcoord_t* tree_mass, __global const nbcoord_t* tree_crit_r2) { int n1 = get_global_id(0) + offset_n1; int stride = points_count; __global const nbcoord_t* rx = y; __global const nbcoord_t* ry = rx + stride; __global const nbcoord_t* rz = rx + 2 * stride; __global const nbcoord_t* vx = rx + 3 * stride; __global const nbcoord_t* vy = rx + 4 * stride; __global const nbcoord_t* vz = rx + 5 * stride; __global nbcoord_t* frx = f; __global nbcoord_t* fry = frx + stride; __global nbcoord_t* frz = frx + 2 * stride; __global nbcoord_t* fvx = frx + 3 * stride; __global nbcoord_t* fvy = frx + 4 * stride; __global nbcoord_t* fvz = frx + 5 * stride; nbcoord_t x1 = rx[n1]; nbcoord_t y1 = ry[n1]; nbcoord_t z1 = rz[n1]; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; int stack_data[MAX_STACK_SIZE] = {}; int stack = 0; int stack_head = stack; stack_data[stack++] = NBODY_HEAP_ROOT_INDEX; while(stack != stack_head) { int curr = stack_data[--stack]; nbcoord_t dx = x1 - tree_cmx[curr]; nbcoord_t dy = y1 - tree_cmy[curr]; nbcoord_t dz = z1 - tree_cmz[curr]; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 > tree_crit_r2[curr]) { if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = tree_mass[curr] / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; res_x -= dx; res_y -= dy; res_z -= dz; } else { int left = left_idx(curr); int rght = rght_idx(curr); if(left < tree_size) { stack_data[stack++] = left; } if(rght < tree_size) { stack_data[stack++] = rght; } } } frx[n1] = vx[n1]; fry[n1] = vy[n1]; frz[n1] = vz[n1]; fvx[n1] = res_x; fvy[n1] = res_y; fvz[n1] = res_z; } 


En la primera versi贸n, el recorrido del 谩rbol comenz贸 en el orden de numeraci贸n de los cuerpos en la matriz original, por lo que los subprocesos vecinos omitieron partes completamente diferentes del 谩rbol, lo que afect贸 negativamente el rendimiento de la memoria cach茅 de la GPU. Por lo tanto, en la segunda realizaci贸n, se aplic贸 un recorrido, comenzando desde la parte superior del 谩rbol, en este caso, los hilos vecinos comienzan a atravesar el 谩rbol a lo largo del mismo camino, porque las copas de los 谩rboles vecinos est谩n cerca y en el espacio. Tambi茅n es importante que elijamos la numeraci贸n en el conjunto de nodos del 谩rbol no desde cero, sino desde uno, en este caso las hojas del 谩rbol se almacenan en la segunda mitad del conjunto, y con el n煤mero de cuerpos igual a la potencia de dos, tendremos el mismo acceso a la memoria por el 铆ndice tn1.

Kernel OpenCL para calcular la fuerza atravesando un 谩rbol (atravesando los nodos de numeraci贸n de un 谩rbol)
 __kernel void ComputeHeapBH(int offset_n1, int points_count, int tree_size, __global const nbcoord_t* y, __global nbcoord_t* f, __global const nbcoord_t* tree_cmx, __global const nbcoord_t* tree_cmy, __global const nbcoord_t* tree_cmz, __global const nbcoord_t* tree_mass, __global const nbcoord_t* tree_crit_r2, __global const int* body_n) { int tree_offset = points_count - 1 + NBODY_HEAP_ROOT_INDEX; int stride = points_count; int tn1 = get_global_id(0) + offset_n1 + tree_offset; int n1 = body_n[tn1]; nbcoord_t x1 = tree_cmx[tn1]; nbcoord_t y1 = tree_cmy[tn1]; nbcoord_t z1 = tree_cmz[tn1]; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; int stack_data[MAX_STACK_SIZE] = {}; int stack = 0; int stack_head = stack; stack_data[stack++] = NBODY_HEAP_ROOT_INDEX; while(stack != stack_head) { int curr = stack_data[--stack]; nbcoord_t dx = x1 - tree_cmx[curr]; nbcoord_t dy = y1 - tree_cmy[curr]; nbcoord_t dz = z1 - tree_cmz[curr]; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 > tree_crit_r2[curr]) { if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = tree_mass[curr] / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; res_x -= dx; res_y -= dy; res_z -= dz; } else { int left = left_idx(curr); int rght = rght_idx(curr); if(left < tree_size) { stack_data[stack++] = left; } if(rght < tree_size) { stack_data[stack++] = rght; } } } __global const nbcoord_t* vx = y + 3 * stride; __global const nbcoord_t* vy = y + 4 * stride; __global const nbcoord_t* vz = y + 5 * stride; __global nbcoord_t* frx = f; __global nbcoord_t* fry = frx + stride; __global nbcoord_t* frz = frx + 2 * stride; __global nbcoord_t* fvx = frx + 3 * stride; __global nbcoord_t* fvy = frx + 4 * stride; __global nbcoord_t* fvz = frx + 5 * stride; frx[n1] = vx[n1]; fry[n1] = vy[n1]; frz[n1] = vz[n1]; fvx[n1] = res_x; fvy[n1] = res_y; fvz[n1] = res_z; } 


Al recorrer el orden de numeraci贸n de los nodos del 谩rbol, obtuvimos un aumento en el rendimiento. Pero esta opci贸n tambi茅n se puede mejorar. La memoria global en la que se encuentran actualmente los nodos del 谩rbol est谩 optimizada para acceso compartido , es decir Los hilos de un grupo deben leer palabras ubicadas en un bloque de memoria. En nuestro caso, el recorrido del 谩rbol comienza a lo largo de los mismos caminos, y solicitamos los mismos datos con todos los hilos del grupo, pero cuanto m谩s profundizamos en el 谩rbol, los caminos de los hilos vecinos divergen cada vez m谩s, y tenemos que solicitar datos diferentes, lo que reduce el rendimiento del subsistema varias veces memoria Pero los nodos de cada sub谩rbol que pertenecen al mismo nivel se encuentran en celdas de memoria relativamente cercanas. Es decirAl atravesar el resto del 谩rbol, los subprocesos vecinos del n煤cleo inform谩tico no acceden a los mismos nodos del 谩rbol, sino que se encuentran en la memoria. Para optimizar este acceso a la memoria, se puede utilizar la memoria de textura. Pero hay un inconveniente. Por el momento, las texturas no admiten trabajar con datos de doble precisi贸n (queremos calcular con precisi贸n). Pero en CUDA hay una funci贸n __hiloint2double que recopila un n煤mero de doble precisi贸n de dos enteros.

C贸digo de solicitud de n煤mero de precisi贸n doble de una textura que almacena enteros
 template<> struct nb1Dfetch<double> { typedef double4 vec4; static __device__ double fetch(cudaTextureObject_t tex, int i) { int2 p(tex1Dfetch<int2>(tex, i)); return __hiloint2double(py, px); } static __device__ vec4 fetch4(cudaTextureObject_t tex, int i) { int ii(2 * i); int4 p1(tex1Dfetch<int4>(tex, ii)); int4 p2(tex1Dfetch<int4>(tex, ii + 1)); vec4 d4 = {__hiloint2double(p1.y, p1.x), __hiloint2double(p1.w, p1.z), __hiloint2double(p2.y, p2.x), __hiloint2double(p2.w, p2.z) }; return d4; } }; 


En este caso, se realizaron dos implementaciones, en una cada elemento del 谩rbol (x, y, z, tree_crit_r2) se solicit贸 de forma independiente, y en la segunda implementaci贸n se combinaron estas solicitudes. La solicitud de la masa del nodo se produce con mucha menos frecuencia, solo si se cumple la condici贸n r2> tree_crit_r2 [curr] , por lo que no tiene sentido combinar esta solicitud con las dem谩s. Otra caracter铆stica 煤til del marco CUDA es la capacidad de controlar la proporci贸n de los tama帽os del cach茅 L1 y el tama帽o de la memoria compartida ( cudaFuncSetCacheConfig ). En el caso del recorrido del 谩rbol, no usamos memoria compartida, por lo que podemos aumentar el cach茅 L1 en detrimento de 茅l.

Kernel CUDA para calcular la fuerza atravesando un 谩rbol (atravesando el orden de numeraci贸n de los nodos del 谩rbol)
 __global__ void kfcompute_heap_bh_tex(int offset_n1, int points_count, int tree_size, nbcoord_t* f, cudaTextureObject_t tree_xyzr, cudaTextureObject_t tree_mass, const int* body_n) { nb1Dfetch<nbcoord_t> tex; int tree_offset = points_count - 1 + NBODY_HEAP_ROOT_INDEX; int stride = points_count; int tn1 = blockDim.x * blockIdx.x + threadIdx.x + offset_n1 + tree_offset; int n1 = body_n[tn1]; nbvec4_t xyzr = tex.fetch4(tree_xyzr, tn1); nbcoord_t x1 = xyzr.x; nbcoord_t y1 = xyzr.y; nbcoord_t z1 = xyzr.z; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; int stack_data[MAX_STACK_SIZE] = {}; int stack = 0; int stack_head = stack; stack_data[stack++] = NBODY_HEAP_ROOT_INDEX; while(stack != stack_head) { int curr = stack_data[--stack]; nbvec4_t xyzr2 = tex.fetch4(tree_xyzr, curr); nbcoord_t dx = x1 - xyzr2.x; nbcoord_t dy = y1 - xyzr2.y; nbcoord_t dz = z1 - xyzr2.z; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 > xyzr2.w) { if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = tex.fetch(tree_mass, curr) / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; res_x -= dx; res_y -= dy; res_z -= dz; } else { int left = nbody_heap_func<int>::left_idx(curr); int rght = nbody_heap_func<int>::rght_idx(curr); if(left < tree_size) { stack_data[stack++] = left; } if(rght < tree_size) { stack_data[stack++] = rght; } } } f[n1 + 3 * stride] = res_x; f[n1 + 4 * stride] = res_y; f[n1 + 5 * stride] = res_z; } 


Un an谩lisis del programa en el perfilador nvprof mostr贸 que incluso cuando se usa memoria de textura para almacenar un 谩rbol, todav铆a hay una carga muy alta en la memoria global.

De hecho, en CUDA, toda la memoria del n煤cleo que se dirige a direcciones 'calculadas' se almacena en la memoria global y, en consecuencia, la pila que se necesita para atravesar el 谩rbol se encuentra en la memoria global y consume una parte significativa del ancho de banda del chip de memoria, porque la pila tiene cada hilo de ejecuci贸n, y hay muchos hilos.

Pero, afortunadamente, ya sabemos c贸mo atravesar un 谩rbol sin usar una pila. Complementando el n煤cleo computacional anterior con las funciones de calcular el siguiente nodo de 谩rbol, obtenemos un nuevo n煤cleo, adem谩s, m谩s compacto.

N煤cleo CUDA para la fuerza inform谩tica al atravesar un 谩rbol sin usar una pila
 __global__ void kfcompute_heap_bh_stackless(int offset_n1, int points_count, int tree_size, nbcoord_t* f, cudaTextureObject_t tree_xyzr, cudaTextureObject_t tree_mass, const int* body_n) { nb1Dfetch<nbcoord_t> tex; int tree_offset = points_count - 1 + NBODY_HEAP_ROOT_INDEX; int stride = points_count; int tn1 = blockDim.x * blockIdx.x + threadIdx.x + offset_n1 + tree_offset; int n1 = body_n[tn1]; nbvec4_t xyzr = tex.fetch4(tree_xyzr, tn1); nbcoord_t x1 = xyzr.x; nbcoord_t y1 = xyzr.y; nbcoord_t z1 = xyzr.z; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; int curr = NBODY_HEAP_ROOT_INDEX; do { nbvec4_t xyzr2 = tex.fetch4(tree_xyzr, curr); nbcoord_t dx = x1 - xyzr2.x; nbcoord_t dy = y1 - xyzr2.y; nbcoord_t dz = z1 - xyzr2.z; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 > xyzr2.w) { if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = tex.fetch(tree_mass, curr) / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; res_x -= dx; res_y -= dy; res_z -= dz; curr = nbody_heap_func<int>::skip_idx(curr); } else { curr = nbody_heap_func<int>::next_up(curr, tree_size); } } while(curr != NBODY_HEAP_ROOT_INDEX); f[n1 + 3 * stride] = res_x; f[n1 + 4 * stride] = res_y; f[n1 + 5 * stride] = res_z; } 


El rendimiento de los n煤cleos que se ejecutan en la GPU depende en gran medida del tama帽o de los bloques en los que dividimos la tarea. Dependiendo de este tama帽o, cu谩ntos registros, memoria local y otros recursos estar谩n disponibles para cada hilo inform谩tico. Tambi茅n debe tener en cuenta que mientras espera el acceso a la memoria en un subproceso, otro subproceso puede realizar c谩lculos en el procesador del sombreador, por lo tanto, con un n煤mero suficiente de subprocesos ejecutados simult谩neamente en un procesador, el tiempo de acceso a la memoria estar谩 oculto detr谩s de los c谩lculos. Por lo tanto, antes de comparar el rendimiento de nuestros n煤cleos, necesitamos calcular el tama帽o de bloque 贸ptimo para cada uno de ellos. Hagamos una comparaci贸n de la mitad disponible de NVidia K80.

Tabla 5. La dependencia de la funci贸n de tiempo de c谩lculo (en segundos) fn GPU N=131072 crit=10
/81632641282565121024
opencl+dense5.772.841.461.131.151.141.141.13
cuda+dense5.442.551.270.960.970.990.99-
opencl+heap+cycle5.885.655.745.965.375.385.355.38
opencl+heap+nested4.543.683.985.255.465.415.485.31
cuda+heap+nested3,622,812.684.264.844.754.84.67
cuda+heap+nested+tex2.61.510.9120.71.851.751.691.61
cuda+heap+nested+tex+stackless2.31.290.7730.50.510.520.520.52


6. ( ) fn GPU N=1M c r i t = 4
Tama帽o de bloque / n煤cleo81632641282565121024
opencl + denso36617989,969,370,369,168,968,0
cuda + denso34616279,660,860,860,759,6-
opencl + mont贸n + ciclo16,218,220,121,221,221,321,221,1
opencl + mont贸n + anidado10,57.636.388.239,959,899.659.66
cuda + mont贸n + anidado8.676.445.395.938.658.618.418.27
cuda + mont贸n + anidado + tex6.383,572,131,443,563.463.303.29
cuda + mont贸n + anidado + tex + sin pila5,783.191,831,211.111.101.111.13

Una situaci贸n dif铆cil, pero, a diferencia de la versi贸n de CPU del recorrido del 谩rbol, est谩 claro que cada paso de optimizaci贸n trae resultados tangibles. La implementaci贸n de 'opencl + heap + cycle' es casi 6 veces m谩s lenta que una soluci贸n exacta con c谩lculo de funci贸n completaf n .Una implementaci贸n de 'opencl + heap + nested', en la que un recorrido del 谩rbol comienza desde los nodos vecinos, 1.4 veces m谩s r谩pido que el anterior, porque La memoria cach茅 se usa mejor. En la implementaci贸n de 'cuda + heap + nested', el cach茅 L1 se aument贸 en detrimento de la memoria compartida, lo que aument贸 la velocidad en 1,4 veces, aunque es posible que en la implementaci贸n de cuda el n煤cleo de compilaci贸n se compile de manera m谩s 贸ptima (en las versiones de 'opencl + dense' y 'cuda Los n煤cleos + densos son id茅nticos, y el rendimiento de la versi贸n de Cuda es ~ 1.2 veces mayor). El aumento m谩s significativo en la velocidad computacional (3.8 veces) se logra cuando el 谩rbol se encuentra en la memoria de textura y se combinan las consultas a los elementos del nodo del 谩rbol. La implementaci贸n con el recorrido del 谩rbol sin usar la pila 'cuda + heap + nested + tex + stackless' es 1.4 veces m谩s r谩pida tambi茅n.en 茅l, todo el ancho de banda del bus de memoria se usa solo para acceder a datos sobre nodos de 谩rbol y no se gasta en la pila. As铆, conc r i t = 10 fue posible lograr la aceleraci贸n dos veces en comparaci贸n con el c谩lculo completo de la funci贸nf n . Pero c r i t = 10 es un valor excesivamente grande del par谩metro; en el gr谩fico del cambio relativo de energ铆a en el sistema versus la relaci贸n de la distancia cr铆tica al nodo del 谩rbol a su radio para la implementaci贸n de la CPU, est谩 claro que los valores m谩s bajosc r i t sin aparente p茅rdida de precisi贸n de la soluci贸n. Intentemos variarc r i t en los tama帽os de bloque 贸ptimos, que determinamos en el paso anterior.




a) N = 128 K
b) N = 1 M
Fig. 11. Dependencia del tiempo de c谩lculo de la funci贸n. f n desde la relaci贸n de la distancia cr铆tica al nodo del 谩rbol a su radio ( c r i t ) para varias implementaciones de GPU Tree Walk

Se puede ver que para peque帽os c r i t todos los m茅todos de c谩lculo de la funci贸nf n para cerrar valores determinados por el tiempo de construcci贸n del 谩rbol kd y la preparaci贸n de datos para la GPU. Adem谩s, el tiempo de construcci贸n de un 谩rbol contribuye significativamente al tiempo total parac r i t4 , entonces este tiempo puede ser descuidado. Es interesante notar que cuandoN = 128 K el rendimiento mejora nuevamente cuando se alcanzac r i t = 1024 , lo m谩s probable es que esto se deba al hecho de que todos los hilos de la GPU rodean los mismos v茅rtices de 谩rbol al mismo tiempo, lo que mejora el uso de cach茅 L1 y elimina por completo la'divergencia de ramificaci贸n'. Tambi茅n se ve que el mejor rendimiento para una implementaci贸n sin una pila (cuda + heap + nested + tex + stackless), supera la versi贸n con la pila aproximadamente1.4 - 1.5 veces. Otras implementaciones son varias veces m谩s lentas. Para consolidar el resultado, calcularemos el tiempo en una GPU con una arquitectura m谩s nueva.

Resultados de lanzamiento en GeForce GTX 1080 Ti (precisi贸n 煤nica)
7. ( ) fn GPU N=1M crit=4
/81632641282565121024
opencl+dense47.823.411.611.5912.812.812.812.8
cuda+dense49.023.811.911.811.711.711.711.7
opencl+heap+cycle7.488.267.737.367.337.277.257.26
opencl+heap+nested1.331.201.412.462.532.492.442.47
cuda+heap+nested1.231.101.312.282.292.282.232.25
cuda+heap+nested+tex0.880.680.6541.61.61.61.61.6
cuda+heap+nested+tex+stackless0.710.470.450.430.430.420.410.395



12. fn ( crit ) GPU


Cuando se usa la GeForce GTX 1080 Ti para el c谩lculo, la diferencia entre las caminatas de 谩rbol con una pila y sin pila es hasta dos veces, siempre que descuidemos el tiempo que lleva construir el 谩rbol. Este hecho nos empuja a garantizar que en el futuro se transfiera a la GPU y se cree un 谩rbol. Una comparaci贸n de las tablas 5-7 muestra que no hay un 煤nico tama帽o de bloque 贸ptimo para un n煤mero diferente de cuerpos, y a煤n m谩s para diferentes arquitecturas de GPU, adem谩s, la diferencia de tiempo de c谩lculo puede alcanzar varias veces, incluso si no tiene en cuenta los valores l铆mite. Por lo tanto, antes de largos c谩lculos, tiene sentido determinar el tama帽o de bloque 贸ptimo para cada configuraci贸n.

Lo principal que hemos logrado es la capacidad de modelar la interacci贸n por pares de poco m谩s de un mill贸n de cuerpos (2 20 ) por un tiempo razonable en uno, no la GPU m谩s nueva. En las GPU m谩s nuevas (Tesla V100), aparentemente, ser谩 posible procesar alrededor de dos millones de cuerpos que interact煤an en un tiempo razonable, porque su rendimiento es aproximadamente cuatro veces mayor que la mitad del Tesla K80. Aunque este n煤mero tambi茅n es incomparable con el n煤mero de estrellas incluso en galaxias enanas, como laPeque帽a Nube de Magallanes, pero, en mi opini贸n, es impresionante.

Conclusi贸n


El uso de m茅todos integrados para resolver ecuaciones diferenciales permite resolver problemas similares con un alto nivel de precisi贸n a un costo de tiempo relativamente peque帽o, y la aplicaci贸n de algoritmos para aproximar la funci贸n de calcular las fuerzas de atracci贸n por pares nos permite procesar vol煤menes colosales de cuerpos que interact煤an. Por lo tanto, a diferencia de los extraterrestres de las "Tareas de tres cuerpos", podemos resolver el problemaN cuerpos, tanto para un peque帽o n煤mero de cuerpos como para grupos de estrellas enteros, aunque a costa de alguna p茅rdida de precisi贸n.

Visualizaci贸n


Para aquellos que han le铆do hasta el final, les dar茅 algunos videos m谩s con visualizaci贸n de la evoluci贸n de los sistemas de los cuerpos.

Simulaci贸n de una colisi贸n de dos galaxias. El n煤mero total de cuerpos es de 60 mil.


Modelando la evoluci贸n de una galaxia que consiste en un mill贸n de estrellas. En el centro hay un cuerpo con una masa del 99% del total. Las part铆culas individuales son casi indistinguibles. Ya m谩s como una ola en una gota de l铆quido. Coloreado seg煤n el m贸dulo de velocidad. Baja velocidad - azul, medio - verde, alto - rojo. Se ve que en el centro la velocidad es mayor, y disminuye gradualmente hasta los bordes, y la m谩s baja en el plano ecuatorial.


Una simulaci贸n m谩s 'din谩mica' de la evoluci贸n de una galaxia de un mill贸n de estrellas. En el centro hay un cuerpo con una masa del 10% del total. El cuerpo central afecta el resto significativamente menos. Primero, las "estrellas" se separan, luego se re煤nen en varios grupos y al final nuevamente forman un gran grupo.

En el proceso de modelado, alrededor del 5% de las "estrellas" abandonaron la regi贸n inicial de manera irrevocable.


En el d茅cimo segundo, se parece mucho a la galaxia de volteretas existente .

El c贸digo del proyecto se puede encontrar en el github .

Source: https://habr.com/ru/post/437014/


All Articles